Kaleidoscope of topological phases with multiple Majorana species

Exactly solvable lattice models for spins and non-interacting fermions provide fascinating examples of topological phases, some of them exhibiting the localized Majorana fermions that feature in proposals for topological quantum computing. The Chern invariant $\nu$ is one important characterization of such phases. Here we look at the square-octagon variant of Kitaev's honeycomb model. It maps to spinful paired fermions and enjoys a rich phase diagram featuring distinct abelian and nonabelian phases with $\nu= 0,\pm1,\pm2,\pm3$ and $ \pm4$. The $\nu=\pm1 $ and $\nu=\pm3$ phases all support localized Majorana modes and are examples of Ising and $SU(2)_2$ anyon theories respectively.

The study of topological quantum states has, in recent years, been boosted by proposals [1][2][3] to use nonabelian anyons for topologically protected quantum computing. Nonabelian anyons in two-dimensional systems were first proposed in the context of strongly correlated systems, with the Moore-Read Fractional Quantum Hall state [4] being the most famous example. Recent years have seen the advent of nonabelian anyons in effectively non-interacting systems. Proposals now exist for their realization in spinless chiral p-wave superconductors [5] and spin lattices [6], at interfaces between an s-wave superconductors and a topological insulator [7] or half metal [8], and on 1-d semiconductor-superconductor heterostructures [9][10][11]. In all of the above examples the nonabelian quasi-particles are realized by isolated, topologically protected Majorana fermions.
A necessary condition for isolated Majorana fermions is particle-hole symmetry of the type found in Bogoliubov-de Gennes (BdG) equations. A general understanding of where to find these conditions has been aided by the general symmetry classification of non-interacting gapped bulk systems [12][13][14][15]. In particular, the symmetry class of 2d chiral p− and f −wave superconductors (class D in [14]) is topologically non-trivial and can host topological bulk and edge-modes. The topological invariant characterizing the ground state in 2d is the Chern number ν, which maps to the set of all integers Z. However,in the chiral p-wave case, for example, only Chern numbers ν = 0 (abelian phase) and ν = ±1 (nonabelian phase) are uncovered from the full set Z. This is to be contrasted with the symmetry class of the Integer quantum Hall effect (class A), that lacks all symmetries. Here the Z invariant is fully realized (in theory) by filling some arbitrary number n of Landau level bands, which will be related to a Chern number ν = ±n [16].
One of our motivations for searching for higher Chern numbers in systems with BdG symmetry is that phases with different odd ν host Majorana fermions of different topological properties. The Chern number dependence of the nonabelian braiding properties is discussed in Kitaev's seminal paper [6]. A summary of the 16-fold way of bulk anyon models -so called because the braiding properties are determined by the value of ν modulo 16-is given in tables 1, 2 and 3 in that pa-per. The cases with odd ν (but not those with even ν) all support localized Majorana modes at the vortices and thus have nonabelian exchange statistics [17]. We will see that the lattice models studied in this paper realizes 9 of these 16 bulk orders, including 4 of the cases with nonabelian anyons, in two parity-doublets ν = ±1 and ν = ±3. The only other example of Chern number ν = ±3 that we are aware of occurs in the work of Mao and co-workers [18], which concerns an f -wave paired system of electrons moving in continuous spacetime. The Chern numbers ν = ±1 and ν = ±3 correspond to anyon models of the Ising and SU (2) 2 type respectively. However, the edge behavior of the system is not fully determined by this classification and will depend on the actual value of ν (not its residue modulo 16). Nevertheless, by counting gapless edge modes for the various phases, we argue that the ν = ±1 edges are described by chiral Ising CFTs (Conformal Field Theories), while the ν = ±3 edges support chiral SU (2) 2 CFTs.
The model we study is the square-octagon (4-8) lattice model discussed by Yang et al. [19]. Like the triangledodecagon (3-12, "Yao-Kivelson") lattice model [21] it is one of the possible generalizations in two dimensions [19] of the original Kitaev honeycomb lattice model [6,[22][23][24][25][26][27][28][29][30]. These exactly solvable 2d Ising-like spin models on trivalent lattices have turned out to be an ideal testing ground for topological properties predicted by the general symmetry arguments. The honeycomb model maps to spinless chiral p-wave system, and exhibits an abelian (ν = 0) as well as nonabelian phases (ν = ±1). Yang and coworkers observed [19] that the squareoctagon variation has two abelian ν = 0 phases as well as the ν = ±1 phases which open up near the critical boundary between the abelian ones. In none of the mentioned models were ground state (vortex-free) sectors with higher Chern numbers encountered. However, since the square-octagon model has a larger unit cell than the model on the honeycomb lattice and can therefore be mapped to a fermionic paired system with an internal pseudo-spin degree of freedom, we suspected that it may support the rich variety of phases found in spin-triplet paired systems, see for example Ref. 20.
We will see below that the the model supports an abundance of abelian and non-abelian phases with Chern num-bers 0, ±1, ±2, ±3, ±4 see Figure 3. Using the fermionization method introduced in Ref. 29, we have also checked directly that this model can be mapped to a fermionic paired system with pairings of type s-, p-, d-, f-g-etc. along with spin-orbit couplings of corresponding winding numbers. In the absence of spin-orbit coupling-for example in spinless pairs-the correspondence between the Chern number and the angular momentum of the pairing is rather simple. A spinless chiral pairing of p-wave type corresponds in the weak-pairing phase to ν = ±1, a chiral f-wave pairing instead to ν = ±3 etc. (In the strong-pairing phase one would in both cases find µ = 0.) The fact that we find phases with Chern numbers ν = ±3 among others concurs generally with the f-wave pairing observed in Ref. 18. However, we stress that in the model considered by us-which as mentioned contains intricate couplings of pseudospin and momentum with a variety of winding numbers-we have not been able to discern any simple relation between dominant pairing and the Chern number.
The outline of the paper is as follows. A general introduction to the spin model will be accompanied by a description of the fermionic solutions, focusing in particular on the vortex free sector, which contains the ground state. We then describe the phase diagram and discuss how the edge spectrum, and indeed the bulk spectrum at phase boundaries reflect the Chern number difference between topological phases. We also provide details of how the contributions of each band in the model combine to give the total Chern number of each phase and describe how these contributions change as the bands cross or touch. A general description of the edge spectra that occur between domains of the various topological phases (including vacuum) is then provided and we give an example where two domains with ν = 0 nevertheless have counter-propagating modes at the interface. We also show that these two abelian ν = 0 phases can be described with toric code like hamiltonians, supported on different sub-lattice structures.
The basic hamiltonian of the Kitaev-type trivalent spin model on the square-octagon (4-8) lattice is defined as where it should be understood that the z-links connect separate squares of the lattice and the x and y-links within the squares are the horizontal and vertical slopes respectively, see Figure 1. We define a basic unit cell of the lattice around the four sites that make up a single square and label the sites inside such a cell by n ∈ 1, 2, 3, 4. Each spin site on the lattice can thus be specified using the position vector q of the unit cell and the index n. A unit cell and fundamental lattice vectors are shown in Figure 1.
Within the model, as in honeycomb and 3-12 lattice models, there are closed loop symmetries that we can generate with overlapping link operators. The simplest of these are the octagonal W (8) and square W (4) plaquette operators. These are defined pictorially in Figure 1. The fact that the hamiltonian commutes with all plaquette operators implies that we may choose energy eigenvectors | n such that W then we say that the state | n carries a square vortex at q. By the vortex-sector we mean the subspace of the system Hilbert space with a particular configuration of vortices. The vortexfree sector for example is the subspace spanned by all eigenvectors such that W = −1 for all q, was shown in [31] to have a Fermi-surface and is thus a spin metal, like certain sectors of the 3-12 lattice variation [32].
In this type of system, it is possible to break time reversal symmetry explicitly, and still preserve the solvable nature of the system, by adding into the hamiltonian three site interactions of the form σ α i σ β j σ γ k . To find these terms we take overlapping products of the adjacent links of the original hamiltonian. These terms are shown in Figure 2. Crucial to obtain-FIG. 2: Time-reversal symmetry breaking terms. There are 8 octagonal (O) terms and 4 square (Q) terms. We construct each term using the overlapping products of two adjacent σ α σ α terms. To fix the sign of each term, we use a clockwise ordering for the adjacent link operators. The resulting signs are noted inside the plaquettes. The O terms are added to the hamiltonian with an overall coupling constant 2κ1 and the Q terms with an overall coupling κ2 ing the rich phase structure shown in Figure 3 is the sign we place in front of each term. Following [19] we use the convention that the individual links are arranged in a clockwise order around each plaquette, which yields the signs shown in the figure. Two overall coupling constants κ 1 for the terms on square plaquettes and κ 2 for the terms on octagonal plaquettes allow us to vary the strength and nature of these T -breaking terms.
To solve the model exactly one converts each vortex sector of the model to a quadratic fermionic system. There are a number of different techniques that do this, each with their own distinct advantages. For the purposes of calculating the fermionic spectra it is easiest to use the Majorana representation and conventions for this model used in [19,31]. It is important to note however that this methodology means we artificially enlarge the Hilbert-space and thus should project back to the physical subspace in order to calculate wavefunctions [6]. Indeed, for the purpose of calculating wavefunctions it is arguably easier to use a Jordan-Wigner type fermionization procedure [27,28]. We have elsewhere given details on a representation that uses a special Jordan-Wigner type string to resolve the fermionic vacuum as Toric code stabilized states [29]. Ground states calculated with this methodology are BCS ground state wave functions over these Toric Code vacua and the full set of excited states can be described in terms of BdG quasi-particles over these ground states.
We now focus on the vortex free sector, using the Majorana fermionic representation of ref. 19. After this fermionization, the Hamiltonian can be written as follows where A † k = (a k,1 , a k,2 , a k,3 , a k,4 ) and the operators a k,n are the Fourier transformed Majorana fermion creation/annihilation operators a k,n = 1 √ N a q,n e ik·q . The indices n label the four sites within a unit cell (see Figure 1), q is the center-of-mass coordinate of the unit cell, and a q,n is the operator which creates or annihilates a Majorana fermion at the site labeled by (q, n). The 4 × 4 matrix M can be de- The Chern invariant can be calculated as where P is the projection onto the negative energy eigenvectors of M . For multi-banded systems one can calculate a Chern integer for each band by using the same formula, but with P replaced by the projection to just that band. Each negative energy band then contributes its own integer to the total Chern invariant. The Chern integers of the completely positive energy bands are just −1 times the analogues for the negative energy bands and therefore contain the same information.
A change in the Chern number (5) indicates a phase transition (the reverse is not always true as Figure 4 illustrates). Changes in Chern number can occur when the positive and negative energy bands touch at one or several points. We can thus explore the phase diagram of the model by mapping out the Chern numbers and the gap as a function of the parameters (J x , J y , J z , κ 1 , κ 2 ). Clearly we can scale one of these parameters to 1 (as long as it was non-zero to start with) and this means we have a 4-dimensional parameter space to deal with. We will present the full structure of the 4D phase diagram in a separate paper [33]. Here, we focus on a two-dimensional slice through parameter space, setting J z = 1, J = J x = J y and κ 1 = κ 2 . In Fig. 3 we show the phases and Chern numbers in this slice.
The information in Figure 3 suggests that we can think of the topological phase transitions as an exchange of the Chern integers between positive and negative energy eigenmodes. Band touching within the negative energy bands also comes with a transfer of Chern integers, but here the total Chern invariant is unchanged. In the model in question the two negative energy bands touch whenever κ = 0, κ = ±J z /2 or κ = ±J/2. These internal 'transitions' are also indicated in Figure 3.
The dispersion relation for the energy E kx,ky of the model can be obtained exactly, but is complicated and we will not write it down here. However, the critical lines shown in Figure 3 have simple analytical descriptions. We find that in all these cases the gap closes linearly and that the number of isolated E kx,ky = 0 points at a given transition dictates how the Chern invariant jumps between phases. For example, the critical lines of positive slope in the diagram, J = 1 2 (4κ ± √ 2(J z + 2κ)) correspond to single zeros of the energy at (k x , k y ) = (0, π). The same is true for the two lines of negative slope, J = 1 2 (−4κ ± √ 2(J z − 2κ)) which have zero energy at the point (k x , k y ) = (π, 0).
As proved in [6], each vortex in the odd Chern number phases is accompanied by a single localized zero energy Majorana fermion excitation (assuming the vortices are wellseparated). These localized Majorana fermions in the different phases with odd ν make the vortices into nonabelian anyons and we obtain here 4 out of the 8 different types nonabelian anyons which occur in Kitaev's 16-fold way for free fermionic bulk topological orders (in 2 pairs of types related by parity). Specifically, the ν = ±1 modes correspond to the so called Ising anyon theory and the ν = ±3 modes correspond to the SU (2) 2 anyon theory. These zero energy excitations give rise to a ground state degeneracy and allow for the possibility of nonabelian statistics and fault tolerant quantum computation [1][2][3]. As far as we know, this is the first time that different types of nonabelian anyons carrying Majorana modes have been found in a lattice model. We find zero energy modes at vortex excitations also in some of the abelian phases of the model, whose topological orders cover 5 of the 8 abelian possibilities. In these cases however, the modes are full Dirac modes (consisting of two Majorana modes) and are thus not topologically protected from local error processes. We will discuss one of these Dirac modes briefly below.
The multitude of topological phases in the model allow one to set up a variety of interesting domain wall scenarios. The general result of Hatsugai [34] always applies here: The net number of left movers minus the number of right movers is equal to the difference in Chern number of the phases on either side of the boundary. In the simplest cases, for example for open boundaries, where the topological medium just ends, we find that this is usually satisfied in the simplest possible way; the edge theory is chiral and has a number of zero modes equal to |ν|. This supports the idea that the open edge theo-ries of the ν = ±1 and ν = ±3 phases are Ising and SU (2) 2 CFTs, the simplest CFTs which satisfy the bulk-boundary correspondence.
For edge theories on domain walls the same is often true but not always. More generally we find that the number and type of zero modes that exist on a domain wall always reflects the critical bulk theory that exists on the phase boundary between the phases on either side of the wall. If one examines the edge spectra in a system with one phase domain below another in such a way that the translational invariance in the x direction is maintained, one finds that the wall contains a band of states that intersect the E = 0 line at value of k x that are consistent with the k x values of the zero modes that occur at the bulk phase transition. If the two phases are separated in the phase diagram by more than one phase boundary, then the edge will reflect all the phase boundaries crossed by an interpolating path. Examples of this behavior are given in Figures 4 and 5.
Careful inspection of the phase diagram reveals a number of situations where phases with the same Chern number touch at the intersections and touchings of critical lines. In these cases, we can have domain wall edge states between phases at the same Chern number. Such domain wall edge states are simply inherited from the boundaries between phases with different Chern numbers that meet at these multicritical points and hence they are stable against perturbations which do not remove the multicritical points. However, some of the multicritical points and corresponding edge states can be easily removed by local perturbations to the Hamiltonian, while others are topologically protected, that is, stable against any local perturbation which preserves particle-hole symmetry. We now give an example of each case. Following Yang and coworkers [19] we call the phase around the κ = 0, |J| < J z / √ 2 domain the A z phase and around the κ = 0, |J| > J z / √ 2 line the A xy phase. It turns out that we can make the tetracritical points that occur where the critical lines of positive slope cross and where the critical lines of negative slope cross vanish by varying J y away from J x , in the process connecting the ν = 0 phases on the left and right to the A xy phase. In contrast, the tetra-critical points between the A xy and A z phases are robust features that cannot be so easily removed. In Figure 4 we show the edge spectrum for two horizontal domain walls between A z phase and the upper A xy phase. We see that on a single edge we have a single band of both left and right moving states that intersects E = 0 at k x = 0 and k x = π. We now argue for the stability of these edge bands. First of all, if the system is wide enough in the y-direction, the edge modes on opposite edges in Fig. 4 decouple exponentially and tunneling between the edges cannot lift the zero modes. Backscattering between the leftmoving and rightmoving zero modes on each edge could in principle lift these modes, but this would violate particle-hole symmetry. This is easiest to see for perturbations which preserve translational symmetry. In this case, because we have E(k) = −E(−k) the modes at k x = 0 and k x = π must remain at zero energy. If translational symmetry is broken (for example by disorder), a slightly more delicate argument is needed. In this case k x is no longer a good quantum number, but if we imagine that the disorder was introduced adiabatically, we may still label the energy eigenstates by the former k x -eigenvalues. The disorder may now make the energy of the k x = 0 and k x = π states nonzero, but in order to preserve particle hole symmetry, it must still preserve the property that every state at positive energy has a partner at negative energy. Therefore, in the continuum limit, where k x takes continuous values and energy is a continuous function of k x , there will always be zero modes on the domain wall, though not necessarily at k x = 0 and k x = π. In fact, there must always be a nonzero even number of zero modes, since zero modes can adiabatically appear and disappear in even numbers only.
Of course k x is discrete for any finite system size and these arguments do not guarantee the existence of modes with exactly zero energy, just increasingly soft modes for growing system size. This is in fact how the system behaves when there is a finite and even number of unit cells in the x-direction. However, if there is an odd number of unit cells in the xdirection, the zero energy modes at k x = 0 and k x = π actually belong to different topological sectors and therefore backscattering between the k x = 0 and k x = π modes would violate conservation of a topologically conserved quantum number. In this situation we can have a robust Majorana zero mode forming on an interface between ν = 0 phases even at finite system size.
It is gratifying that both Abelian phases in the model can be analysed using degenerate perturbation theory. For the ν = 0 phase found in the center of the figure (the isolated z-link limit (J z (J 2 x + J 2 y ) 1/2 ) it was shown in [19] that the low-energy effective system, at the 4th order of perturbation theory, is that of an abelian toric code on a square lattice, To the 4th order the constant terms can be calculated to be where N is number of zlinks or fermions and is thus twice the number of unit cells. In addition, we have found that such a mapping is possible also in the isolated square limit (J z (J 2 x + J 2 y ) 1/2 ). Each isolated square is a miniature spin chain with alternating X and Y links as given by eq. 1. In the absence of a vortex on the plaquette, it has a doubly degenerate ground state. This degenerate subspace acts as an effective spin in much the same way as the ground states of the isolated z-dimers above do. Perturbing away from the fully isolated regime, we can derive, at the 4th order, the effective hamiltonian where r = 2(J 2 x + J 2 y ) 1/2 and E 0 = N (r + . In a suitable basis the operatorsW (8) q are of the formσ y lσ y rσ z uσ z d .
Here the subscript labels indicate left, right, up and down respectively and the operatorsσ operate on the low energy antiperiodic sector of the isolated 4-spin chain, see Fig. 4. In this regime, excitations on the octagonal plaquettes can be interpreted as electric (e) and magnetic (m) particles, in a checkerboard pattern as in ref. [6]. We see then that we have a toric code type phase, but on a different lattice than that found for the isolated z-link limit Further distinctions between these Abelian phases can be made if we examine the system's square vortex excitations. Note first that the spectrum of the 4-spin with a vortex excitation is four fold degenerate. This degeneracy is not lifted when we perturb away from the isolated square limit and so we can infer that in this A xy phase there exists a single massless Dirac fermion mode in the presence of a square-vortex excitation. It is important to note that this degeneracy is dependent on the fact that the chain is in an exact eigenstate of the square plaquette operator and hence it is not topologically protected (it can be broken by local perturbations). The domain wall between the abelian A z and A xy phases at ν = 0 is not the only one which shows counter-propagating zero modes. To illustrate this fact, we include figure 5 which shows that domain walls involving the nonabelian C phase, at ν = 3, exhibit similar phenomena. The C − A z (C − B) boundary spectra shown behave as would be naively expected, with 3 (4) zero modes of the same chirality on each of the respective edges (shown in red and green). However, the C − A xy boundary actually retains the zero modes at k x = 0 and k x = π that comes from either taking the path C − B − A xy or C − A z − A xy through the phase diagram.
In conclusion, in this paper we have seen that the squareoctagon lattice variation of the Kitaev honeycomb model has a startling richness of topological phases. In particular it has two essentially different types of nonabelian anyons at Chern numbers ν = ±1 and ν = ±3 and it realizes 9 of the 16 bulk topological orders allowed by Kitaev's 16-fold way. We demonstrated how the bulk phase diagram of the system can be used to predict the behavior of the edge spectra between different topological domains. We encounter various situations where phases with the same Chern number have a nonzero (but even) number of low energy edge bands due to the interface. An examination of the perturbation theory of the abelian phases reveals that distinct abelian phases can both be mapped to Toric Code type systems that are supported on different sub-lattice structures.
The perhaps most interesting aspect of the model is that it gives an example of a system with Majorana modes that seizes the possibility of higher Chern numbers allowed from general symmetry classifications of Altland and Zirnbauer [12,13]. Compared to the Kitaev honeycomb model that maps to spinless paired fermions, the pseudo fermion spin in the squareoctagon model doubles the number of bands and from this simple doubling of bands we could expect Chern numbers from 0 up to ±2. The real surprise is the occurrence of ν = ±3, ±4 which are beyond the expected doubling. In particular, we have found instances where a single band can contribute with ν = ±3.
A deeper understanding of how in detail the multi-band situation allows for increased Chern numbers is left for future study, but it may be imagined that, given a sufficient number of sites in the unit cell of the lattice, it should be possible to realize all 16 possible bulk topological orders for gapped Majorana systems within the same trivalent spin models. In this respect, some progress has also been made using Bloch Hamiltonians derived from the vortex sectors of the Kitaev honeycomb model and related continuum systems, see for example Ref. [35]. It is also interesting to speculate on the situation for gapped phases in 3d variations [36] of Kitaev type models. The winding number is in this case given by a Wess-Zumino-Witten term [14]. As the number of bands is typically higher, at least four, it seems probable that winding numbers other than ν = 0, ±1 could be found. The general symmetry arguments [14] provide a situation that is almost inverse to that in 2d. Only one of the classes with a Z bulk comes with broken T-symmetry, and this class-AIII-is probably less relevant to the Kitaev type models as it lacks particlehole symmetry. The other two classes with Z bulks, namely DIII and CI, require particle-hole symmetry and time-reversal symmetry. Nevertheless, it remains to be seen whether higher winding numbers can actually be realized, possibly by extending these T-invariant models with, for example, additional 4spin interactions.
G. K. was supported by the Alexander von Humboldt foundation. J. K. is grateful to R. Moessner for comments on the first draft and to R. Roy for discussions. J. K. S. was supported by Science Foundation Ireland Principal Investi-