A picture of me.

Tom Hodson

Maker, Baker Programmer Reformed Physicist RSE@ECMWF

The Kitaev Honeycomb Model

Figure 1: (a) The Kitaev honeycomb model is defined on a honeycomb lattice. The special feature of the honeycomb lattice that makes the model solvable is that each vertex is joined by exactly three bonds, i.e., the lattice is trivalent. One of three labels is assigned to each bond \{x,\;y,\;z\} here represented by colour. (b) After transforming to the Majorana representation we get an emergent gauge degree of freedom u_{jk} = \pm 1 that lives on each bond, the bond variables. These are antisymmetric, u_{jk} = -u_{kj}, so we represent them graphically with arrows on each bond that point in the direction that u_{jk} = +1. (c) The Majorana transformation, discussed later in the main text, can be visualised as breaking each spin into four Majoranas b_i^x,\;b_i^y,\;b_i^z,\;c_i. The b_i^x,\;b_i^y and b_i^z Majoranas then pair along the bonds forming conserved \mathbb{Z}_2 bond operators u_{jk} = \langle i b_i^\alpha b_j^\alpha \rangle. The remaining c_i operators form an effective quadratic Hamiltonian H = \frac{i}{4} \sum_{\langle i,j\rangle_\alpha} 2J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j.

The Spin Hamiltonian

The Kitaev Honeycomb (KH) model is an exactly solvable model of interacting spin\(-1/2\) spins on the vertices of a honeycomb lattice. Each bond in the lattice is assigned a label \(\alpha \in \{ x, y, z\}\) and couples two spins along the \(\alpha\) axis. See fig. 1 for a diagram of the setup.

This gives us the Hamiltonian \[ H = - \sum_{\langle j,k\rangle_\alpha} J^{\alpha}\sigma_j^{\alpha}\sigma_k^{\alpha}, \qquad{(1)}\] where \(\sigma^\alpha_j\) is the \(\alpha\) component of a Pauli matrix acting on site \(j\) and \(\langle j,k\rangle_\alpha\) is a pair of nearest-neighbour indices connected by an \(\alpha\)-bond with exchange coupling \(J^\alpha\). Kitaev introduced this model in his seminal 2006 paper  [1].

The KH can arise as the result of strong spin-orbit couplings in, for example, the transition metal based compounds  [26]. The model is highly frustrated: each spin would like to align along a different direction with each of its three neighbours but this cannot be achieved even classically  [7,8]. This frustration leads the model to have a Quantum Spin Liquid (QSL) ground state, a complex many-body state with a high degree of entanglement but no long-range magnetic order even at zero temperature. While the possibility of a QSL ground state was suggested much earlier  [9], the KH model was the first exactly solvable models of the QSL state. The KH model has a rich ground state phase diagram with gapless and gapped phases, the latter supporting fractionalised quasiparticles with both Abelian and non-Abelian quasiparticle excitations. Anyons have been the subject of much attention because, among other reasons, they can be braided through spacetime to achieve noise tolerant quantum computations  [10]. At finite temperature the KH model undergoes a phase transition to a thermal metal state  [11]. The KH model can be solved exactly via a mapping to Majorana fermions. This mapping yields an extensive number of static \(\mathbb Z_2\) fluxes tied to an emergent gauge field with the remaining fermions are governed by a free fermion hamiltonian.

This section will go over the standard model in detail, first discussing the spin model, then detailing the transformation to a Majorana hamiltonian that allows a full solution while enlarging the Hamiltonian. I will discuss the properties of the emergent gauge fields and the projector. The next section will discuss anyons, topology and the Chern number, using the KH model as an explicit example. I will then discuss the ground state found via Lieb’s theorem as well as work on generalisations of the ground state to other lattices. Finally, I will present the phase diagram.

The Spin Model

Figure 2: A visual introduction to the Kitaev Honeycomb (KH) model. In the spin representation, bond operators k_{ij} = \sigma^\alpha_i \sigma^\alpha_j can be represented by arrows from site i to site j, from these we can form the Hamiltonian and conserved plaquette and loop operators. Note that the bond operators are defined differently in the Majorana representation but the definitions coincide for closed loops of bonds.

As discussed in the introduction, spin hamiltonians like that of the KH model arise in electronic systems as the result the balance of multiple effects  [5]. For instance, in certain transition metal systems with \(d^5\) valence electrons, crystal field and spin-orbit couplings conspire to shift and split the \(d\) orbitals into moments with spin \(j = 1/2\) and \(j = 3/2\). Of these, the bandwidth \(t\) of the \(j= 1/2\) band is small, meaning that even relatively meagre electron correlations (such as those induced by the \(U\) term in the Hubbard model) can lead to the opening of a Mott gap. From there we have a \(j = 1/2\) Mott insulator whose effective spin-spin interactions are again shaped by the lattice geometry and spin-orbit coupling leading some materials to have strong bond-directional Ising-type interactions  [12,13]. In the KH model the bond directionality refers to the fact that the coupling axis \(\alpha\) in terms like \(\sigma_j^{\alpha}\sigma_k^{\alpha}\) is strongly bond dependent.

