Toward Kitaev's sixteenfold way in a honeycomb lattice model

Kitaev's sixteenfold way is a classification of exotic topological orders in which $\mathbb{Z}_2$ gauge theory is coupled to Majorana fermions of Chern number $C$. The $16$ distinct topological orders within this class, depending on $C \, \mathrm{mod} \, 16$, possess a rich variety of Abelian and non-Abelian anyons. We realize more than half of Kitaev's sixteenfold way, corresponding to Chern numbers $0$, $\pm 1$, $\pm 2$, $\pm 3$, $\pm 4$, and $\pm 8$, in an exactly solvable generalization of the Kitaev honeycomb model. For each topological order, we explicitly identify the anyonic excitations and confirm their topological properties. In doing so, we observe that the interplay between lattice symmetry and anyon permutation symmetry may lead to a"weak supersymmetry"in the anyon spectrum. The topological orders in our honeycomb lattice model could be directly relevant for honeycomb Kitaev materials, such as $\alpha$-RuCl$_3$, and would be distinguishable by their specific quantized values of the thermal Hall conductivity.

Kitaev's sixteenfold way is a classification of exotic topological orders in which Z2 gauge theory is coupled to Majorana fermions of Chern number C. The 16 distinct topological orders within this class, depending on C mod 16, possess a rich variety of Abelian and non-Abelian anyons. We realize more than half of Kitaev's sixteenfold way, corresponding to Chern numbers 0, ±1, ±2, ±3, ±4, and ±8, in an exactly solvable generalization of the Kitaev honeycomb model. For each topological order, we explicitly identify the anyonic excitations and confirm their topological properties. In doing so, we observe that the interplay between lattice symmetry and anyon permutation symmetry may lead to a "weak supersymmetry" in the anyon spectrum. The topological orders in our honeycomb lattice model could be directly relevant for honeycomb Kitaev materials, such as α-RuCl3, and would be distinguishable by their specific quantized values of the thermal Hall conductivity.

I. INTRODUCTION
Topological order is an important cornerstone of modern condensed matter physics which facilitates a classification of gapped phases of matter beyond the classical paradigm of spontaneous symmetry breaking [1]. While topologically ordered phases may be fully symmetric and locally featureless, they are characterized by particular patterns of long-range quantum entanglement [2] which manifest in robust global features, such as a topological ground-state degeneracy [3] and a universal correction to the bipartite entanglement entropy [4,5].
Arguably, the most exciting feature of topological order is the fractionalization of fundamental particles into emergent nonlocal quasiparticles. Because of their nonlocal nature, these fractionalized quasiparticles possess unusual "anyonic" particle statistics in two dimensions that is distinct from both bosons and fermions. In particular, moving one anyon around another one ("braiding") may correspond to a nontrivial operation on the underlying quantum state [6]. For Abelian topological orders, these braiding operations act on a single quantum state, while for non-Abelian topological orders, they act on a set of degenerate quantum states within an internal space spanned by the anyons themselves. In addition to their fundamental scientific appeal, such non-Abelian anyons are highly promising from the perspective of topological quantum computation [7].
Each topological order is uniquely characterized by the topological properties of its anyonic quasiparticle excitations: the distinct classes of anyons as well as the fusion and braiding rules between them [8]. To a large extent, anyons generalize the concept of topological defects in classically ordered systems [9]. Indeed, the anyon classes are topologically distinct in the sense that they cannot be locally transformed into each other, while the fusion rules between these classes are analogous to the combination rules between topological defects. Together with the fundamentally quantum braiding rules, these topological properties fully define a given topological order, mathematically described in the language of topological quantum field theory [10].
The simplest and most widely studied topological order is Z 2 gauge theory [6], which gives rise to an entire class of topological orders when coupled to gapped Majorana fermions of Chern number C [8]. This class contains an infinite number of topologically distinct edge theories as the number of chiral Majorana edge modes is given by the Majorana Chern number C itself. Interestingly, however, the bulk topological order is determined by C mod 16, and the infinitely many edge theories thus correspond to only 16 bulk topological orders with distinct topological properties of the bulk anyons.
This classification, commonly known as Kitaev's sixteenfold way [8], contains both Abelian and non-Abelian topological orders, corresponding to even and odd Majorana Chern numbers, respectively. The topological orders of Kitaev's sixteenfold way are relevant for a wide range of topological materials, including fractional quantum Hall systems, topological superconductors, as well as quantum spin liquids. In particular, recent thermal Hall conductivity measurements in the quantum spin liquid candidate α-RuCl 3 [11] indicate a single Majorana edge mode for a range of applied magnetic fields, corresponding to the non-Abelian C = 1 topological order.
In this work, we study an exactly solvable generalization [39] of the Kitaev model that respects all symmetries of the honeycomb lattice and realizes more than half of the topological orders in Kitaev's sixteenfold way, corresponding to Majorana Chern numbers 0, ±1, ±2, ±3, ±4, and ±8. These topological orders contain both Abelian and non-Abelian anyons with a rich variety of fusion and braiding rules, and are experimentally distinguishable by their different quantized values of the thermal Hall conductivity. For each topological order, we use the exact solution of our model to explicitly identify the anyon classes and verify their fusion rules. In some cases, we find that lattice symmetry becomes intertwined with anyon permutation symmetry, corresponding to weak symmetry breaking [8], and gives rise to a "weak supersymmetry" in the excitation spectrum. Since the additional four-spin interactions of our generalized Kitaev model arise naturally from time-reversal-symmetric perturbations [39], in the same way as the three-spin interactions in the original Kitaev model arise from an external magnetic field, we believe that the |C| > 1 topological orders described in this work are likely to be realized in spin-orbit-coupled honeycomb magnets, such as α-RuCl 3 .

