# Critical behavior of dirty parafermionic chains

Akshat Pandey<EMAIL_ADDRESS>Department of Physics, Stanford
University, Stanford, CA 94305, USA  Aditya Cowsik  Department of Physics,
Stanford University, Stanford, CA 94305, USA

###### Abstract

A family of $\mathbb{Z}_{n}$-symmetric non-Hermitian models of Baxter was
shown by Fendley to be exactly solvable via a parafermionic generalization of
the Clifford algebra. We study these models with spatially random couplings,
and obtain several exact results on thermodynamic singularities as the
distributions of couplings are varied. We find that these singularities,
independent of $n$, are identical to those in the random transverse-field
Ising chain; correspondingly the models host infinite-randomness critical
points. Similarities in structure to exact methods for random Ising models, a
strong-disorder renormalization group, and generalizations to other models
with free spectra, are discussed.

Introduction.— Consider an open chain of length $L$ with an $n$-dimensional
Hilbert space on each site. With each site associate analogues of Pauli $Z$
and $X$ matrices, defined as

$Z=\begin{bmatrix}1&&&\\\ &\omega&&\\\ &&\dots&\\\
&&&\omega^{n-1}\end{bmatrix},\quad X=\begin{bmatrix}0&&&1\\\ 1&0&&\\\
&1&\dots&\\\ &&&0\end{bmatrix}$ (1)

respectively. Here $\omega=e^{2\pi i/n}$. The non-Hermitian Hamiltonian

$H=-\sum_{j=1}^{L-1}J_{j}Z_{j}^{\dagger}Z_{j+1}-\sum_{j=1}^{L}h_{j}X_{j}$ (2)

is then a natural $\mathbb{Z}_{n}$-symmetric generalization of the transverse-
field Ising chain (TFIC). It was introduced in relation with the integrable
chiral Potts model by Baxter [1, 2], who demonstrated the remarkable fact that
$H$ has $n^{L}$ eigenvalues which can be written as

$E_{\\{n_{k}\\}}=-\sum_{k=1}^{L}\omega^{n_{k}}\varepsilon_{k},$ (3)

where each $n_{k}\in\\{0\dots n-1\\}$ is chosen independently. The single-
particle energies $\varepsilon_{k}$ are related to the roots of a particular
polynomial whose construction we will discuss. An understanding of the “free
parafermionic” nature of the spectrum was achieved by Fendley [3], who
provided an explicit construction of the parafermionic analogue of raising and
lowering operators associated with each single-particle level. A Jordan-
Wigner-like transformation [4], unlike in the TFIC, does not make this task at
all trivial: the aforementioned parafermionic operators are highly nonlinear
in these Jordan-Wigner parafermions.

While a number of observables in the homogeneous ($J_{j}=J,h_{j}=h$ for all
$j$) version of the model, especially in the vicinity of $h=J$, have been
calculated [5, 6, 7, 8, 9, 10, 11], that the model is in principle exactly
solvable for inhomogeneous couplings has not been exploited. This is the goal
of the present paper.

We have two entwined motivations. First, given the scarcity of theoretical
tools available to analyse strongly disordered quantum critical points in a
controlled manner—the strong-disorder RG being the most versatile but, by its
nature, limited to cases where infinite-randomness scaling is achieved and
errors accumulated in the approach to said scaling are unimportant—we wish to
understand the critical properties of any disordered models on which exact
methods can be brought to bear. As far as we know, the only problem on which
such progress was made is the disordered TFIC, or equivalently the stripy-
random classical 2d Ising model [12, 13, 14, 15, 16]. Second, while the free
parafermionic models have some features in common with the TFIC, the method of
solution is more intricate, so it is not a priori obvious that the free
parafermionic character of the problem is useful in practice in the presence
of disorder. We show, however, that much of the structure of the Ising papers
referred to above carries over in a non-trivial way to the parafermionic
models. These rather attractive relations may be useful for an understanding
of the algebraic structures underlying “free” statistical mechanics in
general, and indeed for the construction of new disordered critical models.

We leverage these results to analyze the singularity structure of the “ground-
state energy” of $H$, and find that disorder washes out the difference between
the critical behaviors among different $n$.

The ground state energy.— The recipe for finding the single-particle energies
$\varepsilon_{k}$ of (3) is as follows [1, 3]. Write the couplings of the
Hamiltonian as
$(b_{1},b_{2},b_{3},\dots,b_{M})\equiv(h_{1},J_{1},h_{2},\dots,h_{L})$, so
that $M=2L-1$. Define a sequence of polynomials $P_{m}(z)$ using the relations

$\displaystyle P_{0}(z)=P_{-1}(z)=1;$ $\displaystyle
P_{m}(z)=P_{m-1}(z)-zb_{m}^{n}P_{m-2}(z),\quad m\in\\{1\dots M\\}.$ (4)

