# Bound states of fractionalized excitations in a modulated Kitaev spin liquid

## Abstract

Fractionalization is a hallmark of spin-liquid behavior; it typically leads to response functions consisting of continua instead of sharp modes. However, microscopic processes can enable the formation of short-distance bound states of fractionalized excitations, despite asymptotic deconfinement. Here we study such bound-state formation for the spin liquid realized in Kitaev’s honeycomb compass model, supplemented by a kekulé distortion of the lattice. Bound states between flux pairs and Majorana fermions form in the Majorana band gaps. We calculate the dynamic spin susceptibility and show that bound states lead to sharp modes in the magnetic response of the spin liquid, with the momentum dependence of the corresponding spectral weight encoding internal symmetry of the bound state. As a byproduct, we also show that isolated fluxes may produce Majorana bound states at exactly zero energy. Generalizations and implications of the results are discussed.

###### pacs:

## I Introduction

Strongly frustrated magnets can host a variety of unconventional states of matter.(1); (2) Particularly fascinating are spin liquids and their descendants. Spin liquids are strongly entangled low-temperature states of local moments which do not display magnetic order and do not break any symmetries of the underlying lattice.(3)

Typically, spin liquids display the phenomenon of fractionalization: Elementary excitations come with quantum numbers which cannot be constructed from electrons or holes, and local operators can create only multiples (often pairs) of such deconfined excitations, such that locally created excitations decay into fractionalized constituents. As a result, measurable response functions, like the dynamic spin structure factor, display multi-particle continua instead of sharp modes. In fact, the absence of sharp modes is often taken as experimental evidence for spin-liquid behavior.(4)

However, this reasoning can be spoiled if fractionalized excitations form short-distance bounds states; this is possible despite long-distance deconfinement. Such bound states may arise from attractive interactions between the fractionalized excitations. Hence, it is conceivable to have a deconfined spin liquid whose dominant response arises from bound states and consists of sharp-mode peaks. This highly unusual situation has been studied relatively little in the theory literature,(5); (6) but may be relevant for real materials. For instance, spin-liquid physics has been proposed to be relevant for underdoped cuprates,(7); (8) but clear-cut evidence for fractionalization is missing, and neutron scattering experiments are typically interpreted in terms of sharp modes.(9); (10); (11)

It is the purpose of this paper to study the formation of bound states of emergent fractionalized spin-liquid constituents in a solvable setting. To this end, we utilize Kitaev’s compass model on the honeycomb lattice which, in its pristine form, realizes a spin liquid.(12) Its elementary excitations are gapless dispersive “matter” Majorana fermions and gapped immobile gauge fluxes (or “visons”). By now, a number of honeycomb-lattice magnets have emerged as candidates for realizing dominant Kitaev interactions,(13) most prominently NaIrO,(14); (15); (16) different polytypes of LiIrO,(16); (17) and -RuCl.(18); (19) These materials display long-range magnetic order at low temperatures, likely due to the presence of additional Heisenberg interactions,(13); (20); (21); (24); (22); (23) but it has been suggested that pressure or doping may be used to suppress magnetic order and to induce spin-liquid behavior.

Here we employ the Kitaev model as a playground for bound-state formation.(6) As well-defined bound states require the existence of energy gaps in the excitation spectrum, we induce dispersion gaps by imposing a kekulé-type modulation (25) of coupling constants in the Kitaev model, i.e., we consider a lattice with a periodic superstructure modulation. We show that bound states between matter fermions and isolated fluxes as well as between matter fermions and flux pairs exist in the resulting spin liquid. These bound states are spatially localized (because the fluxes are localized) and come with different symmetries w.r.t. point-group operations. We calculate the dynamic spin structure factor and show that bound states involving flux pairs are visible spectroscopically; in fact, the sharp-mode response of the bound states may dominate the spectrum despite deconfinement. The momentum dependence of the sharp-mode spectral weight allows one to deduce the spatial symmetries of the bound states. We discuss how to detect deconfinement despite low-energy bound-state formation, and we highlight implications for more realistic models and for carrier-doped spin liquids.