II. LATTICE MODEL
We consider a generalization of the Kitaev spin model on the honeycomb lattice, where the first term is the pure Kitaev model [8] with Ising interactions between the spin components σ α along each α = {x, y, z} bond jk α [see Fig. 1(a)], while the remaining two terms H r with r = 2, 3 contain products of such Ising interactions along paths consisting of r bonds each. If we define jkl αβ to be the path consisting of the two bonds jk α and kl β [see Fig. 1(b)], the second term reads where (αβγ) is a general permutation of (xyz), and (αβγ) is +1 (−1) for even (odd) permutations. Using analogous notation, the third term then takes the form where jklm αβγ and jklm αβα are paths consisting of three bonds each [see Figs. 1(c) and 1(d)]. As it is clear from our construction, the term H r for general r contains (r+1)-spin interactions and thus breaks (preserves) timereversal symmetry for even (odd) r. We remark that the term H 2 was already introduced in Ref. [8] while the term H 3 was first considered in Ref. [39]. It is also important to note that these two terms are respectively generated by time-reversal-breaking and time-reversal-symmetric perturbations on top of the pure Kitaev model.
Remarkably, the generalized Kitaev model in Eq. (1) is exactly solvable in the same way as the original Kitaev model [8]. By expressing each physical spin component as a product of two Majorana fermions, σ α j = ib α j c j , the Hamiltonians H n in Eqs. (2)-(4) become where the Z 2 gauge fields u α jk = −u α kj ≡ ib α j b α k along the bonds jk α are conserved quantities that commute with each other. Therefore, the Hamiltonian H in Eq. (1) describes free fermions coupled to a static Z 2 gauge theory, and c j can be identified as deconfined Majorana fermion ("spinon") degrees of freedom. In terms of these Majorana fermions, each term H r in Eq. (5) corresponds to r-th-neighbor hopping [40]. Also, unlike the gauge fields themselves, the product of the gauge fields around any plaquette p [see Fig. 1(a)] is a gauge-invariant quantity that can be expressed in terms of the physical spins: Thus, W p = ±1 can be identified as static Z 2 gauge flux ("vison") degrees of freedom. While Eq. (5) reduces to a quadratic fermion problem in each flux sector, {W p = ±1}, represented with an appropriate gauge-field configuration, {u α jk = ±1}, it is not immediately clear which flux sector contains the ground state of the physical spin model H. For the pure Kitaev model H 1 , it is guaranteed by Lieb's theorem [41] that the ground state belongs to the 0-flux sector characterized by W p = +1 for all p. However, Lieb's theorem no longer applies if the additional terms H 2 and/or H 3 are included in the spin model. Indeed, it was demonstrated in Ref. [39] that the frustration between H 1 and H 3 can stabilize a wide range of flux sectors as a function of K 3 /K 1 and K 3 /K 1 (see Fig. 2), including the 1-flux sector characterized by W p = −1 for all p, as well as fractional-flux sectors in which a nontrivial fraction of the plaquettes have W p = −1 rather than W p = +1. In these fractional-flux sectors, the plaquettes with W p = −1 form crystalline structures ("vison crystals") that spontaneously break translation symmetry (see Fig. 3).
Assuming an infinitesimally small coupling constant K 2 , we start from the time-reversal-symmetric Hamiltonian H 1 + H 3 [39] and treat the Hamiltonian term H 2 as a time-reversal-breaking perturbation. Due to the finite flux gap, the ground-state flux sectors in Fig. 2  mally small time-reversal-breaking perturbation can have a dramatic effect on their low-energy physics [8].