If all $b_{m}>0$ (we assume this for the rest of the paper), $P_{M}(z)$ has
$L$ positive real zeros $z_{k},k\in\\{1\dots L\\}$. Then
$\varepsilon_{k}=z_{k}^{-1/n}$. We are interested in the state in which all
levels are filled in the $n_{k}=0$ manner, $E_{0}=-\sum_{k}\varepsilon_{k}$.
This is of course the analogue of the ground state energy. Its importance
stems from the fact that it is at the edge of the many-body spectrum in the
complex plane (up to multiplication by $\omega$), and it has singularities as
the distributions of couplings are varied; it is these singularities that we
analyse.

We now show a simple formula for $E_{0}$ in terms of $P_{M}$ which we use for
all subsequent developments. We start by writing

$-E_{0}=\frac{1}{2\pi
i}\int_{C}dz\,z^{-1/n}\frac{P_{M}^{\prime}(z)}{P_{M}(z)},$ (5)

where $C$ is the contour shown in Fig. 1. This converts the zeros $z=z_{k}$ of
$P_{M}$ to poles whose residues contain the requisite factors of
$\varepsilon_{k}$. We can append the contours $C_{U},C_{L}$ since they enclose
no singularities. Their semicircular parts vanish (in the limit of large
radii), so that we can integrate along $C^{\prime}$ which straddles the branch
cut. Integrating by parts, and changing variables to $\Omega=|z|^{-1/n}$, we
obtain

$~{}-E_{0}=\frac{1}{\pi}\sin\frac{\pi}{n}\int_{0}^{\infty}d\Omega\,\ln
P_{M}(-\Omega^{-n}).$ (6)

This formula reflects the free nature of the model: the partition function
factorizes into frequencies which do not talk to each other. For $n=2$,
$E_{0}$ is an appropriate limit of the free energy of a 2d Ising model which
is inhomogeneous in one direction and homogeneous in the other: see the works
by McCoy and Wu [12], and Shankar and Murthy [16]. Indeed the strategy in
these works was to factorize the partition function into wavevectors in the
clean direction (which for us is time). The procedure, however, relied on
manipulations particular to the Ising model (using Pfaffians [17], and Jordan-
Wigner fermions [18], respectively). It is interesting that the weaker notion
of freeness in the models here (and others which we will comment on) is
sufficient to produce the same structures. The connections are also of
practical value, as we will demonstrate presently.

Figure 1: Singularities in the complex $z$ plane of the integrand in
(5)—crosses are poles and the wavy line is a branch cut—and contours involved
in its evaluation.

Singularities of $E_{0}$.— We now study the non-analyticities $E_{0}$ exhibits
as a function of the distributions of couplings. Henceforth all $h_{j}$ will
be set equal to $1$, while $J_{j}$ will be independent random variables drawn
from a distribution $P(J)$. The precise forms of the distributions do not
affect the physics we discuss (as long as tails in the distributions are not
too heavy [19]).

Following the Ising work of Shankar and Murthy [16], we interpret
$P_{M}(-\Omega^{-n})$ as the partition function of a 1d classical Ising model
with random fields. The recursion (4) can be written as

$\begin{bmatrix}P_{m}(-\Omega^{-n})\\\
P_{m-1}(-\Omega^{-n})\end{bmatrix}=\begin{bmatrix}1&b_{m}^{n}\Omega^{-n}\\\
1&0\end{bmatrix}\begin{bmatrix}P_{m-1}(-\Omega^{-n})\\\
P_{m-2}(-\Omega^{-n})\end{bmatrix}.$ (7)

Multiplying the matrices associated with $h_{j+1}(=1)$ and $J_{j}$ yields

$T_{j}(\Omega)=\begin{bmatrix}1+\Omega^{-n}&J_{j}^{n}\Omega^{-n}\\\
1&J_{j}^{n}\Omega^{-n}\end{bmatrix}.$ (8)

