Parafermionic clock models and quantum resonance

TL;DR

This paper investigates the existence of zero-modes in $ Z_N$ parafermionic clock models, revealing conditions under which strong zero-modes can or cannot exist, especially considering prime versus non-prime N and resonance effects.

Contribution

It provides a detailed analysis of zero-mode degeneracies in parafermionic models, highlighting the impact of system size, prime versus non-prime N, and resonance points on zero-mode existence.

Findings

Degeneracies decay exponentially with system size for prime N.

Resonance points lead to band crossings that hinder strong zero-modes.

Non-prime N models generally lack strong zero-modes, except at special anti-resonant points.

Abstract

We explore the $\mathbb{Z}_{N}$ parafermionic clock-model generalisations of the p-wave Majorana wire model. In particular we examine whether zero-mode operators analogous to Majorana zero-modes can be found in these models when one introduces chiral parameters to break time reversal symmetry. The existence of such zero-modes implies $N$-fold degeneracies throughout the energy spectrum. We address the question directly through these degeneracies by characterising the entire energy spectrum using perturbation theory and exact diagonalisation. We find that when $N$ is prime, and the length $L$ of the wire is finite, the spectrum exhibits degeneracies up to a splitting that decays exponentially with system size, for generic values of the chiral parameters. However, at particular parameter values (resonance points), band crossings appear in the unperturbed spectrum that typically result in a splitting of the degeneracy at finite order. We find strong evidence that these preclude the existence of strong zero-modes for generic values of the chiral parameters. In particular we show that in the thermodynamic limit, the resonance points become dense in the chiral parameter space. When $N$ is not prime, the situation is qualitatively different, and degeneracies in the energy spectrum are split at finite order in perturbation theory for generic parameter values, even when the length of the wire $L$ is finite. Exceptions to these general findings can occur at special "anti-resonant" points. Here the evidence points to the existence of strong zero modes and, in the case of the achiral point of the the $N=4$ model, we are able to construct these modes exactly.

Figures

Equations

$$H_J$$

$$H_f$$

$${\sigma }_{i,j}={\omega }^{i-1}{\delta }_{i,j}{\tau }_{i,j}={\delta }_{i+1,j(\text{mod}N)},$$

$$\sigma=\left[\begin{array}[]{ccc}1&0&0\\ 0&\omega&0\\ 0&0&\omega^{2}\end{array}\right],\quad\tau=\left[\begin{array}[]{ccc}0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right]$$

$${\sigma }^N={\tau }^N=I$$

$${\tau }_x{\sigma }_y={\omega }^{{\delta }_{x,y}}{\sigma }_y{\tau }_x$$

$$Q=\prod _{i=1}^L{\tau }_x.$$

$${\gamma }_{2x-1}={\sigma }_x\prod _{j<x}{\tau }_j{\gamma }_{2x}={\omega }^{(N-1)/2}{\sigma }_x\prod _{j\le x}{\tau }_j$$

$${\gamma }_x{\gamma }_y={\omega }^{\mathrm{sgn}(y-x)}{\gamma }_y{\gamma }_x{\gamma }_x^N=1$$

$$H_J$$

$$H_f$$

$$Q={\omega }^{-L(N-1)/2}\prod _{x=1}^L{\gamma }_{2x-1}^{†}{\gamma }_{2x}^{}$$

$$Q|\text{s}{⟩}^S=|\text{s}-1{⟩}^S$$

$$|\text{s}{⟩}_q^S\sim \frac{1}{\sqrt{N}}\sum _{j=0}^{N-1}{\omega }^{-jq}{Q}^j|\text{s}{⟩}^S,$$

$$|\text{s}{⟩}_q^S\equiv \frac{1}{\sqrt{N}}\sum _{j=0}^{N-1}{\sigma }_L^q{Q}^j|\text{s}{⟩}^S.$$

$$|012{⟩}_0^D$$

$$|012{⟩}_1^D$$

$$|012{⟩}_2^D$$

$${\alpha }_x$$

$${\beta }_x$$