III.1. Quadratic Hamiltonians
For each ground-state flux sector in Fig. 2, the gaugefield configuration in Fig. 3 gives rise to a quadratic Majorana problem [see Eq. (5)]. The unit cell of this Majorana problem may consist of n > 1 honeycomb unit cells for two distinct reasons. First, the physical unit cell is enlarged in the fractional-flux sectors because translation symmetry is spontaneously broken. This enlargement is twofold for the 1/2-flux sector, threefold for the 1/3-flux and 2/3-flux sectors, and fourfold for the 1/4-flux and 3/4-flux sectors. Second, if the physical unit cell has an odd number of W p = −1 plaquettes, translation symmetry acts projectively on the Majorana fermions. In this case, the Majorana unit cell, as characterized by the gauge-field configuration, must consist of two physical unit cells. Consequently, the Majorana unit cell has an additional twofold enlargement in all flux sectors except for the 0-flux and 2/3-flux sectors.
For each flux sector, we label the honeycomb sites as j = (r, λ) and the corresponding Majorana fermions as c j = c r,λ , where r is the lattice vector of the Majorana unit cell, and λ = (µ, ν) in terms of the sublattice index µ = A, B and the index ν = 1, . . . , n specifying the particular honeycomb unit cell within the Majorana unit cell. Using this labeling convention, the quadratic Majo- rana Hamiltonian takes the general form where eachH r −r,λ,λ is proportional to the product of the static gauge fields u α jk = ±1 along a path connecting the sites (r, λ) and (r , λ ). Introducing the momentumspace complex fermions where N is the number of honeycomb sites, and arranging them into the 2n-component vector the Hamiltonian in Eq. (7) can then be written as where H q is a 2n × 2n matrix with elements By diagonalizing the matrix H q at each momentum q, one obtains Majorana bands at both positive and negative energies. However, since ψ −q,λ = ψ † q,λ by definition, there is a redundancy in our description, and only the bands with positive energies are physical.

III.2. Projective symmetries
We now discuss the general symmetries of the Majorana problem in Eq. (10). Since the Majorana fermions are fractionalized degrees of freedom, symmetries may act on them projectively [42], i.e., the classical relations between symmetry operations may only be satisfied up to an overall complex phase factor e iϕ . However, as the Majorana fermions are coupled to Z 2 gauge fields, this phase factor must actually be a sign ±1 [43].
We first consider the translation symmetries T 1 and T 2 along the lattice vectors R 1 and R 2 of the physical unit cell (see Fig. 3). Note that the physical unit cell depends on the particular flux sector and may be larger than the original honeycomb unit cell due to spontaneous breaking of translation symmetry in the fractional-flux sectors. If the physical unit cell has no overall Z 2 flux, corresponding to an even number of W p = −1 plaquettes, trans- Conversely, if the physical unit cell has an overall Z 2 flux, corresponding to an odd number of W p = −1 plaquettes, translation symmetry acts projectively on the Majorana fermions, and the two elementary translations anticommute: {T 1 , T 2 } = 0. The eigenvalues of the translations T 1 and T 2 are then no longer compatible quantum numbers for the Majorana fermions. Nevertheless, since [T 2 1 , T 2 ] = 0, one may consider a larger Majorana unit cell spanned by 2R 1 and R 2 , which translates into a smaller Brillouin zone spanned by 1 2 G 1 and G 2 . In fact, the Majorana spectrum is periodic with respect to an even smaller Brillouin zone spanned by 1 2 G 1 and 1 2 G 2 [see Fig. 4(a)] because the residual symmetry T 1 anticommutes with T 2 and hence corresponds to a shift 1 2 G 2 in the Majorana momentum. Thus, one may use a compact Brillouin zone spanned by 1 2 G 1,2 and indicate that each Majorana band has a twofold "translation degeneracy" [see Fig. 4(b)]. Different points in this compact Brillouin zone are labeled by T 2 1 and T 2 2 , while the two degenerate Majorana fermions at a given point are labeled by T 1 and mapped onto each other by T 2 (or vice versa).
We next consider time-reversal symmetry T and inversion symmetry P. Time reversal is an antiunitary operation, {T , i} = 0, and is only a symmetry for K 2 = 0. In each flux sector, it acts on the Majorana fermions as T : c r,(A,ν) → c r,(A,ν) , c r,(B,ν) → −c r,(B,ν) , (12) and hence satisfies T 2 = +1. In contrast to time rever-sal, inversion is a unitary operation, [P, i] = 0, and is a general symmetry of our model. While the action of inversion on the Majorana fermions depends on the given flux sector, it always exchanges the two sublattices A and B, and thus necessarily anticommutes with time reversal: {P, T } = 0. Also, inversion satisfies P 2 = −1 [43] in all flux sectors except for the 3/4-flux sector. For simplicity, we ignore the 3/4-flux sector in the rest of this work and only return to it briefly in Sec. VI.
Finally, the redundancy in our description, corresponding to H −q = −H * q [see Eq. (11)], gives rise to an emergent antiunitary particle-hole symmetry C, which satisfies [C, T ] = 0, [C, P] = 0, and C 2 = +1. We emphasize that particle-hole symmetry is actually an antisymmetry as it anticommutes with the Hamiltonian. While T , P, and C each reverse the fermion momentum q, their two independent products S = T C and R = PC transform the fermions at momentum q among each other: The 2n × 2n transformation matrices are given by where K denotes complex conjugation, I is the n × n unit matrix, and P is an n × n permutation matrix satisfying P · P T = I. The unitary antisymmetry S can thus be identified as sublattice symmetry, while the antiunitary antisymmetry R can be interpreted as an effective momentum-conserving particle-hole symmetry.
Since R is a general antisymmetry of our model, the Hamiltonian matrix H q in Eq. (11) satisfies {R, H q } = 0, which implies that the Majorana spectrum is symmetric around zero energy at each momentum q. Furthermore, in the time-reversal-symmetric limit of K 2 = 0, the antisymmetry S requires {S, H q } = 0 and therefore constrains the Hamiltonian matrix to the form In this limit, the eigendecomposition of the 2n × 2n matrix H q is equivalent to the singular value decomposition of the n × n matrix M q .