We note that bound states between matter fermions and flux excitations were discussed before for the field-induced non-abelian phase of the Kitaev model.(5); (26)

The body of the paper is organized as follows: In Section II we introduce the kekulé-modulated Kitaev model together with its Majorana representation. Section III summarizes the excitation spectrum in the flux-free case. Section IV demonstrates the bound-state formation on the level of Majorana-fermion states while Section V presents the results for bound-state spectroscopy utilizing the dynamical spin susceptibility. Broader implications are discussed in the summary section.

## Ii Kekulé-modulated Kitaev model

### ii.1 Kitaev model

The Kitaev model(12) describes spin-1/2 degrees of freedom at sites of a honeycomb lattice which interact via Ising-like nearest-neighbor exchange interactions . The anisotropy direction in spin space, , is coupled to the bond direction in real space, often dubbed “compass interaction”.(27) Allowing for spatially varying couplings the Hamiltonian reads

(1) |

where are Pauli matrices, and denotes an bond as in Fig. 1. In the homogeneous case , , . For isotropic couplings, , the model then possesses a symmetry of combined real-space and spin rotations.

The Kitaev model is exactly solvable thanks to the existence of an extensive set of conserved quantities, corresponding to fluxes which can be defined for every closed loop on the lattice. For each elementary plaquette of the lattice, the corresponding loop operators can be written as

(2) |

with eigenvalues . For periodic boundary conditions there are, in addition, two “topological” loop operators that wrap around the torus. Flux conservation implies that the Hilbert space of the Kitaev model can be divided into flux sectors, specified by the set of . The ground state is located in the flux-free sector;(12); (28) this remains true in the presence of the kekulé modulation considered below.

### ii.2 Kekulé modulation

The solubility of the Kitaev model relies only on threefold lattice coordination, together with the spin structure of the interactions. Hence, variants of the Kitaev model with vacancies,(29) random or inhomogeneous couplings,(31); (30) or different lattice geometries, both in two(32); (34); (33); (35) and three(36); (37); (38) space dimensions, remain solvable and have been discussed.

Here we impose a modulation of the exchange interactions according to a superstructure on the honeycomb lattice, with the goal to induce spectral gaps. This kekulé modulation, initially discussed in the context of carbon nanotubes,(25) triples the unit cell, but preserves the symmetry of combined real-space and spin rotations.(39)

The spatial pattern of exchange interactions in the kekulé-modulated Kitaev model is shown in Fig. 1. The modulation is parameterized by a ratio of coupling constants, , such that the exchange couplings obey

(3) |

Hence, corresponds to decoupled hexagons, is the original honeycomb lattice, and represents decoupled dimers. We will be mostly interested in the isotropic case, where we choose , i.e., the couplings on the h bonds, as energy unit.

The enlarged unit cell implies that there are two different types of plaquettes and hence two types of flux excitations; for reference we denote a plaquette formed by six h bonds as H plaquette, while a plaquette with three h and three i bonds is dubbed I plaquette. Flux pairs can be created by flipping either a h bond (creating one H and one I flux) or an i bond (creating two I fluxes), see Fig. 1.

### ii.3 Majorana representation

For explicit calculations we utilize Kitaev’s spin representation (12) with four Majorana fermions per site , , and . Defining , the original Hamiltonian in Eq. (1) can be mapped to

(4) |

where , with site on sublattice , and the summation is over all nearest-neighbor bonds. The represent conserved quantities, with eigenvalues , and a given set reduces the Hamiltonian to a bilinear in the Majorana operators:

(5) |

Here is an matrix with elements , with the number of unit cells, and is a vector of Majorana operators on the sublattice. Hence the problem takes the form of non-interacting “matter” Majorana fermions coupled to a static gauge field.