$${\tau }_x=\{\begin{array}{ll}{\alpha }_1& x=1\\ {\alpha }_{x-1}^{†}{\alpha }_x& 1<x<L\\ Q{\alpha }_{L-1}^{†}& x=L\end{array}$$

$$H_J$$

$$H_f$$

$$H_J^q$$

$$H_f^q$$

$$H_{f_e}^q$$

$$H_{f_b}^q$$

$$|\text{d}{⟩}_q^D=\prod _{x=1}^{L-1}({\alpha }_x^{†}{)}^{d_x}|\varnothing {⟩}_q$$

$$|012{⟩}_q^D={\alpha }_2^{†}{\alpha }_3^{†2}|\varnothing {⟩}_q$$

$$\eta |s⟩=|N-s⟩$$

$$\sigma \eta =\eta {\sigma }^{†}\tau \eta =\eta {\tau }^{†}$$

$$RH(\theta ,\varphi )=H(\theta ,-\varphi )R,$$

$$E=-2J\sum _{x=1}^{L-1}\cos (\theta +\frac{2\pi d_x}{N})=\sum _{x=1}^{L-1}{ϵ}_{d_x}$$

$$p=s_L-s_1=\sum _{x=1}^{L-1}{d}_x(\mathrm{mod}N)$$

Full text

There has recently been growing interest in a class of $\mathbb{Z}_{N}$ symmetric one dimensional lattice models known as parafermion chain or quantum clock models Fendley2012 ; Fendley2014 ; Jermyn2014 ; Cobanera2014 ; Alexandradinata2015 ; Iemini2016 . These models generalise the Kitaev wire model Kitaev2001 which exhibits localised unpaired Majorana zero-modes at each end. The recent surge of interest is inspired in part by proposals for their physical realisation and their potential application to universal topological quantum computationClarke2013 ; Alicea2015 ; Mong2014 , something which is not possible with Majorana zero-modes.

Clock-like systems of this type have been studied earlierHowes1983 ; Gehlen1985 , and much is known about the exactly solvable chiral Potts models that occur at special values of the coupling constants (see e.g. Refs Au-Yang1997, ; Perk2015, ). From the perspective of topological quantum computation however, there are a number of interesting open problems surrounding the topological degeneracies and potential zero-modes of the models. In particular, one may ask whether zero-modes exist which result in topological degeneracies throughout the energy spectrum Jermyn2014 ; Alexandradinata2015 . In this context one speaks of strong vs. weak zero-modes (see e.g. Ref. Alicea2016, ) and the question is important because degeneracies that exist at energies above the ground-state potentially allow for topologically fault-tolerant quantum devices at higher temperaturesHuse2013 . This area is of course also interesting on a fundamental level, as it addresses if and when decoupled/free zero-mode quasi-particles can exist in complex systems that are otherwise interacting.

We study this question here by direct perturbative and numerical calculations of the energy spectra of the models. This approach can be seen as complementary to the iterative position space constructions where (see Ref. Fendley2012, ; JackKemp2017, for example), a strong zero-mode on a chain of length $L$ is defined as an operator $\Psi$ such that $[H,\Psi]\rightarrow 0$ as $L\rightarrow\infty$, with finite-size corrections whose expectation values vanish exponentially with increasing system size. For a strong zero-mode it is also required that $\Psi^{N}$ is proportional to the identity operator for some integer $N$, and that $\Psi$ does not commute with some operator $Q$ which generates a discrete symmetry, so $[H,Q]=0$ but $[\Psi,Q]\neq 0$.

In the thermodynamic limit this definition implies topological degeneracies throughout the spectrum, as the action of $\Psi$ cycles through degenerate states with different eigenvalues with respect to $Q$. Although directly probing the energy spectrum lacks the spatial resolution of the aforementioned iterative approach, it has a number of advantages for this type of study. In particular, the ability to resolve the degeneracy at any energy scale means that it can be straightforward to identify weaker variants of the zero-mode. For example we will see that although there is often no strong zero-mode, the topological degeneracy is still preserved up to energies far above those of the ground state manifold. Similarly, we find scenarios where the degeneracy is preserved at all energies, but only between subsets of the discrete symmetry sectors. In this latter case it becomes straightforward, within these sub-spaces, to formally define hidden zero-modes that have a lower $\mathbb{Z}_{N}$ symmetry.

For the Kitaev wire, it is well known that the topological degeneracy is present throughout the spectrum. This model has strong Majorana zero-modes on the edges which provide each state with even fermion number a (nearly) degenerate partner at odd fermion number and vice versa. Of course the Kitaev wire is equivalent to an Ising spin chain through Jordan-Wigner transformation and its spectrum is exactly solvable.

The behaviour for generic $\mathbb{Z}_{N}$-parafermion models is more complicated and depends on the chiral angle $\theta$ (see Eq. 1). In Ref. Fendley2012, , it was noted that the $\mathbb{Z}_{3}$ model does not have strong zero-modes at the time reversal invariant

points ($\theta=n\frac{2\pi}{3}$ and $\theta=n\frac{2\pi}{3}+\frac{\pi}{3},\quad\forall n=0,1,2$). However an iteratively constructed perturbative expansion for an exact zero-mode operator was given for regimes where $\theta$ takes generic values (sufficiently far from the time reversal invariant points). The radius of convergence of this expansion was conjectured to be dependent on $\theta$.

In Ref. Jermyn2014, , the stability of zero-modes in the $\mathbb{Z}_{3}$ model was explored further using perturbation theory, diagonalisation of approximate hamiltonians, and DMRG. It was shown explicitly that the breakdown of zero-modes at the time reversal invariant points is due to domain-wall tunnelling processes which lead to power law splitting of the degeneracy between excited states. Evidence was then shown for the suppression of these processes for generic values of $\theta$, leading to restoration of apparent exponential splitting and “zero-mode revival”.

A central result of this manuscript is to show that there is in fact power law splitting of the topological degeneracy (and hence no strong zero-modes) in large regions of the model’s parameter space. This can be explained in terms of perturbative effects occurring at resonant crossings between bands labelled by different types and numbers of generalised domain-wall excitations. Sufficiently far from these crossings, the degeneracy may be recovered and this is observed for finite systems. However as the system size increases, so too does the number of resonant band crossings, spoiling the degeneracy over ever larger areas of parameter space. Our analysis suggests that in the limit $L\rightarrow\infty$, strong zero-modes do not exist for generic values of the chiral parameter $\theta$, at any value of $N$ (although there are at least isolated $\theta$-values that allow for strong zero-modes for $N=3$ and guaranteed exactly solvable points for $N=4$). This also implies that the radius of convergence of the expansion for the zero-mode given in Ref. Fendley2012, in $\theta$-space is generically zero.

To show the power law splitting, we use degenerate perturbation methods like those employed in Ref. Kells2015, for the interacting Kitaev chain. In a similar way we find evidence that the topological degeneracy is broken only at order $\sim L$, provided that the un-perturbed system does not contain any band crossings. However, because these crossings shift and become more dense as we make $L$ longer, we cannot say that strong zero-modes exist in the $L\rightarrow\infty$ limit. Therefore while topological degeneracies can persist at finite size, and also for many low lying energy levels in the large $L$ limit, they do not hold generally at finite energy density. It will be outlined elsewhere Kellsinprep why similar scaling applies to the interacting Kitaev chain.

We consider the model at any value of $N$ and note that there is an important difference between the cases where $N$ is prime and those where it is not. In particular we note that if $N$ is composite, certain bands with different number and type of domain-walls are degenerate for all values of the chiral parameter. This leads to a uniform power law splitting of the topological degeneracy and hence no $\mathbb{Z}_{N}$ parafermionic zero-modes occur for generic $\theta$ when $N$ is composite. Of course the actual splitting varies with $\theta$ and exceptional values of $\theta$ where the degeneracy is preserved can occur.

In the course of our study several special ‘anti-resonant’ cases arise. One of these occurs in the $N=4$ clock hamiltonian, as defined in Ref. Clarke2013, and closely related to the Ashkin-Teller model (see e.g. Ref. Kohmoto1981, ). We show that this model can be rewritten as a frustrated spin-$1/2$ ladder in such a way that the couplings along the rungs vanish precisely at the achiral resonant points. Here the system is exactly solvable and can be mapped to a $\mathbb{Z}_{4}$ time-reversal invariant Majorana model, i.e. two uncoupled Majorana chains. At these points, the model also displays exact $\mathbb{Z}_{4}$ parafermionic zero-modes which can be expressed in terms of the Majorana modes of the two decoupled chains.

Another interesting special case is the region around the $\theta=\pi/6$ chiral point of the $N=3$ model, where an apparent resonance is effectively cancelled by a precise match up of contributing domain-wall excitation energies. This region was pinpointed in earlier studies of the problemJermyn2014 as being a promising candidate for a strong zero-mode. Although our results show the special region around the $\pi/6$ point in parameter space again becomes vanishingly small in the limit $L\rightarrow\infty$, our analytic and numerical results suggest that strong zero-modes do exist at the point itself. We note furthermore that when the universal degeneracy breaks down in the nearby parameter space, it typically does so at very high energies.

An outline of the paper is as follows. In section I, we present the model and rewrite it using a domain-wall picture which is useful for the subsequent analysis. In section II we discuss the unperturbed system and characterise the resonant crossings between bands. In section III we analyse the effect of perturbations at the resonant crossings and also far away from any crossings. In section IV we discuss the results of our exact diagonalisation analysis and compare them to the aforementioned perturbative analysis. In section V the recovery of zero-modes at anti-resonant points, the exactly solvable nature of the $N=4$ model at achiral points and hidden zero-modes are discussed. Finally we wrap up with the conclusions and outlook in section VI.

I The Model

I.1 Clock Spins and Parafermions

We consider the $\mathbb{Z}_{N}$ quantum clock hamiltonian as given in Refs: Clarke2013, ; Jermyn2014, and written $H=H_{J}+H_{f}$, where

[TABLE]

$N$ is the number of ’hour-marks’ on the clock at each site. The operators $\sigma$ and $\tau$ are $N\times N$ matrices defined by

[TABLE]

with $\omega=e^{\frac{2i\pi}{N}}$. The operator $\sigma$ measures the clock value at a given site. The operator $\tau$ reduces the clock value on a given site by one. For $N=2$, $\sigma$ and $\tau$ are the $\sigma^{z}$ and $\sigma^{x}$ Pauli operators respectively. For $N=3$, they take the form

[TABLE]

Note that for $N=2$, this model is the transverse field Ising spin chain. (In this case the angles $\theta$ and $\phi$ are redundant and can be removed by dividing $J$ and $f$ by real factors $\cos(\theta)$ and $\cos(\phi)$). When $N>2$ and $\theta$ is not a multiple of $\frac{2\pi}{N}$, the model is called the chiral clock model since spatial-parity and time reversal symmetries are brokenFendley2012 . For $N=3$ this model is equivalent to the clock models analysed in Ref. Fendley2012, , while for higher $N$, it is a special case of the general $\mathbb{Z}_{N}$ models of in Ref. Fendley2012, . We consider the regime with $f\ll J$ which is the ordered phase. This can be mapped exactly to the topologically ordered phase of the corresponding parafermion chain model. The $\phi$ angle does not have much of an effect on many features and we will often just consider the $\phi=0$ case. There is one special feature at $\phi=0$, which is the presence of dihedral symmetry – see Sec. I.4

The operators $\tau$ and $\sigma$ are $N^{\text{th}}$ roots of the identity

[TABLE]

and satisfy

[TABLE]

Using this, one may see that the hamiltonian possesses a global symmetry: the operator $Q$ which moves the clock back one click at all sites,

[TABLE]

The eigenvalues of Q are $\omega^{q}$ with $q=0,1,...N-1$. We call the different eigenspaces of the $Q$ operator the $q$-sectors and label each with its respective $q$ value.

Using a non-local transformation due to Fradkin and KadanoffFradkin80 and analogous to the Jordan–Wigner transformation, we can rewrite the spin model in terms of parafermionic variables. Parafermionic operators are defined as follows,

[TABLE]

These satisfy the relations

[TABLE]

The factor of $\omega^{(N-1)/2}$ in $\gamma_{2x}$ is necessary for the second equation. One now finds that

[TABLE]

and

[TABLE]

In these variables, $Q$ is a topological quantum number analogous to the fermionic parity $(-1)^{F}$. In fact $\gamma_{x}Q=\omega Q\gamma_{x}$ and we see that $Q$ will commute only with operators constructed from products of $N$ parafermion operators. $Q$ is called $N$-ality, or sometimes “parafermionic parity”. When $f=0$, the hamiltonian does not contain the operators $\gamma_{1}$ and $\gamma_{N}$ and these operators are then strong edge zero-modes, giving rise to an $N$-fold degeneracy in the entire spectrum. Of course in the spin model we can understand this degeneracy as a trivial case of spontaneous breaking of the $\mathbb{Z}_{N}$ spin symmetry since there are no spin flip terms in the hamiltonian at $f=0$. Note that at $f=0$, the operators $\sigma_{x}$ will in fact act as zero-modes for all $x$. However, for $x$ away from the edges, these operators are non-local in the parafermion language. Also, as $f$ is taken to nonzero values, only the edge zero-modes have a chance to persist as strong zero-modes. In the rest of the paper, we will focus on determining to what extent the degeneracies of the model persist when $f\neq 0$.

I.2 $Q$-eigenstates

The hamiltonian (1) is written in the position space clock basis, where states are written $|\,\textbf{s}\rangle^{S}=|\,s_{1},s_{2},...,s_{L}\rangle^{S}$, and $s_{x}\in\{0,1,...,N-1\}$ is the clock value on site $x$, or in other words $\sigma_{x}|\,s_{1},s_{2},...,s_{L}\rangle^{S}=\omega^{s_{x}}|\,s_{1},s_{2},...,s_{L}\rangle^{S}$. The superscript $S$ explicitly indicates position space.

These states are not eigenstates of the $Q$ operator (5), but are super-positions of eigenstates from different $q$-sectors. The action of $Q$ on $|\,\textbf{s}\rangle^{S}$ changes the value of the clock at each site by $-1\ (\text{mod}\ N)$.

[TABLE]

For each state $|\,\textbf{s}\rangle^{S}$, we can write down a state in a given $q$-sector by performing a discrete Fourier transform as in

[TABLE]

which is an eigenstate of $Q$, with eigenvalue $\omega^{q}$. For what follows it is necessary to choose a convention for the global phase of these states. Here it will be convenient to fix the phase according to the value of the clock on the last site. The component of the state where the last site has its clock set to zero will be taken to have a positive real coefficient. Since $\sigma_{L}$ measures the clock value at site $L$, this can be achieved using the definition

[TABLE]

I.3 Domain-wall picture

We now introduce a domain-wall representation where one focuses on the differences between clock values on neighbouring sites, rather than the clock values themselves. This picture allows us to directly write down the hamiltonian in each $q$-sector. We find that by doing this, and as a consequence of the phase convention we chose in the previous section, there are only two terms in the hamiltonian which depend on the particular $q$-sector. Better still, both of these act at the right end of the chain. This makes the domain-wall representation very convenient for numerical calculations and a natural setting for perturbative expansions.

Each $q$-sector is spanned by $(\mathbb{C}^{N})^{\otimes(L-1)}$ states which we label $|\,\textbf{d}\rangle^{D}_{q}=|\,d_{1},d_{2},...,d_{L-1}\rangle$ where $d_{x}\in\{0,1,...,N-1\}$. $|\,\textbf{d}\rangle^{D}_{q}$ is a sum over position spin space states like in (10), where s satisfies $d_{x}=s_{x+1}-s_{x}\ (\text{mod}\ N)$, for $1\leq x\leq L-1$. We illustrate this with an explicit example for a state with $N=3$ and $L=4$ for all $q$.

[TABLE]

We now proceed to write the hamiltonian (1) in this representation. First we define operators $\alpha$ and $\beta$ given by

[TABLE]

These operators have the same matrix representations as operators $\tau$ and $\sigma$ respectively, except that they act on the domain-wall basis. In terms of these operators $\tau_{x}$ can be written as

[TABLE]

and the hamiltonian (1) is now given by

[TABLE]

Note that in this form all $q$-dependence appears only in $H_{f}$ on terms $-fe^{i\phi}Q\alpha_{L-1}^{\dagger}$ and $-fe^{-i\phi}\alpha_{L-1}Q^{\dagger}$ which act on the end of the chain.

Since $Q$ commutes with $H$ and with $\alpha^{\phantom{\dagger}}_{L-1}$ and $\alpha^{\dagger}_{L-1}$, this hamiltonian is block diagonal, and we can write the hamiltonian in each $q$-sector by replacing $Q$ with the appropriate eigenvalue,

[TABLE]

For later sections it will be useful to decompose $H_{f}$ further into terms that act on the edge and those containing bulk terms. This gives $H_{f}^{q}=H_{f_{e}}^{q}+H_{f_{b}}^{q}$ where

[TABLE]

We can now interpret $(\alpha^{\dagger}_{x})^{k}$ as creating a domain-wall of type $k$ at site $x$. For arbitrary $|\,\textbf{d}\rangle_{q}^{D}$ we may then write

[TABLE]

where $|\,\emptyset\rangle_{q}=|\,0,0...0\rangle_{q}^{D}$ is the state without domain-walls and with the appropriate eigenvalue of $Q$. As an example we write

[TABLE]

While this way of writing the states is reminiscent of second quantization, one should remember that the ”vacuum states” $|\,\emptyset\rangle_{q}^{D}$ are not annihilated by the domain-wall annihilation operators and also that $[\alpha^{\phantom{\dagger}}_{x},\alpha^{\dagger}_{y}]=0$, so the domain-wall creation operators are not the usual creation operators of bosonic or fermionic Fock space. There have been attempts to formulate a Fock space for parafermions, notably the work of Cobanera and Ortiz Cobanera2014 .

I.4 Dihedral symmetries

When the angle $\phi$ is a multiple of $\frac{2\pi}{N}$, the models we study here have an enlarged group of symmetries, isomorphic to $D_{N}$, the symmetry group of a regular polygon with $N$ sides. The $\mathbb{Z}_{N}$ generated by $Q$ is included in this $D_{N}$ as the group of rotations of the polygon. The extension of the symmetry group from $\mathbb{Z}_{N}$ to $D_{N}$ is not generally present in the models of Ref. Fendley2012, for $N>3$ (even at $\phi=0$). Since $D_{N}$ is a non-Abelian group, its presence causes degeneracies in the spectrum. In particular, as long as $\phi=0$, any eigenstate of $Q$ with eigenvalue $\omega^{q}$ such that $\omega^{q}\neq\omega^{-q}$ comes with a partner at the same energy whose $Q$ eigenvalue is $\omega^{-q}$.

We now present the $D_{N}$ symmetry explicitly111Our presentation of the $D_{N}$ symmetry group here is in accordance with the results given by Motruk, Berg and PollmannMotruk2013 and generalizes some of their findings.. First of all recall that the dihedral group $D_{N}$ is generated by two elements $a$ and $b$ subject to the relations $a^{N}=b^{2}=(ab)^{2}=e$ (where $e$ is the group identity). Here, $a$ can be viewed as a rotation over an angle $\frac{2\pi}{N}$ and $b$ as a reflection acting on the regular $N$-gon. The group has a total of $2N$ elements, the rotations $a^{j}$ and the reflections $a^{j}b$, with $j\in\{0,...,N-1\}$. We now introduce an operator $R$, such that $R^{2}=(QR)^{2}=I$, which together with $Q$ (from (5)) generates a $D_{N}$ symmetry. Next we define $\eta$ to be the operator that transforms each clock value as

[TABLE]

with $\eta|\,0\rangle=|\,0\rangle$. Note that $\eta$ corresponds to the complex conjugation of the eigenvalues of $\sigma$. We see that $\eta^{2}=1$ and that the following exchange relations with $\tau$ and $\sigma$ hold

[TABLE]

We can now define $R_{0}=\prod_{i}\eta_{i}$. This satisfies $(R_{0})^{2}=I$, $(R_{0}Q)^{2}=1$ and $R_{0}H(\theta,\phi)=H(-\theta,-\phi)R_{0}$. We note that $R_{0}$ and $Q$ generate a $D_{N}$ group together, but $R_{0}$ is only a symmetry if $\theta=\phi=0$. Combining $R_{0}$ with the operator $F$ which flips the clock state, that is, which switches the clock state at position $i$ on the chain with the clock state in position $L-i$, we obtain an operator $R=R_{0}F$ which is a symmetry for all $\theta$ as long as $\phi=0$. Explicitly, we define $F|\,s_{1},s_{2},...,s_{L}\rangle^{S}=|\,s_{L},...,s_{2},s_{1}\rangle^{S}$ and $R=R_{0}F$. We then find that $R^{2}=(RQ)^{2}=I$ and

[TABLE]

which in particular gives $[H(\theta,\phi=0),R]=0$. We see that we have $D_{N}$-symmetry at $\phi=0$ and in fact, $\mathbb{Z}_{2}\times D_{N}$ symmetry if $\phi=\theta=0$, as we have $F$ as an additional commuting symmetry when $\theta=0$.

To see how this leads to degeneracies, note that $QR=RQ^{-1}$ which means that if $|\,\psi\rangle$ is an eigenstate of $H$ and of $Q$ with eigenvalues $(E,\omega^{q})$ then $R|\,\psi\rangle$ is an eigenstate of $H$ and $Q$ with eigenvalues $(E,\omega^{-q})$. As a result all eigenstates of $Q$ with eigenvalues $\omega^{q}\neq\omega^{-q}$ are at least doubly degenerate. Only states with real eigenvalues of $Q$ have a chance of being nondegenerate. Generically only states within a given irreducible representation of $D_{N}$ will be degenerate (barring the existence of some further symmetry but we don’t observe this). Since $D_{N}$ has only one dimensional and two dimensional irreducible representations for all $N$ (see for instance Ref. james_and_liebeck, ) we see that the $D_{N}$ symmetry only causes doublets and not higher multiplets of degenerate states (or bands) to appear. We find that the $D_{N}$ doublets are split when $\phi\neq 0$, so it is not possible to extend the $D_{N}$ symmetry to this regime.

II Unperturbed system

The regime of interest for topologically non-trivial phases in these systems is where $f\ll J$. In this work we will explore this regime by performing perturbative expansions in $H_{f}^{q}$. Before doing this however, we will spend some time exploring the unperturbed system, when $f=0$ which will serve to guide our later analysis.

II.1 Unperturbed energy levels

The unperturbed system $H_{J}$, is diagonal in both representations discussed in the previous section. In the domain-wall picture, the energy of a state $|\,\textbf{d}\rangle_{q}^{D}$ has the simple form

[TABLE]

where we recall that $d_{x}=s_{x+1}-s_{x}$. We also define $\epsilon_{j}\equiv-2J\cos(\theta+\frac{2\pi j}{N})$ for compactness. There is no $q$ dependence in the unperturbed hamiltonian which means there is an exact $N$ fold degeneracy for the entire spectrum, consistent with the presence of the zero-modes at $f=0$.

Figure 1 shows the energy spectra of the unperturbed system plotted against $\theta$ for $N=2,3$ and $4$. Despite the fact that these spectra are for very small chains, they capture the important features.

There are degenerate bands of states, with each band labeled by the number of each type of domain wall present. In the absence of $H_{f}$, the domain walls do not disperse and hence there is a large degeneracy associated with different choice of the domain wall positions. A further useful label for the bands, which is also indicated with colour in the figures, is defined by

[TABLE]

We can think of $p$ as the total number of domain-walls modulo $N$. This quantity, which we call total domain-wall angle, is crucial to understanding the perturbative analysis discussed in later sections. The underlying symmetry $P$ of $H_{J}$ can be written as

[TABLE]

and has eigenvalues $\omega^{p}$.

At $\theta=0$, the $N$ fold degenerate ground state is ’ferromagnetic’ and spanned by states with all neighbouring clocks aligned. The ground state in each q-sector is written $|\,00...0\rangle^{D}_{q}$. At $\theta=\frac{2\pi k}{N}$, for nonzero integer $k$, the spectrum is identical, but ground states are spanned by states with all neighbouring clocks $k$ clock positions apart, written $|\,kk...k\rangle^{D}_{q}$. It is interesting to note that for $N$ odd, the spectrum is not symmetric around $E=0$, except at $\theta=\frac{k\pi}{2N}$, with $k$ odd, where the system also has its “superintegrable” point for $N=3$ Fendley2012 .

As noted, the energy of a band in the unperturbed system depends only on the number and type of domain-walls, but not their location. If we count the trivial domain walls where adjacent clock variables are equal, then every pair of neighbouring sites gives rise to a domain wall and we have $L-1$ domain walls of $N$ distinct types. The maximum number of unique energy levels for a chain of length $L$ is then the number of ways of distributing the $L-1$ domain-walls among $N$ numbered boxes. This is the number of compositions of $L-1$ with $N$ parts, denoted $c_{S}(L-1,N)$ Eger2013 . $S$ here is the set of possible values the parts can take, so here $S=\{0,...,L-1\}$. The $c_{S}(L-1,N)$ are extended binomial coefficients Eger2013 and satisfy the recurrence relation

[TABLE]

which can be used to calculate these efficiently for large systems.

Each composition is an ordered tuple of non-negative integers $(p_{0},p_{1},...,p_{N-1})$ that satisfy

[TABLE]

We use these tuples to label the bands of states of the unperturbed system. Here $p_{j}$ is the number of domain-walls of type $j$ in the band. The energy of a band in the unperturbed system described by $\vec{p}$ can be written

[TABLE]

II.2 Mapping of resonance points

We define a resonance point to be a place where two bands of the unperturbed system cross and become degenerate. At and in the vicinity of resonance points all crossing bands must be considered when performing perturbative expansions. It is then of interest to ask for which $\theta$ values these resonance points occur, and whether they are dense on the $\theta$-axis in the limit $L\rightarrow\infty$. In this section we derive expressions for the $\theta$ values at which resonance points occur and following this show that for $N>2$ they are dense on the $\theta$-axis.

Let us consider two bands, labelled by tuples $\vec{a}$ and $\vec{b}$ where $a_{j}$ gives the number of domain-walls of type $j$ for the first band, etc. We write $\vec{c}=\vec{a}-\vec{b}$, and note that $\vec{c}$ satisfies the following constraints:

[TABLE]

The quantity $n_{\vec{c}}$ is useful as it is related to the order in perturbation theory necessary to connect states in different bands at a resonance point.

From (27), we see that the bands labelled by $\vec{a}$ and $\vec{b}$, are degenerate precisely when

[TABLE]

Assuming that both sides of this equation are nonzero, we find that the bands are degenerate at $\theta$ values which satisfy

[TABLE]

and hence the bands cross twice over the full range of possible $\theta$ values. Alternatively, it may be that both sides of (29) are zero. In this case the bands are degenerate for all $\theta$ and we can write (still with $\omega=e^{2\pi i/N}$)

[TABLE]

When $N=3$, it is easy to see that this can never be satisfied with $\sum_{j}c_{j}=0$ and $\vec{a}\neq\vec{b}$. Hence there are no bands which are everywhere degenerate for $N=3$. However for $N=4$, there are nontrivial difference vectors $\vec{c}$ which satisfy equation (31) and the constraints given in Eq. (II.2); the vectors $\vec{c}=(p,-p,p,-p)$ give solutions for all integer $p$ such that $4|p|\leq 2(L-1)$. This leads to the appearance at $N=4$ of precise degeneracies between bands, at all $\theta$, that comprise of a different number and type of domain-walls. Figure 2 shows the distribution of the bands (in green) that are everywhere degenerate (e.g. bands with different domain-wall composition lie on top of each other everywhere) and the resonance points and for a $L=6$ system with $N=4$. A contrasting plot for $N=3$, Figure 3, shows that bands with different number and type of domain-walls typically have different energies. The exception in this case of course being at the resonant crossings. These examples reflect the general situation for any $N$, with everywhere-degenerate bands appearing always when $N$ is composite and never when $N$ is prime.

Fully characterising the degenerate bands comes down to a characterisation of all linear relations of $N^{\rm th}$ roots of unity with integer coefficients. This is an interesting problem in number theory – some recent results and references can be found in Ref. Steinberger2008, . The essential difference between prime and composite $N$ can be easily seen however. First of all, if $N$ is prime, then all nontrivial linear combinations of the $N^{\rm th}$ roots of unity with integer coefficients satisfying equation (31) have constant coefficients. This means $c_{j}=p$ for some $p\in\mathbb{Z}$ for all $j$ and we see that the constraint $\sum_{j}c_{j}=0$ is never satisfied. Hence when $N$ is prime there are never any bands that are degenerate for all $\theta$. We now show why this is the case for $N$ prime, but not for $N$ composite. For $N$ prime the cyclotomic polynomial $1+x+...+x^{N}$ is the minimal polynomial for the primitive roots $\omega^{j}$ (for $j\neq 0$). If there were a linear combination of the powers of $\omega$ with non-constant coefficients, such that $\sum_{j=0}^{N-1}d_{j}\omega^{j}=0$ then we could obtain a nontrivial linear combination which omits $\omega^{N-1}$, namely $\sum_{j=0}^{N-2}(d_{j}-d_{N-1})\omega^{j}=0$. However, this would show that $\omega$ is the root of the polynomial $\sum_{j=0}^{N-2}(d_{j}-d_{N-1})x^{j}$ which is of lower degree than the minimal polynomial leading to a contradiction. For composite $N$, we can always choose integer factors $e$ and $f$ of $N$, and find non-trivial coefficients $\{d_{j}\}$ which satisfy $\sum_{j=0}^{e-1}d_{j}\sum_{k=0}^{f-1}\omega^{ek+j}=0$ and $\sum_{j=0}^{e-1}d_{j}=0$. These satisfy the constraint that the sum of coefficients vanish and that these combinations all equal zero since $\sum_{k=0}^{f-1}\omega^{ek+j}=0$ forall $0\leq j\leq e$. Hence we have everywhere-degenerate bands for all composite $N$.

II.3 Resonance points are dense on the $\theta$ axis as $L\rightarrow\infty$

We can now deal with the question of whether resonance points occur arbitrarily near any particular value of $\theta$. This is clearly not the case for finite $L$, but as $L$ grows, ever more resonance points appear and in the limit $L\rightarrow\infty$ they do form a dense set on the $\theta$-axis, for all $N>2$. To see this we first note that the two solutions to Eq. (30) are in fact given by

[TABLE]

It is clear that we can make the linear combination $\sum_{j}c_{j}\omega^{j}$ take any fixed argument $\alpha$ to arbitrary accuracy if we are allowed to take the coefficients $c_{j}$ as large as we want. However a complication is once again the condition that $\sum_{j}c_{j}=0$. We can circumvent this as follows. First of all note that we can choose any two powers of $\omega$, for example $\omega$ and $\omega^{0}=1$ as a basis for the complex plane over the reals. We can then take a point $z$ with $\arg(z)=\alpha$ in the complex plane and write $z=a+b\omega$ with $a,b\in\mathbb{R}$. The coefficients $a,b\in\mathbb{R}$ can be arbitrarily well approximated by rationals $a^{\prime},b^{\prime}\in\mathbb{Q}$, in particular, for any chosen $\epsilon$, we can make sure that $|\arg(a^{\prime}+b^{\prime}\omega)-\alpha|<\epsilon$. We can now write $a^{\prime}=\frac{a_{1}}{a_{2}}$ and $b^{\prime}=\frac{b_{1}}{b_{2}}$ for integers $a_{1},a_{2},b_{1},b_{2}$ and consider $Z=Na_{2}b_{2}(\frac{a_{1}}{a_{2}}+\frac{b_{1}}{b_{2}}\omega)-(a_{1}b_{2}+a_{2}b_{1})(\sum_{j=0}^{N-1}\omega^{j})$. We note that $Z\neq 0$ and $Z$ is of the form $\sum_{j}c_{j}\omega^{j}$ with $c_{j}\in\mathbb{Z}$ and $\sum_{j}c_{j}=0$ and importantly that $\arg(Z)=\arg(a^{\prime}+b^{\prime}\omega)$. This means we can indeed get within any $\epsilon$ of our chosen argument $\alpha$ using a judicious choice of coefficients $c_{j}$. Hence the resonance points are dense on the $\theta$ axis in the limit $L\rightarrow\infty$.

As an illustration of how the $\theta$ axis is eventually covered by resonance points, we show all resonance points in an $N=3$ chain of length $L=30$ in Fig. 4. This figure shows that while the resonance points are eventually dense on the axis, there are notable gaps without any resonances even at large finite sizes, surrounding the most prominent resonances. Also higher order resonance points tend to appear at high energy, meaning that if only states with energy below a given bound are considered, there may be appreciable regions of the $\theta$-axis where no resonance effects are observed.

III Perturbation theory

In the last section we explored in detail the appearance and characterisation of resonance points in the unperturbed spectrum. The reason we are so interested in resonance points is that the behaviour of the perturbative series depend very much on whether one is on- or off-resonance. We show that at the off-resonant points $q$-sector dependent contributions only occur at an order of perturbation theory that scales with the length of the system. On the other hand, when we perform the perturbative expansions at resonant points (where bands of different domain-wall type and number cross), we typically see that there are $q$-sector dependent terms which split the degeneracy at an order that does not scale with $L$. In this section we examine these pertubative series at both resonant and off-resonant points.

Our approach is to employ the Raleigh-Schrodinger degenerate perturbative methods pioneered by Kato and Bloch Kato1949 ; Bloch1958 ; Bloch1958b ; Messiah1962 ; Lowdin1962 . In these approaches one generates effective hamiltonians acting within a degenerate band of the unperturbed system, which can be diagonalised to approximate the true eigenvalues of the full system. We outline the method in detail in appendix A. Here we simply state the effective hamiltonian $H^{\text{eff}(3)}_{q}$ for a Bloch expansion in the $Q=q$ sector for a band $b$, to order $f^{3}$

[TABLE]

In this expression, $P_{0}$ is the projector onto the band $b$, $Q_{0}=1-P_{0}$ is its complement, while $E_{0}$ is the unperturbed energy of the band, as given by Eq. (27).

Recall the perturbing hamiltonian in terms of $\alpha$ operators (16) is written as

[TABLE]

There are only two terms of this operator that actually depend on $q$, the terms proportional to $\alpha^{\phantom{\dagger}}_{L-1}$ and $\alpha^{\dagger}_{L-1}$. Only terms in the perturbative expansion which contain these terms can lead to an energy splitting between $q$-sectors. It is also important to note that both these terms break total domain-wall angle, defined in (24).

III.1 Off resonant $\theta$

The main argument we make is that, in the absence of resonance points, the degeneracy between energy eigenvalues for different q-sectors is exact to an order of the perturbing parameter that depends on the length of the wire. This is provided of course the respective perturbative expansions converge, which is only the case for sufficiently small perturbing parameter $f$ and sufficiently far away from resonant points.

Our strategy for the off-resonant scenario is as follows. We show explicitly that for all bands and all values of N there can be no q-dependence in perturbative corrections up to the third order. We do this explicitly in this section up to second order and discuss generally why this is the case for third order, with additional details provided in appendix B. For higher orders we have been unable to show explicitly the lack of q-dependence but argue why this should also be the case. We also provide numerical evidence that there is no such q-dependence up to eighth order for some finite systems, and show how these results agree to machine precision with exact diagonalisation calculations. For brevity some of the more technical aspects of these calculations have been moved to the appendix.

III.1.1 First order corrections

In the domain-wall picture it is very easy to show that there can be no $q$-dependent first order corrections for any band $b$. According to Eq. (III), to get the energy levels of the band to first order in $f$, we must diagonalise the operator

[TABLE]

However, it is not possible to connect two states from band $b$ via either of the two $q$-dependent terms appearing in $H^{q}_{f_{e}}$. These terms change the total domain-wall angle and thus connect to states in different bands, orthogonal to the original band $b$. Concretely, if $|\,i\rangle$ and $|\,j\rangle$ are two states in band $b$, then

[TABLE]

Hence all potentially $q$-dependent matrix elements of $H^{\text{eff}(1)}_{q}$ are actually zero, meaning $H^{\text{eff}(1)}_{q}$ is independent of $q$ and the eigenvalues of $H$ are independent of $q$ to first order in $f$.

III.1.2 Second order corrections

We now show that there is no $q$-dependence in the energy to second order. According to Eq. (III), to obtain the energy levels of any band $b$ to second order in $f$, we must diagonalise the operator

[TABLE]

We already know that $H^{\text{eff}(1)}_{q}$ does not depend on $q$, so we can concentrate on potentially $q$-dependent matrix elements of the second order part of $H^{\text{eff}(2)}_{q}$. Consider the matrix elements between two states $|\,i\rangle=|\,\textbf{d}^{(i)}\rangle_{q}^{D}$ and $|\,j\rangle=|\,\textbf{d}^{(j)}\rangle_{q}^{D}$, in band $b$. We can see from the fact that the total domain-wall angle of both states is the same, that all $q$-dependent contributions to the matrix elements of $H^{\text{eff}(2)}_{q}$ connecting these states must involve either $\alpha^{\dagger}_{L-1}$ and $\alpha^{\phantom{\dagger}}_{1}$, $\alpha^{\phantom{\dagger}}_{L-1}$ and $\alpha^{\dagger}_{1}$ or potentially $\alpha^{\dagger}_{L-1}$ and $\alpha^{\phantom{\dagger}}_{L-1}$. The contributions involving both $\alpha^{\dagger}_{L-1}$ and $\alpha^{\phantom{\dagger}}_{L-1}$ are independent of $q$ since the accompanying factors of $\omega^{q}$ and $\bar{\omega}^{q}$ cancel any $q$-dependence. The contributions from other terms add to zero which we will now show. We show this for $\alpha_{1}$ and $\alpha_{L-1}^{\dagger}$, but the same is true for $\alpha_{1}^{\dagger}$ and $\alpha_{L-1}$ (since this is equivalent to exchanging $|\,i\rangle$ and $|\,j\rangle$).

The contributions involved are

[TABLE]

where $|\,k\rangle=|\,\textbf{d}^{(k)}\rangle_{q}^{D}$ and $|\,l\rangle=|\,\textbf{d}^{(l)}\rangle_{q}^{D}$ run through the eigenstates of $H_{0}$, which do not lie in band $b$. $E_{k}$ and $E_{l}$ are the unperturbed energies of $|\,k\rangle$ and $|\,l\rangle$ respectively. Both the sum over $k$ and the sum over $l$ contain at most a single nonzero term. These nonzero terms occur in pairs where the energy denominators have equal magnitude but opposite signs, leading to exact cancellation (since all nonzero matrix elements of $\alpha_{x}$ are equal to $1$).

In cases where there are nonzero terms, all the bulk domain-walls in states $|\,i\rangle,|\,j\rangle,|\,k\rangle$ and $|\,l\rangle$ must be equal. This means $d^{(i)}_{x}=d^{(j)}_{x}=d_{x}^{(k)}=d_{x}^{(l)}$ for all $x\in\{2,...,L-2\}$. For the domain-wall excitations on the ends, the following conditions must hold

[TABLE]

It then follows that

[TABLE]

and

[TABLE]

Summing these we see that $E_{0}-E_{k}=-(E_{0}-E_{l})$ and the two nonzero matrix elements between $\langle i\,|$ and $|\,j\rangle$ cancel.

III.1.3 Third and higher order corrections

As one might expect, dealing with third and higher order corrections is challenging due to the ever increasing combinations and permutations of terms that must be considered. In these cases it is instructive to perform the perturbative expansions numerically. From these calculations, we observe an energy splitting between $q$-sectors in the third order corrections at generic values of $\theta$. However this splitting is incrementally reduced by corrections at higher orders. This leads us to conjecture that in the absence of resonance points and for sufficiently small $f$, the maximum splitting possible scales as $f^{L}$. We proceed by first showing some results from these numerical calculations and then discussing why we might intuitively expect this behaviour. In appendix B we show that when using the alternative perturbative expansion of SoliverezSoliverez1969 there are no longer any $q$-dependent third order corrections. We also demonstrate how cancellations between terms in the perturbative series at third and higher orders can occur.

We use the perturbative expansions of Bloch in our numerical calculations (see appendix A for details). While all expansions are equivalent in the large order limit, Bloch’s expansions contain fewer terms, making them simpler and more efficient to compute. At modest system size, we can compare the perturbative results to the exact energy spectra obtained by numerical diagonalisation. For example, figure 5 shows numerical results for an $N=3$ system at non-resonant values of $\theta$. It is clear that the results from the numerical perturbation theory calculations agree with the exact results up to corrections which decrease with the order $n$, whenever they are not too small to be numerically resolved. In this particular example the corrections appear to scale as $f^{n+1}$. Similar convergent behaviour is observed for all other bands. At and in the vicinity of resonance points though it is necessary to start from a subspace containing all bands which cross at the resonance point.

Figure 6 shows the maximal difference between the perturbative energy estimates for states in the $q=0$ and $q=1$ sectors of the first excited band of the same $N=3$ system. We note that the difference is nonzero from the third order onwards, but decreases with increasing order $n\geq 3$ and is always small enough to be explained as an error of the approximation. That is, the energy splitting between the $q=0$ and $q=1$ sectors, calculated from $H^{\text{eff}(n)}$, is less than $f^{n}$. In general we expect that the maximum difference between $q$-dependent matrix elements of the effective hamiltonians decay as $f^{n}$. How these matrix elements effect the energy splitting between $q$-sectors is not completely clear. However we may argue that these corrections should in general decay faster than $f^{n}$, since the number of $q$-dependent matrix elements is relatively small and the dimension of the effective hamiltonians grows with $L$. This behaviour is observed for all bands as long as $\theta$ is not resonant and $f$ is sufficiently small. For short chains, when $n\geq L$ a non-vanishing splitting between $q$-sectors can be observed which is always less than $f^{L}$.

Showing this behaviour directly in the degenerate perturbation series for all cases is a formidable task. Appendix B shows some specific examples and arguments in this direction. However we can argue intuitively why we expect this to be true. For a process to lead to $q$-dependent energy splitting it must include at least one $q$-dependent term acting on the end of the chain. Since both these terms change the total domain-wall angle, such a process must also include a term acting at the start of the chain to connect back to the original band. This means that at orders lower than $L$ any process that contributes will contain disconnected sets of operators acting on opposite ends of the chain. We may intuitively expect that contributions to the energy arising from operators acting on spatially disconnected parts of the chain should vanish. General arguments in this direction can be made starting from extensive energy scaling: If a system consists of independent subsystems the total energy should equal the sum of the energies of the subsystems and one would expect cross terms in the perturbation series to vanish, or at least the energy per site due to such terms should vanish in the thermodynamic limit $L\rightarrow\infty$ (see e.g. Ref. Shavitt2009, ). These considerations are closely related to the existence of linked-cluster-theorems for degenerate systems where non-local terms corresponding to un-linked diagrams cancelLowdin1951 ; Brandow1967 ; Tang2013 . In section III.1.4 we provide numerical evidence that this is indeed the case by excluding disconnected terms from the pertubative expansions.

Further general arguments that point in this direction can be made based on the work of Fendley in Ref. Fendley2012, . Here for the $\mathbb{Z}_{3}$ model at generic values of $\theta$ an asymptotic expansion of a parafermionic zero-mode operator in powers of $f$ (or $f/J$) is well defined.222Here, “generic values” refers to a condition requiring the invertibility of a set of $\theta$-dependent matrices which appear in the construction of terms in the expansion. We expect that this is equivalent to the requirement that $\theta$ is not a resonance point. For finite systems, the arguments given by Fendley should indeed imply that the splitting between $q$-sectors vanishes to order $f^{L}$ at generic $\theta$ and for small enough values of $f$.

III.1.4 Excluding disconnected terms

It is possible to perform the perturbative expansions numerically in such a way that processes with operators acting on both ends are explicitly excluded. While this is not a well defined pertubative expansion, we find that amazingly it agrees exactly with the exact results up to an error which decreases with the order of the perturbative expansion $n$, exactly as with the regular perturbative expansion (figure 5). By construction there is no $q$-dependence appearing at any order. Of course since this method cannot account for processes containing operators acting on both ends, we cannot observe the splitting due to processes connecting both ends of the system. This is not an issue though when the order $n$ is less than the system size $L$.

The method proceeds by calculating three different effective hamiltonians using different modifications to the perturbing hamiltonian. We label these effective hamiltonians $H^{\text{left}(n)}_{q}$, $H^{\text{right}(n)}_{q}$ and $H^{\text{neither}(n)}_{q}$ where the modifications consist of excluding terms acting only on the leftmost site, only on the rightmost site and only on the leftmost or the rightmost sites respectively. Once we have these effective hamiltonians we combine them to get the final effective hamiltonian as

[TABLE]

III.2 Resonant $\theta$

We now turn our attention to the resonance points themselves and see that the behaviour here is qualitatively different to that found at off-resonant points. At resonant points, bands of states with different domain-wall configurations become degenerate making it possible for local processes to split the degeneracy between $q$-sectors. We begin by looking at first order processes, particularly those connecting the two lowest lying bands at $\theta=0$ for the $N=3$ case. Here we write down explicit expressions for the matrix elements of the effective hamiltonian. We can then diagonalise these effective hamiltonians analytically to get expressions for the energy splitting between $q$-sectors. We then look at other higher order resonance points and using numerical perturbative expansions we pick out the order at which the splitting occurs. We find that the splittings observed are consistent with the exact diagonalisation results discussed in section IV and with the characterisation of resonance points given in section II.2. We then show that not all resonance points lead to an energy splitting between $q$-sectors. Particular examples of this are the resonance points at $\theta=0$ for $N=4$ and at $\theta=\frac{\pi}{6}$ for $N=3$. Both of these are discussed more in section V. Note that for convenience, we set $\phi=0$ for most of the rest of this section.

III.2.1 First order resonance points

We first discuss resonance points where the splitting between $q$-sectors is first order in $f$. This can only happen when states from two bands that cross can be connected directly by the $q$-dependent terms in $H^{q}_{f}$, i.e. by $f\omega^{q}\alpha^{\dagger}_{L-1}$ and its conjugate. Since these terms can only change a single domain-wall, it is necessary that there is a degeneracy between consecutive single domain-wall energies, $\epsilon_{d}=\epsilon_{d\pm 1}$ for some $d$. In Fig. 7, we illustrate where degeneracies of domain-wall energies occur for systems at $\theta=0$ and $\theta=\frac{\pi}{N}$ for $3\leq N\leq 6$. We note that degeneracies between consecutive domain-wall types occur for all $N$ at $\theta=\frac{\pi}{N}$, whereas at $\theta=0$ we see that adjacent degenerate wall types occur only for $N$ odd. More generally making use of our earlier characterisation of resonance points in section II.2, we observe that first order processes can only take place at resonance points where $n_{\vec{c}}=2$. There is an additional constraint that the two non-zero elements of $\vec{c}$ must appear consecutively, where we consider $\vec{c}$ to be periodic. From plots 3 and 2 we see that these only appear when $\theta$ is a multiple of $\frac{\pi}{N}$.

We now look in detail at the resonance point between the first and second excited bands of the $N=3$ model at $\theta=0$. These are the bands just above the ground state in figure 1. Each of these bands contains a single domain-wall excitation.

A useful shorthand for what follows will be to employ a labelling of the states where only domain-walls with $d\neq 0$ and their locations are indicated, e.g.

[TABLE]

Using this shorthand, we label states in these bands $|\,1_{x}\rangle_{q}=\alpha_{x}^{\dagger\;1}|\,\emptyset\rangle_{q}$ and $|\,2_{x}\rangle=\alpha_{x}^{\dagger\;2}|\,\emptyset\rangle_{q}$. The only $q$-dependent first order processes that connect these states occur at the end. The relevant matrix elements are (setting $\phi=0$)

[TABLE]

Ordering the basis elements as:

[TABLE]

we see that the effective hamiltonian for the two bands can be written as a $2L-2$ Toeplitz matrix:

[TABLE]

which can be diagonalised with a generalised Discrete-Fourier-Transform:333We give the results below for general $N$ because these are also the energy splittings for resonances between the lowest bands at $\theta=\frac{\pi}{N}$, for arbitrary $N$.

[TABLE]

where $L^{\prime}=2L-2$. This gives a spectrum of the form

[TABLE]

where $k\in[1,...,L^{\prime}]$. The difference in energies between two states $\Delta E=E(k_{1})_{q_{1}}-E(k_{2})_{q_{2}}$ is therefore

[TABLE]

The relevant comparison to make for the $q$-dependent splitting is between states where $k_{1}=\pm k_{2}$, as these states would have the same energy when $\Delta q=q_{1}-q_{2}$ is taken to zero. In these cases, for large $L$ we find that $\Delta E\propto\frac{f}{L^{2}}$ if $k_{1}\ll L$ while $\Delta E\propto\frac{f}{L}$ if $k_{1}\approx L/2$. We also note that there is no splitting if $\Delta q=N$. In this case the states $|\,k_{1}\rangle_{q_{1}}$ and $|\,-k_{1}\rangle_{q_{2}}$ remain exactly degenerate as a consequence of the dihedral symmetry (which is present since we set $\phi=0$). At $N=3$, this degeneracy is observed when $q_{1}=1,q_{2}=2$. Of course it is also clear that there can be no splitting between sectors in this case as the effective matrix (37) for $q=1$ is the Hermitian conjugate of the effective matrix for $q=2$.

For higher energy bands at $\theta=0$, there are many more bands crossing and each of these bands has a larger numbers of states (we may for example consider bands with states $|\,1_{x},1_{y}\rangle$, $|\,1_{x},2_{y}\rangle$, $|\,2_{x},2_{y}\rangle$). Many of these bands can still be connected via first order processes like the ones we just demonstrated. The result is that the splitting between $q$-sectors here also scales with $f$. The factor will in general be smaller than in (40), since the effective matrices will have larger dimension and thus the $q$-dependent matrix elements have less impact on the eigenvalues.

Similar results for this first order crossing where found by Jermyn et. al. in ref. Jermyn2014, , although using their notation the effective matrices are not written directly in the each $q$-sector and by projecting to sectors with constant domain-wall number, it is not possible to observe higher order resonance points at $\theta\neq 0$.

III.2.2 Higher order resonance points

Resonance points can also introduce $q$-dependent splitting at higher orders in $f$. These can occur between bands which differ only at a single domain wall, e.g. the bands containing states $|\,d_{x}\rangle$ and $|\,(d+m)_{x}\rangle$ are in resonance if $\epsilon_{d}(\theta)=\epsilon_{d\pm m}(\theta)$ which can can potentially cause splitting at order $l=\min(m,N-m)$. However, most higher order resonances are not of the simple type just described. Instead they result from combinations of domain-wall energies that become degenerate. Using the notation for bands that was introduced in section II.2, we describe a resonance point between two bands $\vec{a}$ and $\vec{b}$ by a vector $\vec{c}=\vec{a}-\vec{b}$. We then write $n_{\vec{c}}$ for the Hamming distance between $\vec{a}$ and $\vec{b}$. If the bands are in resonance, then the order of perturbation theory at which $q$-dependent splitting appears, $m_{\mathrm{split}}$, must satisfy

[TABLE]

To see this, note that the application of $H_{f_{e}}^{q}$ to the band with domain-wall vector $\vec{a}$ changes the Hamming distance of the state to $\vec{b}$ by $2$, while application of $H_{f_{b}}^{q}$ changes the Hamming distance by at most $4$. At least one application of $H_{f_{e}}^{q}$ is necessary to introduce $q$-dependence. Hence to get from $\vec{a}$ to $\vec{b}$ we need to operate with $H_{q}^{f}$ at least $1+(n_{\vec{c}}-2)/4$ times. This gives a lower bound for the order at which $q$-dependent splitting can occur at such a resonance point.

As a specific example consider the resonance point between the fifth and sixth excited bands of the $N=3$ model. This resonance is marked in figure 3 and has $n_{\vec{c}}=6$. The lower (upper) band from the point of view of increasing $\theta$ contains states of type $|\,1_{x},1_{y}\rangle$ ($|\,2_{x}2_{y}2_{z}\rangle$). Using the notation introduced in section II.2, these are characterised by the tuples $\vec{a}=(L-3,2,0)$ and $\vec{b}=(L-4,0,3)$ respectively. The domain-wall energies at this resonance point satisfy $\epsilon_{0}+2\epsilon_{1}=3\epsilon_{2}$. Using Eq. (30), we find the exact $\theta$ for this crossing is $\theta=\arctan(\frac{\sqrt{3}}{5})$. From the domain-wall content of these bands, we see that the lowest order on which these can in principle be connected with $H_{f}^{q}$ is the second order. However, the $q$-dependent second order terms cancel and we find $q$-dependent splitting at third order. Figure 8 shows results of numerical perturbative expansions at this point for a range of $f$ values. We see that there is splitting in the energy between $q$-sectors appearing at third order. Unlike the behaviour at the off-resonant points (figure 6) where the spltting is reduced at higher orders, here it remains. This agrees with observations from the exact numerics shown in figure 9 and discussed in section IV.

IV Full diagonalisation of small systems

For small systems we can examine the maximum splitting between $q$-sectors over a range of values of $\theta$ and $f$ by fully diagonalising the hamiltonian. We can do this directly in each $q$-sector using the domain-wall picture discussed in section I.3. Figure 9 shows a plot of the maximum splittings between the $q=0$ and $q=1$ sectors for an $N=3$ system with $L=9$ (note that the $q=1$ and $q=2$ sectors are exactly degenerate due to the dihedral symmetry described in section I.4). In each $q$-sector the energy levels have been ordered and we have plotted the maximal absolute difference between eigenvalues at the same position in this ordering.

It is clear from figure 9 that the largest splitting between sectors occurs around the $\theta=0$ point. This can be attributed to the first order processes occurring at resonance points here, which are discussed in detail in section III.2. For larger $f$, we see other branches of splittings emerge. These can be attributed to third, fourth and fifth order processes occurring at the resonance points indicated in figure 3. While there are many resonance points at $\theta=\frac{\pi}{6}$, these do not result in any splitting between the sectors. This is because the bands crossing here have the same total domain-wall angle. This is examined in more detail in section V.1. Away from resonance points, we see that there is very little splitting (especially at low $f$). This is consistent with the perturbation theory results for off-resonant $\theta$ in section III.1. Some splitting does appear at larger $f$ even away from resonances, but this can be ascribed to finite size effects which allow for $q$-dependent perturbative processes at order $f^{L}$.

V Anti-resonance and hidden zero-modes in the $\mathbb{Z}_{3}$ and $\mathbb{Z}_{4}$ models

In previous sections we outlined why for $N$-prime, finite chain length $L$, sufficiently small $f$ and at generic values of the chiral parameter $\theta$ we observe that the energy splitting between $q$-sectors decays exponentially as the length of the chain is increased. We were also able to show that the regions where this $f^{L}$ splitting is observed become vanishingly small as $L$ gets larger, and hence there is typically no strong zero-mode in the thermodynamic limit. We also showed that strong zero-modes do not typically exist anywhere if $N$ is composite, owing to $q$-dependent terms appearing at very low orders between bands that are everywhere degenerate in the un-perturbed limit.

In this section we outline counter-examples where these general observations do not hold. The first case is for the region around $\theta=\frac{\pi}{6}$ of the $N=3$ model which was previously noted in Ref. Jermyn2014, as being the most likely region to support strong zero-modes. We show below that, despite the presence of prominent resonance points at $\theta=\frac{\pi}{6}$ (see Figure 1), these do not result in any $q$-dependent power-law energy splitting. This fact, coupled to the knowledge that $N$ in this case is prime, allows us to argue that a region $\delta\theta$ about this point is also degenerate to order $f^{L}$, but that the size of this stable region decreases as $\delta\theta\sim 1/\sqrt{3}L$. We go on to show that there are in fact additional values of $\theta$ where the resonance points do not result in any $q$-dependent power-law splitting. However these are less prominent and have a smaller surrounding stable region than for $\theta=\frac{\pi}{6}$.

The second case is the $\theta=\frac{n\pi}{2}$ points of the $N=4$ model. These $\theta$ values are also prominent resonant points (see Figure 1), but are actually exactly solvable; mappable to two uncoupled Majorana chains. Although the case for strong-modes in this scenario is non-perturbative, because the $N=4$ system is composite, we actually see that the degeneracy away from this point is broken at low orders of $f$ and so the parafermionic strong zero-modes can only exist exactly at the special point i.e $\delta\theta=0$. We make the case however that there exist hidden topological zero-modes between some of the q-sectors and write down what these modes look like in position space using an iterative approach similar to that employed in Ref. Fendley2012, .

V.1 The $\mathbb{Z}_{3}$ system in the vicinity of $\theta=\pi/6$

The $\theta=\frac{\pi}{6}$ point at $N=3$ is unique for odd $N$ in that the single domain-wall energies at this point are equally spaced, see Fig. 10. This property is quite powerful as it is one of the main ingredients in showing the exact solvability and $N$-fold degeneracy of a class of $\mathbb{Z}_{N}$ models studied in Ref. Gehlen1985, . In particular, the so called super-integrable point of the $N=3$ model also occurs at $\theta=\frac{\pi}{6}$, but in addition has $\phi=\frac{\pi}{6}$.

On a perturbative level this precise match-up between domain-wall excitation energies is also very important because it brings together states (in energy) that can only be connected by total domain-wall angle (see (24)) preserving terms. This means that the quadratic bulk terms $H_{f_{b}}^{q}$ in (19) can quite easily map between states in the two bands (even on the first order) but that terms that break this total domain-wall angle (those from $H_{f_{e}}^{q}$ from (19)) must occur in total domain-wall angle preserving pairs. As a consequence the qualitative behaviour of the PT at this point is identical to that of the off-resonant PT discussed in section III.1.

As additional evidence for this claim we note that the effective hamiltonians generated within the degenerate subspace are $q$-independent up to the third order, when using Soliverez’s Soliverez1969 symmetrized perturbative expansion. Further numerical evidence is provided in Figure 11 where we note that, identical to the off-resonant cases, the $q$-dependency becomes progressively smaller (in a manner consistent with a topological degeneracy) as we include more orders in the expansion.

For finite $L$, there is a substantial region around the $\theta=\frac{\pi}{6}$ point which can contain exact zero-modes, since the resonance points which break the degeneracy only approach $\theta=\frac{\pi}{6}$ slowly with increasing $L$, see figures 4 and 9. The resonances near $\theta=\frac{\pi}{6}$ also occur at high energy (in the middle of the spectrum) and at high orders in $f$. It is therefor likely that even in the thermodynamic limit, there will be many states, even states with finite energy density, for which the topological degeneracy is preserved to extremely good approximation. We should also mention here that the perturbation theory and exact diagonalisation give the same qualitative results regardless of value chosen for the other chiral parameter $\phi$. This suggests that strong zero-modes exist all along the $\theta=\frac{\pi}{6}$ line. Note that their existence is not dependent on the super-integrability at $\theta=\pi/6$ as this occurs only at the point $(\theta,\phi)=(\frac{\pi}{6},\frac{\pi}{6})$.

We can find other so-called anti-resonance points where bands with the same total domain-wall angle cross and thus no $q$-dependent splitting results. This is illustrated in figure 4 where the anti-resonance points are shown in red. These points though are not as prominent as the point at $\theta=\frac{\pi}{6}$ because: 1) the bands are not as well separated and thus perturbative arguments are only valid for very small $f$ and 2) the distance to surrounding resonance points is much smaller meaning the stable region in the vicinity of these points is also much smaller.