III.3. Generic Majorana nodes
In terms of the low-energy physics, the gapless nodes of the momentum-space Majorana spectrum are of particular interest. Because of the antisymmetry R, a generic nodal momentum Q has two zero-energy fermions that correspond to distinct Majorana bands of opposite energies. By projecting onto these two low-energy Majorana bands around q = Q, one then obtains an effective lowenergy theory of the given Majorana node.
For simplicity, we start our discussion from the timereversal-symmetric limit of K 2 = 0. Since the Hamiltonian matrix H Q takes the form of Eq. (15), we can choose the low-energy subspace to be spanned by two fermions located on the two respective sublattices A and B, where u Q (v Q ) is the left (right) eigenvector of the matrix M Q corresponding to zero eigenvalue [44]. If we project onto these two low-energy fermions, the antisymmetries S and R are represented with the 2 × 2 matriceŝ and the most general Hamiltonian matrix anticommuting with bothŜ andR takes the form where τ 1,2,3 are the Pauli matrices. Since there are two independent real coefficients, β 1 (q) and β 2 (q), that must vanish at the nodal momentum Q itself, the generic nodal structures in two dimensions are point nodes. Expanding β 1,2 (q) up to linear order in δq ≡ q − Q = (δq x , δq y ), these point nodes are generically Dirac nodes with linear dispersions [45]. Also, by considering the complex phase of β 1 (q) + iβ 2 (q) at q = Q + (cos ϑ, sin ϑ)δq as a function of ϑ, one can assign a winding number W Q to each Dirac node, which is generically given by with sgn x ≡ x/|x|. Therefore, each Dirac node is a stable U(1) vortex protected by time-reversal symmetry or, equivalently, by sublattice symmetry. If we then break time-reversal symmetry with an infinitesimally small K 2 = 0, the Hamiltonian matrixĤ q still anticommutes withR but no longer withŜ. Thus, its most general form readŝ where the third coefficient may be expanded up to linear order in K 2 such that β 3 (q) = m q K 2 . Since there are three independent real coefficients, nodes can only emerge as a result of fine tuning, and the spectrum is generically gapped. In particular, at each Dirac node of the K 2 = 0 limit, the Hamiltonian matrix becomeŝ and the Dirac node at momentum Q is thus gapped out by a fermion mass term ∝ m Q K 2 .