Consider now a classical Ising model defined by the Hamiltonian $-\beta
H_{I}=\sum_{j}[\mathcal{E}_{j}+\mathcal{J}(\sigma_{j}\sigma_{j+1}-1)+\mathcal{H}_{j}\sigma_{j}+\mathcal{H}^{\prime}_{j}\sigma_{j+1}]$.
The transfer matrix that assembles these Boltzmann weights is clearly $T_{j}$
if we identify
$\mathcal{E}_{j}=(1/2)\ln[J_{j}^{n}\Omega^{-n}(1+\Omega^{-n})]$,
$\mathcal{J}=\mathcal{H}_{j}=(1/4)\ln(1+\Omega^{-n})$ and
$\mathcal{H}^{\prime}_{j}=(1/4)\ln(1+\Omega^{-n})-(1/2)\ln(J_{j}^{n}\Omega^{-n})$.
The fields can be redistributed among sites so that all
$\mathcal{H}^{\prime}_{j}$ vanish, and
$\mathcal{H}_{j}=-(1/2)\ln(1+\Omega^{-n})+(1/2)\ln(J_{j}^{n}\Omega^{-n})$—we
have also multiplied the field by a minus sign for later convenience. This
puts the Hamiltonian into a more standard random-field Ising form. This sort
of “gauge” freedom in the transfer matrix (namely the ability to modify
couplings as long as $\mathcal{J}$, and all
$\mathcal{H}_{j}+\mathcal{H}^{\prime}_{j-1}$, up to a $j$-independent sign,
remain unchanged—equivalently to insert factors of $1=AA^{-1}$ for appropriate
matrices $A$) is useful to make a number of connections between quantities
most naturally expressed as products of different families of random matrices
that ultimately yield the same polynomials. Still different factorizations
will be exploited in the next two sections.

Returning to the problem at hand, to find the ground state energy density
$E_{0}/L$, we have to integrate a quantity of the form
$f(\Omega)=(1/L)\ln\bra{v}\prod_{j}T_{j}(\Omega)\ket{w}$, where the boundary
vectors $v$ and $w$ can be read off $\eqref{eq:recursion}$, and the product
goes down in $j$ rightward. Each $f(\Omega)$ is a 1d random-field Ising
model’s free energy density, as we have just explained. These are self-
averaging quantities in the $L\to\infty$ limit, and so is $E_{0}/L$.
Singularities in $E_{0}/L$ are inherited from the approach to singularities in
$f(\Omega)$ in the zero-temperature limit, which translates to $\Omega\to 0$.
Shankar and Murthy used results from Refs. [20, 21] to extract these
singularities, and our calculation now closely follows theirs. Let us suppose
$P(J)$ has support for $J\in[J_{1},J_{2}]$. We expect Griffiths-type
singularities when $J_{1}$ crosses $1$ and when $J_{2}$ crosses $1$. These
correspond to points when the $\Omega\to 0$ classical chain goes from having
all up-pointing fields to having arbitrarily long regions where the field
points down, and vice versa. Now let us specialize to
$P(J)=[\delta(J-J_{1})+\delta(J-J_{2})]/2$, $J_{1}\ll 1$ and $J_{2}=1+t_{G}/n$
with $0<t_{G}\ll 1$, to focus on the vicinity of the latter of these Griffiths
singularities, viz. the paramagnetic one. In this regime (from [20, 16])
$f(\Omega)$ has a singular part $e^{-4\alpha\mathcal{J}},$ with the exponent
$\alpha$ defined by $\overline{e^{2\alpha\mathcal{H}}}=1$, for $\Omega$ small
enough that not all $\mathcal{H}$ are the same sign. (Overlines denote
disorder averages.) The latter relation gives $\alpha\approx\ln
2/(t_{G}-\Omega^{n})$, and so we are to integrate $e^{-4\alpha\mathcal{J}}$ up
to $\Omega<t_{G}^{1/n}$. We get that the singular part of $E_{0}/L$ is
proportional to

$\int_{0}^{t_{G}^{1/n}}d\Omega\,\,\Omega^{n\ln 2/(t_{G}-\Omega^{n})}\sim
t_{G}^{c/t_{G}},$ (9)

where $c$ is a constant. There will be an identical singularity as $J_{1}$
approaches $1$ from below.

The other interesting point in the classical chain is the one where the
average field switches direction, $\overline{\mathcal{H}}=0$. This corresponds
to the ferromagnet-paramagnet transition in the quantum chain, which occurs
when $\overline{\ln h}=\overline{\ln J}$. Let us now consider the behaviour in
the vicinity of this point. Shankar and Murthy spliced together results from
Refs. [20, 21] to argue that all the singularities of $f$ in the vicinity of
$\overline{\mathcal{H}}=0$ are captured by the function
$f\sim-\overline{\mathcal{H}}\coth(2\alpha\mathcal{J})$ where $\alpha$ is
defined as above. Again assuming all $h=1$ and defining $t=-(n/2)\overline{\ln
J}$, $\alpha\sim-\overline{\mathcal{H}}\sim t+\Omega^{n}$, where we have
dropped, and will continue to drop, a number of constants. The energy then has
the same singularity as

$\int_{0}^{\epsilon}d\Omega\,(t+\Omega^{n})\coth[(t+\Omega^{n})\ln\Omega^{-n}].$
(10)