To find the locations of the anti-resonances analyticallly, we first write the single domain wall energies $\epsilon_{i}$ for $i>0$ in the form $\epsilon_{i}=\epsilon_{0}+m_{i}\delta$ where $\{m_{i}\}$ are integers and $\delta$ is a real number. Using equation (27) we find that at a resonance point described by vector $\vec{c}$ we have $\sum_{i}c_{i}\epsilon_{i}=0$ which implies that

[TABLE]

The total domain-wall angle (24) of a band described by vector $\vec{a}$ is $p_{\vec{a}}=\sum_{i}ia_{i}\ (\text{mod}\ N)$. At resonance points where both bands have the same domain-wall angle

[TABLE]

must hold. While (42) and (43) are valid for all $N$ when $N=3$, they lead to the straightforward condition that $m_{1}+m_{2}=0\ (\text{mod}\ 3)$. Using the expression for the single domain-wall energies from (23) we find that for a given $\{m_{i}\}$ satisfying this constraint the resonance point can be found at $\theta=\arctan\left(\frac{\sqrt{3}(m_{1}-m_{2})}{m_{1}+m_{2}}\right)$ (assuming the system is sufficiently large).

For example, the simplest solution $(m_{1},m_{2})=(2,1)$ corresponds to the $\theta=\frac{\pi}{6}$ point. Another solution $(m_{1},m_{2})=(5,4)$ gives $\theta=\arctan\left(\frac{1}{3\sqrt{3}}\right)$ which corresponds to the line of $n_{\vec{c}}=10$ resonance points on the left side of figure 3. From the exact numerics shown in figure 9 we see that there is indeed no apparent splitting in the vicinity of this point, while there is for the other line of $n_{\vec{c}}=10$ resonance points on the right side of the figure.

