2 Background
The Kitaev Honeycomb Model
The Spin Hamiltonian
This section introduces the Kitaev honeycomb (KH) model. The 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 that bond couple 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 [2–6]. The model is highly frustrated: each spin would like to align along a different direction with each of its three neighbours. This cannot be achieved even classically [7,8]. This frustration leads the 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 one of the first concrete 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. We will discuss the properties of the emergent gauge fields and the projector. We 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 we will look at the phase diagram.
The next section will discuss anyons, topology and the Chern number, using the Kitaev model as an explicit example of these topics.
The Spin Model
As discussed in the introduction, spin hamiltonians like that of the Kitaev 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 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 Kitaev 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 Falikov-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.
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 Falikov-Kimball model, the KH model has a extensive number of conserved quantities, the fluxes. As with the FK model it will make sense to 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 picture.
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 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 approach 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}\) emphasis 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 = i\), this relation is not determined by the group properties of \(G_{3,0}\). 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 [16] for more discussion of this group theoretic view.
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 $^x ^y ^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|= |$ 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. Working through the algebra we see that the \(D\) operator corresponds to the fermion parity \(D = -(2n_f - 1)(2n_g - 1)\) where \(n_f,\; n_g\) are the number operators. 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 and overall fermion parity \(\prod_i D_i\). This tells us that any arbitrary state can be made to have non-zero overlap with the physical subspace via the addition or removal of a single fermion. This implies that in the thermodynamic limit the projection step is not generally necessary to extract physical results, see [17] or appendix A.5 for more details.
We can now rewrite the spin hamiltonian in Majorana form with 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’re 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 \[\tilde{H}_u = \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
\[ \tilde{H}_u = \sum_m \epsilon_m (n_m - \tfrac{1}{2}).\]
The ground state \(|n_m = 0\rangle\) of the many body system at fixed \(\{u_{ij}\}\) is
\[E_{u,0} = -\frac{1}{2}\sum_m \epsilon_m \]
and we can construct any state from a particular choice of \(n_m = 0,1\).
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)}\]
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\) so form 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. 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.
The ground state of the Kitaev Honeycomb model has flux configuration \(\{\phi_i = +1\}\). This can be proven via Lieb’s theorem [19] which gives the 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 but the connection is quite instructive to understanding the Kitaev Model and its generalisations.
Lieb discusses a model of mobile electrons
\[H = \sum_{ij} t_{ij} c^\dagger_i c_j\]
where the hopping term \(t_{ij} = |t_{ij}|\exp(i\theta_{ij})\) incorporate Aharanhov-Bohm [20] (AB) phases \(\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 the \(u_{ij} = \pm 1\) correspond to AB phases \(\theta_{ij} = 0,\pi\) along each bond. Stokes’ theorem tells us that the magnetic flux \(\tilde{\phi}_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 &= i^n \prod_{\mathcal{P}_i} \exp{(i\theta_{jk}})\\ &= i^n \exp \left( i \sum_{\mathcal{P}_i} \theta_{jk} \right)\\ &= i^n \exp \left( i \oint_{\mathcal{P}_i} \vec{A} \cdot d\vec{l} \right) \end{aligned}\qquad{(5)}\]
Thus we can interpret the fluxes \(\phi_i\) as the exponential of magnetic fluxes of some fictitious gauge field \(\vec{A}\) and the bond operators as \(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. The fact that fluxes can be written as products of bond operators and composed is a consequence of (the exponential of) Stokes’ theorem. The additional phase factors of \(i^n\) can be incorporated as additional \(\tfrac{\pi}{2}\) phases but then make little difference when all the plaquettes are hexagons or have an even number of sides. However if the lattice contains odd plaquettes, as in the Yao-Kivelson model [22], the complex fluxes that appear are a sign that chiral symmetry has been broken.
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.
A final but important point to mention is that 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 are. 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 this would not have much relevance to physical realisations of the Kitaev model that will likely have a planar geometry with open boundary conditions. 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

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 three dimensions 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 [24]. Put another way, the system supports quasiparticles with a definite location in space and a 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.
First we realise that in two dimensions, swapping identical particles twice is not topologically equivalent to the identity, see fig. 5. Instead it corresponds to encircling one particle around the other. This means we can in general pick up any complex phase \(e^{i\theta}\), hence the name any-ons upon exchange. These are known as Abelian anyons because complex multiplication commutes and hence the group of braiding operations forms and 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 us around within the ground state manifold and they need not commute. This leads to non-Abelian anyons. As Kitaev points 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 of relevant even for the planar case.

In condensed matter systems, the existence of anyonic excitations automatically implies the system has topological ground state degeneracy on the torus [25] and indeed anyons and topology are intimately linked [26–28]. Originally a concept used to describe complex vector bundles in algebraic topology [29], the Chern number has found use in physics as a powerful tool diagnostic tool for topological systems. Kitaev showed that vortices in the KH model are Abelian when the Chern number is even and non-Abelian when the Chern number is odd. In the odd case the non-Abelian statistics of the vortices arise due to unpaired Majorana modes that are bound to them.
Recently, topological systems have gained interest because of proposals to use their ground state degeneracy to implement both passively fault tolerant and actively stabilised quantum computations [30–32].
Ground State Phases
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 [30]. 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 that 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.
At finite temperatures, recent work has shown that the KH model undergoes a transition to a thermal metal phase.
Summary
We have seen that…
Summary The Kitaev Honeycomb model is remarkable because it combines three key properties. First, the form of the Hamiltonian 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 Kitaev Honeycomb model with small corrections [5,34]. Second, its ground state is the canonical example of the long sought after quantum spin liquid state. Its excitations are anyons, particles that can only exist in two dimensions 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 Kitaev Honeycomb Model, like the Falikov-Kimball model of chapter 1, 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.
“dynamical two spin correlation functions are identically zero beyond nearest neighbor separation in the Kitaev Model” [35]
- RuCl_3 is the classic QSL candidate material
- really follows the Kitaev-Heisenberg model
- experimental probes include inelastic neutron scattering, Raman scattering
Next Section: Disorder and Localisation