Stability of zero modes in parafermion chains
Abstract
Onedimensional topological phases can host localized zeroenergy modes that enable highfidelity storage and manipulation of quantum information. Majorana fermion chains support a classic example of such a phase, having zero modes that guarantee twofold degeneracy in all eigenstates up to exponentially small finitesize corrections. Chains of ‘parafermions’—generalized Majorana fermions—also support localized zero modes, but, curiously, only under much more restricted circumstances. We shed light on the enigmatic zeromode stability in parafermion chains by analytically and numerically studying the spectrum and developing an intuitive physical picture in terms of domainwall dynamics. Specifically, we show that even if the system resides in a gapped topological phase with an exponentially accurate groundstate degeneracy, higherenergy states can exhibit a splitting that scales as a power law with system size—categorically ruling out exact localized zero modes. The transition to powerlaw behavior is described by critical behavior appearing exclusively within excited states.
I Introduction
Topological quantum computing represents a promising and conceptually elegant route to scalable quantum computation.Kitaev (2003); Nayak et al. (2008) Underlying this approach are topological phases of matter that harbor emergent particles known as nonAbelian anyons. Such particles exhibit two defining features: they generate a groundstate degeneracy that scales exponentially with the number of anyons present in the system,Moore and Read (1991); Nayak and Wilczek (1996) and braiding the anyons around one another noncommutatively rotates the system’s quantum state within the groundstate manifold.Moore and Seiberg (1989); Witten (1989) The advantage that topological quantum computation offers over more traditional quantum computing schemes is that information is encoded and processed nonlocally in the braiding history of nonAbelian anyons. Local environmental perturbations that ordinarily cause decoherence are thereby expected to be relatively benign.
One conceptually simple realization of nonAbelian anyons are quasiparticles (or defects) that bind exponentially localized, topologically protected zeroenergy modes.Read and Green (2000); Kitaev (2001) These modes are described by operators with appreciable weight on some finite length scale and that commute with the system’s Hamiltonian up to exponentially small corrections , where is the separation between adjacent anyons. The localized character of the zero modes ensures welldefined braiding relations for the anyons that bind them, while the fact that they carry no energy guarantees the groundstate degeneracy necessary for nonAbelian statistics. It is worth emphasizing, however, that zero modes so defined make an extremely strong statement about the system’s spectrum: they imply an exponentially accurate degeneracy not just for ground states, but in fact for all eigenstates.
The Kitaev chain Kitaev (2001) provides an illuminating example. The Hamiltonian reads
(1) 
where denotes a Hermitian Majorana fermion operator satisfying the commutation relation . As Fig. 1 illustrates, the couplings and favor competing Majorana dimerization patterns. In the special case with , the outermost Majorana operators completely decouple from the rest of the system, commute with , and thus form exact localized zero modes. Consequently, every eigenstate assumes at least twofold degeneracy. Turning on finite preserves the zero modes—and with it the degeneracy in the spectrum—provided that .Kitaev (2001) Rather than localizing to one site the zero modes then simply decay exponentially into the bulk on the scale of the correlation length (which diverges at ). The survival of localized zero modes indicates that the chain resides in the same topologically nontrivial phase for any . Throughout this phase the ends of the chain behave as nonAbelian anyons whose nontrivial exchange statistics can be meaningfully harvested in networks.Alicea et al. (2011); Clarke et al. (2011); Halperin et al. (2012); Bonderson (2013) Braiding this type of anyon, however, enables only rather limited (i.e., nonuniversal) faulttolerant quantum information processing.Bravyi and Kitaev (2002)
In the pursuit of nonAbelian anyons with greater utility for quantum computation, a variation of the Kitaev chain due to FendleyFendley (2012) has proven influential. The Hamiltonian for this ‘parafermion chain’ (which Sec. II discusses in depth) is given by
(2) 
Here denotes a parafermion operator satisfying a generalization of the Majorana fermion algebra:
(3) 
Figure 1 still illustrates the structure of the couplings, which we assume are real and nonnegative throughout.
With the outermost operators drop out from the Hamiltonian—precisely as in the Kitaev chain—and represent localized parafermion zero modes
Similarities with the Majorana case, however, largely end here. Most strikingly, there is strong evidence that localized zero modes disappear entirely upon introducing arbitrarily small !Fendley (2012) Such dramatic behavior defies intuition given that for the system resides in a gapped topological phaseFendley (2012); Motruk et al. (2013); Bondesan and Quella (2013) where one would naturally expect to yield only perturbative effects. Stable localized zero modes were instead found only in a ‘chiral’ deformation of the Hamiltonian wherein with nonzero [more precisely, Fendley constructed localized zero modes when ].
Understanding the stark differences from the Kitaev chain and diagnosing implications for quantum information applications seem particularly pressing given the growing literature devoted to realizing parafermion zero modes (see, e.g., Refs. Barkeshli and Qi, 2012; Clarke et al., 2013a; Lindner et al., 2012; Cheng, 2012; Vaezi, 2013; Barkeshli et al., 2013a; Barkeshli and Qi, 2013; Barkeshli et al., 2013b, c; Hastings et al., 2013; Oreg et al., 2014; Burrello et al., 2013; Mong et al., 2014; Klinovaja and Loss, 2014, 2013; Barkeshli et al., 2014; Klinovaja et al., 2014; Clarke et al., 2013b; Zhang and Kane, 2014; Orth et al., 2014; Li et al., 2014; Tsvelik, 2014). Despite all this work, the curious state of affairs regarding the stability of zero modes in parafermion chains has remained largely unexplored. The purpose of this paper therefore is to explain the generic absence of localized zero modes in Eq. (2) as well as their resurrection in the chiral case Fendley (2012).
Based on various complementary analytical and numerical methods, our work paints the following picture: In the topological phase with the ground states, as expected, remain threefold degenerate up to corrections that vanish exponentially as one approaches the thermodynamic limit. Surprisingly, however, even the lowestlying excited states that would otherwise be exactly degenerate at exhibit a powerlaw splitting with system size for , implying the destruction of localized zero modes. We show that the onset of powerlaw splitting can be understood via domainwall tunneling processes that simply have no analogue in the Kitaev chain. We further demonstrate that chirally deforming the Hamiltonian frustrates these domainwall tunneling events, eventually restoring exponential splitting of the states (at least in part of the spectrum) consistent with zeromode revival.
One noteworthy implication of our work is that the disappearance of localized zero modes should not be conflated with a demise of nonAbelian anyons.
On the contrary, throughout the topological phase exhibited by Eq. (2) the parafermion chain still allows one to demonstrate nonAbelian statistics since the allimportant groundstate degeneracy persists modulo exponentially small corrections.
For such cases it should be possible to define weaker zero mode operators that are localized and commute with a projected Hamiltonian.
For completeness we note that parafermion braiding, while carrying some advantages over the Majorana case, remains nonuniversal.Clarke et al. (2013a); Lindner et al. (2012); Hastings et al. (2013) One can, nevertheless, leverage parafermionic systems to generate new twodimensional phases that do permit universal topological quantum computation.Mong et al. (2014) Interestingly, similar physics can even appear in local bosonic twodimensional systems.Barkeshli et al. (2014).
To flesh out the above results, we begin in Sec. II by describing basic properties of the parafermion chain model—in particular explaining a nonlocal mapping to ‘spins’ of the threestate Potts model—and introduce the criteria used for evaluating the existence of zero modes. Sections III and IV then explore the ground states and excited states of the Hamiltonian using perturbation theory, exact diagonalization of a truncated Hilbert space model, and density matrix renormalization group (DMRG) simulations. Finally, we conclude in Sec. V by highlighting additional implications and extensions of this study.
Ii The Model and ZeroMode Criterion
The most general parafermion chain Hamiltonian studied in this paper is given by
(4) 
where again satisfies the properties in Eq. (3). Note that in total the system consists of parafermion sites (to define a sensible Hilbert space this number is necessarily even). Without loss of generality we will restrict the chiral phase appearing in the second term to the range , since symmetry relates Hamiltonians with and to Eq. (4).
As noted above, in the limit the existence of zero modes at the ends of the chain is obvious since then for any . To appreciate the implications of these zero modes it is useful to define a ‘triality’ operator
(5) 
akin to the total fermion parity in the Kitaev chain. Since , admits eigenvalues , , or , where
(6) 
For any choice of couplings commutes with the Hamiltonian. Crucially, however, the zeromode operators and do not commute with —they cycle the triality by . It follows that the entire spectrum can be grouped into triplets of energy eigenstates with trialities , , and that are exactly degenerate at . The spectrum of the Kitaev chain in the analogous limit consists of degenerate doublets with opposite fermion parity.
Our goal is to explore the fate of these localized zeromode operators in the generic situation with . We now define the precise criteria used in evaluating whether or not exact edge zero modes exist. Finitesize effects at nonzero generically split the exact degeneracy between different triality states except with finetuning. Let denote the system’s energies, where labels the triality and indexes the levels such that . The existence of exponentially localized zero modes implies that holds for every with some length scale . An exceedingly useful corollary is that the existence of such modes may be categorically ruled out in large swaths of parameter space merely by demonstrating subexponential (e.g., powerlaw) splitting of a single triplet of energy levels . Demonstrating the presence of exact zero modes, by contrast, poses a much more difficult problem, as doing so requires proving a global property of all energy levels. In this paper we will content ourselves with identifying regimes where zero modes are definitely absent.
One can obtain a great deal of intuition by exploiting an exact mapping between the parafermion chain Hamiltonian in Eq. (4) and the chiral threestate Potts model. This mapping is analogous to that between the Kitaev chain and the transversefield Ising model, and accordingly is implemented with a variation of the JordanWigner transformation introduced by Fradkin and KadanoffFradkin and Kadanoff (1980). In particular, writing
(7) 
decomposes the parafermions via strings of Potts model ‘spin’ operators that satisfy
(8) 
All other commutators amongst and vanish, a straightforward result of Eq. (3).
Under this nonlocal change of variables the Hamiltonian becomes
(9) 
which defines a local bosonic model whose states are much easier to analyze than those in the parafermion chain. We stress, however, that both models are equivalent and exhibit precisely the same spectra. Equation (9) can be obtained from the anisotropic limit of a wellknown classical lattice model. With and this is the ferromagnetic threestate Potts model. When , this is typically called the chiral clock or sometimes the chiral Potts model.
Using a basis of eigenstates, denoted , the Potts operators on a given site can be represented as
(10) 
We then have
(11) 
from which it follows that cycles the ‘spin’ measured by . Notice that and represent a straightforward generalization of anticommuting Pauli matrices in the Ising model. The triality defined earlier is simply the generator of a global symmetry exhibited by the Potts Hamiltonian, winding every spin:
(12) 
Hereafter we find it much more illuminating to work in the Potts model representation. Thus it is worth translating into Potts language the consequences of localized zero modes in the parafermion chain. For simplicity, consider and suppose for the moment that —where parafermion zero modes definitely exist. In this case Eq. (9) reduces to the ferromagnetic Potts model with a vanishing ‘transverse field’. The Hamiltonian has three brokensymmetry ground states that transform into one another under the action of in Eq. (12); we label these by , , and . One can of course also define a basis of Schrödingercatlike ground states with definite triality via
(13) 
The zeromode operators and cycle the system amongst the degenerate triplet defined in Eqs. (13) [to see this recall Eq. (7)]. Similar conclusions hold for the excited states, which one can fruitfully view in terms of domain walls (i.e., kinks and antikinks) separating different ferromagnetically ordered regions. All such excited states may also be grouped into degenerate triplets of triality eigenstates that transform into one another under the action of the localized parafermion zeromode operators. Analogous results apply to the chiral case with despite the fact that chirality can nontrivially rearrange the spectrum (see the next section).
Restoring nonzero lifts the degeneracy among the groundstate and excited triplets via a transparent physical mechanism. Namely, the term in Eq. (9) creates mobile domain walls that can tunnel the system between states related by a global transformation. Repeated action of the term can, for instance, send , , and , thereby splitting the threefold groundstate degeneracy in the ferromagnetic case. The question we address in the following sections is thus: how do the degeneracies split?
Iii Perturbative Regime
With the insights of the previous section we now distill our earlier criterion for the existence of localized zero modes to the following question: Do domainwall hopping processes preserve the degeneracy amongst all triplets of energy eigenstates, up to corrections that decay to zero exponentially with system size? If the answer is ‘no’ then localized zero modes in the parafermion chain are ruled out. In this section we address this question by studying both the splitting of ground states and of lowlying excited states in the limit . Section IV then explores the complementary regime where is of order one.
As described above the entire spectrum consists of triplets of exactly degenerate states when . With , degenerate perturbation theory allows one to quantify the effect of nonzero on these triplets. In what follows we discuss the ferromagnetic model with from this perturbative perspective, then attack the chiral case with , and finally address the antiferromagnetic limit . In both the ferromagnetic and antiferromagnetic limits we will show that localized zero modes certainly do not exist even for arbitrarily small but finite , consistent with Ref. Fendley, 2012. Rigorous conclusions are harder to obtain for the chiral case, though in our analysis we definitively rule out zero modes over ranges of and in the vicinity of and ; we argue that outside of these regimes chirality restores zero modes in the perturbative regime. Figure 2 summarizes our results for this section.
iii.1 Ferromagnetic limit,
Figure 3 (left side) illustrates the spectrum for the ferromagnetic case. Let us now discuss these levels in greater detail. There are three ground states , and corresponding to all spins uniformly polarized along one of three possible directions. As in the Ising model every excited state can be viewed in terms of domain walls between differently polarized regions, although here six flavors exist rather than two because of the larger groundstate degeneracy. The energy cost of a domain wall comes from the single link connecting different values of the spins, giving
(14) 
The lowest excited levels of Fig. 3 are singledomainwall states. The six flavors are denoted as , , , and , where for example indicates the presence of a domain wall where spins on the left are in ground state and on the right are in . Each flavor can sit at any of positions, so a total of such states exist—all exactly degenerate at . There exists a larger set of degenerate twodomainwall states corresponding to , etc., and so on up the spectrum. Crucially, the three ground states and each set of degenerate domainwall states form subspaces that are energetically wellseparated from one another by a gap .
Including infinitesimal enables domainwall creation, annihilation, and hopping. The splitting of the levels within a given subspace can then be computed using standard degenerate perturbation theory methods. We are particularly interested in processes that take one eigenstate to another related by a global transformation, e.g., or . These are precisely the processes that produce finitesize splitting of the degeneracy encoded by any localized zero modes.
The splitting of the groundstate energies due to nonzero follows from standard arguments. These imply that mixing between the three ground states is suppressed exponentially in system size. Namely, the ‘cheapest’ way to evolve from, say, to is to create an domain wall at one of the end chain—leaving the system in an excited state, tunnel the wall over to the opposite end, and annihilate the domain wall to reenter the groundstate manifold. One can visualize the process graphically via Fig. 4. Here the horizontal axis denotes the position along the chain, the vertical axis represents the perturbation step (which can also be interpreted as time), and the diagonal line indicates the domainwall trajectory. Such a process necessitates exiting the groundstate subspace for a macroscopic number of steps, resulting in a groundstate splitting
(15) 
In short, an energy barrier prevents efficient mixing of these states through local spin flips. Thus the groundstate degeneracy encoded by the exact zero modes of the limit survives the introduction of nonzero , up to exponentially small corrections in system size.
Qualitatively different behavior arises in the singledomainwall sector. In fact here the energy barrier that suppresses groundstate mixing disappears altogether so that the dominant contribution to the splitting that we seek arises already at first order in degenerate perturbation theory. To see this consider the process illustrated in Fig. a which takes without leaving the singlewall subspace. This process proceeds by first hopping the domain wall all the way to the right end of the chain. If we were dealing with an Ising model, then the only way to remain in the original subspace would be to subsequently tunnel the domain wall leftward—which fails to accomplish the domainwall change required. The Potts chain, however, offers an alternative: the rightmost spin can wind, converting to . This newly formed domain wall can then tunnel all the way to the chain’s left end where it can similarly convert to . Moving the domain wall back to its original location completes the process.
The splitting amongst triplets of singledomainwall states resulting from such processes certainly depends subexponentially on system size. Indeed, the energy denominators responsible for the exponential dependence in Eq. (15) are entirely absent since the system never leaves the degenerate singlewall subspace (again, all the action takes place at first order in perturbation theory). We can explore the splitting more quantitatively by examining the Hamiltonian projected onto the singledomainwall subspace, which resembles a tightbinding model for six flavors of particles reflecting the domainwall types; for details see Appendix B. If one were to ignore the fact that these particles are domain walls and give them periodic boundary conditions, then the energy would be
(16) 
where labels the momentum. These energy levels are trivially sixfold degenerate since the domainwall flavor is conserved with periodic boundary conditions. With open boundaries, however, the system’s edges mediate backscattering processes that do convert domainwall flavors into one another as in Fig. 4. Since the backscattering happens within a continuum of extended levels, powerlaw splitting generically arises for all triplets of singlewall states even with arbitrarily small nonzero .
We can compute analytically the powerlaw splitting arising with open boundaries in the singledomainwall projection. The key observation is that there is only one way for domain walls to change flavor at each end, e.g., to at the right boundary. As a consequence, we can ‘unfold’ the system into a periodic chain of size with a single flavor. In this unfolded picture each flavor of the original model corresponds to a region of size so that a domain wall bouncing from an edge is replaced by a particle moving into a different region. Cycling through all six flavors corresponds to going around the periodic chain.
We can then obtain the spectrum using Fourier analysis; details of the calculation appear in Appendix B. The solution exhibits an interesting form that can be anticipated from previously known results for the scaling limit of the threestate Potts model: the amplitude for changing flavor when a kink scatters off the boundary is of magnitude 1, while the amplitude for not changing (i.e., just bouncing back as is) vanishes Chim (1995). Indeed we find that the eigenstates take the form of simple right or leftmoving plane waves that propagate unreflected around the unfolded periodic chain. More quantitatively, the results are as follows. Labeling the energies of a single domainwall triplet as for integer , we define
(17) 
In the large limit this becomes
(18) 
Near the band edges a powerlaw of the form emerges; elsewhere a scaling takes over. We verified the splittings extracted above by numerically simulating the projected Hamiltonian with open ends; very large system sizes can be readily simulated since the truncated Hilbertspace dimension scales only linearly with . Figure 5 (red circles) illustrates our simulation results. The vertical axis denotes the natural log of the splitting of the three lowestenergy domainwall eigenstates while the horizontal axis represents . Clear splitting indeed appears as shown by the solid line.
The powerlaw behavior captured here immediately precludes the existence of exponentially localized zero modes in the ferromagnetic limit for any small but finite . We stress that this result is quite robust as processes that involve neglected subspaces are higherorder in and in the perturbative regime cannot remove the powerlaw splitting arising at first order.
iii.2 Chiral case,
Having ruled out exact localized zero modes in the ferromagnetic case, we turn now to the chiral regime where lies between 0 and noninclusive. It is instructive to first discuss the spectrum illustrated schematically on the right side of Fig. 3. For any in this range remain the unique exact ground states. Chirality does, however, alter the excited states in an important way. More precisely, the , , and domain walls now carry energy that differs from the energy carried by their mirror counterparts , , and . Explicitly,
(19) 
The singledomainwall sector therefore splits into two degenerate subspaces, the twowall sector splits into three (, and all yield different energy), etc. Note that Fig. 3 illustrates the levels only for small chiral phases ; larger reorders the excited states as we discuss below.
Chirality imparts only minor quantitative effects on the groundstate splitting induced by nonzero . Our arguments from the previous subsection indeed carry over straightforwardly since a finite gap to the first excited state persists whenever . Exponentially small splitting is thus again guaranteed over a range of , the primary difference being that the denominator in Eq. (15) should be replaced by the gap .
For the excited states, by contrast, chirality yields more dramatic consequences. This originates from the reduced degeneracy of the domainwall eigenstates relative to the ferromagnetic case, which tends to suppress the tunneling processes that previously led to powerlaw splitting.
It is simplest to examine small chiral phases where the system remains close to the ferromagnetic point; values of near the antiferromagnetic point will be discussed separately below. Following Sec. III.1 we can capture the leading effects of on the lowlying excited states by projecting the Hamiltonian onto the singledomainwall subspace. The effective Hamiltonian (see Appendix B) again resembles a tightbinding model for six flavors of particles, half of which now experience different onsite energies resulting in two types of excitation branches. With periodic boundary conditions one obtains band energies
(20) 
that exhibit threefold degeneracy for each momentum . The upper branch represents the bands formed by states while the lower branch similarly corresponds to .
Our aim is to now understand how these branches mix at the ends of a system with open boundaries. As Fig. 6 illustrates, there are three distinct regimes distinguished by the degree to which the bands overlap. Implications for zero modes depend sensitively on the ratio of to the energy difference . It is thus often convenient to utilize the ratio
(21) 
when analyzing the spectrum.
In the limit , the bands exhibit no overlap as shown in Fig. 6(a). The hybridization between the upper and lower branches at the edges of the chain thus can be treated perturbatively. If we revisit the process sketched in Fig. a, a clear qualitative difference from the ferromagnetic case appears. Now, when the domain wall (originating from the upper branch) tunnels to the right end of the chain, conversion to necessitates moving to the energetically wellseparated lower band. Returning the domain wall to the left end—where it can transition back to a wall in the upper branch—requires the system to remain in the lowerenergy space for a macroscopic number of steps. The key physics is that chirality precludes resonantly switching domainwall flavors at the system’s boundaries. Thus such processes produce exponential rather than powerlaw splitting among all triplets of singlewall eigenstates. Numerical simulations of the projected Hamiltonian confirm this picture. The yellow diamonds in Fig. 5 illustrate the natural log of the splitting for the lowest three singledomain states versus system size assuming ; the splitting is wellcaptured by an exponential.
The obvious indications for the absence of localized zero modes that we uncovered in the ferromagnetic limit thus disappear upon introducing a small chiral phase in the regime. This finding is in harmony with the construction in Ref. Fendley, 2012 of edge zero modes at and with sufficiently small. Our results here do not prove, however, that zero modes persist to finite , as subexponential splitting could arise from various other sources including tunneling within multidomainwall subspaces and higherorder processes in perturbation theory.
Further observations do nevertheless suggest that localized zero modes indeed persist. First, we have studied numerically the effective Hamiltonian projected onto the twodomainwall subspace and again found only exponential splitting when . Second, it is possible to address higherorder processes rather efficiently by developing a set of Feynmandiagramlike rules governing domain wall dynamics. Appendix A describes the methodology. There we briefly sketch a calculation indicating that through second order in the splitting amongst singlewall states remains exponential. This is a result of the cancellation of contributions from various diagrams related by symmetries of the allowed tunneling processes. Though we have not proven that the cancellation occurs at all orders, the symmetries by which it works at low orders suggest that this is the case.
Lowering brings the upper and lower branches closer to one another in energy until they eventually begin to overlap at a critical . Indirect gapless interband excitations are then permitted as shown in Fig. 6(b). The critical value follows from the condition that , yielding . This value corresponds to an dependent critical chiral phase
(22) 
With open boundaries the chain’s edges provide the momentum transfer needed here to resonantly scatter domain walls in the upper branch into walls in the lower branch. Processes like that of Fig. a consequently yield powerlaw splitting among the triplets at the bottom of the upper band and top of the lower band.
With the bands overlap in a finite energy window [see Fig. 6(c)] within which all triplets similarly admit powerlaw splitting. This picture is supported by our numerics in Fig. 5 for , which reveal splitting that is exponential for the lowest excited triplet (squares) yet powerlaw for triplets at the middle of the lower branch (diamonds). In the perturbative regime, we can therefore conclusively rule out localized zero modes not just in the ferromagnetic limit, but in fact along the finite interval .
Mapping the problem to an ‘unfolded’ system as described in the previous subsection and Appendix B provides an enlightening perspective. The projected singlewall Hamiltonian for the open chain again effectively describes a particle hopping in a periodic system of length . Due to chirality, however, the unfolded Hamiltonian now includes a squarewell potential that alternates between and every sites. For simplicity let us consider small chiralities and focus on states near the bottom of the upper and lower branches so that the particle’s kinetic energy can be approximated with a freeparticle dispersion . (Analogous arguments follow for states near the top of the bands.) Given sufficient kinetic energy to overcome the squarewell barrier, the particle will efficiently sample the entire extended periodic chain. Translating back to Potts model language, this means that domain walls can resonantly convert between different flavors at the edges of the system—resulting in powerlaw splitting at such energies. With insufficient kinetic energy the particle will remain predominantly within the squarewell minima, decaying evanescently into the regions corresponding to squarewell maxima. The evanescent decay translates into inefficient domain wall conversion and a splitting of triplets that decays exponentially with system size.
Suppose next that the chiral phase is slightly below so that the system resides close to the antiferromagnetic point. Because in this regime, the ordering of domainwall subspaces changes dramatically compared to the cases examined above. At the lowestenergy subspace consists of states while the secondlowest subspace arises from the twodomain wall states . By contrast the states , that we often invoked earlier lie much higher up in the spectrum so that conversion processes sketched in Fig. a are now irrelevant.
We can, however, identity a new series of domainwall processes that yield powerlaw splitting over a range of values. As an example, one can convert to an domain wall related by a global transformation by utilizing the lowestlying twowall subspace. This is accomplished by moving the original domain wall to the system’s right end and then inserting a second domain wall, i.e., . Upon hopping both domain walls to the other end one can remove the leftmost wall, sending . The remaining domain wall, having switched flavors, finally returns to its original location. Such processes generically split the degeneracy among singlewall triplets in the lowest excited subspace.
To quantify the induced splitting it is once again useful to project onto the relevant subspaces. For periodic boundary conditions on the resulting effective Hamiltonian, the singlewall states broaden into a band with energies given in Eq. (20). The twowall states are similar, although the energies for the latter are nontrivial since domain walls interact. Introducing open boundaries enables transitions between these subspaces at the system’s edges, which can supply any momentum necessary. Powerlaw splitting arises for triplets in the lower band if they overlap in energy with the upper band. The critical at which powerlaw behavior first appears follows from , where the right side approximates the energy needed to insert an extra domain wall into a zeromomentum singlewall state. (Note that this estimate neglects domain wall interactions; this is reasonable since they are dilute.) We then obtain
(23) 
Interestingly, the dependence is remarkably similar to Eq. (22) despite the rather different processes involved. In the interval at least one excited triplet exhibits powerlaw splitting in the perturbative regime, ruling out localized zero modes here too.
iii.3 Scaling behavior
In the previous subsections we gained quite a bit of mileage by studying an effective Hamiltonian projected onto the singledomainwall subspace. This truncated model is expected to capture the dominant contribution to the splitting amongst triplets of lowlying excited states in the limit . Through analytical arguments corroborated by numerics we concluded that chirality induces a sharp change in the level spacing among triplets of excited states.
More precisely, at we showed that all singledomainwall triplets exhibit powerlaw splitting with system size for arbitrary nonzero . Introducing small nonzero immediately restores exponential splitting for a range of the lowest and highestenergy triplets but preserves powerlaw splitting for all intermediate singlewall states. Further increasing shrinks the window of powerlawsplit states until exponential splitting fully returns beyond a critical value given in Eq. (22)—presumably restoring exact localized zero modes. It is tempting to view this behavior as indicating chiralitytuned critical behavior in the excitedstate spectrum, despite the fact that no singular behavior arises in the groundstate sector. We now provide evidence for such a scenario by showing that our numerical results exhibit data collapse consistent with a simple scaling ansatz.
Provided , the singledomainwall splitting of a given triplet should depend only on from Eq. (21) and the effective system size for domain walls. Let us consider a particular triplet that transitions from exponential to powerlaw splitting at some (statedependent) critical value . We stress that is distinct from discussed in the previous subsection; the former applies to just one particular multiplet. Under rescaling of the length by we postulate the following scaling form,
(24) 
where denotes the deviation from and are critical exponents to be determined. Setting then allows us to write
(25) 
for some function .
As a concrete example consider the lowestlying excited triplet, for which .
Our results for the ferromagnetic case imply that in the perturbative regime while is some nonzero constant.
iii.4 Antiferromagnetic limit,
Finally, we consider the antiferromagnetic limit where . One can anticipate that this case is rather special from the excitation energy quoted in Eq. (19); in particular, with half of the domainwall flavors cost no energy, resulting in a macroscopic groundstate degeneracy in the limit. Here localized zero modes disappear for any nonzero in the perturbative regime, just as in the ferromagnetic case (and for similar reasons). Indeed, the construction of edge zero modes in Ref. Fendley, 2012 works only for , and so fails at the antiferromagnetic point.
To demonstrate this result directly we need only examine the groundstate manifold perturbed by infinitesimal . It is convenient here to perform a gauge transformation , which effectively transforms the Potts Hamiltonian to one with . In the transformed model the term is minimized by any state for which no two adjacent spins align with one another (hence the ‘antiferromagnetic’ nomenclature). Thus the limit supports ground states. Consider, for instance, a ground state of the form . As usual we are interested in tunneling processes that connect states related by a global transformation, i.e., , since such events split the degeneracy encoded by the exact zero modes present at . One possible pathway is for repeated action of the operators to first wind all of the even spins, sending , followed by all odd spins, sending . The system remains in the groundstate manifold throughout such processes, so that a splitting amongst groundstate triplets arises at first order in degenerate perturbation theory.
This splitting is thus expected to scale as a powerlaw with system size—similar to what we observed previously for the ferromagnetic case—precluding localized zero modes even for arbitrarily small nonzero as claimed. Interference effects that might alter this conclusion are absent since all nonzero matrix elements in the groundstate manifold are unity. Higherorder processes that require exiting the groundstate manifold are also negligible in the limit considered here.
There is a complementary way of arguing for powerlaw splitting in the antiferromagnetic case. In the perturbative regime one can obtain an effective lowenergy Hamiltonian for the chain by projecting the operator appearing in the term onto the groundstate manifold. (The term of course projects trivially.) To do so it is useful to define bond operators that measure whether the spin at a given site is wound ‘clockwise’ or ‘counterclockwise’ relative to its neighbor. More specifically, we can introduce Pauli matrices such that with the groundstate projector. As an example for the spin pair while for the mirrored configuration . With this notation one can see that flips both and if they have opposite signs, and otherwise projects trivially. More formally, we obtain
(26) 
Thus the bulk maps onto an XY model
Iv NonPerturbative Regime
Both the analytical arguments and numerical diagonalization exploited in Sec. III break down when becomes of order unity. Here we complement our earlier perturbative analysis using extensive DMRG simulations. To understand the nature of the spectrum in the nonperturbative regime, we computed the splitting among the ground states and firstexcited triplet for a wide range of and over which DMRG exhibits good convergence. The simulations turn out to perform particularly well at larger .White (1992); Schollwöck (2005); McCulloch (2007) System sizes were taken to be at least to minimize finitesize effects but sufficiently small to avoid the splitting falling below machine precision; the maximum considered varies from to depending on the Hamiltonian parameters. More details of our simulations are described in Appendix C.
To quantify the splittings obtained from the DMRG, our data were fit to an exponential form . The color scales in Figs. 8 and 8 show the optimal value versus (radius) and (angle) for the groundstate and first excited splittings, respectively. The ground states exhibit an exponential splitting characterized by a nonzero throughout the regime of convergence, which notably does not include the antiferromagnetic limit because of its gaplessness. As in the perturbative regime, however, richer physics arises for the first excited triplet. In the ferromagnetic limit we obtain , indeed consistent with the powerlaw splitting found in the perturbative regime. For the chiral cases, by contrast, takes on finite values indicating exponential splitting with system size, also in agreement with our perturbative analysis.
We showed earlier that powerlaw splitting for excited states survives over a finite range of chiral phases . However, accessing these levels within DMRG rapidly becomes prohibitive as one increases the energy. Thus the finite transitions are difficult to capture numerically. Instead we crudely estimate the global zeromode stability regions by naively extrapolating our results for the critical values obtained for the perturbative regime in Sec. III.2. Figure 8 displays the results of our extrapolation. The shaded regions indicate the parameter space over which zero modes are expected to survive based on our perturbative criterion. While of course this scheme is not quantitatively reliable, we do expect to capture the qualitative trends. [One feature that is not expected to be robust is the ‘accidental’ symmetry of Fig. 8 with respect to sending , as is clear from the numerics in Fig. 8. This pathological property should be absent even in the perturbative regime upon including interactions between domain walls, which were neglected in our estimate of .]
It is interesting to discuss these findings in relation to Ref. Fendley, 2012, where using an iterative method localized zero modes were constructed explicitly in the limit where the control parameter
(27) 
is much less than one. In this case corrections to the zero modes could be arranged in a series that clearly decays exponentially into the bulk of the chain. Remarkably, our extrapolation in Fig. 8 follows the lobelike form of the control parameter —strongly suggesting that we have correctly identified the essential physics that determines the robustness of exact zero modes in parafermion chains.
V Conclusions
In this paper we employed a variety of techniques to diagnose the puzzling stability of local parafermion zero modes—which differs markedly from the Majorana case—identified by Ref. Fendley, 2012. Viewing the physics from the mathematically equivalent lens of the chiral threestate Potts model proved particularly illuminating. In Potts language physically intuitive domainwall hopping/conversion processes produce powerlaw splitting among lowlying excited states in the ferromagnetic and antiferromagnetic limits of the model, ensuring the generic demise of localized parafermion zero modes in both cases. Chirally deforming the Hamiltonian tends to suppress these processes, and for sufficiently large deviations from the ferromagnetic/antiferromagnetic limit restores exponential splitting of all lowlying excited states that we examined. We speculate that the restoration of exponential splitting coincides with the emergence of localized zero modes; although a general proof remains unavailable such a scenario gels nicely with the results from Ref. Fendley, 2012.
We also showed that the transition from powerlaw to exponential behavior reflects a subtle type of chiralitytuned quantum criticality that, interestingly, emerges only at energies above the ground states. Indeed, except in the antiferromagnetic limit the ground states remain degenerate up to exponentially small corrections in all cases and do not exhibit any singular behavior as a function of chirality.
Although we focused on parafermion chains for simplicity, many of our results extend straightforwardly to systems (where similar stability issues ariseFendley (2012)). For instance, it is clear that here too domainwall tunneling events generically preclude exponentially localized modes in the ferromagnetic and antiferromagnetic limits. Our findings also naturally explain the comparative robustness of localized Majorana zero modes in the Kitaev chain, since the relevant domainwall processes that we invoked for the Potts model have absolutely no analogue in that context.
Since we motivated this work from the vantage point of topological quantum computation, it is worth reemphasizing that the existence of localized zero modes as defined here is most certainly overkill for this application. Harnessing nonAbelian statistics in parafermion chains (and related systems) merely requires degenerate ground states of a topological phase to within exponential accuracy; this weaker condition appears much more broadly as noted above, even when exact zero mode operators are definitely absent. Systems supporting exact zero modes do, nevertheless provide the appealing possibility of performing topological quantum computation at finite energy density—a possibility first proposed in the framework of manybody localizationHuse et al. (2013) (see also Ref. Bahri et al., 2013). We suggest that parafermion chains offer particularly interesting platforms to explore in this regard. Apart from chirality, localization via quenched disorder should provide another means of suppressing the domainwall processes that produced powerlaw splitting in the ferromagnetic limit. It would be quite interesting to perform largescale exactdiagonalization studies to explore this scenario further in future work.
Acknowledgements.
We are indebted to David Clarke, Nate Lindner and Lesik Motrunich for illuminating conversations, as well as to Miles Stoudenmire for invaluable discussions on numerics. We also acknowledge funding from the NSF through grants DMR1341822 (A. J. & J. A.) and DMR/MPS1006549 (P. F.); the Alfred P. Sloan Foundation (J. A.); the Sherman Fairchild Foundation (R. M.); the Caltech Institute for Quantum Information and Matter, an NSF Physics Frontiers Center with support of the Gordon and Betty Moore Foundation; the Caltech Summer Undergraduate Research Fellowship program along with partial support from the family of Jean J. Dixon (A. J.); and the Walter Burke Institute for Theoretical Physics at Caltech.Appendix A Path cancellations at higher orders in perturbation theory
Section III performed a perturbative analysis that incorporated the lowestorder processes that produced a splitting among eigenstates with different trialities. For excited states it is possible to capture very similar transformations to that of Fig. a in far fewer perturbation steps if one considers terms higher order in . One thus might worry that such processes, while parametrically smaller in , nevertheless produce the dominant scaling with system size. In the chiral case with , for instance, if powerlaw splitting were to appear at higher orders then that would immediately imply a further reduction in stability window for localized zero modes. However, our DMRG results—which are nonperturbative and hence include all orders—suggest that this is not the case since our simulations captured exponential splitting among the lowest excited triplet for any nonzero chirality. In this appendix we provide more support for this conclusion by examining higherorder perturbative computations.
While cumbersome, higherorder calculations can be greatly facilitated by using of a diagrammatic representation of domainwall events similar to Fig. a. For this purpose it is useful to attach an upwardpointing arrow to lines indicating domain walls of the type , , and and a downward arrow for the other three flavors. The following vertex rules then arise:

The boundary of the system acts as a source and sink for domain walls. This is the oneline vertex.

Two lines may be created or destroyed at a point if and only if they have opposite arrow directions. This is the twoline vertex.

A threeline vertex may exist if and only if all of the arrows are incoming or outgoing.

In a given perturbation step each line may either remain fixed or move one spatial unit.
Since all matrix elements in the perturbation theory are either 0 or 1 the assignment of weights to diagrams conforming to these rules is straightforward. Any step that moves the system out of the original subspace gets penalized by a factor of where is the magnitude of the energy change incurred. The sign arises if the system moves to a higherenergy state, while the sign arises if a lowerenergy state results.
Consider first the diagram in Fig. 9 depicting a secondorder process wherein . By themselves, events like this yield powerlaw splitting with , even in the chiral case. There is, however, a compensating process shown in Fig. 9 where . Since here the system enters a higherenergy sector these events produce the exact opposite contribution . The existence of such cancellations is the main message of this appendix. Preliminary calculations at third order point to a similar outcome, and we expect that they are a generic feature of higherorder perturbation theory in certain regions of parameter space. Again, this conclusion is supported by our DMRG results at intermediate .
Appendix B Projection domainwall Hamiltonians
In the limit one can to a good approximation neglect tunneling processes that mix sectors with different numbers of domain walls. Formally this is achieved by projecting in Eq. (9) onto a subspace with fixed domain wall number. The procedure is generally straightforward though some care is necessary at the boundaries. In the onewall sector the resulting effective Hamiltonian admits a particularly clean block form:
(28) 
with the identity matrix,
(29) 
As noted in Sec. III.1, because of the form of the boundary terms described by and it is possible to recast in terms an enlarged periodic chain of size . In this mapping the domain wall becomes a particle living on sites labeled by position , with the domainwall flavor corresponding to . We will take to represent , , and , and to represent , , and , respectively. Let us now specialize to the ferromagnetic limit, . Denoting the position modulo by , it is straightforward to verify that the eigenstates of the projected Hamiltonian are given by
(30) 
with eigenvalues
(31) 
These energies are the same as those in Eq. (16) up to a constant; however, because of the enlarged system size in this description the allowed momenta are now . Note that the wavefunctions specified in Eq. (30) are also triality eigenstates with eigenvalue
(32) 
[The action of simply maps .]
For extracting the splittings of interest, it is convenient to label the energies and momenta according to both the triplet to which they belong and their triality , i.e., as and . We need only specify that and require that . Using Eqs. (17) and (31), we then find
(33) 
which in the large limit converges to Eq. (18).
For our numerical simulations of the truncated singlewall model, we explicitly split the spectrum into the three triality sectors by simultaneously diagonalizing the triality operator and the effective Hamiltonian. In practice this may be done by adding the two operators together with random coefficients, determining the diagonalizing unitary operator for the combination, and then using it to diagonalize the two operators individually. For almost all choices of random numbers this diagonalizes both operators. Pathological cases may be handled separately, or by rerunning the calculation.
To handle the twowall sector, a Python+Cython+NumPy code was written to explicitly enumerate all states of interest and compute the relevant matrix elements of between them. In either sector the effective Hamiltonians can be diagonalized numerically via standard routines to obtain splittings for any triplet of singledomainwall states.
Appendix C DMRG methods
All DMRG computations in this paper utilized the Developer Branch of ITensor (http://itensor.org) commit 475352f76c6209db865ea4405cb86f665f40fae5. A control file, a Hamiltonian file, and a model file were created based on existing ITensor code and with the assistance of Miles Stoudenmire, one of the authors of ITensor. These files extended ITensor to perform calculations on the Potts model, which is not a native function of the code. The Eigen C++ library version 3.0 was also used in the main control file for diagonalizing the states that ITensor produced. This generates excited states with less groundstate overlap than DMRG alone gives.
The results from ITensor were verified for small system sizes with the opensource Quantum Chains package (https://github.com/adamjermyn/QuantumChains/) written by the first author. Any of the code used in this paper is available upon request from the first author.
Footnotes
 To distinguish from parafermions in conformal field theory, which are related but distinct, these zero modes are sometimes referred to as ‘parafendleyons’.
 For an explicit construction see A. Alexandradinata, C. Fang, N. Regnault, and B. A. Bernevig, unpublished.
 In writing the Potts representation of the parafermion chain, we dropped a factor of in the term; this factor can always be absorbed into the chiral phase .
 Of course in the extreme limit we have . Note also that the scaling forms provided here assume . For larger transverse fields mixing between other sectors matters, and hence becomes a function of three parameters—, and .
 We are grateful to David Clarke for pointing out to us this XY model mapping.
References
 A. Y. Kitaev, Ann. Phys. 303, 2 (2003).
 C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
 G. Moore and N. Read, Nucl. Phys. B 360, 362 (1991).
 C. Nayak and F. Wilczek, Nucl. Phys. B 479, 529 (1996).
 G. Moore and N. Seiberg, Commun. Math. Phys. 123, 177 (1989).
 E. Witten, Commun. Math. Phys. 121, 351 (1989).
 N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
 A. Y. Kitaev, Sov. Phys.–Uspeki 44, 131 (2001).
 J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nat. Phys. 7, 412 (2011).
 D. J. Clarke, J. D. Sau, and S. Tewari, Phys. Rev. B 84, 035120 (2011).
 B. I. Halperin, Y. Oreg, A. Stern, G. Refael, J. Alicea, and F. von Oppen, Phys. Rev. B 85, 144501 (2012).
 P. Bonderson, Phys. Rev. B 87, 035113 (2013).
 S. Bravyi and A. Kitaev, Annals of Physics 298, 210 (2002).
 P. Fendley, J. Stat. Mech. 2012, 11020 (2012).
 To distinguish from parafermions in conformal field theory, which are related but distinct, these zero modes are sometimes referred to as ‘parafendleyons’.
 J. Motruk, E. Berg, A. M. Turner, and F. Pollmann, Phys. Rev. B 88, 085115 (2013).
 R. Bondesan and T. Quella, Journal of Statistical Mechanics: Theory and Experiment 2013, P10024 (2013).
 M. Barkeshli and X.L. Qi, Phys. Rev. X 2, 031013 (2012).
 D. J. Clarke, J. Alicea, and K. Shtengel, Nature Commun. 4, 1348 (2013a).
 N. H. Lindner, E. Berg, G. Refael, and A. Stern, Phys. Rev. X 2, 041002 (2012).
 M. Cheng, Phys. Rev. B 86, 195126 (2012).
 A. Vaezi, Phys. Rev. B 87, 035132 (2013).
 M. Barkeshli, C.M. Jian, and X.L. Qi, Phys. Rev. B 87, 045130 (2013a).
 M. Barkeshli and X.L. Qi, “Synthetic Topological Qubits in Conventional Bilayer Quantum Hall Systems,” (2013), unpublished, 1302.2673 .
 M. Barkeshli, C.M. Jian, and X.L. Qi, Phys. Rev. B 88, 241103 (2013b).
 M. Barkeshli, C.M. Jian, and X.L. Qi, Phys. Rev. B 88, 235103 (2013c).
 M. B. Hastings, C. Nayak, and Z. Wang, Phys. Rev. B 87, 165421 (2013).
 Y. Oreg, E. Sela, and A. Stern, Phys. Rev. B 89, 115402 (2014).
 M. Burrello, B. van Heck, and E. Cobanera, Phys. Rev. B 87, 195422 (2013).
 R. S. K. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, P. Fendley, C. Nayak, Y. Oreg, A. Stern, E. Berg, K. Shtengel, and M. P. A. Fisher, Phys. Rev. X 4, 011036 (2014).
 J. Klinovaja and D. Loss, Phys. Rev. Lett. 112, 246403 (2014).
 J. Klinovaja and D. Loss, “Timereversal invariant parafermions in interacting Rashba nanowires,” (2013), unpublished, arXiv:1312.1998 .
 M. Barkeshli, Y. Oreg, and X.L. Qi, “Experimental proposal to detect topological ground state degeneracy,” (2014), unpublished, arXiv:1401.3750 .
 J. Klinovaja, A. Yacoby, and D. Loss, “Kramers pairs of Majorana fermions and parafermions in fractional topological insulators,” (2014), unpublished, arXiv:1403.4125 .
 D. J. Clarke, J. Alicea, and K. Shtengel, “Exotic circuit elements from zeromodes in hybrid superconductor/quantum Hall systems,” (2013b), unpublished, arXiv:1312.6123 .
 F. Zhang and C. L. Kane, “Timereversalinvariant fractional Josephson effect,” (2014), unpublished, arXiv:1404.1072 .
 C. P. Orth, R. P. Tiwari, T. Meng, and T. L. Schmidt, “NonAbelian parafermions in timereversal invariant interacting helical systems,” (2014), unpublished, arXiv:1405.4353 .
 W. Li, S. Yang, H.H. Tu, and M. Cheng, “Criticality in translationinvariant parafermion chains,” (2014), unpublished, arXiv:1407.3790 .
 A. M. Tsvelik, “ parafermion zero modes without Fractional Quantum Hall effect,” (2014), unpublished, arXiv:1407.4002 .
 For an explicit construction see A. Alexandradinata, C. Fang, N. Regnault, and B. A. Bernevig, unpublished.
 M. Barkeshli, H.C. Jiang, R. Thomale, and X.L. Qi, “Generalized Kitaev models and slave genons,” (2014), unpublished, arXiv:1405.1780 .
 E. Fradkin and L. P. Kadanoff, Nucl. Phys. B 170, 1 (1980).
 In writing the Potts representation of the parafermion chain, we dropped a factor of in the term; this factor can always be absorbed into the chiral phase .
 L. Chim, J. Phys. A 28, 7039 (1995).
 Of course in the extreme limit we have . Note also that the scaling forms provided here assume . For larger transverse fields mixing between other sectors matters, and hence becomes a function of three parameters—, and .
 We are grateful to David Clarke for pointing out to us this XY model mapping.
 S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
 I. P. McCulloch, J. Stat. Mech. P10014 (2007).
 D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013).
 Y. Bahri, R. Vosk, E. Altman, and A. Vishwanath, “Localization and topology protected quantum coherence at the edge of ’hot’ matter,” (2013), unpublished, arXiv:1307.4092 .