V.2 The $\mathbb{Z}_{4}$ system

Let us first show how the $N=4$ model can be written as a spin-ladder that decouples into two spin-$\frac{1}{2}$ chains when $\theta=\frac{n\pi}{2}$ and $\phi=m\pi$ for integer $n,m$. At these special points the resulting model can be solved exactly (for all $f$). To see how this comes about, we first introduce the local unitary transformation $\xi$ which permutes the third and fourth clock states. We write this as an operator on $\mathbb{C}^{4}=\mathbb{C}_{u}^{2}\otimes\mathbb{C}_{d}^{2}$ as follows

[TABLE]

where we have introduced two copies of the Pauli matrices (up and down). The operator $\xi$ satisfies the following properties,

[TABLE]

Thus after applying the unitary operator $\Xi=\prod_{i=1}^{L}\xi_{i}$, the $N=4$ hamiltonian takes the following form,

[TABLE]

Setting $\theta=\frac{n\pi}{2}$ and $\phi=m\pi$, we see that many of the terms vanish and we end up with two decoupled chains. For $\theta=n\pi$, these are transverse Ising chains on the legs of the ladder, whereas for $\theta=n\pi+\frac{\pi}{2}$ they are zig-zag chains with alternating ferromagnetic and anti-ferromagnetic couplings. This structure is shown in figure 12.