IV.1. Definition and numerical results
Since the Majorana spectrum is symmetric around zero energy at each momentum q and generically gapped for K 2 = 0 in each ground-state flux sector, diagonalizing the quadratic Hamiltonian in Eq. (10) gives equal numbers of Majorana bands at strictly positive and strictly negative energies. Therefore, one can define a ground-state Chern number of the Majorana fermions by summing the Chern numbers of all the negative-energy Majorana bands. As it was argued in Ref. [8], the low-energy physics of the corresponding topological order is completely determined by this Majorana Chern number C.
Mathematically, the eigendecomposition of the 2n×2n matrix H q in Eq. (11) gives 2n eigenvalues ε q,κ and 2n corresponding eigenvectors w q,κ with κ = 1, . . . , 2n at each Majorana momentum q. The Chern number of each Majorana band, labeled with κ, is then obtained as where the Berry curvature at momentum q is given by in terms of the corresponding Berry connection We note that the cross product of two vectors is a scalar in two dimensions. Also, since the integral in Eq. (23) is defined over the conventional Brillouin zone in Fig. 4(a), it must be multiplied by 2 when calculated over the compact Brillouin zone in Fig. 4(b). Using Eqs. (23)-(25), the Majorana Chern number can then be numerically computed in each flux sector via where the bands κ = 1, . . . , 2n are arranged by increasing energy eigenvalues so that the summation is over all negative-energy Majorana bands. The numerical results for the Majorana Chern numbers are summarized in Fig. 5. For most of the flux sectors in Fig. 2, the Chern number only depends on K 2 and is otherwise the same throughout the entire flux sector. In contrast, for the 2/3-flux sector, there are two disconnected ("upper" and "lower") phases with distinct Chern numbers, while for the 1-flux sector, there are several phases with distinct Chern numbers that are separated by topological transitions as a function of K 3 /K 1 . We note that K 3 /K 1 is an irrelevant parameter in the 1-flux sector as symmetry-related pairs of K 3 interactions [see Fig. 1(c)] give rise to equivalent Majorana hopping terms with a perfect destructive interference between them [39].

IV.2. Analytical understanding
By studying the phase transitions between the various phases in Fig. 5, we can also understand their Majorana Chern numbers analytically. If the Majorana spectrum is gapped for K 2 = 0, the Chern number C vanishes due to time-reversal symmetry and is robust against an infinitesimally small K 2 = 0. Thus, the 1/3-flux phase, the 1/2-flux phase, and the "lower" 2/3-flux phase of Fig. 2 are all characterized by C = 0. If the Majorana spectrum has gapless nodes for K 2 = 0, these nodes are all gapped out by an infinitesimally small K 2 = 0, and the resulting gapped phases have opposite Chern numbers ±C for opposite signs of K 2 . Since a change in the Chern number is always connected to a closing gap, the Chern number C at K 2 > 0 can be understood as a sum of contributions from the various nodes at K 2 = 0.
For K 2 = 0, the contribution to the Chern number from the given Dirac node,Ĉ Q , is the Chern number of the negative-energy band in the low-energy theory. For the Hamiltonian matrix in Eq. (27), this quantity can be calculated by means of a standard formula [8]: where d(q) ≡ β(q)/|β(q)|. Geometrically,Ĉ Q is simply the number of "skyrmions" in the vector field d(q). Since the vector-field configuration in Eq. (28) corresponds to a half skyrmion or meron (see Fig. 6), the contribution of the given Dirac node to the Chern number becomeŝ  Fig. 7 for the 0-flux phase, the 1/4-flux phase, the "upper" 2/3-flux phase, and the Dirac 1-flux phase of Fig. 2. For K 2 > 0, the resulting total Chern numbers, C = QĈ Q , are 1, 2, 3, and 2, respectively. In general, Dirac nodes emerge in pairs related by inversion symmetry. Since the two nodes in any pair have opposite mass coefficients m Q as well as opposite winding numbers W Q , each pair contributes ±1 to the Chern number. If translation symmetry acts projectively on the Majorana fermions, corresponding to an overall Z 2 flux in the physical unit cell, the total Chern number is then necessarily even as these contributions ±1 come in identical pairs due to translation degeneracy.

