1 Introduction

The concept of free fermions is fundamental to the celebrated exact solution of the two-dimensional Ising model in zero magnetic field [1,2,3] and its one-dimensional quantum counterpart [4]. Much later, it was revealed [5] that in solving the Ising model Onsager initially diagonalised the associated row transfer matrix by hand: first for strip width \(L=2\), then \(L=3\), and so on. Eventually, by the \(L=6\) case, he confirmed that the \(2^6=64\) transfer matrix eigenvalues were all of the form \(\exp (\pm \gamma _1 \pm \cdots \pm \gamma _6)\). This observation suggested an underlying product algebra which, in turn, led Onsager to the mathematical structure underlying his original derivation. Insights obtained from the Ising model – specifically the exact results – played a central role in the development of the theory of phase transitions and critical phenomena. The concept of free fermions also underpins the exact solution of many other fundamental models, with the XY chain playing an early key role [6].

Fig. 1
figure 1

A free parafermion spectrum in the complex plane, for \(N=3\), \(L=4\). The black dots are the \(N^L\) energy eigenstates. The spectrum is built up by starting at zero and adding each parafermion multiplied by a power of the root of unity \(\omega = \exp (2 \pi \textrm{i} / N)\), as per Eq. (3). Here the values of \(\epsilon _j\) are chosen arbitrarily and for most realistic values the paths will overlap each other, but will have the same essential branching structure. Algebraic expressions are shown for some example states, including the ground state \(-\epsilon _1-\epsilon _2-\epsilon _3-\epsilon _4\)


Fig. 2
figure 2

Free parafermion quasi-energy levels for \(N=2\), 3, 4, 5. Each diagram represents a single energy eigenstate of the Hamiltonian. For each \(j\in \{1,\ldots ,L\}\), one power of \(\omega = \exp (2 \pi \textrm{i} / N)\) is chosen. In other words, each quasi-energy is subject to a Fermi “exclusion circle”. This gives \(N^L\) possible choices, determining all the energy eigenvalues


In 1989, Baxter [7, 8] derived a simple spin chain Hamiltonian from a 2D classical model known as the \(\tau _2\) model, a model which was also essential to the solution of the integrable chiral Potts model [9, 10]. As we shall see below, the Hamiltonian of this spin chain takes a similar form to that of the quantum Ising chain, but the spins have N allowed states, with Baxter’s Hamiltonian reducing to the Ising model for \(N=2\). Baxter found that the energy spectrum of this more general model decomposes into a sum of independent terms involving powers of \(\omega\), an Nth root of unity [7].

In 2012, Fendley [11] studied the edge modes of a chiral Hermitian Z(N) spin chain by rewriting the Hamiltonian in terms of parafermions. Soon after, he realised Baxter’s non-Hermitian model has algebraic properties which generalise those of the fermion representation of the Ising model, and other free fermionic models. This lead to a rigorous derivation of the spectrum, charges and the corresponding generalised Clifford algebra in 2014 [12]. This was the pivotal observation of the previously elusive free parafermions.

In this article we briefly outline related developments for free parafermions, along with further recent progress.

2 The free parafermion model

Baxter’s free parafermion model is defined by the simple Hamiltonian

\(\begin{aligned} H=-\sum \limits _{j=1}^{L-1} Z_j^\dagger Z_{j+1} - \lambda \sum \limits _{j=1}^L X_j, \end{aligned}\)
(1)

with N states per site. X and Z are N-state generalisations of the Pauli matrices given by the \(N\times N\) matrices

\(\begin{aligned} Z = \left( \begin{array}{ccccc} 1 &{} 0 &{} 0 &{} \ldots &{} 0 \\ 0 &{} \omega &{} 0 &{} \ldots &{} 0 \\ 0 &{} 0 &{} \omega ^2 &{} \ldots &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \ldots &{} \omega ^{N-1} \\ \end{array}\right) ,\quad X = \left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} \ldots &{} 0 &{} 1 \\ 1 &{} 0 &{} 0 &{} \ldots &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} \ldots &{} 0 &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \ldots &{} 1 &{} 0 \\ \end{array}\right) , \end{aligned}\)
(2)