When the system can be written as decoupled spin-$\frac{1}{2}$ chains (for $\theta=\frac{n\pi}{2}$ and $\phi=m\pi$), there is an exact four fold degeneracy and we expect there to be exact parafermionic zero-modes. Here we explore the case where $\theta=0$ and write expressions for these in terms of the fermionic zero-modes of the decoupled chains. Analogous constructions are possible for other special points. Starting with the zero-modes for the two Ising chains the resulting parafermionic zero-modes:

[TABLE]

where the operators $\alpha^{u/d}_{L/R}$ are Majorana zero-modes acting at the left and right hand edge of the upper and lower Ising chains. The operator

[TABLE]

exchanges the states between upper and lower chains and satisfies

[TABLE]

Using the anti-commutation relations of the fermionic zero-modes (note that zero-modes on different chains commute), and the fact that

[TABLE]

where $P_{d}=\prod_{i}\sigma_{i,d}^{x}$ and $P_{u}=\prod_{i}\sigma_{i,u}^{x}$ are the $\mathbb{Z}_{2}$ symmetry operators (or fermionic parity operators) of the upper and lower chains, one may check directly that the $\alpha_{L}$ and $\alpha_{R}$ given above satisfy the correct parafermionic relations (6), in this case $\alpha_{L}^{4}=\alpha_{R}^{4}=1$ and $\alpha_{L}\alpha_{R}=-i\alpha_{R}\alpha_{L}$.