with any $\epsilon<1$ (the non-analyticity comes from the integral near zero).
Shankar and Murthy’s observations on an integral of this form with $n=2$ in
fact generalize to all $n$ [16]: to obtain the most important piece of the
singularity we can drop the $\Omega^{n}$ additive pieces, so that the integral
becomes $\int_{0}^{\epsilon}d\Omega\,t\coth(nt\ln\Omega^{-1})$. Up to a
trivial scaling by $n$, this is the same function of $t$ for all $n$. It is
also infinitely differentiable yet non-analytic at $t=0$ [12, 16].

The energy $E_{0}$ therefore contains three infinitely differentiable
singularities as the distributions of $h$ and $J$ are tuned across each other,
of the same form as those in the Ising ($n=2$) case [12, 16]. This should be
contrasted with the clean Hamiltonians, in which one finds a usual critical
exponent $\alpha=1-2/n$ [5].

Boundary magnetization.— While introducing magnetic field ($Z$ or
$Z^{\dagger}$) type terms in the bulk destroys the free-parafermion
solvability of the model, a field on the boundary does not. To be explicit,
the Hamiltonian

$H_{\mathfrak{H}}=H-\mathfrak{H}Z_{L}^{\dagger}\quad(\mathfrak{H}>0)$ (11)

is just as solvable as $H$. This is because $Z_{L}^{\dagger}$ obeys the same
algebra with $X_{L}$ as two adjacent terms in $H$, e.g. $X_{j}$ and
$Z^{\dagger}_{j-1}Z_{j}$. One simply introduces $b_{M+1}=\mathfrak{H}$ and
works with the polynomial

$P_{M+1}(z)=P_{M}(z)-z\mathfrak{H}^{n}P_{M-1}(z).$ (12)

Note that this has the same degree as $P_{M}(z)$—the number of single-particle
levels remains the same. The energy is given by (6) with $M$ replaced by
$M+1$. We define the boundary “magnetization” in the corresponding state in
the presence of the boundary field as
$\mathfrak{M}(\mathfrak{H})=-dE_{0}(\mathfrak{H})/d\mathfrak{H}$ which,
following some trivial manipulations, can be written as

$\mathfrak{M}(\mathfrak{H})=\frac{1}{\pi}\sin\frac{\pi}{n}\int_{0}^{\infty}d\Omega\,\frac{n\mathfrak{H}^{n-1}}{\Omega^{n}r_{M}+\mathfrak{H}^{n}},$
(13)

where $r_{m}=P_{m}(-\Omega^{-n})/P_{m-1}(-\Omega^{-n})$ is a random variable
that depends on all the $J_{j}$ that $P_{m}$ contains. Unless the chain is
clean it fluctuates as a function of $m$—see (7)—even as $m\to\infty$, which
is the limit of interest. However, in this limit it will approach a stationary
distribution [12]. We will be particularly interested in the spontaneous
magnetization, $\lim_{\mathfrak{H}\to 0^{+}}\mathfrak{M}(\mathfrak{H})$.

At this point it is convenient to write our matrix recursion somewhat
differently using the random-field Ising picture discussed in the previous
section. It is easy to verify that if we define the matrix (cf. (8))

$\tilde{T}_{j}(\Omega)=\begin{bmatrix}\Omega^{n}+1&J_{j}^{n}\Omega^{n/2}\\\
\Omega^{n/2}&J_{j}^{n}\end{bmatrix}$ (14)

then

$\left(\prod_{j=1}^{L-1}\tilde{T}_{j}\right)\begin{bmatrix}\Omega^{n}+1\\\
\Omega^{n/2}\end{bmatrix}=\Omega^{nL}\begin{bmatrix}P_{M}(-\Omega^{-n})\\\
\Omega^{-n/2}P_{M-1}(-\Omega^{-n})\end{bmatrix}.$ (15)

Therefore, if we define $x=r_{M}\Omega^{n/2}$, and suppose that in the
$M\to\infty$ limit the law of $x$ is $\nu(x)$, then the disorder-averaged
boundary magnetization is given by

$\overline{\mathfrak{M}(\mathfrak{H})}=\frac{1}{\pi}\sin\frac{\pi}{n}\int_{0}^{\infty}d\Omega\int_{0}^{\infty}dx\,\nu(x)\frac{n\mathfrak{H}^{n-1}}{\Omega^{n/2}x+\mathfrak{H}^{n}}.$
(16)

This rearrangement allows a connection to the work of McCoy and Wu [12], who
studied the distribution $\nu(x)$ for the form of random matrix shown in (14),
with $n=2$. This was accomplished for a power law distribution of $J$’s,
particularly $\mu(J^{n})=P(J)/(nJ^{n-1})=N\lambda_{0}^{-N}J^{n(N-1)}$ for
$0<J^{n}<\lambda_{0}$ and $0$ otherwise. The distribution $\mu(\lambda)$ is
assumed to be narrow around $\lambda_{0}$, i.e. $N\gg 1$. The tuning parameter
for the transition is therefore $\lambda_{0}$; the phase transition in
magnetization takes place when $t=(1/N-\ln\lambda_{0})/2=0$. We work with the
same distribution. McCoy then used this $\nu(x)$ to study boundary
magnetization [14], manipulating Pfaffians for a 2d Ising lattice with an
extra column to get to a formula similar to (16).