with \(\omega = \exp (2\pi \textrm{i}/N)\) a root of unity. The Hamiltonian is non-Hermitian and has complex eigenvalues. Baxter found that these eigenvalues have a simple form that generalises free fermions:

\(\begin{aligned} E=\sum \limits _{k=1}^L \omega ^{s_k} \epsilon _k . \end{aligned}\)
(3)

Here \(\epsilon _k\) are quasienergies which depend on \(\lambda\), and the set of integers \(s_k \in \{0,\ldots ,N-1\}\) label each state. The N choices for each \(s_k\) generalise the two choices (\(+1\) and \(-1\), corresponding to particles and holes) in a free fermion model. This gives the model its name of free parafermions, although this was only coined later by Fendley [12]. Figure 1 shows an example of a free parafermionic energy spectrum in the complex plane, subject to the generalised exclusion rule (Fig. 2).

2.1 \(\mathcal{P}\mathcal{T}\) symmetry

Although the Hamiltonian is not Hermitian, it is parity-time (\(\mathcal{P}\mathcal{T}\))-symmetric, with the action of the operators \(\mathcal {P}\) and \(\mathcal {T}\) such that

\(\begin{aligned} \mathcal {P} Z_j \mathcal {P} = Z_{L+1-j} \,,{} & {} \;\mathcal {P} X_j \mathcal {P} = X_{L+1-j},\end{aligned}\)
(4)

\(\begin{aligned} \quad \mathcal {T} Z_j \mathcal {T} = Z_j^\dagger \,,{} & {} \;\mathcal {T} X_j\mathcal {T} = X_j. \end{aligned}\)
(5)

Bender and Boettcher showed that \(\mathcal{P}\mathcal{T}\)-symmetric models may have real spectra [13] if the symmetry is unbroken. They can be given unitary time evolution with the appropriate choice of metric, even with a broken symmetry and complex eigenvalues [14]. This has lead to a great deal of activity including many experiments with \(\mathcal{P}\mathcal{T}\)-symmetric systems, many of which are detailed in the extensive review by Ashida et al. [15]. In the case of the free parafermion model, \(\mathcal{P}\mathcal{T}\) symmetry is always broken, and the spectrum appears in complex conjugate pairs, with the conjugation corresponding to the action of the symmetry operator.

For a finite diagonalisable system, the metric

\(\begin{aligned} G = \sum \limits _i |L_i \rangle \langle L_i |\end{aligned}\)
(6)

is sufficient, where \(|L_i\rangle\) is the ith left eigenstate. An inner product between two states \(|\Psi _1 \rangle\) and \(|\Psi _2 \rangle\) is then evaluated as \(\langle \Psi _1 |G |\Psi _2 \rangle\). The expectation values of an eigenstate are easily expressed with this metric: for a right eigenstate \(|R\rangle\) and corresponding left eigenstate \(|L\rangle\) and operator A, the expectation value is

\(\begin{aligned} \langle \langle A \rangle \rangle = \langle L |A |R \rangle . \end{aligned}\)
(7)

Because the model is non-Hermitian, it is possible that the Hamiltonian is non-diagonalisable for some values of its parameters, known as exceptional points. These points have interesting physical properties and are discussed in detail in Ashida’s review [15]. The behaviour of the metric as the system passes through an exceptional point is the subject of current research [16, 17]. At least in the case of the free parafermion model, the exceptional points are isolated, and Eq. (6) is usable for most calculations. Recent work on the model’s exceptional points is discussed in Section 5.

3 From fermions to parafermions

3.1 The quantum Ising chain and fermions

The widely studied Hamiltonian of the quantum Ising model on a chain of length L is defined in terms of Pauli matrices by

\(\begin{aligned} H_\mathrm{Ising}= -\sum \limits _{j=1}^{L-1} \sigma ^z_j \sigma ^z_{j+1} - \lambda \sum \limits _{j=1}^L \sigma ^x_j, \end{aligned}\)
(8)

where the subscripts indicate on which lattice site an operator is acting. The spin-spin interaction strength is scaled to unity for simplicity. For convenience we have defined the model with open boundary conditions in contrast to period boundary conditions for which \(\sigma _{L+1}^z=\sigma _1^z\). Hamiltonian (8) possesses a Z(2) symmetry, \(\left[H_\mathrm{Ising},(-1)^F\right]=0\), with