In the spin hamiltonian eq. 1, we can already tease out a set of conserved fluxes that will be key to the model’s solution. These fluxes are the expectations of Wilson loop operators

\[\hat{W}_p = \prod_{\langle i,j\rangle_\alpha\; \in\; p} \sigma_i^{\alpha}\sigma_j^{\alpha},\]

the products of bonds winding around a closed path \(p\) on the lattice. These operators commute with the Hamiltonian and so have no time dynamics. The winding direction does not matter so long as it is fixed. By convention, we will always use clockwise. Each closed path on the lattice is associated with a flux. The number of conserved quantities grows linearly with system size and is thus extensive. This is a common property for exactly solvable systems and can be compared to the heavy electrons present in the Falicov-Kimball model. The square of two loop operators is one so any contractible loop can be expressed as a product of loops around plaquettes of the lattice, as in fig. 3. For the honeycomb lattice, the plaquettes are the hexagons. The expectations of \(\hat{W}_p\) through each plaquette, the fluxes, are therefore enough to describe the whole flux sector. We will focus on these fluxes, denoting them by \(\phi_i\). Once we have made the mapping to the Majorana Hamiltonian, I will explain how these fluxes can be connected to an emergent \(B\) field which makes their interpretation as fluxes clear.

Figure 3: In the KH model, Wilson loop operators \hat{W}_p = \prod_{\langle i,j\rangle_\alpha\; \in\; p} \sigma_i^{\alpha}\sigma_j^{\alpha} can be composed via multiplication to produce arbitrary contractible loops. As a consequence, we need only to keep track of the value of the flux through each plaquette \phi_i. This relationship between the u_{ij} around a region and the fluxes inside it is evocative of Stokes’ theorem from classical electromagnetism. Indeed, it turns out to be closely related as we shall see later.

It is worth noting in passing that the effective Hamiltonian for many Kitaev materials incorporates a contribution from an isotropic Heisenberg term \(\sum_{i,j} \vec{\sigma}_i\cdot\vec{\sigma}_j\). This is referred to as the Heisenberg-Kitaev model  [14]. Materials for which the Kitaev term dominates are generally known as Kitaev Materials. See  [5] for a full discussion of Kitaev Materials.

As with the Falicov-Kimball model, the KH model has an extensive number of conserved quantities, the fluxes. So again we will work in the simultaneous eigenbasis of the fluxes and the Hamiltonian so that we can treat the fluxes like a classical degree of freedom. This is part of what makes the model tractable. We will find that the ground state of the model corresponds to some particular choice of fluxes. We will refer to local excitations away from the flux ground state as vortices. In order to fully solve the model however, we must first move to a Majorana representation.

The Majorana Model

Majorana fermions are something like ‘half of a complex fermion’ and are their own antiparticle. From a set of \(N\) fermionic creation \(f_i^\dagger\) and anhilation \(f_i\) operators, we can construct \(2N\) Majorana operators \(c_m\). We can do this construction in multiple ways subject to only mild constraints required to keep the overall commutations relations correct  [1]. Majorana operators square to one, but otherwise have standard fermionic anti-commutation relations.

\(N\) spins can be mapped to \(N\) fermions with the well-known Jordan-Wigner transformation and indeed this approach can be used to solve the Kitaev model  [15]. Here I will introduce the method Kitaev used in the original paper as this forms the basis for the results that will be presented in this thesis. Rather than mapping to \(N\) fermions, Kitaev maps to \(4N\) Majoranas, effectively \(2N\) fermions. In contrast to the Jordan-Wigner transformation which makes fermions out of strings of spin operators in order to correctly produce fermionic commutation relations, the Kitaev transformation maps each spin locally to four Majoranas. The downside is that this enlarges the Hilbert space from \(2^N\) to \(4^N\). We will have to employ a projector \(\hat{P}\) to come back down to the physical Hilbert space later. As everything is local, I will drop the site indices \(ijk\) in expressions that refer to only a single site.

The mapping is defined in terms of four Majoranas per site \(b_i^x,\;b_i^y,\;b_i^z,\;c_i\) such that

\[\tilde{\sigma}^x = i b^x c,\; \tilde{\sigma}^y = i b^y c,\; \tilde{\sigma}^z = i b^z c.\qquad{(2)}\]