IV.2.2. Line nodes
While the generic Majorana nodal structures at K 2 = 0 are point nodes, these point nodes seem to coexist with line nodes in the three Fermi 1-flux phases of Fig. 2. As expected, these accidental line nodes are unstable against generic further-neighbor Majorana hopping terms that respect the projective symmetries of the system. In fact, the accidental "phases" with line nodes can be understood as phase transitions between two generic phases with point nodes as a function of a fifth-neighbor hop-  ping amplitude K 5 . While each line node is gapped out into six Dirac nodes at the same momenta for K 5 > 0 and K 5 < 0, the winding number of each Dirac node changes sign at K 5 = 0 (see Fig. 8). Thus, Dirac nodes with opposite winding numbers must be connected by line nodes at the phase transition so that they can exchange their winding numbers with each other.
Instead of constructing a low-energy theory for each accidental line node, it is then more natural to include an infinitesimally small K 5 = 0 and consider the low-energy theories of the resulting Dirac nodes. Their contributions to the Chern number at K 2 = 0, as given by Eq. (30), are illustrated in Fig. 8 for all three Fermi 1-flux phases and for both signs of K 5 . For the first Fermi phase, the Dirac nodes corresponding to each line node have a vanishing net contribution to the Chern number for both K 5 > 0 and K 5 < 0. Thus, we can deduce that the gapped phase at K 5 = 0 and K 2 > 0 is adiabatically connected to that obtained from the Dirac phase at K 2 > 0 and that its total Chern number is C = 2. Conversely, for the second and the third Fermi phases, the Dirac nodes corresponding to each line node have opposite net contributions to the Chern number for K 5 > 0 and K 5 < 0. Thus, for K 5 = 0, there are two possible gapped phases at K 2 > 0 with total Chern numbers 8 and −4, respectively. To discriminate between these two scenarios, we consider the multicritical point at K 2 = 0 and K 3 = 1 2 K 1 .

IV.2.3. Multicritical point
The second and the third K 2 = 0 and K 3 = 1 2 K 1 [46] where the accidental line node shrinks to a single momentum Q (see Fig. 9). Remarkably, at this momentum Q, the 4 × 4 Hamiltonian matrix in Eq. (11) takes the exact general form where "⊗" denotes the Kronecker product, while τ 0,1,2,3 and η 0,1,2,3 are the Pauli matrices acting on the µ = A, B and ν = 1, 2 degrees of freedom, respectively. Since the two terms in H Q commute, we can use an appropriate canonical transformation to recast it in the simpler form This matrix has zero eigenvalues, corresponding to phase transitions, along the lines K 2 = ±(K 3 − 1 2 K 1 ) in parameter space, while it vanishes identically at the intersection of these lines, i.e., at the multicritical point. In the following, we determine the Chern numbers of the gapped phases around the multicritical point by following the two-step scheme shown in Fig. 5(g).
First, we fix K 3 = 1 2 K 1 and construct the low-energy theory of the multicritical point node at momentum Q for K 2 K 1 . By expanding H q in Eq. (11) up to first order in both K 2 and δq ≡ q − Q, and projecting onto the basis vectors in Eq. (32), we obtain Therefore, the low-energy theory contains a pair of Dirac nodes corresponding toη 3 = ±1. Since the two respective Dirac nodes have mass coefficients m Q = ±2 √ 3 and winding numbers W Q = ±1, they contributeĈ Q = 1 to the Chern number at K 2 > 0. Together with the con-tributionsĈ Q1 =Ĉ Q2 = 1/2 from the two Dirac nodes at momenta Q 1,2 (see Fig. 9), the total Chern number of the phase at K 3 = 1 2 K 1 and K 2 > 0 is then where the additional factor of 2 comes from translation degeneracy in the 1-flux sector.
Next, we fix a particular value of K 2 > 0 and consider the phase transition at K 3 = 1 2 K 1 + K 2 . At this phase transition, the low-energy subspace at the critical momentum Q is spanned by theτ 3 = −1 basis vectors in Eq. (32). By expanding H q in Eq. (11) up to first order in δK ≡ K 3 − 1 2 K 1 − K 2 and second order in δq ≡ q − Q, and projecting onto these basis vectors, we obtain where θ and χ are positive numbers. Therefore, the lowenergy theory of the phase transition is a quadratic point node at momentum Q. The corresponding change in the Chern number, δĈ Q , across the phase transition is the difference between the Chern numbers of the negativeenergy bands at δK > 0 and δK < 0. Since Eq. (35) assumes the general form of Eq. (27), these Chern numbers are given by Eq. (29) in terms of the respective vector fields d(q) plotted in Fig. 10. For δK < 0, the vectorfield configuration is topologically trivial, and the Chern number is thusĈ Q,− = 0. For δK > 0, the vector-field configuration corresponds to a double skyrmion, and the Chern number is thusĈ Q,+ = 2. Remembering translation degeneracy, the change in the total Chern number across the phase transition is then and the total Chern number of the gapped phase next to the third Fermi phase is C = 4 + 4 = 8. Since the phase transition at K 3 = 1 2 K 1 −K 2 is governed by an analogous low-energy theory, the gapped phase next to the second Fermi phase also has a total Chern number C = 8.