\(\begin{aligned} (-1)^F = \prod \limits _{j=1}^L \sigma ^x_j. \end{aligned}\)
(9)

The spin operators can be rewritten in terms of fermionic operators using the Jordan-Wigner transformation. To be specific, we define the fermionic operators

\(\begin{aligned} \phi _{2j-1}=\prod \limits _{k=1}^{j-1}\sigma ^x_k \sigma ^z_j, \quad \phi _{2j}=\prod \limits _{k=1}^{j-1}\sigma ^x_k \sigma ^y_j. \end{aligned}\)
(10)

As fermionic operators, these operators satisfy the Clifford algebra and anticommute with each other. Moreover, they are Hermitian and square to the identity, so they are actually Majorana fermionic operators.

One can check that bilinears in neighbouring fermionic operators results in \(\sigma ^z_j \sigma ^z_{j+1}\) and \(\sigma ^x_j\) for some j. Then it is possible to rewrite the Hamiltonian (8) as

\(\begin{aligned} H=\textrm{i}\sum \limits _{j=1}^{2L-1} t_j \phi _j \phi _{j+1}, \end{aligned}\)
(11)

with \(t_{2j-1}=\lambda\) and \(t_{2j}=1\). The key here is that Hamiltonians expressed as bilinears in fermionic operators can be diagonalised exactly by choosing an appropriate basis, of which an early example is shown in Appendix A of Ref. [6]. The total energy is then the sum of individual fermions without any interactions in this new basis. These fermions are free and therefore such Hamiltonians are free fermionic.

3.2 The Z(N) spin chain and parafermions

As mentioned in Section 2, the free parafermion model is defined using the operators X and Z which generalise the Pauli matrices. A generalisation of the Y operator may also be defined as

\(\begin{aligned} Y_{i,j} =\omega ^{(N-1)/2}\left(\delta _{N,1}+\omega ^{N-i}\delta _{i+1,j}\right)=\omega ^{(N-1)/2} \left( \begin{array}{ccccc} 0 &{} \omega ^{N-1} &{} 0 &{} \dots &{} 0 \\ 0 &{} 0 &{} \omega ^{N-2} &{} \dots &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \dots &{} \omega \\ 1 &{} 0 &{} 0 &{} \dots &{} 0 \end{array}\right) _{i,j}, \end{aligned}\)
(12)

These operators satisfy the commutation relations

\(\begin{aligned} XY=\omega YX,\quad YZ=\omega ZY, \quad ZX=\omega XZ, \end{aligned}\)
(13)

and have the property that their N-th power is equal to the identity operator:

\(\begin{aligned} A^N=\mathbb {I}, \quad A=X,Y,Z. \end{aligned}\)
(14)

This implies \(A^{-1}=A^\dagger =A^{N-1}\) for \(A=X,Y,Z\). Also, products of two different generalised Pauli matrices in the correct order gives \(\omega ^{(N-1)/2}\) times the conjugate transpose of the remaining one:

\(\begin{aligned} XY= \omega ^{(N-1)/2}Z^\dagger ,\quad YZ=\omega ^{(N-1)/2} X^\dagger ,\quad ZX=\omega ^{(N-1)/2}Y^\dagger . \end{aligned}\)
(15)

Such matrices have a long history and were defined by Sylvester in 1882 as generalisations of quaternions [18]. Like quaternions, they satisfy a (generalised) Clifford algebra defined by the above commutation relations. This leads to a set of parafermion generators, or simply “parafermions", which have appeared in various forms and were studied by Yamazaki [19] and Morris [20] in the context of the generalised Clifford algebra. For other early related references to parafermions the reader is referred to Jaffe and Pedrocchi [21], who investigated the relation with reflection positivity, and also Fradkin and Kadanoff [22] who introduced the Z(N) generalisation of the Jordan-Wigner transformation.

The free parafermion model (1) can be obtained from the Ising chain by directly substituting the generalised operators. It is straightforward to check that (1) is invariant under the Z(N) symmetry operator

\(\begin{aligned} \omega ^\mathcal {P}\equiv \prod _{j=1}^L X_j, \end{aligned}\)
(16)