The tildes on the spin operators \(\tilde{\sigma_i}^\alpha\) emphasise that they live in this new extended Hilbert space and are only equivalent to the original spin operators after applying a projector \(\hat{P}\). The form of the projection operator can be understood in a few ways. From a group-theoretic perspective, before projection, the operators \(\{\tilde{\sigma}^x, \tilde{\sigma}^y, \tilde{\sigma}^z\}\) form a representation of the gamma group \(G_{3,0}\). The gamma groups \(G_{p,q}\) have \(p\) generators that square to the identity and \(q\) that square (roughly) to \(-1\). The generators otherwise obey standard anticommutation relations. The well-known gamma matrices \(\{\gamma^0, \gamma^1, \gamma^2, \gamma^3\}\) represent \(G_{1,3}\) the quaternions \(G_{0,3}\) and the Pauli matrices \(G_{3,0}\).

The Pauli matrices, however, have the additional property that the chiral element \(\sigma^x \sigma^y \sigma^z = \pm i\) is not fully determined by the group properties of \(G_{3,0}\), but it is equal to \(i\) in the Pauli matrices. Therefore, to fully reproduce the algebra of the Pauli matrices, we must project into the subspace where \(\tilde{\sigma}^x \tilde{\sigma}^y \tilde{\sigma}^z = +i\). The chiral element of the gamma matrices for instance \(\gamma_5 = i\gamma^0 \gamma^1 \gamma^2 \gamma^3\) is of central importance in quantum field theory. See ref.  [16] for more discussion of this group theoretic view.

So the projector must project onto the subspace where \(\tilde \sigma^x \tilde \sigma^y \tilde \sigma^z = i\). If we work this through, we find that in general \(\tilde \sigma^x \tilde \sigma^y \tilde\sigma^z = iD\) where \(D = b^x b^y b^z c\) must be the identity for every site. In other words, we can only work with physical states \(|\phi\rangle\) that satisfy \(D_i|\phi\rangle = |\phi\rangle\) for all sites \(i\). From this we construct an on-site projector \(P_i = \frac{1 + D_i}{2}\) and the overall projector is simply \(P = \prod_i P_i\).

Another way to see what this is doing physically is to explicitly construct the two intermediate fermionic operators \(f\) and \(g\) that give rise to these four Majoranas. Denoting a fermion state by \(|n_f, n_g\rangle\) the Hilbert space is the set \(\{|00\rangle,|01\rangle,|10\rangle,|11\rangle\}\). We can map these to Majoranas with, for example, this definition

\[\begin{aligned} b^x = (f + f^\dagger),\;\;& b^y = -i(f - f^\dagger),\\ b^z = (g + g^\dagger),\;\;& c = -i(g - g^\dagger). \end{aligned}\]

Working through the algebra, we see that the operator \(D = b^x b^y b^z c\) is equal to the fermion parity \(D = -(2n_f - 1)(2n_g - 1) = \pm1\) where \(n_f,\; n_g\) are the number operators. So setting \(D = 1\) everywhere is equivalent to restricting to the states \(\{|01\rangle\) and \(|10\rangle\}\), though we could equally well have used \(D = -1\).

Expanding the product \(\prod_i P_i\) out, we find that the projector corresponds to a symmetrisation over \(\{u_{ij}\}\) states within a flux sector and overall fermion parity \(\prod_i D_i\), see ref.  [17] or appendix A.6 for the full derivation. The significance of this is that an arbitrary many-body state can be made to have non-zero overlap with the physical subspace via the addition or removal of just a single fermion. This implies that, in the thermodynamic limit, the projection step is not generally necessary to extract physical results

We can now rewrite the spin hamiltonian in Majorana form with the caveat that they are only strictly equivalent after projection. The Ising interactions \(\sigma_j^{\alpha}\sigma_k^{\alpha}\) decouple into the form \(-i (i b^\alpha_i b^\alpha_j) c_i c_j\). We factor out the bond operators \(\hat{u}_{ij} = i b^\alpha_i b^\alpha_j\) which are Hermitian and, remarkably, commute with the Hamiltonian and each other.

\[\begin{aligned} \tilde{H} &= - \sum_{\langle i,j\rangle_\alpha} J^{\alpha}\tilde{\sigma}_i^{\alpha}\tilde{\sigma}_j^{\alpha}\\ &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} \hat{u}_{ij} \hat{c}_i \hat{c}_j. \end{aligned}\]

The bond operators \(\hat{u}_{ij}\) square to one so have eigenvalues \(\pm1\). As they are conserved we will work in their eigenbasis and take off the hats in the Hamiltonian.

\[\begin{aligned} H &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j. \end{aligned}\qquad{(3)}\]

The Fermion Problem

We now have a quadratic Hamiltonian, eq. 3, coupled to a classical field \(u_{ij}\). What follows is relatively standard theory for quadratic Hamiltonians  [18].