For every flux sector, described by a set of set of , the Majorana hopping model can be diagonalized on finite-size lattices using standard techniques, and we refer the reader to the literature.(12); (29); (31) The result is a model of free canonical fermions representing the matter degrees of freedom, with

(6) |

where the mode energies are non-negative. For the flux-free unmodulated isotropic Kitaev model, momentum takes the role of the quantum number , and the resemble the positive-energy part of the spectrum of graphene, with Dirac points at momenta .

## Iii Phase diagram and Majorana spectrum in the flux-free sector

We now turn to the properties of the kekulé-modulated Kitaev model. To set the stage, we discuss the spectrum and dispersion of the matter Majorana fermions in the flux-free sector. We allow for coupling anisotropies, and consider the case . In general, there are now three positive-energy bands of matter Majorana fermions due to the tripled unit cell.

The effect of a weak superstructure can be considered perturbatively. In the isotropic case, the kekulé modulation couples the two Dirac points which are gapped out, hence the spectrum is gapped for any . In contrast, the anisotropic Kitaev model displays gapless Dirac cones which are shifted in momentum space for(12) and . In these cases, the kekulé modulation does not couple the Dirac points, such that extended gapless phases remain for . We have computed the phase boundaries, and the resulting phase diagram is shown in Fig. 2. Remarkably, the gapless phases extend all the way to and .

In the following we focus on the isotropic case, where we show matter Majorana bandstructures for selected values in Fig. 3 and an overview of the spectrum in Fig. 4. As noted, all bands are gapped for . For the three Majorana bands overlap, whereas an inter-band gap opens for . In the limit of , the spectrum is that of decoupled hexagons, with three discrete local modes, two of which are degenerate, see Fig. 4. In contrast, the limit corresponds to isolated dimers, with a threefold degenerate local mode. Interestingly, at leading order in the spectrum corresponds to that of a Kagome lattice.

We note that the transition from to is accompanied by a band inversion, such that the band structure for is characterized by a non-trivial topological index.(40) This implies the existence of topologically protected matter Majorana edge states for .

## Iv Bound-state formation

Next we turn our attention to the physics in excited flux sectors. Here we demonstrate the existence of spatially localized states of matter Majorana fermions bound to static fluxes which act as “impurities”. These impurity bound states arise in one of the gaps of the Majorana spectrum.

### iv.1 Flux pair

We start with flux configurations consisting of a pair of fluxes on adjacent plaquettes; those are relevant for the zero-temperature spin response.(41); (42)

We have computed the matter Majorana spectrum by diagonalizing the hopping problem (5) on finite-size lattices with periodic boundary conditions. Bound-state formation is signalled by the existence of (one or more) isolated mode energies in energy windows corresponding to the band gaps of the flux-free case. For a single flipped bond, the same information can be obtained by analyzing the poles of a suitable T matrix.

We find matter Majorana bound states for any , and their energies are summarized in Fig. 4. For flipping a single h bond generates one bound state close to zero energy; in addition two (almost degenerate) bound states appear for within the upper band gap. In the limit these three bound states correspond to the states of an isolated hexagon threaded by a flux, with energies , , . No bound states are obtained by flipping a i bond. In contrast, for bound states are generated by flipping either an h or an i bond; in each case a single bound state emerges below the lower band edge.

The wavefunctions of the bound states are portrayed in Fig. 5. They can be classified according to their signature under reflection at a axis perpendicular to the flipped bond. (We recall that these wavefunctions are not gauge-invariant.) For the low-energy bound state is even, Fig. 5(a). The same applies to the lower of the elevated-energy bound states emerging for , Fig. 5(b), whereas the higher on is odd, Fig. 5(c). In all cases, the main weight is located on the H plaquette sharing the flipped bond. For the bound state obtained from an h flip is even, with the main weight on the I plaquette. Flipping an i bond leaves two mirror symmetries intact (along the bond axis and perpendicular to it), and the bound state is even under the first and odd under the second. Its main weight is equally distributed on the two I plaquettes (not shown).