V.2.1 Stability and hidden zero-modes

As $N=4$ is composite, there are bands which are everywhere degenerate and this should allow them to be split in a $q$-dependent fashion for generic $\theta$. The special points above are an exception to this, but once we move away from these we expect to rapidly loose the strong zero-mode protection, see section II. This is clear from figure 13, where we show the maximal energy splitting between $q$-sectors for all energy levels as a function of $\theta$, for a small finite system. Away from the $\theta=0$ there is generically splitting between all $q$-sectors which differ by $1$ or $3$. Figure 14 shows the splitting as a function of $f$ at each order of perturbation theory for two such everywhere degenerate bands of a seven site system. This figure clearly demonstrates that there are third order processes which split the degeneracy between $q$-sectors and that this splitting does not disappear as higher order terms are added.

For $q$-sectors which differ by $2$, there is no splitting observed at generic $\theta$. The absence of splitting between the $q=1$ and $q=3$ sectors is explained by the dihedral symmetry discussed in section I.4 (which is present since we set $\phi=0$). However between the $q=0$ and $q=2$ sectors there is splitting around the most prominent resonance points (cf. Fig. 2), but behaviour consistent with exact degeneracy away from these points.

The fact that that the sectors with $\Delta q=2$ are not generally split points to the survival of order two zero-mode operators related to $\alpha_{L}^{2}$ and $\alpha_{R}^{2}$. To see this note that whenever we have a $\mathbb{Z}_{N}$ parafermionic zero-mode such that $N$ is even, then the operator $\Phi=\Psi^{(N/2)}$ is a zero-mode with $\Phi^{2}=1$ and $\Phi Q=-Q\Phi$. Considering the two sides of the chain, we also note that $\Phi_{L}\Phi_{R}=(-1)^{N/2}\Phi_{R}\Phi_{L}$, so $\Phi$ is a fermionic mode if $N/2$ is odd. If $N/2$ is even, as is the case here then $\Phi$ is bosonic, but its existence nevertheless requires order two degeneracy throughout the spectrum, as it anti-commutes with $Q$. At $\theta=0$ we can immediately write $\Phi_{L}=(\alpha_{L})^{2}=\alpha_{L}^{u}\alpha_{L}^{d}$ and $\Phi_{R}=$ for any $f$. Away from $\theta=0$, the parafermionic zero-modes no longer exist but the mode $\Phi$ may survive. We can attempt to construct it using an iterative technique similar to that used in Ref. Fendley2014, which yields a power series expansion in the parameter $f/J$. We write