V.1. Anyon classes and fusion rules
Since the gapped phases described in the previous section are all topologically ordered, they have anyonic excitations characterized by particular fusion and braiding properties. Given that gapped Majorana fermions with a total Chern number C are coupled to a Z 2 gauge theory, the topological classification is understood in terms of Kitaev's sixteenfold way [8]. While all phases with different Chern numbers C are topologically distinct theories with different numbers of chiral Majorana edge modes, and are experimentally distinguishable by their thermal Hall conductivities, κ xy = πCT /12, at temperature T , there are only 16 distinct classes in terms of their bulk anyon properties, determined by C mod 16.
For each phase, we can use the exact solution of our lattice model to explicitly identify the topologically distinct anyon classes of the corresponding theory and to verify the expected fusion rules between them [8]. In general, the Majorana fermions are identified as the fermion excitations , while the gauge fluxes are related to the various classes of vortex excitations. If the Chern number C is odd, there is only one vortex class σ, and each flux excitation with respect to the ground-state flux sector corresponds to such a vortex excitation σ. If the Chern number C is even, there are two topologically distinct vortex classes denoted by e and m for C mod 4 = 0 and by a andā for C mod 4 = 2. Therefore, a flux excitation at any given plaquette may correspond to either of the two vortex classes. Since the two vortex classes differ by a fermion in each case, as indicated by the fusion rules the vortex classes corresponding to the various plaquettes can be mapped out by considering the fermion parity of the ground state within the flux sector (p, p ) that contains two flux excitations at a general plaquette p and at a far-away reference plaquette p . If the ground-state fermion parities of the flux sectors (p 1 , p ) and (p 2 , p ) are identical (opposite), the flux excitations at the plaquettes p 1 and p 2 correspond to identical (distinct) vortex excitations. For each phase with even C, the resulting map of the vortex classes is depicted in Fig. 11. We note that the two vortex classes are related by an anyon permutation symmetry and that the same maps are thus equally valid with e ↔ m and a ↔ā.
In terms of the anyon fusion rules, the fermions have similar properties in all phases. Indeed, the general fusion rule × = 1 indicates that two fermion excitations fuse into a topologically trivial excitation. In contrast, there are three distinct scenarios for the fusion rules between two identical vortices [8]. If the Chern number C is odd, the fusion rule σ × σ = 1 + indicates that the vortices are non-Abelian anyons as a pair of them has a degenerate internal space. The two states in this internal space correspond to two fusion channels into a trivial excitation and a fermion excitation, respectively. In the lattice model, the internal space manifests as a zero-energy fermion in any flux sector containing two flux excitations far away from each other. Conversely, if the Chern number C is even, the vortices are Abelian anyons with only one fusion channel. If C mod 4 = 0, the fusion rules e × e = m × m = 1 (38) indicate that two identical vortices fuse into a trivial excitation, whereas if C mod 4 = 2, the fusion rules a × a =ā ×ā = indicate that two identical vortices fuse into a fermion excitation. In the lattice model, these fusion rules are reflected in the ground-state fermion parity of a flux sector containing two far-away flux excitations that correspond to the same vortex class. For C mod 4 = 0 (C mod 4 = 2), this fermion parity is even (odd) with respect to the overall ground state of the model. We finally note that the anyon braiding rules are even more specific to the topological order than the anyon fusion rules. To explicitly check these braiding rules in our lattice model, we would need to calculate the complex hopping matrix elements as one flux excitation is moved around another one in such a way that the distance between the two flux excitations is much larger than the correlation length at each step. Unfortunately, for most of our phases, the correlation length exceeds the maximal system size for which such a calculation would be feasible at all. Nevertheless, for all phases with sufficiently small correlation lengths, the results of such a calculation are in agreement with the braiding rules in Ref. [8].