McCoy and Wu’s results show, mutatis mutandis, that deviations from a clean
model occur only for $|t|,\Omega^{n/2}\lesssim N^{-2}$, motivating a change of
variables $\delta\sim N^{2}t,\phi\sim N^{2}\Omega^{n/2},$ and correspondingly
$\mathfrak{h}=N^{2/n}\mathfrak{H}$. In the $\delta,\phi\sim 1$ regime, another
change of variables $x=\lambda_{0}^{1/2}\phi e^{-q}/2$ is useful. The
distribution of $q$ was determined by McCoy and Wu to be
$U(q)\equiv(\phi/2)^{\delta}\exp\\{-\delta
q-e^{q}-\phi^{2}e^{-q}/4\\}/(2K_{\delta}(\phi))$ where $K_{\delta}$ is a
modified Bessel function of the second kind [12]. Throughout this discussion
we are cavalier about factors that do not matter for the form of the
singularity. We have that $N^{2/n}\overline{\mathfrak{M}(\mathfrak{H})}$ goes
like

$\int_{0}^{\infty}d\phi\int_{-\infty}^{\infty}dq\,\phi^{2/n-1}\left(\frac{\phi}{2}\right)^{\delta}\frac{e^{-\delta
q-e^{q}-\phi^{2}e^{-q}/4}\mathfrak{h}^{n-1}}{K_{\delta}(\phi)(\phi^{2}e^{-q}+\mathfrak{h}^{n})}.$
(17)

In the $\mathfrak{h}\to 0$ limit the $\phi^{2}$ term in the exponent can be
dropped. Making a further change of variables $\phi=\mathfrak{h}^{n/2}\alpha$,
and plugging in the small-$\delta$ asymptotics of $K_{\delta}$, we have
(dropping factors of order 1):

$\int_{0}^{\infty}d\alpha\int_{-\infty}^{\infty}dq\frac{(\mathfrak{h}^{n/2}\alpha/2)^{\delta+|\delta|}}{\Gamma(|\delta|)}\frac{\alpha^{2/n-1}}{\alpha^{2}e^{-q}+1}e^{-\delta
q-e^{q}}$ (18)

so that all $\mathfrak{h}$ dependence is now in the factor
$(\mathfrak{h}^{n/2}\alpha/2)^{\delta+|\delta|}$. When $\delta>0$ this factor
vanishes as $\mathfrak{h}\to 0$, so there is no magnetization in the
disordered phase. For $\delta<0$, the factor is $1$, and the rest of the
integral can be easily performed to yield constant factors times
$\Gamma(1/n+|\delta|)/\Gamma(|\delta|)\sim\Gamma(1/n)|\delta|$. Therefore

$\lim_{\mathfrak{H}\to
0^{+}}\overline{\mathfrak{M}(\mathfrak{H})}\sim\Theta(-t)|t|^{\beta_{\text{boundary}}};\quad\beta_{\text{boundary}}=1.$
(19)

We see that the boundary magnetization exponent is also independent of $n$.

Strong-disorder RG.— While the focus of this paper has been on exact results,
it is worth adding that our formalism allows us to make contact rather
naturally with Fisher’s strong-disorder RG [19]. Our comments in fact
generalize those of Fisher’s Appendix E, in which he explains the connection
(for the TFIC) between his results and Shankar and Murthy’s. Here one
introduces yet another matrix factorization of the polynomial
$P_{M}(\Omega^{-n})$. With matrices

$\mathcal{T}^{(h)}(h)=\begin{bmatrix}\Omega^{n/2}/4-h^{n}\Omega^{-n/2}&-1/2\\\
1/2&-\Omega^{-n/2}\end{bmatrix},$ (20)

and

$\mathcal{T}^{(J)}(J)=\begin{bmatrix}\Omega^{-n/2}&1/2\\\
-1/2&J^{n}\Omega^{-n/2}-\Omega^{n/2}/4\end{bmatrix},$ (21)

one can observe that

$P_{2L-1}(\Omega^{-n})=\begin{bmatrix}\Omega^{-n/2}&1/2\end{bmatrix}\mathcal{T}^{(h)}(h_{L})\mathcal{T}^{(J)}(J_{L-1})\dots{\\\
}\dots\mathcal{T}^{(J)}(J_{1})\mathcal{T}^{(h)}(h_{1})\begin{bmatrix}1\\\
-\Omega^{n/2}/2\end{bmatrix}.$ (22)