[TABLE]

and require that

[TABLE]

where $H_{0}=H_{J}/J$ and $V=H_{f}/f$. We can now determine $\Phi^{(n)}$ from $\Phi^{(n-1)}$, determining the entire series from $\Phi^{(0)}$, by requiring that $[\Phi^{(m)},H_{0}]=-[\Phi^{(m-1)},V]$. Taking $\Phi^{(0)}=\sigma^{z}_{1,u}\sigma^{z}_{1,d}$ (which is clearly correct at $\theta=0$), we obtain, to first order in $\frac{f}{J}$,

[TABLE]

The number of terms on the right hand side grows fast with the order of approximation. For example $\Phi^{(2)}$ involves terms with products of the operators $\sigma^{x}_{u/d,1/2}$ and $\sigma^{z}_{u/d,1/2/3}$ and we spare the reader the explicit expression. We do not know whether the operator $\Phi$ is generically normalizable. However, it is clear that even at first order, there are singularities in the expansion. We note that the $\cos(2\theta)$ in the numerator of $\Phi^{(1)}$ becomes zero precisely at the resonance point $\theta=\pi/4$ where there is indeed $q$-dependent splitting of the degeneracy at first order in $f$ (this is prominently visible in Fig.13). We expect that, similarly to the parafermionic mode at $N=3$, the zero-mode $\Phi$ at $N=4$ has an expansion to arbitrarily high orders at generic $\theta$, which may work well in finite systems but which will have a radius of convergence equal to zero at all $\theta$ in the $L\rightarrow\infty$ limit.

VI Conclusions and outlook

In this work we looked into the presence of strong zero-modes in $\mathbb{Z}_{N}$ parafermionic chain models. We did this through a detailed study of topological degeneracies in the entire energy spectra of these models, a necessary condition for strong zero-modes. When $N>2$ is prime, and in particular when $N=3$, the introduction of a chiral parameter $\theta$ can lead to regions where the necessary degeneracy exists is of order $f^{L}$. However, in the thermodynamic limit, resonance points where bands of the unperturbed model cross become dense on the $\theta$-axis and break the degeneracy at a lower order, at least in states which are nearby in energy. However not all resonance points result in splitting of the energy and we find so-called ’anti-resonance’ points where no such splitting occurs. These are found where bands with the same total domain-wall angle cross. A prominent example of this for $N=3$ is at $\theta=\frac{\pi}{6}$. While there are many other ’anti-resonance’ points for $N=3$, the point at $\theta=\frac{\pi}{6}$ is particularly interesting because the energy difference between bands here remains of order $J$ as $L\to\infty$. This suggests that strong zero-modes could persist at this point. It would be interesting to further explore the existence and structure of these modes using iterative methods similar to those employed in Refs. Fendley2012, ; JackKemp2017, and numerical methods similar to those discussed in Ref. Kells2015b, . For $N=4$ ’anti-resonance’ points are found at the achiral points. Here the model is exactly solvable and expressions for the parafermionic zero-modes were derived.

Our approach includes the introduction of a basis of domain-wall states which are also eigenstates of the topological symmetry $Q$. This enables us to write the hamiltonian directly in each $q$-sector and to isolate the $q$-dependent terms on the end of the chain. We used this to show directly that there can be no first or second order processes splitting the degeneracy between $q$-sectors at off-resonant points for general $N$, and none at third order for $N=3$. Using numerical perturbation theory methods we were able to show for accessible prime $N>2$ (notably $N=3$) that at off-resonant points, there is no splitting between $q$-sectors on orders less than the chain length. We were also able to numerically pick out the order at which the degeneracy is split at particular resonance points, obtaining results consistent with a lower bound for the order of splitting, given in Eq. (41). We showed that by excluding the unlinked processes responsible for $q$-dependent splitting in the pertubative expansions, the perturbative approximations still converge to the exact values. While it is not clear a priori that this is a valid perturbative expansion, it nevertheless shows that the unlinked processes cancel to the order of approximation that we have calculated, and we expect this to continue at higher orders.

We also showed that when $N$ is not prime, there is qualitatively different behaviour resulting from pairs of bands which remain degenerate for all values of $\theta$. This causes energy splitting at fixed order over the entire $\theta$ axis (bar some exceptional points). Hence there are generically no strong parafermionic zero-modes when $N$ is not prime. However at special values of $N$ and $\theta$, strong parafermionic zero-modes do occur. In particular this is the case at $N=4$ and $\theta=\frac{m\pi}{2}$. These points are special as the model can be rewritten as a pair of Ising chains and solved exactly for all $f$ and $J$. The $N=4$ model also has interesting features away from the special $\theta$ values. For example, despite its lack of parafermionic zero-modes it does feature a bosonic zero-mode which again has convergence issues at resonance points, but only at those resonance points which involve pairs of bands which are not everywhere degenerate.

Another aspect of this work is the analysis of symmetries from section I.4. When $\theta$ is a multiple of $\frac{\pi}{N}$ and $\phi=0$, the models have both a conjugation (time reversal) and a spatial parity symmetry. For general $\theta$, but still at $\phi=0$, these are no longer symmetries individually but their product is still a symmetry, which does not commute with $Q$. This gives us the non-Abelian dihedral group as a symmetry. From the representation theory of this group we can determine that (at $\phi=0$), the $q$ and $N-q$ sectors are exactly degenerate for each $N$.

A natural direction for further work is to characterize the energy scales below which degeneracies can persist in the thermodynamic limit as a function of $\theta$. A look at Fig. 4 which shows the resonance points for an $L=30$ chain shows clearly that the lowest energy at which a resonance point can be found varies considerably with $\theta$. In fact even at this relatively large size, there are still considerable windows without any resonance points around $\theta=0$, $\theta=\frac{\pi}{6}$, and a number of other prominent resonance points. We have shown that the window around the $\theta=\frac{\pi}{6}$ point closes as $\delta\theta\sim\frac{1}{\sqrt{3}L}$ and that the order of the processes that can cause splitting at the closest resonance point increases linearly with $L$. It would be interesting to see if this is a particular feature of the $\theta=\frac{\pi}{6}$ point or if similar behaviour is observed around other anti-resonance points.

Another natural question is how widely the influence of a particular resonance extends along the $\theta$-axis. In other words, how far from the resonance point can we expect appreciable $q$-dependent splitting and how does this length scale on the $\theta$-axis scale with $L$? In Ref. Jermyn2014, the scaling behaviour of the energy splitting for the lowest states of the first excited band around $\theta=0$ was explored numerically. It was found that the width of the region of influence around this particular resonance point decays rapidly with $L$ which suggests that similar behaviour might be observed at other resonance points.

Appendix A Degenerate perturbation theory formalism

In this appendix we outline some details of the Raleigh-Schrödinger degenerate perturbative expansions that are used in this work. For further details and full derivations, consult Refs Kato1949, ; Bloch1958, ; Bloch1958b, ; Messiah1962, ; Lowdin1962, .

We consider hamiltonians of the form $H=H_{0}+\lambda V$, where $H_{0}$ is the unperturbed part, and $\lambda V$ the perturbation. We focus on an eigenspace of of $H_{0}$ with energy $E_{0}$. The projector onto this eigenspace is written $P_{0}$ and we write $Q_{0}=1-P_{0}$ for its complement. We also define the operator $a=E_{0}-H_{0}$ which gives the difference between the unperturbed energy of the chosen eigenspace and the unperturbed energy of the state on which we apply it.