with \((\omega ^\mathcal {P})^N=\mathbb {I}\) following from (14). Following Fendley’s work [12], we can rewrite this Hamiltonian (1) in terms of parafermionic operators \(\psi _j\) acting on a lattice of length 2L:

\(\begin{aligned} \psi _{2j-1}=\prod \limits _{k=1}^{j-1}X_k Z_j, \quad \psi _{2j}=\omega ^{-1}\prod \limits _{k=1}^{j-1}X_k Y_j^\dagger . \end{aligned}\)
(17)

These operators satisfy \(\psi _j^N = 1\), and the \(\omega\)-commutation relation

\(\begin{aligned} \psi _a \psi _b=\omega \psi _b \psi _a, \quad a
(18)

This can be verified with (14) and (13). In fact, Fendley’s analysis extrends to any nonuniform couplings \(t_j\). These operators are parafermions and generalise Majorana fermions to the Z(N) case. It is not hard to check that the Hamiltonian (1) is a sum of bilinears in parafermionic operators:

\(\begin{aligned} H=-\omega ^{-\frac{N-1}{2}}\sum \limits _{j=1}^{2L-1} t_j \psi _j^\dagger \psi _{j+1}, \end{aligned}\)
(19)

with \(t_{2j-1}=\lambda\) and \(t_{2j}=1\). Starting with this form, Fendley [12] performs a series of linear transformations which determine the free spectrum and eigenvectors. Although H is expressed as a sum of bilinears of parafermions, this does not guarantee that it has a free spectrum, unlike the fermionic case [23]. A general condition for free spectra is given by the exchange algebra discussed in Section 4.

3.3 Spectrum and other physical quantities

Baxter’s method [7, 8] to calculate the spectrum of (1) can be summarised as follows. Analogously to the Ising case, we interlace the coupling constants 1 and \(\lambda\) to define the \(2L\times 2L\) symmetric matrix B:

\(\begin{aligned} B= \left( \begin{array}{cccccc} 0 &{} \lambda ^{N/2} &{} \\ \lambda ^{N/2} &{} 0 &{} 1 \\ &{} 1 &{} 0 &{} \ \lambda ^{N/2} \\ &{} &{} \ddots &{} \ddots &{} \ddots \\ &{} &{} &{} 1 &{} 0 &{} \lambda ^{N/2} \\ &{} &{} &{} &{} \lambda ^{N/2} &{} 0 \end{array}\right) , \end{aligned}\)
(20)

where the powers \(\lambda ^2\) found in the Ising case are replaced by \(\lambda ^{N/2}\). The off-diagonal entries are the coefficients \(t_j^{N/2}\) in (19), so the matrix B represents the Hamiltonian in the parafermion basis. The next step is calculating the quasienergies \(\epsilon _k\) as the solutions of the equation

\(\begin{aligned} \det {\epsilon ^{N/2} \mathbb {I}-B}=0. \end{aligned}\)
(21)

This equation is of degree L in \(\epsilon ^{N}\). From the Z(N) symmetry, we know the eigenvalues are N-fold degenerate, thus only L independent \(\epsilon _k\)’s are left. This implies the form of the free spectrum, as given in (3). Baxter’s argument is based on analogies to the \(N=2\) Ising case and does not give any information beyond the eigenvalues. Fendley’s derivation of the free spectrum [12] also determines the eigenvectors in terms of the parafermion operators, as well as other aspects such as conserved charges and higher Hamiltonians.

4 The exchange algebra

Fendley’s solution [12] makes use of the fact that the terms in the Hamiltonian obey a simple exchange algebra. He uses a similar algebra in related work [24] on free fermions. These are both forms of a more general algebra identified by Alcaraz and Pimenta [25], which is satisfied by a large class of free parafermionic models. This algebra is defined for a general Hamiltonian of the form

\(\begin{aligned} H = \sum \limits _{i=1}^M h_i \,. \end{aligned}\)
(22)

The M generators \(h_i\) satisfy the algebra

\(\begin{aligned} h_i h_{i+m}{} & {} = \omega \, h_{i+m}h_i \quad \textrm{for}\quad 1\le m\le p, \nonumber \\ \left[ h_i, h_j\right]{} & {} = 0 \qquad \qquad ~ \, {\textrm{for}}\quad |i-j |> p, \end{aligned}\)
(23)

and the closure relation \(h_i^N = \lambda _i^N\), where \(\omega\) and \(\lambda _i\) are complex numbers. The parameter p is a positive integer determining the range of the “interaction” between terms within which they \(\omega\)-commute, i.e.,  \(p > 1\) is a multispin model. Alcaraz and Pimenta show that any Hamiltonian obeying this algebra is integrable and has a free parafermion spectrum of the form (3), for some values of the parafermion quasienergies \(\epsilon _j\). The Baxter Hamiltonian (1) satisfies the algebra with \(\omega = \exp (2 \pi \textrm{i}/N)\) and \(M=2L-1\) the total number of terms in the Hamiltonian.

This algebra has been used to explore generalised free parafermion models, including a class of multispin XY-type models with spectra composed of combinations of free parafermions [26, 27]. It has also been used to develop an efficient numerical method for calculating the mass gaps associated with these free particle modes [28]. Most recently, it has been shown [29] how to build standard quantum Ising chains with inhomogeneous couplings which have the same spectra as the new family of free fermionic quantum spin chains with multispin interactions.

Minami [30] has provided a general list of models composed of the X, Y and Z operators which satisfy generalised Onsager algebras.

4.1 Polynomial expressions

Alcaraz and Pimenta [25] show that the quasienergies \(\epsilon _j\) for a general Hamiltonian satisfying the algebra (23) are given by the roots \(z_i\) of a polynomial \(P_M^{(p)}(z)\), with \(\epsilon _i = z_i^{-1/N}\). This provides a practical way to determine the quasienergies for a general model satisfying the algebra, although it is less efficient than diagonalising (20). There are \(\bar{M} = \lfloor \frac{M+p}{p+1}\rfloor\) quasienergies, and the polynomial is determined by the recursion relation

\(\begin{aligned} P_M^{(p)}(z) = \sum \limits _{l=0}^{\bar{M}} C_M(l) z^l, \end{aligned}\)
(24)

which satisfies a recurrence relation for \(M\ge\) 1

\(\begin{aligned} P_M^{(p)}(z) = P^{(p)}_{M-1}(z) - z\lambda _M^N P_{M-(p+1)}^{(p)}(z), \end{aligned}\)
(25)

and an initial condition \(P^{(p)}_M(z) = 1\) for \(M \le 0\).

5 Other developments

Before turning to other developments, we remark that Hamiltonians obeying relations similar to the exchange algebra with \(\omega =-1\) have been shown [31] to have explicit free fermionic spectra. More specifically, frustration graphs are drawn for Hamiltonians representing their local Hamiltonian commutation relations. If the resulting graph of a Hamiltonian is free of “even-hole” or “claw” structures, then a free fermion spectrum can be constructed. It would be worthwhile attempting this approach with Z(N) clock models.

Much is known, including zero mode criticality and low energy CFT [32], about the free parafermion model when Hermitian conjugate (h.c.) terms are added to the Hamiltonian. For example, including the h.c. terms for \(N=3\) results in the well known 3-state Potts model. So far various critical properties have been calculated for the free parafermion Hamiltonian (1). These include the bulk ground state energy and critical exponents [33] and certain correlations [34]. Insights from the free parafermionic structure of the Z(N) spin chain were also used to study the \(\tau _2\) model [35,36,37].

It should be stressed that the free parafermion model has only been solved for open boundary conditions. An unexpected numerical observation is that the ground state energy depends on the boundary conditions in the bulk (\(L\rightarrow \infty\)) limit [38]. This is a characteristic property of non-Hermitian systems, likely being an example of a non-Hermitian skin effect [39], which can only occur in a non-Hermitian system and is caused by macroscopic occupation of boundary states.

By extending \(\lambda\) to complex values, it has most recently been shown that a series of exceptional points appears where two quasienergies become degenerate [40]. This leads to a macroscopic number of degeneracies in the eigenvectors of the full Hamiltonian. Exceptional points are particular to non-Hermitian systems and have many interesting properties [15]. These exceptional points are also seen to appear in the \(N=2\) Ising case for complex values of \(\lambda\).