For the computation of the spin structure factor (see below) we will also need the flux gap, i.e., the energy cost of flipping a single bond. The corresponding numerical results are shown in Fig. 6.

### iv.2 Single flux

To complete the physical picture, we now analyze the fermion spectrum in the presence of isolated fluxes. Given that periodic boundary conditions require the total number of plaquette fluxes to be even, i.e., the product of all from Eq. (2) to be , we study a configuration with two fluxes placed at a maximum distance of . Such configurations correspond to flipping a chain of bonds (“Dirac string”) as shown in Fig. 7.

Interestingly we find for that only isolated H fluxes generate bound states. The configurations shown in Fig. 7 have one H and one I flux, and we find a total of three bound states as in the flux-pair case. Their energetics is interesting: For the energies are again , , . For finite the system preserves the symmetry of the isolated hexagon, and consequently the effect of finite on the bound states is equivalent to a renormalization. As a result, the two upper bound states are degenerate, and the lower bound state is exactly at zero energy. Hence, an isolated H flux in the kekulé-modulated Kitaev model with provides a means to generate zero-energy Majorana bound states. Concerning the (gauge-dependent) wavefunctions, we see that lowest bound state one is now odd under reflection w.r.t. the Dirac-string axis. For a configuration with two isolated H fluxes, we find a total of six bound states (not shown).

In contrast, for we observe that isolated I fluxes generate one bound state each whereas H fluxes do not induce bound states. Hence, for a configuration with two distant I fluxes we find two bound states (not shown).

## V Bound-state spectroscopy

So far, we have demonstrated the presence of bound states formed by matter Majorana fermions and gauge fluxes or pairs thereof. The latter bound states are directly visible in spectroscopic probes, and we show this by determining the dynamic spin structure factor for the kekulé-modulated Kitaev model.

### v.1 Dynamical spin correlations

Dynamical spin correlations in the Kitaev model(41) have been explicitly calculated in Ref. (42); (5). Consider the zero-temperature spin correlation function

(7) |

where is the many-body ground state. The application of a operator creates a flux pair in the plaquettes which involve the bond emanating from site . This leads to the dynamical rearrangement of matter fermions in the modified gauge field akin to a quantum quench. The spin correlator can be expressed purely in terms of matter fermions in the ground-state flux sector, subject to a perturbation .(41); (42) For instance, the off-site correlator reads

(8) | ||||

where is non-zero only if and are nearest neighbors connected by an bond, i.e., vanishes beyond nearest neighbors. Further, is the matter-fermion ground state of in the flux-free sector.(43) runs over all matter-fermion states in the two-flux sector, i.e., the eigenstates of . and are the corresponding many-body energies. The prefactor depending on the spin component. Physically, Eq. (8) expresses that a spin-flip excitation decays into a matter Majorana fermion and a flux pair.

### v.2 Sources of sharp-mode peaks

Analyzing Eq. (8) shows that the results crucially depend on whether the ground states in the two flux sectors have (i) the same or (ii) different total fermion parity.(42); (5) In case (i) non-vanishing contributions to arise from excited states with an odd number of matter-fermion excitations while in case (ii) an even number is required.(43) As a result, there are two distinct sources for sharp-mode peaks in the dynamical spin response.(42); (5) In case (i), the leading contributions come from excited single-fermion states. Hence, reflects the matter-fermion density of states in the two-flux sector, and sharp peaks are caused by energetically isolated single-fermion states – these are precisely the bound states discussed in Section IV above. The peak energy in is then given by the sum of the flux gap and the fermion bound-state energy. In case (ii), the first contribution is from the zero-fermion state in the two-flux flux sector, and hence always yields a delta peak, located at the flux-gap energy. The next contributions come from two-fermion states, such that the presence of a single fermionic bound state does not produce a sharp-mode peak.