We use the version of the formalism developed by BlochBloch1958 . Here, the eigenvalues of $H$ which reduce to $E_{0}$ as $\lambda\rightarrow 0$ are approximated to order $\lambda^{n}$ by the eigenvalues of an effective hamiltonian

[TABLE]

Note that this acts non-trivially only on the chosen eigenspace of $H_{0}$. The $U^{(j)}$ are given by $U^{(0)}=P_{0}$ and for $j>0$,

[TABLE]

where $\sum_{\vec{k}}^{(j)}$ is a sum over all sets of non-negative integers $k_{1},k_{2},...,k_{j}$ satisfying the conditions:

[TABLE]

and the $S^{k}$ are given by

[TABLE]

All these expressions are slightly modified from earlier work by KatoKato1949 who showed that the expansion converges absolutely when $||\lambda V||<\Delta_{0}(E_{0})/2$, where $\Delta_{0}(E_{0})$ is the distance from $E_{0}$ to the nearest eigenvalue $E\neq E_{0}$ of $H_{0}$.

The $\mathcal{U}^{(j)}$ can be efficiently calculated using $\mathcal{U}^{(0)}=P_{0}$ and the recurrence relation

[TABLE]

To fourth order, the explicit effective hamiltonian is

[TABLE]

Up to order $\lambda^{n}$, the eigenvectors of $H^{\text{eff}(n)}$ in the $E_{0}$ eigenspace of $H_{0}$ are the projections of the corresponding eigenvectors of $H$ to this eigenspace. While the eigenvalues of $H$ are orthogonal, their projection onto the $E_{0}$ eigenspace of $H_{0}$ need not be. This is reflected by the fact that $H^{\text{eff}(n)}$ is typically not Hermitian for $n\geq 3$.

Appendix B Vanishing of third and higher order corrections in off-resonant perturbation theory

Here we look at the third and higher order pertubative expansions at off-resonant points. We first talk about how at orders $n\geq 2$ the choice of expansion used becomes important. Following this we show a specific example of a third order processes that leads to $q$-dependent matrix elements when using Bloch’s expansion. Next we show that any $q$-dependent third order processes supported at non-overlapping positions cancel. Finally we consider the third order processes at overlapping positions and show that using Soliverez’s expansion these cancel in the specific case mentioned earlier, and using numerics cancel in general for $N=3$.

B.1 Different pertubative expansions

An additional complicating issue, which appears in degenerate perturbation theory, is that there is no canonical choice of perturbation series. In fact there are multiple perturbation series which converge to the same limit (when they converge), but which can differ from each other at any finite order by higher order contributions. For example, in our numerical examinations we use the expansion by BlochBloch1958 which involves diagonalising effective hamiltonians which are completely $q$-independent to second order in $f$ but which do depend on $q$ at order $f^{3}$ and onward. On the other hand, the effective hamiltonians in the expansion by KatoKato1949 (of which Bloch’s is a modification) actually have $q$-dependent terms already at order $f^{2}$, despite the fact that the two expansions converge to the same limit. Another expansion by SoliverezSoliverez1969 also has eigenvalues which converge to the same limit, but has effective hamiltonians which are independent of $q$ to order $f^{3}$. This means in particular that there is no $q$-dependent energy splitting up to order $f^{3}$, despite the fact that the Kato and Bloch effective hamiltonians already depend on $q$ at this order. Unfortunately the hamiltonians in Soliverez’ expansion also become $q$-dependent from the next order up. While it may be possible to devise a perturbation scheme tailor made for this problem which has an explicitly $q$-independent effective hamiltonian up to order $f^{L-1}$ we will not pursue this here.

B.2 Specific example

First of all, let us show explicitly that the effective hamiltonian $H^{\text{eff}(3)}$ has nonzero $q$-dependent matrix elements. For small $\theta\neq 0$ we consider the first band above the ground state consisting of states $|\,1_{x}\rangle_{q}$. Within this band, $H^{\text{eff}(3)}$ connects the states $|\,1_{2}\rangle$ and $|\,1_{L-1}\rangle$, which host domain-walls on the second and last sites, respectively. The terms in $H^{\text{eff}(3)}$ which connect these states must involve the operators $\alpha_{1}^{\phantom{\dagger}}$, $\alpha^{\dagger}_{1}\alpha_{2}^{\phantom{\dagger}}$ and $\alpha^{\dagger}_{L-1}$ from $H^{q}_{f}$. These can potentially appear in the third order terms of $H^{\text{eff}(3)}$ in six possible permutations. Applying each of these permutations we find, up to multiplication by an overall factor of $\omega^{q}f^{3}e^{3i\phi}$,

[TABLE]

We now see that the processes on the third and fourth lines here do not actually lead to terms in $H^{\text{eff}(3)}$. This is because $|\,1_{1}\rangle$ is in the same band $b$ as the the initial and final states and is thus projected out by the $Q_{0}$ operator. In terms of the domain-wall energies $\epsilon_{d}$ in (23) we then find

[TABLE]

which shows that $H^{\text{eff}(3)}$ has at least one element that depends on $q$. We will comment on more general analytical calculation of $q$-dependent matrix elements of effective hamiltonians in the next subsection.

B.3 General arguments for the absence of splitting

All of the perturbative schemes mentioned in section B.1 involve the evaluation of the matrix elements of a characteristic set of operators at each order. For example, at third order, the following operators play a role

[TABLE]

At order $n$, we have operators $O_{p_{1}...p_{n-1}}$ containing $n-1$ factors $(E_{0}-H_{0})^{-p_{i}}$ with $p_{i}\geq 0$ and $\sum{p_{i}}=n-1$. The matrix elements of these operators are built from the various terms in $H^{q}_{f}$. Each occurrence of $H^{q}_{f}$ may be replaced by either $\alpha_{1}$, $\alpha_{L-1}$ or $\alpha^{\dagger}_{x}\alpha_{x+1}$ or their complex conjugates, with the appropriate prefactor from Eq. (33). We can ignore the prefactors for now and note that the $\alpha_{x}$ have matrix elements which can only take the values [math] and $1$.

B.3.1 Non-overlaping positions

If we take three terms from $H^{q}_{f}$ which are supported at non-overlapping positions (for example $\alpha_{1}$, $\alpha^{\dagger}_{L-1}$ and $\alpha^{\dagger}_{x}\alpha_{x+1}$ with $2\leq x\leq L-3$) then it is easy to calculate the corresponding contributions to the matrix elements of all three operators above. Applying each of the three operators causes a change in the energy, let’s call these energy differences $a$, $b$ and $c$, so for instance we may have the situation that applying $\alpha_{1}$ to a state with energy $E$ gives a state with energy $E+a$, subsequently applying $\alpha^{\dagger}_{x}\alpha_{x+1}$ to the state will give an energy $E+a+b$ and finally applying $\alpha^{\dagger}_{L-1}$ gives energy $E+a+b+c$. We can only get a nonzero matrix element if the final state is in the same band as the initial state, which means we must require $a+b+c=0$. Assuming that $a\neq 0$ and $a+b\neq 0$ this set of operators will only give a contribution to $O_{11}$ and this will be equal to $\frac{1}{a(a+b)}$ up to the appropriate factors from Eq. (33). However the triple of operators corresponding to the energy changes by $a$, $b$ and $c$ can occur in $6$ different orders. If $a$,$b$ and $c$ (and hence also $a+b$, $a+c$ and $b+c$) are all nonzero, there are $6$ corresponding contributions to any matrix element of $O_{11}$. These all cancel, since

[TABLE]

which vanishes, since we know that $a+b+c=0$. Similar identities make many contributions to the matrix elements of the operators $O_{11...1}$ which occur at higher orders vanish.

If one of the energy differences $a$, $b$ or $c$ equals zero, the contribution to the corresponding matrix element of $O_{11}$ no longer vanishes, but can in principle be cancelled by contributions to the matrix elements of $O_{02}$ or $O_{20}$. For example if $a=0$ then necessarily also $b+c=0$ and there are only two ways to produce a non-vanishing contribution to $O_{11}$, by adding energy according to the following steps,

[TABLE]

These give a total contribution proportional to $\frac{2}{b^{2}}$ (or equally $\frac{2}{c^{2}}$). The other orders of adding the energy now yield contributions to $O_{02}$ and $O_{20}$. We find that these contributions are also proportional to $\frac{2}{b^{2}}$ for $O_{02}$ (from $a+b+c$ and $a+c+b$) and again proportional to $\frac{2}{b^{2}}$ for $O_{20}$ (from $b+c+a$ and $c+b+a$), all with the same proportionality constants from Eq. (33). Similarly, if $b=0$ or $c=0$, we find that all operators get contributions proportional to $\frac{2}{a^{2}}$. It is clear that the contributions from $O_{11}$, $O_{02}$ and $O_{20}$ can all cancel each other, depending on how these operators appear in the effective hamiltonian. The Bloch hamiltonian $H^{\text{eff}(3)}$ contains the combination $O_{11}-O_{20}$ and does not contain $O_{02}$. We see that the contributions we have been considering here all vanish. The same is true for Soliverez’ third order hamiltonian,

[TABLE]

Again, the contributions from the $O_{p_{1}p_{2}}$ cancel.

B.3.2 Overlaping positions

We see that, in order to have non-vanishing contributions, we must take operators from $H^{q}_{f}$ which have overlapping support, that is, at least two of the operators form a “chain”. An examples is the triple $\alpha_{1}$, $\alpha^{\dagger}_{1}\alpha^{\phantom{\dagger}}_{2}$, $\alpha_{L}$ which we analysed in section B. In such cases the energy differences which appear depend on the order of application of the operators and we showed already that $\langle 1_{2}\,|H^{\text{eff}(3)}|\,1_{L-1}\rangle=(\epsilon_{0}-\epsilon_{1})^{-2}\neq 0$. We also noted during the calculation that in this case $\langle 1_{2}\,|O_{20}|\,1_{L-1}\rangle=0$. We can also calculate $\langle 1_{2}\,|O_{02}|\,1_{L-1}\rangle$, which appears in Soliverez’ hamiltonian. We see that the third and fourth permutations shown in Eq. (B.2) give

[TABLE]

The contribution from $O_{02}$ to Soliverez’ hamiltonian thus precisely cancels the contribution from $O_{11}$ and this potentially $q$-dependent matrix element is actually zero. One may similarly show that all further potentially $q$-dependent matrix elements of Soliverez’ third order hamiltonian are zero, in all bands (assuming $L\geq 3$). This shows in particular that there is no $Q$-dependent energy splitting in any band at order $f^{3}$, as long as $\theta$ is not at a resonance point.

We should note that all results we have presented about the absence of $Q$-dependent energy splitting at first, second and third order in $f$ are still valid if the coefficients $f$ and $J$ are allowed to vary along the chain (assuming they remain bounded away from zero). This is easy to see in retrospect. All results on vanishing were obtained by showing cancellation of contributions from sets of operators taken from of $H^{q}_{f}$ acting in different orders but at fixed positions along the chain. The cancellations are never dependent on any relationship between coupling constants at different sites or links of the chain.