The operation of the RG comprises repeated multiplications of these matrices,
three at a time. The decomposition above has the nice property that,
approximately,

$\mathcal{T}^{(h)}(h_{1})\mathcal{T}^{(J)}(J)\mathcal{T}^{(h)}(h_{2})\propto\mathcal{T}^{(h)}\left(\frac{h_{1}h_{2}}{J}\right),$
(23)

with corrections that go to zero as $J$ becomes big and $\Omega$ small. An
identical relation holds with $h$’s and $J$’s interchanged. This means that if
$J$ is a bond with transverse fields $h_{1,2}\ll J$ next to it, we can replace
the two spins and the bond with a new spin with transverse field
$h_{1}h_{2}/J$, so long as we are interested in small single-particle energies
compared to these scales (or equivalently the low-frequency contribution to
$E_{0}$, as in (6)). This is exactly the result one would derive in a more
pedestrian approach via second-order perturbation theory.

This implies immediately that strong-disorder RG at criticality leads to an
identical infinite-randomness fixed point to the TFIC with, e.g., the same
tunnelling exponent $\psi=1/2$ [19]. (This fact could have been deduced by
simply inspecting $\varepsilon_{k}=z_{k}^{-1/n}$.) Moreover the RG
trajectories are identical decimation-by-decimation for a family of models
labelled by $n$. This should be exploited to reproduce the results of this
paper, but also to calculate distributions of correlation functions [19, 22],
which we do not see a way to attack exactly. Strong-disorder RG often renders
infinite-randomness critical points more tractable than their clean
counterparts; this is certainly such an example, since little is known about
correlation functions in the clean chain. Also, the RG rule (23) works just as
well for complex couplings $\\{J_{j},h_{j}\\}$, with the correctness of the RG
now controlled by broad distributions of their magnitudes.

Disguised free fermions.— There is another class of one-dimensional
Hamiltonians introduced by Fendley [23] which are free fermionic in the sense
of having a spectrum of the form (3) with $n=2$, despite not being quadratic
in Jordan-Wigner fermions. See also generalizations in Refs. [24, 25, 26, 27].
These models are solved using polynomials which obey recursion relations
similar to (4). A disordered version of Fendley’s model was studied by Alcaraz
et al. [28], numerically and using strong-disorder RG (formulated as
perturbation theory). A number of infinite-randomness critical points were
found, all of the Ising type. The matrices which generate these recursions are
larger than $2\times 2$, so the theory is naively more complicated than that
used here, which relied on the consideration of a single random variable of
the form $P_{M}/P_{M-1}$ [12, 20, 21]. Nevertheless the random couplings
appear in these other polynomials in similar combinations [29], which ought to
make some progress possible. It would also be nice to reverse-engineer free
models whose single-particle energies are related to other interesting classes
of products of random matrices.

Discussion.— Let us make two more general comments on subjects which greatly
interest us.

It is remarkable that the critical behaviour of the disordered free
parafermionic models is independent of $n$. This strange superuniversal
behaviour—identical critical exponents for a number of different theories with
clearly different symmetries and therefore naively different “operator
content”—seems to exist for other families of infinite-randomness fixed points
in one spatial dimension [30], as well as in higher dimensions [31]. (Other
universality classes are known, however [32, 33, 34, 35, 36, 37].) We hope
that access to exact results that demonstrate superuniversality of infinite-
randomness fixed points will be useful for progress in understanding this
phenomenon.

A particularly transparent solution of the 2d square-lattice Ising model with
fully inhomogeneous couplings is via a mapping to a dimer problem which is
solvable via Pfaffians [38, 17]. The relation between this solution and free
fermions can be made explicit [39]. What the analogues of the 2d Ising
model/dimers/Pfaffians are for parafermions is an open question. We believe
our Eq. (6) will be of value in answering this. The link to the so-called
$\tau_{2}$ model [2, 40, 41, 42], effectively associated with a particular
frequency of the quantum chain, is tantalizing but not quite what we want. It
is at least worth fully understanding the relations between these various
models for $n=2$.

Acknowledgments.— We are grateful to V. Calvera, P. Fendley, J.A. Hoyos, D.A.
Huse, C. Murthy and N. O’Dea for discussions. A.P. and A.C. were supported by
Stanford Graduate Fellowships.