For the kekulé-modulated Kitaev model with isotropic couplings, we find for that the parity of the two-flux state matches that of the flux-free state both for flipped h and i bonds, i.e., case (i) is realized. In contrast, for this only applies to flipped h bonds, while flipping an i bond generates the parity mismatch of case (ii).

### v.3 Results

We now discuss the dynamic spin structure factor at momentum ,

(9) |

for the kekulé-modulated Kitaev model. We calculate this in a few-mode approximation which has been shown to yield highly accurate results for the standard Kitaev model. The explicit calculation in case (ii) requires to perform a local gauge transformation, and we refer the reader to the literature for details.(5)

Numerical results for at zero temperature are shown in Fig. 8. Two striking features are apparent: First, multiple sharp-mode peaks show up, with some of them even occuring on top of a continuum contribution (despite the fact that bound states are located in a band gap). This can be rationalized as follows: For all peaks in arise from the bound states discussed in Section IV.1, with all bound states occurring for flipped h bonds, Fig. 4. Importantly, involves contributions both from flipped h and flipped i bonds. Since the corresponding flux gaps are different, Fig. 6, the two contributions come with continua which are shifted by a different flux “offset”. As a result, a bound-state peak from the h-bond contribution may energetically overlap with the continuum from the i-bond contribution – this is what happens in Fig. 8(a,b). For the spin structure factor contains two sharp-mode peaks of different origin, Fig. 8(c). For a flipped h bond we obtain a bound-state peak, located very close to the lower edge of the continuum, while a flipped i bond leads to the parity mismatch and hence a zero-particle peak at the i flux gap (but no bound-state peak).

Second, there is little momentum dependence in . For the continuum contributions this is similar to the standard Kitaev model where it results from fractionalization and the fact that spin correlators vanish beyond nearest neighbors. In addition, the bound-state peaks are non-dispersive which is a simple consequence of the fluxes (and hence the bound states) being localized.

The weights of the sharp-mode peaks are generically momentum-dependent, as shown in Fig. 9. First, the weights can be large: For the peaks account for roughly of the total spectral weight in . Hence, the peaks can dominate the response which also becomes clear from the insets in Fig. 8. Second, the momentum dependence is determined by the internal symmetries of the bound state. This is nicely seen in Fig. 9(a) where the two bound states which are even under reflection have a finite weight for all momenta, whereas the odd bound state has vanishing weight at .

We can reasonably expect that the response at small non-zero temperatures has the same qualitative properties: The response function will involve a trace over initial states with both fermions and fluxes, but in all cases sizeable matrix elements are only obtained if the intermediate state contains, in addition to the excitation effectively created by , the same excitations as the initial state .(31) This implies that all thermal contributions to the low-temperature spin correlator will display the same bound-state peaks, and the main effect of temperature is that of thermal broadening.

## Vi Discussion and outlook

In this paper we have studied spin excitations of a kekulé-modulated Kitaev model. This model provides an explicit and exactly solvable example for a deconfined phase with fractionalized excitation whose dynamic response is nevertheless dominated by sharp-mode peaks which arise via bound-state formation.(6) We have provided detailed results for the dynamic spin structure factor which we have rationalized via an analysis of spatially localized Majorana fermion states in excited flux sectors.

### vi.1 Beyond the Kitaev model

We now turn to a broader view, and start with physics beyond the solvable Kitaev limit. The gapless spin liquid of the standard Kitaev model has been shown to be stable against small perturbations (like a Heisenberg interaction).(13); (20) Two qualitative modifications occur: The flux excitations (visons) become mobile, and spin operator acquires an additional decay channel into two Majorana fermions.(45) Upon perturbing the kekulé-modulated Kitaev model we expect – by continuity – that bound states between vison pairs and Majorana fermions continue to exist, but these bound states will themselves become dispersive. Their spatial structure will become more complicated – in particular, it involves also configurations with the two visons not being on adjacent plaquettes – but the bound states will continue to contribute sharp modes in the dynamic spin structure factor. Hence, the phenomenon studied in this paper is robust. In addition, it is also conceivable that bound states between pairs of matter Majorana fermions exist, but this requires a separate analysis which is beyond the scope of this paper. If the Majorana spectrum does not display spectral gaps, then true bound states can only occur above the upper band edge; alternatively quasi-bound states may exist inside the spectrum. However, neither of the two happen in the standard Kitaev model, .