Because of the antisymmetry \(J^{\alpha} u_{ij}\), the eigenvalues of eq. 3 come in pairs \(\pm \epsilon_m\). We organise the eigenmodes of \(H\) into pairs, such that \(b_m\) and \(b_m'\) have energies \(\epsilon_m\) and \(-\epsilon_m\). The transformation \(Q\)

\[(c_1, c_2... c_{2N}) Q = (b_1, b_1', b_2, b_2' ... b_{N}, b_{N}'),\]

puts the Hamiltonian into normal mode form

\[H = \frac{i}{2} \sum_m \epsilon_m b_m b_m'.\]

The determinant of \(Q\) appears when evaluating the projector explicitly, otherwise, the \(b_m\) are merely an intermediate step. From them, we form fermionic operators

\[ f_i = \tfrac{1}{2} (b_m + ib_m'),\]

with their associated number operators \(n_i = f^\dagger_i f_i\). These let us write the Hamiltonian neatly as

\[ H = \sum_m \epsilon_m (n_m - \tfrac{1}{2}).\]

The energy of the ground state \(|n_m = 0\rangle\) of the many-body system at fixed \(\{u_{ij}\}\) is

\[E_{0} = -\frac{1}{2}\sum_m \epsilon_m, \]

and we can construct any state from a particular choice of \(n_m = 0,1\). If we only care about the ground state energy \(E_{0}\), it is possible to skip forming the fermionic operators. The eigenvalues obtained directly from diagonalising \(J^{\alpha} u_{ij}\) come in \(\pm \epsilon_m\) pairs. We can take half the absolute value of the set to recover \(\sum_m \epsilon_m\) directly.

An Emergent Gauge Field

We have transformed the spin Hamiltonian into a Majorana hamiltonian \(H = i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j\) describing the dynamics of a classical field \(u_{ij}\) and Majoranas \(c_i\). It is natural to ask how the classical field \(u_{ij}\) relates to the fluxes of the original spin model. We can evaluate the fluxes \(\phi_i\) in terms of the bond operators

\[\phi_i = \prod_{\langle j,k\rangle \in \mathcal{P}_i} i u_{jk}.\qquad{(4)}\]

Figure 4: A honeycomb lattice with edges in grey, along with its dual, the triangle lattice in red. The vertices of the dual lattice are the faces of the original lattice and, hence, are the locations of the vortices. (Left) The action of the gauge operator D_j at a vertex is to flip the value of the three u_{jk} variables (black lines) surrounding site j. The corresponding edges of the dual lattice (red lines) form a closed triangle. (Middle) Composing multiple adjacent D_j operators produces a large, closed dual loop or multiple disconnected dual loops. Dual loops are not directed like Wilson loops. (Right) A non-contractable loop which cannot be produced by composing D_j operators. All three operators can be thought of as the action of a vortex-vortex pair that is created, one of them is transported around the loop, and then the two annihilate again. Note that every plaquette has an even number of u_{ij}s flipped on its edge. Therefore, all retain the same flux \phi_i.

In addition, the bond operators form a highly degenerate description of the system. The operators \(D_i = b^x_i b^y_i b^z_i c_i\) commute with \(H\) forming a set of local symmetries. The action of \(D_i\) on a state is to flip the values of the three \(u_{ij}\) bonds that connect to site \(i\). This changes the bond configuration \(\{u_{ij}\}\) but leaves the flux configuration \(\{\phi_i\}\) unchanged. Physically, we interpret \(u_{ij}\) as a gauge field with a high degree of degeneracy and \(\{D_i\}\) as the set of gauge symmetries. Products of the gauge symmetries correspond to closed loops on the dual lattice, see fig. 4. The Majorana bond operators \(u_{ij}\) are an emergent, classical, \(\mathbb{Z}_2\) gauge field! The flux configuration \(\{\phi_i\}\) is what encodes physical information about the system without all the gauge degeneracy.

The ground state of the KH model is the flux configuration where all fluxes are one \(\{\phi_i = +1\; \forall \; i\}\). This can be proven via Lieb’s theorem  [19] which gives the lowest energy magnetic flux configuration for a system of electrons hopping in a magnetic field. Kitaev remarks in his original paper that he was not initially aware of the relevance of Lieb’s 1994 result. This is not surprising because at first glance the two models seem quite different. Yet the connection is quite instructive for understanding the KH and its generalisations.

Lieb discussed a model of mobile electrons

\[H = \sum_{ij} t_{ij} c^\dagger_i c_j,\]

where the hopping terms \(t_{ij} = |t_{ij}|\exp(i\theta_{ij})\) incorporate Aharanhov-Bohm (AB) phases  [20] \(\theta_{ij}\). The AB phases model the effect of a slowly varying magnetic field on the electrons through the integral of the magnetic vector potential

\[\theta_{ij} = \int_i^j \vec{A} \cdot d\vec{l},\]

a Peierls substitution  [21]. If we map the Majorana form of the Kitaev model to Lieb’s model, we see that our \(t_{ij} = i J^\alpha u_{ij}\). The \(i u_{ij} = \pm i\) correspond to AB phases \(\theta_{ij} = \pi/2\) or \(3\pi/2\) along each bond.

Figure 5: Lieb’s theorem gives the ground state flux configuration for even sided plaquettes in systems with at least one translationally invariant direction  [19]. These labels correspond to the Kitaev fluxes \phi_i rather than the magnetic fluxes Q_i of Lieb’s original paper (\phi_i = \exp{iQ_i}). Other work has extended Lieb’s theorem numerically to arbitrary plaquettes  [22–25]. The additional twofold degeneracy of the \pm i, \mp i terms is a consequence of the odd sided plaquettes breaking chiral symmetry. Chiral symmetry is spontaneously broken in the ground state  [23].

Stokes’ theorem tells us that the magnetic flux \(Q_m\) through a surface \(S\) is related to the line integral of \(\vec{A}\) along the boundary of the surface \(\partial S\), which in the discrete case reduces to the sum of the AB phases along the path. We can thus rewrite eq. 4 as

\[\begin{aligned} \phi_i &= \prod_{\mathcal{P}_i} i u_{jk}\\ &= \prod_{\mathcal{P}_i} \exp{(i\theta_{jk}})\\ &= \exp \left( i \sum_{\mathcal{P}_i} \theta_{jk} \right)\\ &= \exp \left( i \oint_{\mathcal{P}_i} \vec{A} \cdot d\vec{l} \right)\\ &= \exp \left( i Q_i \right). \end{aligned}\qquad{(5)}\]

Thus, we can interpret the fluxes \(\phi_i\) as the exponential of magnetic fluxes \(Q_m\) of some fictitious gauge field \(\vec{A}\) and the bond operators as \(i u_{ij} = \exp i \int_i^j \vec{A} \cdot d\vec{l}\). In this analogy to classical electromagnetism, the sets \(\{u_{ij}\}\) that correspond to the same \(\{\phi_i\}\) are all gauge equivalent as we have already seen via other means. The fact that fluxes can be written as products of bond operators and composed is a consequence of eq. 5. If the lattice contains odd plaquettes, as in the Yao-Kivelson model  [26], the complex fluxes that appear are a sign that chiral symmetry has been broken.

In full, Lieb’s theorem states that the ground state has magnetic flux \(Q_i = \sum_{\mathcal{P}_i}\theta_{ij} = \pi \; (\mathrm{mod} \;2\pi)\) for plaquettes with \(0 \; (\mathrm{mod}\;4)\) sides and \(0 \; (\mathrm{mod}\;2\pi)\) for plaquettes with \(2 \; (\mathrm{mod}\;4)\) sides. In terms of our fluxes, this means \(\phi = -1\) for squares, \(\phi = 1\) for hexagons and so on.

While Lieb’s theorem is restricted to bipartite lattices with translational symmetry, other works have shown numerically that it tends to hold for more general lattices too  [2225]. From this we find that the generalisation to odd sided plaquettes is similar but with an additional chiral symmetry. \(\phi = \pm i\) for plaquettes with \(1 \; (\mathrm{mod}\;4)\) sides and \(\mp i\) for those with \(3 \; (\mathrm{mod}\;4)\) sides. Overall we can write

\[\phi = -(\pm i)^{n_{\mathrm{sides}}}.\qquad{(6)}\]

Later I will present numerical evidence that this rule continues to hold for general amorphous lattices.

Understanding \(u_{ij}\) as a gauge field provides another way to understand the action of the projector. The local projector \(P_i = \frac{1 + D_i}{2}\) applied to a state constructs a superposition of the original state and the gauge equivalent state linked to it by flipping the three \(u_{ij}\) around site \(i\). The overall projector \(P = \prod_i P_i\) can thus be understood as a symmetrisation over all gauge equivalent states, removing the gauge degeneracy introduced by the mapping from spins to Majoranas.

Figure 6: Wilson loops that wind the major or minor diameters of the torus measure flux winding through the hole of the doughnut/torus (red) or through the filling (green). If they made doughnuts which had both had a jam filling and a hole, this analogy would be a lot easier to make  [27].

A final but important point to mention is that the local fluxes \(\phi_i\) are not quite all there is. We’ve seen that products of \(\phi_i\) can be used to construct the flux associated with arbitrary contractible loops. On the plane contractible loops are all there is. However, on the torus we can construct two global fluxes \(\Phi_x\) and \(\Phi_y\) which correspond to paths tracing the major and minor axes. The four sectors spanned by the \(\pm1\) values of these fluxes are gapped away from one another, but only by virtual tunnelling processes so the gap decays exponentially with system size  [1]. Physically \(\Phi_x\) and \(\Phi_y\) could be thought of as measuring the flux that threads through the hole of the doughnut. In general, surfaces with genus \(g\) have \(g\) ‘handles’ and \(2g\) of these global fluxes. At first glance, it may seem that these fluxes would not have much relevance to physical realisations of the Kitaev model which are likely to have a planar geometry. However, these fluxes are closely linked to topology and the existence of anyonic quasiparticle excitations in the model, which we will discuss next.

Anyons, Topology and the Chern number

Figure 7: Worldlines of particles in 2D can become tangled or braided with one another.

To discuss different ground state phases of the KH model, we must first review the topic of anyons and topology. The standard argument for the existence of fermions and bosons goes like this: the quantum state of a system must pick up a factor of \(\pm1\) if two identical particles are swapped. Only \(\pm1\) are allowed since swapping twice must correspond to the identity. This argument works in 3D for states without topological degeneracy, which seems to be true of the real world, but condensed matter systems are subject to no such constraints.

In gapped condensed matter systems, all equal time correlators decay exponentially with distance  [28]. Put another way, gapped systems support quasiparticles with a definite location in space and finite extent. As such it is meaningful to consider what would happen to the overall quantum state if we were to adiabatically carry out a series of swaps as described above. This is known as braiding. Recently, braiding in topological systems has attracted interest because of proposals to use ground state degeneracy to implement both passively fault tolerant and actively stabilised quantum computations  [2931].

First we realise that, in 2D, swapping identical particles twice is not topologically equivalent to the identity, see fig. 7. Instead it corresponds to encircling one particle around the other. This means we can in general pick up any complex phase \(e^{i\theta}\) upon exchange, hence the name any-ons. Those that pick up a complex phase are known as Abelian anyons because complex multiplication commutes and hence the group of braiding operations on them forms an Abelian group.

The KH model has a topologically degenerate ground state with sectors labelled by the values of the topological fluxes \((\Phi_x\), \(\Phi_y)\). Consider the operation in which a quasiparticle pair is created from the ground state, transported around one of the non-contractible loops, and then annihilated together, call them \(\mathcal{T}_{x}\) and \(\mathcal{T}_{y}\). These operations move the system around within the ground state manifold and they need not commute. This leads to non-Abelian anyons. As Kitaev pointed out, these operations are not specific to the torus: the operation \(\mathcal{T}_{x}\mathcal{T}_{y}\mathcal{T}_{x}^{-1}\mathcal{T}_{y}^{-1}\) corresponds to an operation in which none of the particles crosses the torus, rather one simply winds around the other. Hence, these effects are relevant even for the planar case!

In condensed matter systems, the existence of anyonic excitations automatically implies that the system has topological ground state degeneracy on the torus  [32] and indeed anyons and topology are intimately linked  [3335]. The Chern number \(\nu\), was a concept originally used to describe complex vector bundles in algebraic topology  [36]. It has now found use in physics as a powerful diagnostic tool for topological systems. Kitaev showed that there is a full classification of anyonic statistics in terms of \(\nu\;(\mathrm{mod}\;16)\), but relevant to us is that vortices in the KH model are Abelian when the Chern number is even and non-Abelian when the Chern number is odd. Non-Abelian statistics of the vortices in the KH model arise due to unpaired Majorana modes that are bound to them.

Ground State Phases

Figure 8: Setting the energy scale of the KH model with the constraint that J_x + J_y + J_z = 1 yields a triangular phase diagram where each of the corners represents J_\alpha = 1. For each corner \alpha the region |J_\alpha > |J_\beta| + |J_\gamma| supports a gapped non-Abelian phase equivalent to that of the Toric code  [29,37]. The point around equal coupling J_x = J_y = J_z, the B phase, is gapless. The B phase is known as a Majorana metal and on the honeycomb lattice it has a Dirac cone dispersion similar to that of graphene.

Setting the overall energy scale with the constraint \(J_x + J_y + J_z = 1\) yields a triangular phase diagram. In each of the corners one of the spin-coupling directions dominates, \(|J_\alpha > |J_\beta| + |J_\gamma|\), yielding three equivalent \(A_\alpha\) phases while the central triangle around \(J_x = J_y = J_z\) is called the B phase. Both phases support two kinds of quasiparticles, fermions and \(\mathbb{Z}_2\)-vortices. In the A phases, the vortices have bosonic statistics with respect to themselves but act like fermions with respect to the fermions, hence they are Abelian anyons. This phase has the same anyonic structure as the Toric code  [29]. The B phase can be described as a semi-metal of the Majorana fermions  [5]. Since the B phase is gapless, the quasiparticles aren’t localised and so don’t have braiding statistics.

An external magnetic can be used to break chiral symmetry. The lowest order term which breaks chiral symmetry but retains the solvability of the model is the three spin term \[ \sum_{(i,j,k)} \sigma_i^{\alpha} \sigma_j^{\beta} \sigma_k^{\gamma}, \] where the sum \((i,j,k)\) runs over consecutive indices around plaquettes. The addition of this to the spin model leads to two bond terms in the corresponding Majorana model. The effect of breaking chiral symmetry is to open a gap in the B phase. The vortices of the gapped B phase are non-Abelian anyons. This phase has the same anyonic exchange statistics as \(p_x + ip_y\) superconductor  [38], the Moore-Read state for the \(\nu = 5/2\) fractional quantum Hall state  [39] and many other systems  [4044]. Collectively these systems have attracted interest as possible physical realisations for quantum computers whose operations are based on braiding operations.

At finite temperatures, recent work has shown that the KH model undergoes a transition to a thermal metal phase. Vortex disorder causes the fermion gap to fill up and the DOS has a characteristic logarithmic divergence at zero energy which can be understood from random matrix theory  [11].

To surmise, the KH model is remarkable because it combines three key properties. First, the form of the Hamiltonian can plausibly be realised by a real material. Candidate materials, such as \(\alpha\mathrm{-RuCl}_3\), are known to have sufficiently strong spin-orbit coupling and the correct lattice structure to behave according to the KH model with small corrections  [5,45]. Second, its ground state is the canonical example of the long sought after QSL state, its dynamical spin-spin correlation functions are zero beyond nearest neighbour separation  [46]. Its excitations are anyons, particles that can only exist in 2D that break the normal fermion/boson dichotomy. Third, and perhaps most importantly, this model is a rare many-body interacting quantum system that can be treated analytically. It is exactly solvable. We can explicitly write down its many-body ground states in terms of single particle states  [1]. The solubility of the KH model, like the FK model, comes about because the model has extensively many conserved degrees of freedom. These conserved quantities can be factored out as classical degrees of freedom, leaving behind a non-interacting quantum model that is easy to solve.

Next Section: Disorder and Localisation


A. Kitaev, Anyons in an Exactly Solved Model and Beyond, Annals of Physics 321, 2 (2006).
G. Jackeli and G. Khaliullin, Mott Insulators in the Strong Spin-Orbit Coupling Limit: From Heisenberg to a Quantum Compass and Kitaev Models, Physical Review Letters 102, 017205 (2009).
M. Hermanns, I. Kimchi, and J. Knolle, Physics of the Kitaev Model: Fractionalization, Dynamic Correlations, and Material Connections, Annual Review of Condensed Matter Physics 9, 17 (2018).
S. M. Winter, A. A. Tsirlin, M. Daghofer, J. van den Brink, Y. Singh, P. Gegenwart, and R. Valentí, Models and Materials for Generalized Kitaev Magnetism, J. Phys.: Condens. Matter 29, 493002 (2017).
S. Trebst and C. Hickey, Kitaev Materials, Physics Reports 950, 1 (2022).
H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and Realization of Kitaev Quantum Spin Liquids, Nat Rev Phys 1, 4 (2019).
S. Chandra, K. Ramola, and D. Dhar, Classical Heisenberg Spins on a Hexagonal Lattice with Kitaev Couplings, Phys. Rev. E 82, 031113 (2010).
E. Sela, H.-C. Jiang, M. H. Gerlach, and S. Trebst, Order-by-Disorder and Spin-Orbital Liquids in a Distorted Heisenberg-Kitaev Model, Phys. Rev. B 90, 035113 (2014).
P. W. Anderson, Resonating Valence Bonds: A New Kind of Insulator?, Materials Research Bulletin 8, 153 (1973).
M. Freedman, A. Kitaev, M. Larsen, and Z. Wang, Topological Quantum Computation, Bull. Amer. Math. Soc. 40, 31 (2003).
C. N. Self, J. Knolle, S. Iblisdir, and J. K. Pachos, Thermally Induced Metallic Phase in a Gapped Quantum Spin Liquid - a Monte Carlo Study of the Kitaev Model with Parity Projection, Phys. Rev. B 99, 045142 (2019).
G. Khaliullin, Orbital Order and Fluctuations in Mott Insulators, Progress of Theoretical Physics Supplement 160, 155 (2005).
J. Chaloupka, G. Jackeli, and G. Khaliullin, Kitaev-Heisenberg Model on a Honeycomb Lattice: Possible Exotic Phases in Iridium Oxides A₂IrO₃, Phys. Rev. Lett. 105, 027204 (2010).
M. Petitjean, Chirality of Dirac Spinors Revisited, Symmetry 12, 4 (2020).
F. L. Pedrocchi, S. Chesi, and D. Loss, Physical solutions of the Kitaev honeycomb model, Phys. Rev. B 84, 165414 (2011).
J. Blaizot, G. Ripka, J. W. Negele, H. Orland, and G. E. Brown, Quantum Theory of Finite Systems and Quantum Many‐Particle Systems, Physics Today 41, 106 (1988).
E. H. Lieb, Flux Phase of the Half-Filled Band, Physical Review Letters 73, 2158 (1994).
Y. Aharonov and D. Bohm, Significance of Electromagnetic Potentials in the Quantum Theory, Phys. Rev. 115, 485 (1959).
R. Peierls, Zur Theorie des Diamagnetismus von Leitungselektronen, Z. Physik 80, 763 (1933).
T. Eschmann, P. A. Mishchenko, K. O’Brien, T. A. Bojesen, Y. Kato, M. Hermanns, Y. Motome, and S. Trebst, Thermodynamic Classification of Three-Dimensional Kitaev Spin Liquids, Phys. Rev. B 102, 075125 (2020).
H. Yao, S.-C. Zhang, and S. A. Kivelson, Algebraic Spin Liquid in an Exactly Solvable Spin Model, Phys. Rev. Lett. 102, 217202 (2009).
T. Eschmann, P. A. Mishchenko, T. A. Bojesen, Y. Kato, M. Hermanns, Y. Motome, and S. Trebst, Thermodynamics of a Gauge-Frustrated Kitaev Spin Liquid, Phys. Rev. Research 1, 032011 (2019).
V. Peri, S. Ok, S. S. Tsirkin, T. Neupert, G. Baskaran, M. Greiter, R. Moessner, and R. Thomale, Non-Abelian Chiral Spin Liquid on a Simple Non-Archimedean Lattice, Phys. Rev. B 101, 041114 (2020).
H. Yao and S. A. Kivelson, An Exact Chiral Spin Liquid with Non-Abelian Anyons, Phys. Rev. Lett. 99, 247203 (2007).
M. B. Hastings, Lieb-Schultz-Mattis in Higher Dimensions, Phys. Rev. B 69, 104431 (2004).
A. Yu. Kitaev, Fault-Tolerant Quantum Computation by Anyons, Annals of Physics 303, 2 (2003).
D. Poulin, Stabilizer Formalism for Operator Quantum Error Correction, Phys. Rev. Lett. 95, 230504 (2005).
M. B. Hastings and J. Haah, Dynamically Generated Logical Qubits, Quantum 5, 564 (2021).
T. Einarsson, Fractional Statistics on a Torus, Phys. Rev. Lett. 64, 1995 (1990).
M. Oshikawa, Y. B. Kim, K. Shtengel, C. Nayak, and S. Tewari, Topological Degeneracy of Non-Abelian States for Dummies, Annals of Physics 322, 1477 (2007).
S. B. Chung, H. Yao, T. L. Hughes, and E.-A. Kim, Topological Quantum Phase Transition in an Exactly Solvable Model of a Chiral Spin Liquid at Finite Temperature, Phys. Rev. B 81, 060403 (2010).
H. Yao, S.-C. Zhang, and S. A. Kivelson, Algebraic Spin Liquid in an Exactly Solvable Spin Model, Phys. Rev. Lett. 102, 217202 (2009).
S. Chern, Characteristic Classes of Hermitian Manifolds, Annals of Mathematics 47, 85 (1946).
A. Y. Kitaev, Quantum Communication, Computing, and Measurement, in Proceedings of the 3rd International Conference of Quantum Communication and Measurement, New York: Plenum (1997).
G. Moore and N. Read, Nonabelions in the Fractional Quantum Hall Effect, Nuclear Physics B 360, 362 (1991).
J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Non-Abelian Statistics and Topological Quantum Information Processing in 1d Wire Networks, Nature Phys 7, 5 (2011).
R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana Fermions and a Topological Phase Transition in Semiconductor-Superconductor Heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
Y. Oreg, G. Refael, and F. von Oppen, Helical Liquids and Majorana Bound States in Quantum Wires, Phys. Rev. Lett. 105, 177002 (2010).
J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Generic New Platform for Topological Quantum Computation Using Semiconductor Heterostructures, Phys. Rev. Lett. 104, 040502 (2010).
A. Banerjee et al., Proximate Kitaev Quantum Spin Liquid Behaviour in {\Alpha}-RuCl$_3$, Nature Mater 15, 733 (2016).
G. Baskaran, S. Mandal, and R. Shankar, Exact Results for Spin Dynamics and Fractionalization in the Kitaev Model, Phys. Rev. Lett. 98, 247201 (2007).