## References

  * Baxter [1989a] R. J. Baxter, A simple solvable ZN Hamiltonian, Physics Letters A 140, 155 (1989a).
  * Baxter [1989b] R. J. Baxter, Superintegrable chiral Potts model: Thermodynamic properties, an “inverse” model, and a simple associated Hamiltonian, Journal of Statistical Physics 57, 1 (1989b).
  * Fendley [2014] P. Fendley, Free parafermions, Journal of Physics A Mathematical General 47, 075001 (2014), arXiv:1310.6049 [cond-mat.stat-mech] .
  * Fradkin and Kadanoff [1980] E. Fradkin and L. P. Kadanoff, Disorder variables and para-fermions in two-dimensional statistical mechanics, Nuclear Physics B 170, 1 (1980).
  * Alcaraz  _et al._ [2017] F. C. Alcaraz, M. T. Batchelor, and Z.-Z. Liu, Energy spectrum and critical exponents of the free parafermion Z N spin chain, Journal of Physics A Mathematical General 50, 16LT03 (2017), arXiv:1612.02617 [cond-mat.stat-mech] .
  * Alcaraz and Batchelor [2018] F. C. Alcaraz and M. T. Batchelor, Anomalous bulk behavior in the free parafermion Z (N ) spin chain, Phys. Rev. E 97, 062118 (2018), arXiv:1802.04453 [cond-mat.stat-mech] .
  * Liu  _et al._ [2019] Z.-Z. Liu, R. A. Henry, M. T. Batchelor, and H.-Q. Zhou, Some ground-state expectation values for the free parafermion Z(N) spin chain, Journal of Statistical Mechanics: Theory and Experiment 12, 124002 (2019), arXiv:1905.08384 [cond-mat.stat-mech] .
  * Alcaraz and Pimenta [2020a] F. C. Alcaraz and R. A. Pimenta, Free fermionic and parafermionic quantum spin chains with multispin interactions, Phys. Rev. B 102, 121101 (2020a), arXiv:2005.14622 [cond-mat.stat-mech] .
  * Alcaraz and Pimenta [2020b] F. C. Alcaraz and R. A. Pimenta, Integrable quantum spin chains with free fermionic and parafermionic spectrum, Phys. Rev. B 102, 235170 (2020b).
  * Alcaraz  _et al._ [2021] F. C. Alcaraz, J. A. Hoyos, and R. A. Pimenta, Powerful method to evaluate the mass gaps of free-particle quantum critical systems, Phys. Rev. B 104, 174206 (2021), arXiv:2109.01938 [cond-mat.stat-mech] .
  * Henry and Batchelor [2023] R. Henry and M. Batchelor, Exceptional points in the Baxter-Fendley free parafermion model, SciPost Physics 15, 016 (2023), arXiv:2301.11031 [cond-mat.stat-mech] .
  * McCoy and Wu [1968] B. M. McCoy and T. T. Wu, Theory of a Two-Dimensional Ising Model with Random Impurities. I. Thermodynamics, Physical Review 176, 631 (1968).
  * McCoy and Wu [1969] B. M. McCoy and T. T. Wu, Theory of a Two-Dimensional Ising Model with Random Impurities. II. Spin Correlation Functions, Physical Review 188, 982 (1969).
  * McCoy [1969] B. M. McCoy, Theory of a Two-Dimensional Ising Model with Random Impurities. III. Boundary Effects, Physical Review 188, 1014 (1969).
  * McCoy [1970] B. M. McCoy, Theory of a Two-Dimensional Ising Model with Random Impurities. IV. Generalizations, Phys. Rev. B 2, 2795 (1970).
  * Shankar and Murthy [1987] R. Shankar and G. Murthy, Nearest-neighbor frustrated random-bond model in d=2: Some exact results, Phys. Rev. B 36, 536 (1987).
  * McCoy and Wu [1973] B. M. McCoy and T. T. Wu,  _The two-dimensional Ising model_ (Harvard University Press, 1973).
  * Schultz  _et al._ [1964] T. D. Schultz, D. C. Mattis, and E. H. Lieb, Two-Dimensional Ising Model as a Soluble Problem of Many Fermions, Reviews of Modern Physics 36, 856 (1964).
  * Fisher [1995] D. S. Fisher, Critical behavior of random transverse-field Ising spin chains, Phys. Rev. B 51, 6411 (1995).
  * Derrida and Hilhorst [1983] B. Derrida and H. J. Hilhorst, Singular behaviour of certain infinite products of random 2 × 2 matrices, Journal of Physics A Mathematical General 16, 2641 (1983).
  * Nieuwenhuizen and Luck [1986] T. M. Nieuwenhuizen and J. M. Luck, Exactly soluble random field Ising models in one dimension, Journal of Physics A Mathematical General 19, 1207 (1986).
  * Fisher and Young [1998] D. S. Fisher and A. P. Young, Distributions of gaps and end-to-end correlations in random transverse-field Ising spin chains, Phys. Rev. B 58, 9131 (1998), arXiv:cond-mat/9802246 [cond-mat] .
  * Fendley [2019] P. Fendley, Free fermions in disguise, Journal of Physics A Mathematical General 52, 335002 (2019), arXiv:1901.08078 [cond-mat.stat-mech] .
  * Elman  _et al._ [2021] S. J. Elman, A. Chapman, and S. T. Flammia, Free Fermions Behind the Disguise, Communications in Mathematical Physics 388, 969 (2021), arXiv:2012.07857 [quant-ph] .
  * Chapman  _et al._ [2023] A. Chapman, S. J. Elman, and R. L. Mann, A Unified Graph-Theoretic Framework for Free-Fermion Solvability, arXiv e-prints , arXiv:2305.15625 (2023), arXiv:2305.15625 [quant-ph] .
  * Fendley and Pozsgay [2023] P. Fendley and B. Pozsgay, Free fermions with no Jordan-Wigner transformation, arXiv e-prints , arXiv:2310.19897 (2023), arXiv:2310.19897 [cond-mat.stat-mech] .
  * Pozsgay [2024] B. Pozsgay, Quantum circuits with free fermions in disguise, arXiv e-prints , arXiv:2402.02984 (2024), arXiv:2402.02984 [quant-ph] .
  * Alcaraz  _et al._ [2023] F. C. Alcaraz, J. A. Hoyos, and R. A. Pimenta, Random free-fermion quantum spin chain with multispin interactions, Phys. Rev. B 108, 214413 (2023), arXiv:2308.16249 [cond-mat.dis-nn] .
  * de Calan  _et al._ [1985] C. de Calan, J. M. Luck, T. M. Nieuwenhuizen, and D. Petritis, On the distribution of a random variable occurring in 1D disordered systems, Journal of Physics A Mathematical General 18, 501 (1985).
  * Senthil and Majumdar [1996] T. Senthil and S. N. Majumdar, Critical Properties of Random Quantum Potts and Clock Models, Phys. Rev. Lett. 76, 3001 (1996), arXiv:cond-mat/9511121 [cond-mat] .
  * Kang  _et al._ [2020] B. Kang, S. A. Parameswaran, A. C. Potter, R. Vasseur, and S. Gazit, Superuniversality from disorder at two-dimensional topological phase transitions, Phys. Rev. B 102, 224204 (2020), arXiv:2008.09617 [cond-mat.str-el] .
  * Hyman and Yang [1997] R. A. Hyman and K. Yang, Impurity Driven Phase Transition in the Antiferromagnetic Spin-1 Chain, Phys. Rev. Lett. 78, 1783 (1997), arXiv:cond-mat/9608039 [cond-mat] .
  * Monthus  _et al._ [1997] C. Monthus, O. Golinelli, and T. Jolicœur, Percolation Transition in the Random Antiferromagnetic Spin-1 Chain, Phys. Rev. Lett. 79, 3254 (1997), arXiv:cond-mat/9705231 [cond-mat.str-el] .
  * Refael  _et al._ [2002] G. Refael, S. Kehrein, and D. S. Fisher, Spin reduction transition in spin-3/2 random Heisenberg chains, Phys. Rev. B 66, 060402 (2002), arXiv:cond-mat/0111295 [cond-mat.dis-nn] .
  * Damle and Huse [2002] K. Damle and D. A. Huse, Permutation-Symmetric Multicritical Points in Random Antiferromagnetic Spin Chains, Phys. Rev. Lett. 89, 277203 (2002), arXiv:cond-mat/0207244 [cond-mat.str-el] .
  * Hoyos and Miranda [2004] J. A. Hoyos and E. Miranda, Random antiferromagnetic SU (N) spin chains, Phys. Rev. B 70, 180401 (2004), arXiv:cond-mat/0406130 [cond-mat.str-el] .
  * Quito  _et al._ [2016] V. L. Quito, J. A. Hoyos, and E. Miranda, Random SU(2)-symmetric spin-S chains, Phys. Rev. B 94, 064405 (2016), arXiv:1512.04542 [cond-mat.dis-nn] .
  * Kasteleyn [1963] P. W. Kasteleyn, Dimer Statistics and Phase Transitions, Journal of Mathematical Physics 4, 287 (1963).
  * Samuel [1980] S. Samuel, The use of anticommuting variable integrals in statistical mechanics. I. The computation of partition functions, Journal of Mathematical Physics 21, 2806 (1980).
  * Baxter [2014] R. J. Baxter, The $\tau$2 model and parafermions, Journal of Physics A Mathematical General 47, 315001 (2014), arXiv:1310.7074 [cond-mat.stat-mech] .
  * Au-Yang and Perk [2014] H. Au-Yang and J. H. H. Perk, Parafermions in the $\tau$2 model, Journal of Physics A Mathematical General 47, 315002 (2014), arXiv:1402.0061 [math-ph] .
  * Au-Yang and Perk [2016] H. Au-Yang and J. H. H. Perk, Parafermions in the tau-2 model II, arXiv e-prints , arXiv:1606.06319 (2016), arXiv:1606.06319 [math-ph] .