V.2. Weak symmetry breaking and supersymmetry
For each map of vortex classes in Fig. 11, it is instructive to consider the interplay between the unbroken lattice symmetries in the given flux sector and the relevant anyon permutation symmetry (e ↔ m or a ↔ā). Interestingly, for the 1/2-flux phase and the "lower" 2/3-flux phase of Fig. 2, certain lattice symmetries become intertwined with the anyon permutation symmetry e ↔ m in the sense that they map plaquettes corresponding to the two vortex classes e and m onto each other. In each case, one such lattice symmetry is a twofold rotation around a z bond separating two W p = −1 plaquettes.
This kind of interplay, commonly known as weak symmetry breaking [8], was first discussed for the spatially anisotropic gapped phase ("A phase") of the original Kitaev model, where the anyon permutation symmetry e ↔ m is intertwined with translation symmetry. While there is no symmetry breaking in the conventional sense as all ground-state correlations are fully symmetric, there is a symmetry breaking in the topological properties of the anyonic excitations as symmetries map topologically distinct anyons onto each other.
Even more remarkably, for both the 1/2-flux phase and the "lower" 2/3-flux phase of Fig. 2, there are certain plaquettes that do not correspond to any particular vortex class e or m (see Fig. 11). Since these plaquettes are mapped onto themselves by a lattice symmetry that is intertwined with the anyon permutation symmetry e ↔ m, a vortex excitation at such a plaquette has a degenerate internal space consisting of two states that correspond to the two vortex classes e and m. In the lattice model, this degenerate internal space manifests as a zero-energy fermion in any flux sector that contains a flux excitation at such a plaquette. The creation and annihilation operators of the zero-energy fermion can then be identified as the generators of a fermionic symmetry that is reminiscent of supersymmetry [47].
We emphasize that this fermionic symmetry is not a supersymmetry in the conventional sense because it relates the vortex excitations e and m that are both bosons in terms of their self statistics [8]. Nevertheless, it may be interpreted as a generalized "weak supersymmetry" because it relates topologically distinct anyonic excitations that are different in their mutual statistics with respect to each other. We also note that such a fermionic symmetry does not manifest in the A phase of the original Kitaev model because translation symmetry does not have a fixed point and does not map any plaquette onto itself. However, we expect it to be a generic feature of symmetry-enriched topological order whenever an anyon permutation symmetry is intertwined with a point-group symmetry, such as a rotation or a reflection.

VI. DISCUSSION
In this work, we have realized a wide range of distinct topological orders in an exactly solvable spin model on the honeycomb lattice. Each of these topological orders is a Z 2 gauge theory coupled to Majorana fermions with a total Chern number C. Given their respective Majorana Chern numbers 0, ±1, ±2, ±3, ±4, and ±8, these topological orders correspond to more than half of Kitaev's sixteenfold way. In particular, the C = ±3 phases realize non-Abelian topological orders that are distinct from the C = ±1 phases of the Kitaev honeycomb model both in the number of Majorana edge modes and in the braiding properties of the non-Abelian anyons. Also, the C = ±8 and C = ±4 phases realize Abelian topological orders in which the gauge fluxes have fermionic and semionic particle statistics, respectively [8]. Since our model emerges naturally from the Kitaev honeycomb model in the presence of perturbations [39], we expect that its topological orders are likely to be realized in spin-orbit-coupled honeycomb magnets, such as α-RuCl 3 .
From an experimental perspective, the key signature of each C = 0 topological order is a specific quantized value of the thermal Hall conductivity, κ xy = πCT /12, at any temperature T below the bulk energy gap. This quantized value is directly proportional to the chiral central charge, c = C/2, of the edge theory, whose integer (fractional) values correspond to Abelian (non-Abelian) bulk topological orders. In the fractional-flux sectors, the topological order is also accompanied by a spontaneous breaking of translation symmetry which, in the presence of any spin-lattice coupling, gives rise to a periodic lattice distortion and can thus be picked up with nuclear magnetic resonance or elastic x-ray scattering. Moreover, the spontaneous breaking of discrete translation leads to a finite-temperature phase transition [39] that is readily observable in the specific heat.
Conceptually, the generalization of the Kitaev honeycomb model in this work facilitates a convenient bandstructure engineering for Majorana fermions coupled to a Z 2 gauge theory. In many ways, the resulting topological phases are analogous to those studied in the context of noninteracting electrons. For example, the C = 0 topological orders in this work correspond to Chern insulators as the Majorana fermions form topologically nontrivial bands with finite Chern numbers. In the future, it would be interesting to engineer other topological band structures for the Majorana fermions and thereby realize Majorana analogs of other noninteracting topological phases, such as topological insulators. Alternatively, it could be worth completing the realization of Kitaev's sixteenfold way by engineering Majorana band structures with total Chern numbers ±5, ±6, and ±7.
We finally note that, compared to the other flux sectors in this work, the Majorana fermions have a different behavior in the 3/4-flux sector because inversion symmetry acts differently on them. In the presence of time-reversal symmetry, the entire Majorana spectrum is twofold de-generate, and the Dirac nodes thus come in pairs at each nodal momentum. If time-reversal symmetry is broken, these Dirac nodes then expand into line nodes instead of gapping out. While the line nodes are stable in the noninteracting Majorana theory, they are expected to have instabilities against Majorana interactions [48]. The study of these instabilities and the resulting topological phases is the subject of ongoing further work.