For spin liquids with other fractionalization schemes similar bound state formation, e.g., spinon-vison or spinon-spinon, may occur depending on microscopic details, but explicit examples have not been studied to our knowledge. We note that bound-state formation is in fact expected upon approaching a confinement transition: Gauge-field fluctuations produce bound states on all energy scales on the confined side of the transition, and a large confinement length on the deconfined side implies the existence of low-energy bound states. This will studied in future work.

Finally, bound-state formation is relevant beyond insulators. For instance, in doped Mott insulators with deconfinement one may expect charge-neutral spin- spinons and spinless charge- holons as elementary excitations. Bound-state formation can yield objects with different quantum numbers. In particular, the idea of a fractionalized Fermi liquid (FL) phase(46) realized in a doped one-band Mott insulator implies the existence of charge- spin- particles which can be understood as bound states of spinons and holons.(8)

### vi.2 Detecting confinement

From an experimental point of view, the response from bound states in a deconfined phase is hard to distinguish from the response of a conventional, non-fractionalized phase. In the following we discuss a few aspects of this. First, a deconfined phase with bound states can be expected to show sharp peaks and continua in response functions, see Fig. 8. Notably, such response may also arise from a conventional phase where single-particle peaks and multi-particle continua can coexist. A careful analysis of quantum numbers (e.g. via the response to a magnetic field) may help to distinguish the two cases. Second, the combination of multiple observables typically provides additional information. For instance, for the Kitaev model it has been shown that Raman scattering mainly probes the spectrum in the ground-state flux sector(47) where no bound states occur in the case studied in this paper. Third, a convincing proof of deconfinement may require to show the existence of emergent gauge-field excitations, for instance via flux trapping in suitable geometries.(48)

###### Acknowledgements.

We thank L. Fritz, I. Göthel, J. Knolle, R. Moessner, S. Rachel, S. Ray, and F. Zschocke for discussions and collaborations on related work. This research was supported by the DFG through SFB 1143 and GRK 1621.### References

- A. P. Ramirez, Ann. Rev. Mat. Sci. 24, 453 (1994).
- O. A. Starykh, Rep. Prog. Phys. 78, 052502 (2015).
- L. Balents, Nature 464, 199 (2010).
- T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature 492, 406 (2012).
- J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moessner, Phys. Rev. B 92, 115127 (2015).
- The formation of probe-induced bound states in modifications of the Kitaev model has been noticed in Ref. (30) for the special case of triaxial strain.
- P. W. Anderson, Science 235, 1196 (1987).
- M. Punk, A. Allais, and S. Sachdev, PNAS 112, 9552 (2015); S. Sachdev and D. Chowdhury, Prog. Theor. Exp. Phys 12C102 (2016).
- R. J. Birgeneau, C. Stock, J. M. Tranquada, and K. Yamada, J. Phys. Soc. Jpn. 75, 111003 (2006).
- J. M. Tranquada, in: Handbook of High-Temperature Superconductivity: Theory and Experiment, J. R. Schrieffer, J. S. Brooks, eds. (Springer, Amsterdam, 2007).
- M. Vojta, Adv. Phys. 58, 699 (2009).
- A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
- J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 105, 027204 (2010).
- Y. Singh and P. Gegenwart, Phys. Rev. B 82, 064412 (2010).
- X. Liu, T. Berlijn, W.-G. Yin, W. Ku, A. Tsvelik, Y.-J. Kim, H. Gretarsson, Y. Singh, P. Gegenwart, and J.P. Hill, Phys. Rev. B 83, 220403 (2011).
- Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Phys. Rev. Lett. 108, 127203 (2012).
- T. Takayama, A. Kato, R. Dinnebier, J. Nuss, H. Kono, L. S. I. Veiga, G. Fabbris, D. Haskel, and H. Takagi, Phys. Rev. Lett. 114, 077202 (2015).
- K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, Phys. Rev. B 90, 041112 (2014).
- J. A. Sears, M. Songvilay, K. W. Plumb, J. P. Clancy, Y. Qiu, Y. Zhao, D. Parshall, and Y.-J. Kim, Phys. Rev. B 91, 144420 (2015).
- J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 110, 097204 (2013).
- J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Phys. Rev. Lett. 112, 077204 (2014).
- J. Reuther, R. Thomale, and S. Rachel, Phys. Rev. B 90, 100405(R) (2014).
- Y. Sizyuk, C. Price, P. Wölfle, and N. B. Perkins, Phys. Rev. B 90, 155126 (2014).
- I. Kimchi, R. Coldea, and A. Vishwanath, Phys. Rev. B 91, 245134 (2015).
- C. Chamon, Phys. Rev. B62, 2806 (2000).
- V. Lahtinen, A. W. W. Ludwig, and S. Trebst, Phys. Rev. B 89, 085121 (2014).
- For a comprehensive review on compass and Kitaev models see: Z. Nussinov and J. van den Brink, Rev. Mod. Phys. 87, 1 (2015) and preprint arXiv:1303.5922.
- E. H. Lieb, Phys. Rev. Lett. 73, 2158 (1994).
- A. J. Willans, J. T. Chalker, and R. Moessner, Phys. Rev. Lett. 104, 237203 (2010).
- S. Rachel, L. Fritz, and M. Vojta, Phys. Rev. Lett. 116, 167201 (2016).
- F. Zschocke and M. Vojta, Phys. Rev. B 92, 014403 (2015).
- H. Yao and S. A. Kivelson, Phys. Rev. Lett. 99, 247203 (2007).
- H. Yao, S.-C. Zhang, and S. A. Kivelson, Phys. Rev. Lett. 102, 217202 (2009).
- G. Baskaran, G. Santhosh, and R. Shankar, preprint arXiv:0908.1614
- M. Kamfor, S. Dusuel, J. Vidal, and K. P. Schmidt, J. Stat. Mech. P08010 (2010).
- T. Si and Y. Yu, Nucl. Phys. B 803, 428 (2008).
- S. Mandal and N. Surendran, Phys. Rev. B 79, 024426 (2009).
- M. Hermanns and S. Trebst, Phys. Rev. B 89, 235102 (2014).
- A different type of superstructure in the Kitaev model, which involves a non-trivial distribution of , , and bonds (but no modulation of their strength) has been discussed in E. Quinn, S. Bhattacharjee, and R. Moessner, Phys. Rev. B91, 134419 (2015).
- L. H. Wu and X. Hu, Sci. Rep. 6, 24347 (2016).
- G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev. Lett. 98, 247201 (2007).
- J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moessner, Phys. Rev. Lett. 112, 207203 (2014).
- Note that we ignore the fermion parity condition for the ground state(44); (31) because we are interested in the thermodynamic limit where its effects vanish.
- F. L. Pedrocchi, S. Chesi, and D. Loss, Phys. Rev. B84, 165414 (2011).
- X.-Y. Song, Y.-Z. You, and L. Balents Phys. Rev. Lett. 117, 037209 (2016).
- T. Senthil, S. Sachdev, and M. Vojta, Phys. Rev. Lett. 90, 216403 (2003); T. Senthil, M. Vojta, and S. Sachdev, Phys. Rev. B 69, 035111 (2004).
- J. Knolle, G.-W. Chern, D. L. Kovrizhin, R. Moessner, and N. B. Perkins, Phys. Rev. Lett. 113, 187201 (2014).
- T. Senthil and M. P. A. Fisher, Phys. Rev. Lett. 86, 292 (2001).