Vacancy-induced low-energy density of states in the Kitaev spin liquid
Wen-Han Kao, Johannes Knolle, Gábor B. Halász, Roderich Moessner, Natalia B. Perkins
VVacancy-induced low-energy density of states in the Kitaev spin liquid
Wen-Han Kao, Johannes Knolle,
2, 3, 4
G´abor B. Hal´asz, Roderich Moessner, and Natalia B. Perkins ∗ School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA Department of Physics TQM, Technische Universit¨at M¨unchen, James-Franck-Straße 1, D-85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), 80799 Munich, Germany Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA Max-Planck-Institut fr Physik komplexer Systeme, 01187 Dresden, Germany (Dated: July 24, 2020)The Kitaev honeycomb model has attracted significant attention due to its exactly solvable spin-liquid groundstate with fractionalized Majorana excitations and its possible materialization in magnetic Mott insulators withstrong spin-orbit couplings. Recently, the 5d-electron compound H LiIr O has shown to be a strong candidatefor Kitaev physics considering the absence of any signs of a long-range ordered magnetic state. In this work, wedemonstrate that a finite density of random vacancies in the Kitaev model gives rise to a striking pileup of low-energy Majorana eigenmodes and reproduces the apparent power-law upturn in the specific heat measurementsof H LiIr O . Physically, the vacancies can originate from various sources such as missing magnetic momentsor the presence of non-magnetic impurities (true vacancies), or from local weak couplings of magnetic momentsdue to strong but rare bond randomness (quasivacancies). We show numerically that the vacancy effect is readilydetectable even at low vacancy concentrations and that it is not very sensitive neither to nature of vacancies norto different flux backgrounds. We also study the response of the site-diluted Kitaev spin liquid to the three-spininteraction term, which breaks time-reversal symmetry and imitates an external magnetic field. We propose afield-induced flux-sector transition where the ground state becomes flux free for larger fields, resulting in a clearsuppression of the low temperature specific heat. Finally, we discuss the effect of dangling Majorana fermionsin the case of true vacancies and show that their coupling to an applied magnetic field via the Zeeman interactioncan also account for the scaling behavior in the high-field limit observed in H LiIr O . I. INTRODUCTION
Various types of disorder in quantum spin liquids (QSLs)have recently attracted a lot of attention from both experimen-tal and theoretical points of view [1–15]. There are three mainreasons for this interest. First, some level of disorder in var-ious forms of dislocations, vacancies, impurities, and bonddisorder is inevitable in real materials. Second, disorder cansignificantly affect the low-energy properties of these systems.In particular, if the system is close to a QSL state, quencheddisorder on top of the quantum disordered strongly correlatedspin state of a QSL can give rise to diverse and often puz-zling behaviors [8, 12, 13, 16]. Namely, sometimes disorderis detrimental for the QSL state since it either localizes theresonating spin singlets or induces competing glassy statesinstead of entangled ones [17–20]. However, in some othercases, e.g., in classical spin ice materials, certain forms ofdisorder can instead enhance the quantum dynamics of spinsthroughout the system and generate a QSL with long-rangeentanglement [6, 16, 21]. Third, given that the properties ofQSLs are difficult to detect directly because such states lackany local order parameter, much additional information canbe obtained by studying the distinctive responses to local per-turbations, such as static defects, dislocations, and magneticor non-magnetic impurities. In particular, these perturbationsmay nucleate exotic excitations characteristic to the spin liq-uid under consideration [1, 2].Of specific interest is the role of disorder in the materials ∗ [email protected] that have been suggested to be potential candidates [22–26]to realize the Kitaev QSL [27]. In a flurry of recent experi-ments on the honeycomb ruthenium chloride α -RuCl , it wasshown that both bond disorder and stacking disorder are notnegligible [28–32]. Perhaps, disorder also plays a crucial rolefor a potential proximity of Ag LiIr O to a Kitaev QSL state[33]. However, perhaps the most remarkable and intriguingconsequences of disorder have been observed in a presump-tive quantum spin liquid state of the hydrogen intercalatediridate H LiIr O [8]: (i) the specific heat displays a low-temperature divergence of C/T ∝ T − / ; (ii) only a smallfraction of the total magnetic entropy is released at these lowtemperatures; and (iii) there is a non-vanishing contributiondown to the lowest temperature in the NMR rate /T and analmost flat Knight shift. All of these observations signal thepresence of abundant low-energy density of states (DOS) re-lated to magnetic excitations. However, despite the presenceof dominant Kitaev exchange, this phenomenology is at oddswith the thermodynamics of the pure Kitaev QSL [34–36],which has a vanishing specific heat C/T ∝ T and a signifi-cant release of half of its total entropy at low T .Motivated by these experimental findings, some of us haverecently considered a minimal model of a bond disordered Ki-taev QSL [11] that can account for these salient experimentalobservations in H LiIr O [8]. However, in order to recoverthe low temperature scaling of the specific heat, the Kitaev-like model of Ref. 11 assumed a somewhat ad-hoc form ofbinary bond disorder and invoked a random flux backgroundeven at very low temperatures.In this work, we show that a finite density of vacancies inthe Kitaev model induces a pileup of the low-energy DOS, N ( E ) , and a low temperature divergence of the specific heat, a r X i v : . [ c ond - m a t . s t r- e l ] J u l C/T . These results are consistent with an algebraic diver-gence with an exponent around ν = 1 / over a broad rangeof low energies/temperatures. As a finite density of static,randomly-located vacancies is always present in the two di-mensional limit of layered materials, i.e., similar to the caseof graphene [37–41], our simple model provides a natural ex-planation for the experimental observations in H LiIr O [8].We also carefully treat the flux background and show thatthe energy of a random-vacancy configuration is minimizedwhen a flux is bound to each vacancy. Hence, we resolvethe problem of a random flux background since now the fluxconfiguration is determined by the vacancy distribution. Fi-nally, we show numerically that the vacancy effect for the low-energy DOS is robust with respect to the addition of bond-randomness or a random flux background. II. THE MODEL
The extended Kitaev honeycomb model is defined in termsof localized spin 1/2 degrees of freedom that are coupled in abond-anisotropic manner on the honeycomb lattice [27]. Thespin Hamiltonian reads H = − (cid:88) (cid:104) ij (cid:105) J (cid:104) ij (cid:105) α ˆ σ αi ˆ σ αj − κ (cid:88) (cid:104)(cid:104) ik (cid:105)(cid:105) ˆ σ αi ˆ σ βj ˆ σ γk , (1)where ˆ σ αi denotes Pauli spin operators with α = x, y, z and (cid:104) ij (cid:105) α labels the nearest-neighbor sites i and j along an α -type bond. The second term is the three-spin interaction withstrength κ ∼ h x h y h z J that imitates an external magnetic fieldand breaks time-reversal symmetry while preserving the exactsolvability [27]. By rewriting each spin operator in terms offour Majorana fermions, ˆ σ αi = i ˆ b αi ˆ c i , and defining the linkoperators ˆ u ij = i ˆ b αi ˆ b αj , the Hamiltonian takes the form H = i (cid:88) (cid:104) ij (cid:105) J (cid:104) ij (cid:105) α ˆ u (cid:104) ij (cid:105) α ˆ c i ˆ c j + iκ (cid:88) (cid:104)(cid:104) ik (cid:105)(cid:105) ˆ u (cid:104) ij (cid:105) α ˆ u (cid:104) kj (cid:105) β ˆ c i ˆ c k . (2)The solvability of the Kitaev model relies on the extensivenumber of conserved fluxes defined on each hexagonal pla-quette, ˆ W p = ˆ σ x ˆ σ y ˆ σ z ˆ σ x ˆ σ y ˆ σ z = (cid:81) (cid:104) ij (cid:105)∈ p ˆ u (cid:104) ij (cid:105) α , which canblock-diagonalize the Hamiltonian (1) into flux sectors sincefluxes commute with each other, [ ˆ W p , ˆ W p (cid:48) ] = 0 , and with theHamiltonian, [ ˆ W p , H ] = 0 . Both the flux operators ˆ W p andthe link operators ˆ u (cid:104) ij (cid:105) α have eigenvalues ± . Once the linkvariable is specified for each bond, the physically relevant fluxsector is determined, and the Hamiltonian can be solved ex-actly as a tight-binding model of Majorana fermions. For thepure Kitaev model ( κ = 0 ), it has been proven that the ground-state sector is flux free [42], i.e., the fluxes have eigenvalues W p = +1 for all plaquettes.Since the honeycomb lattice is bipartite, each unit cell con-tains two sublattice sites labeled as A and B . Thus, for a sys-tem with N unit cells, there are N lattice sites and N hexag-onal plaquettes. Under periodic boundary conditions (PBC),the fluxes can be excited only in pairs since flipping one linkvariable in the zero-flux sector results in W p = − on both sides of that link. As a result, there exists a global constraintfor the fluxes, (cid:81) p W p = 1 , such that the number of inde-pendent fluxes is reduced by one. Since two additional fluxdegrees of freedom are introduced by the toric topology, thetotal number of different flux sectors is then N +1 .With this decomposition, the spin degrees of freedom inthe original Hamiltonian are now fractionalized into itinerantMajorana fermions and static Z gauge fluxes. In a given fluxsector, the gauge can be fixed and all link variables can bespecified, leading to a Hamiltonian of non-interacting Majo-rana matter fermions, H = i (cid:0) c A c B (cid:1) (cid:18) F M − M T − D (cid:19) (cid:18) c A c B (cid:19) , (3)where the hopping amplitudes between sites on sublattice Aand sublattice B are M ij = J α ˆ u (cid:104) ij (cid:105) α , while the entries in thediagonal blocks F and − D contain hopping amplitudes be-tween sites on the same sublattice. Two adjacent Majoranafermions from the same unit cell can be combined into a com-plex matter fermion, (cid:40) ˆ f = (ˆ c A + i ˆ c B ) / f † = (ˆ c A − i ˆ c B ) / (4)such that the Hamiltonian can be written in a Bogoliubov − de-Gennes form and diagonalized in the standard way, H = 12 (cid:0) f f † (cid:1) (cid:18) ˜ h ∆∆ † − ˜ h T (cid:19) (cid:18) f † f (cid:19) = (cid:88) n (cid:15) n ( a † n a n −
12 ) , (5)where ˜ h = ( M + M T ) + i ( F − D ) and ∆ = ( M T − M ) + i ( F + D ) . Therefore, for a given flux configuration, thefermionic ground-state energy reads E ( { u ij } ) = − (cid:80) n (cid:15) n ,and the global density of states is given by N ( E ) = (cid:88) n δ ( E − (cid:15) n ) . (6) III. KITAEV MODEL WITH RANDOM VACANCIESA. Vacancies, quasivacancies, and fluxes
A vacancy is usually a simple absence of an atom at a givensite. However, in the present work, we will use this term moregenerally since it can also correspond to a non-magnetic im-purity or a magnetic moment that is weakly connected to itsneighbors due to extremely strong but relatively rare bond ran-domness. To distinguish it from the simple absence of an atom(which we call a true vacancy), we will refer to the latter typeof defect as a quasivacancy .In order to introduce randomly distributed vacancies intothe Kitaev honeycomb model (2), we first consider the time-reversal symmetric case at κ = 0 . In this case, the second J’ u = -1(a) Quasivacancy (b) Quasilocalized mode (c) Bound-flux sector
FIG. 1. (Color online) (a) A pair of quasivacancies introduced on different sublattices. The blue dashed lines represent reduced couplings J (cid:48) → , while the red circles depict the real-space wave functions of the resulting quasivacancy modes. (b) The real-space wave functions ofthe zero-energy quasilocalized modes introduced by the vacancies. (c) Bound-flux sector. By flipping a string of link variables from u = +1 to u = − , a pair of fluxes can be attached to each pair of vacancies in order to minimize the energy of the system. Note that, for J (cid:48) → , thethree flux degrees of freedom around the vacancy site effectively merge into one. term in Eq. (2) is absent, while the first term can be written as H = i (cid:88) (cid:104) ij (cid:105) i,j ∈ P J (cid:104) ij (cid:105) α ˆ u (cid:104) ij (cid:105) α ˆ c i ˆ c j + i (cid:88) (cid:104) kl (cid:105) k ∈ V ,l ∈ P J (cid:48)(cid:104) kl (cid:105) α ˆ u (cid:104) kl (cid:105) α ˆ c k ˆ c l , (7)where P denotes the subset of normal lattice sites and V de-notes the subset of vacancy sites. We consider a compensatedcase with equal numbers of vacancies on the two sublattices ofthe honeycomb lattice. By taking the limit of J (cid:48) α (cid:28) J α , sitesbelonging to V behave as quasivacancies. For such a quasiva-cancy, a Majorana fermion ˆ c remains on the vacancy site, butits nearest-neighbor hopping amplitudes are removed in thelimit of J (cid:48) α → . Therefore, the way we diagonalize the pureKitaev Hamiltonian remains valid, even though the number offlux degrees of freedom is effectively reduced.For each vacancy with J (cid:48) α → , there are three hexagonalplaquettes around the vacancy site, resulting in = 8 distinctflux sectors labeled by W , , = ± . However, the MajoranaHamiltonian in each flux sector is completely determined by W v = W W W = ± , which corresponds to a large va-cancy plaquette merged from the three individual plaquettesaround the vacancy. Since the remaining two degrees of free-dom do not affect the Majorana DOS and can even be partially“gauged away” in the case of true vacancies [43], we ignorethem in the rest of this work and characterize the vacancy witha single vacancy flux ˆ W v = ˆ W ˆ W ˆ W . Consequently, fora system with N v isolated vacancies, the number of flux de-grees of freedom is effectively reduced by N v . With periodicboundaries, the global constraint for the fluxes then becomes (cid:81) p W h,p (cid:81) q W v,q = 1 , where W h correspond to bulk hexag-onal plaquettes. Since there are two physically distinct fluxoperators ˆ W h and ˆ W v , introducing a flux pair falls into oneof three situations: both fluxes on hexagonal plaquettes, bothfluxes on vacancy plaquettes, and one flux on each kind ofplaquette. We will show that, in order to minimize the totalenergy, fluxes must be bound to the vacancy plaquettes. B. Quasilocalized eigenmodes and flux binding
When considering vacancies in the Kitaev model, many re-sults are completely analogous to those in graphene [37–41].For example, it was found for both systems that introducing avacancy leads to a zero-energy eigenmode with a quasilocal-ized wavefunction on the other sublattice around the vacancysite [1, 2, 37, 38].Indeed, if we consider only the nearest-neighbor couplingsand ignore the possible flux excitations in the Kitaev model,then there is a one-to-one mapping between the Majoranahopping in the Kitaev model and the fermionic hopping ingraphene (see Fig. 1 (b)).It was shown by Willans et al. through analytical calcula-tions that a single vacancy binds a flux in the gapped phase ofthe Kitaev model ( J x , J y (cid:28) J z ) [1, 2]. In the gapless phase(including the isotropic point J x = J y = J z ), the flux-bindingeffect can be verified numerically. Practically, we considertwo vacancies, one on an A sublattice site and the other on aB sublattice site (see Fig. 1 (a)), which are separated by a dis-tance ∼ L/ , where L is the linear dimension of the system,and then calculate the energy difference between the bound-flux (Fig. 1 (c)) and the zero-flux sectors (Fig. 1 (a)): E bind = E bound − E zero . (8)In Fig. 2 (a), different system sizes up to L = 120 are con-sidered and, by extrapolation, we show that the flux-bindingenergy converges to E bind = − . J , which is consistentwith the previous result [1]. In the same way, we also calculatethe binding energies for non-zero J (cid:48) and show (see Fig. 2 (b))that the flux-binding effect remains for J (cid:48) /J < . , whichcan be compared to another previous result with a slightlydifferent setup [4]. The flux-binding effect for non-zero J (cid:48) indicates that we can extend the quasivacancy picture to abond-disordered model where the vacancy sites are not truly (a)(b) FIG. 2. (Color online) Binding energy of a flux to a vacancy. (a)The flux-binding energy converges to E bind = − . J in thethermodynamic limit by extrapolation. System sizes are classified by L mod 3 since they have different rates of convergence. The bindingenergies are calculated for finite-size systems up to L = 120 with J (cid:48) = 0 . (b) The J (cid:48) dependence of the flux-binding energy. Byinterpolation, the critical value of J (cid:48) is estimated as . J , whichsets the upper limit of the flux-binding effect. removed or replaced by non-magnetic ions but the couplingstrengths are strongly suppressed by structural disorder.Note that, in this setting, the boundary conditions play acrucial role in the flux binding effect. For open boundary con-ditions, it is possible to create a single flux on each vacancyplaquette so that the energy is minimized. However, for peri-odic boundary conditions, the fluxes must be created in pairs.For a system with only one vacancy, there exists only one va-cancy plaquette to bind the flux and thus the other flux mustbe bound to a hexagonal plaquette. This arrangement resultsin a higher total energy since the flux excitation energy on ahexagonal plaquette is . J [27] which can not be fullycompensated by E bind . Therefore, the ground-state sector isstill flux free. In contrast, when the system contains an evennumber of isolated vacancies, it is possible to bind fluxes toall vacancy plaquettes, as we describe in the following. C. Bound-flux sector
In order to minimize the total energy of a random-vacancyconfiguration, fluxes must be introduced and bound to eachvacancy plaquette. The most unbiased approach is to apply aMarkov chain Monte Carlo simulation that samples the fluxconfigurations at low temperatures [34, 35]. However, thisapproach is not easily realized for large systems with disor-der since the tight-binding Hamiltonian ( L × L matrix)must be diagonalized for each update. Thus, here we discusshow to generate the low-energy bound-flux sector in whicheach vacancy has a flux attached to it. Note that when two ormore vacancies are connected to each other and thus the dis-cussion in Sec. III A is no longer valid, our approach may notprovide the actual ground-state flux sector. However, in thedilute limit, most vacancies are isolated and the ground stateis well approximated with the bound-flux sector.In all our calculations, we use periodic boundary condi-tions because open boundaries lead to additional zero-energyeigenmodes and make it difficult to examine the direct con-sequences of adding vacancies into a finite-size system. Asdiscussed in the previous subsection, the consequence of PBCis that fluxes always appear in pairs. When a link variable u ij is flipped, two fluxes are introduced on adjoining hexagons.Therefore, when we flip a link variable on the edge of a va-cancy plaquette, one flux emerges on this vacancy plaquetteand the other one emerges on a neighboring hexagonal pla-quette. The former changes the energy by ∆ E v = E bind < ,while the latter changes the energy by ∆ E f > .Numerical studies of the single-vacancy effect show that theenergy decrease of a flux binding to a vacancy cannot com-pensate the energy increase by the other flux [2]. Indeed, atthe isotropic point, ∆ E v converges to − . J , while ∆ E f converges to . J . In contrast, for a system with two iso-lated vacancies, it is possible to generate a flux pair on hexag-onal plaquettes, and then propagate the two single fluxes in-dividually until they bind to the vacancies, lowering the totalenergy by approximately | ∆ E v | compared to the zero-fluxcase (see Fig. 3 (b)).However, finding the ground-state flux configuration formore than two vacancies in terms of link variables is a non-trivial task because fluxes may create or annihilate during theflipping of link variables. Instead, we propose an algorithmthat generates a vacancy configuration along with the boundfluxes. First, we randomly choose a position on the lattice andplace a vacancy pair. The pair contains one vacancy on anA-sublattice site and the other one on a nearby B-sublatticesite, as shown in Fig. 3(a). One shared link of the two va-cancy plaquettes is flipped such that each plaquette binds aflux. Second, we randomly move one vacancy in the direc-tions of the two primitive vectors. When a vacancy migrates,the link variables are flipped along the same path, keeping theflux bound to the vacancy plaquette (Fig. 3 (b)). Thus, insteadof propagating fluxes to minimize the energy of a given va-cancy configuration, as is shown in Fig. 1 (c), we propagatecomposites of one vacancy and one flux to generate a randomconfiguration of vacancies with bound fluxes.Note that this method does not give the ground-state flux (a)(b) FIG. 3. (Color online) The recipe for creating vacancies with boundfluxes. (a) Randomly place one pair of vacancies on the lattice. Flipone shared link variable (thick black line) such that two fluxes arecreated and attached to the two vacancy plaquettes. (b) Randomlymove one of the vacancies in the pair, and flip a string of link vari-ables along the path, such that the fluxes are always bound to thevacancies. After the migration, place another pair of vacancies withfluxes and repeat the same process. sector for all vacancy configurations. When two or more va-cancies are connected after migration, the binding fluxes mayannihilate each other such that no flux binds to the merged va-cancy plaquette. Nevertheless, this method works relativelywell for a low density of vacancies since connected vacanciesare then rare. In the following sections, we denote the finite-flux sector generated by this method as the bound-flux sector . D. Time-reversal symmetry broken case ( κ (cid:54) = 0 ) The effective three-spin interaction term in Eq. (1) breakstime-reversal symmetry and changes the energetics of themodel by simultaneously gapping out the fermionic spectrumand introducing localized zero-energy Majorana modes in thepresence of isolated fluxes [27]. Here we test the disorder-averaged total energy of the system in the bound-flux andzero-flux sectors for various κ . By comparing the energiesof the two flux sectors, a ground-state transition from thebound-flux to the zero-flux sector is observed. This transi-tion is shown in Fig. 4, where we plot the difference betweenthe two energies, each obtained by averaging over 4000 disor- FIG. 4. (Color online) Ground-state transition between the bound-flux and the zero-flux sectors. Below (above) κ = 0 . , the bound-flux sector has lower (higher) energy than the zero-flux sector. Theresults are averaged over 4000 disordered samples ( L = 40 ) with2% quasivacancies ( J (cid:48) = 0 . ).FIG. 5. (Color online) Density of states in the bound-flux sector with2% true vacancies ( J (cid:48) = 0 ). The red dashed line shows the analyt-ical density of states of the pure Kitaev honeycomb model withoutfluxes and vacancies. The spectral weight is transferred towards thelow-energy region from the high-energy region near the Van Hovesingularity. The inset shows the data in base-10 logarithmic scalesand demonstrates the N ( E ) ∼ E − / behavior at low energies. dered samples with 2% quasivacancies at J (cid:48) = 0 . . Whilethe bound-flux sector is lower in energy for ≤ κ ≤ . , thezero-flux sector becomes energetically favorable for κ ≥ . . IV. DENSITY OF STATES AND SPECIFIC HEAT
In this section, we discuss how the presence of vacanciesresults in a low-temperature divergence in the specific heat,
C/T , which might be related to the recent experimental obser-vation by Kitagawa et al. on the Kitaev spin-liquid candidateH LiIr O [8]. We consider both the case of true vacancies with J (cid:48) = 0 and the case of quasivacancies with J (cid:48) > .At finite temperatures, both itinerant Majorana fermionsand fluxes contribute to the specific heat and ln 2 in the ther-mal entropy [34, 35]. However, in H LiIr O , it is reportedthat the low-energy excitations release only - of ln 2 en-tropy at 5K, implying that the flux degrees of freedom mightbe frozen. Therefore, we assume that the temperature depen-dence of the specific heat is solely due to the thermal occupa-tion of itinerant Majorana fermions in static flux sectors, C ( T ) = (cid:90) E · N ( E ) · ∂n F ( E, T ) ∂T d E = (cid:88) n (cid:16) (cid:15) n T (cid:17) e (cid:15) n /T ( e (cid:15) n /T + 1) , (9)where n F ( E, T ) = ( e E/T + 1) − is the Fermi function, andwe use the definition of the density of states in Eq. (6) to reachthe final expression.As before, we first consider the time-reversal symmetriccase with κ = 0 in Secs. IV A and IV B, and then discussthe effect of a finite κ in Sec. IV C. A. True vacancies
We introduce a certain amount of vacancies into the Majo-rana problem, as described by Eq. (7), at the isotropic point( J x = J y = J z ≡ J = 1 ) and see how the fermionic specificheat and density of states are affected. The vacancy concen-tration n v is defined as the number of vacancies ( N v ) dividedby the number of lattice sites ( L ). Half of the vacanciesare on the A sublattice and the other half are on the B sublat-tice. For example, in a system with L = 20 and n v = 2% ,we randomly put 8 vacancies on A sites and 8 vacancies onB sites. When generating the random vacancy configurations,we avoid removing the same site twice or more. Thus, fora given concentration, each disorder realization contains thesame amount of vacancies. In the remainder of the paper, theresults are presented for systems with linear dimension be-tween L = 20 and L = 40 on a torus, and all the data areaveraged over to disorder realizations.First, we demonstrate the pileup of low-energy states in asystem with 2 % true vacancies. The density of states in thebound-flux sector, averaged over 4000 realizations ( L = 40 ),is shown in Fig. 5. Since the low concentration of vacanciesacts like a weak disorder on top of the Kitaev spin liquid, theoverall behavior of the density of states is similar to the ana-lytical result for the pure Kitaev model in the zero-flux sector,except for the low-energy region. Note that the states at ex-actly zero energy are removed from the density of states, sothat the pileup in the low-energy region is exclusively fromstates with small but non-zero energies. The upturn is clearlyseen in the log-log plot (inset of Fig. 5), and it fits well to apower-law form N ( E ) ∼ E − ν with ν ≈ / .Next, the fermionic specific heat (9) is calculated fromthe density of states and shown in Fig. 6. As expected, thepileup of low-energy states approximated by the power law (a)(b) FIG. 6. (Color online) (a) Specific heat
C/T for L = 40 systems cal-culated from the fermionic density of states. Compared to the pureKitaev model without fluxes and vacancies, adding 2% true vacan-cies ( J (cid:48) = 0) leads to a clear upturn in both the bound-flux and thezero-flux sectors. The inset shows the data in logarithmic scales. (b)The system-size dependence of C/T . The low-temperature upturnbecomes more robust for larger systems due to the reduction of thefinite-size effect. The curves are averaged over 4000-10000 disorderrealizations, depending on the system size. N ( E ) ∼ E − ν brings about a similar behavior in the specificheat, C/T ∝ T − ν with an exponent close to 0.5. This power-law behavior is consistent with what has been found in theKitaev spin-liquid candidate H LiIr O [8], where the authorsmentioned that an approximately 2 % density of magnetic im-purities or vacancies in the material may be responsible forthe magnetization and specific heat results. Thus, even with-out the presence of bond randomness [11], the low-energyfermionic states produced by the vacancies in the ground-stateflux sector can give rise to a similar power-law behavior.In Fig. 6 (b), the size dependence of this power-law upturnis presented. Because of the gapless nature of the Kitaev spinliquid at the isotropic point, the finite-size effects are consid-erable at low temperatures. We found that L = 40 is a rea-sonable choice in practice such that the finite energy of thevacancy states can be extended to the scale of − . In therest of this paper except Sec. V, systems with L = 40 will be (a) Bound-flux sector (b) Zero-flux sector (c) Random-flux sector FIG. 7. (Color online) The low-energy density of states and fermionic specific heat with different concentrations of true vacancies ( J (cid:48) = 0) inthe (a) bound-flux, (b) zero-flux, and (c) random-flux sectors. All the three cases show a clear upturn in both the density of states and C/T . used for all calculations.In addition to the bound-flux sector, we also considered thezero-flux and random-flux sectors for the same set of random-vacancy configurations. The corresponding densities of statesand fermionic specific heats are presented in Fig. 7 for vari-ous vacancy concentrations. The pileup of vacancy-inducedstates (Fig. 7, upper row) and the corresponding upturn in thespecific heat (Fig. 6 (a) and Fig. 7, lower row) appears in allthe flux sectors, indicating that the effect of vacancies plays amajor role in the low-energy region. Besides, we also see thatthe upturn power is slightly dependent on the vacancy con-centration. In the bound-flux and zero-flux sectors, densitiesof states with different n v start splitting below a characteristicenergy scale, which is similar to the tight-binding model ofgraphene with compensated vacancies [40]. B. Quasivacancies
In order to study the effect of quasivacancies introducedin the Kitaev model, we computed the density of states forthe model (7) with different coupling strengths J (cid:48) α = J (cid:48) upto 0.05 in the bound-flux (Fig. 8 (a)), zero-flux (Fig. 8 (b))and random-flux (Fig. 8 (c)) sectors. Note that, based on theenergetic analysis shown in Fig. 2 (b), the bound-flux sector islower in energy than the zero-flux and the random-flux sectorsfor all vacancy couplings J (cid:48) ≤ . .In the limit of J (cid:48) → , the resonant peak around zeroenergy is present in all three cases, indicating that the low-energy physics is governed by vacancies rather than fluxes. Inthe bound-flux sector, a finite value of J (cid:48) leads to a couplingof the quasivacancy mode to its surroundings and results in alarger width of the zero-energy peak. In the zero-flux sector, however, a tiny energy gap opens when increasing the mag-nitude of J (cid:48) . This phenomenon comes from the hybridiza-tion of the two different zero-energy modes corresponding tothe same quasivacancy. When switching on the coupling J (cid:48) ,the quasivacancy mode (Fig. 1(a)) begins to hybridize withthe zero-energy quasilocalized mode (Fig. 1(b)), leading to asplitting of the corresponding energy levels. The formation ofthis pseudo-gap is also reported in site-diluted graphene [38].The hybridization picture of low-energy modes will becomeeven more clear when we open a larger bulk gap by addingthe time-reversal symmetry breaking term to the system. C. Three-spin interaction
In the previous sections we have demonstrated how thevacancy-induced low-energy states give rise to upturns in boththe fermionic DOS and the specific heat
C/T , similar to theexperimental findings in H LiIr O [8]. Now we turn our at-tention to the remaining question from the experiment: howdoes this low-energy upturn get suppressed in the presence ofan external magnetic field? As discussed above, the three-spin interaction with strength κ ∼ h x h y h z J in Eq. (1) rep-resents the leading-order perturbation effect of the Zeemanterm [27]. This interaction breaks time-reversal symmetryand introduces zero-energy Majorana modes in the presenceof fluxes. Therefore, the presence of vacancies in conjunc-tion with the flux-binding effect provides a natural scenariofor creating Majorana zero modes at low temperatures.By introducing the three-spin term into the pure Kitaev hon-eycomb model, the gapless spin liquid becomes gapped dueto the next-nearest-neighbor hopping of Majorana fermions.This effect is clearly seen in Fig. 9, where the fermonic den- (a) Bound flux(b) Zero flux(c) Random flux FIG. 8. (Color online) Density of states for various quasivacancycoupling strengths J (cid:48) with a concentration n v = 2% of quasivacan-cies. The low-energy states induced by the quasivacancies are accu-mulated around the Dirac point of the Majorana spectrum. The zero-energy peak is suppressed for larger J (cid:48) due to the stronger couplingto the bulk. For the zero-flux sector, a pseudo-gap gradually formsby the hybridization of the quasivacancy mode and the quasilocalizedmode around the same quasivacancy. sity of states is shown for different three-spin couplings κ inboth the bound-flux and the zero-flux sectors for a 2% con-centration of vacancies (recall that the bound-flux sector isthe ground state flux sector only for ≤ κ ≤ . ). Fig. 9 alsoclearly shows that there is a pronounced difference betweenthe two flux sectors in the number of resonant peaks insidethe bulk energy gap. In the zero-flux sector, two broad peaksappear inside the energy gap and move away from E = 0 with increasing κ . In contrast, one additional resonant peakis present around E = 0 in the bound-flux sector whose po-sition is independent of κ . Note that, due to its finite width,this central peak contains a number of eigenmodes with smallbut nonzero energy, resulting in measurable signatures in ther-modynamic quantities such as the specific heat at low temper-atures. The presence or absence of this peak along with theflux-sector transition shown in Fig. 4 plays a crucial role inunderstanding the C/T results. A cartoon picture of the eigen- mode hybridization leading to different numbers of peaks inthe two flux sectors will be discussed in Sec. V.The localized nature of these vacancy-induced eigenmodescan be illustrated by the inverse participation ratio (IPR). Thisquantity is defined as P n = (cid:88) i | φ n,i | , (10)where the index n labels the eigenmode wave function φ n,i and the index i labels the lattice site. In Fig. 9 the IPR for eacheigenmode is shown by red dots. For a delocalized mode, theIPR scales roughly as ∼ /N in a system with N sites sincethe wavefunction is spread out uniformly over the entire lat-tice. This behavior is precisely what we see for the fermionicbulk modes. However, for the in-gap modes introduced by thevacancies, the IPR is significantly larger since the wave func-tion is confined to a small portion of the lattice. Similarly tographene, when κ = 0 , each vacancy leads to a zero-energyeigenmode with a quasilocalized wave function on the othersublattice around the vacancy site. For a single vacancy, thiswave function can be written in an analytical form [37, 38]: Ψ( x, y ) ∼ e i K (cid:48) · r x + iy + e i K · r x − iy . (11)The wavevectors K and K (cid:48) denote the two different Diracpoints on the corner of the first Brillouin zone. While the an-alytical form of the quasilocalized wave function is no longeravailable for a finite density of vacancies, we can still relatethe enhanced values of the IPR for the in-gap states in ournumerical calculations to the /r decay of the quasilocalizedwavefunctions Ψ( x, y ) . The IPR is particularly large for the E = 0 mode in the bound-flux sector.In Fig. 10, the temperature dependence of C/T is presentedfor various values of κ . According to the flux-sector transitionshown in Fig. 4, the ground-state flux sector is the bound-fluxsector for κ < . and the zero-flux sector for κ ≥ . .As discussed in the previous sections, the upturn of the curvefor κ = 0 can be extended to very low temperatures for largesystem sizes, and it is comparable to the experimental resultfor H LiIr O without magnetic field. For κ = 0 . , a smallnumber of localized modes appear in the DOS (see Fig. 9 (a)),giving rise to a steeper upturn in C/T . Still, there is no sup-pression at the lowest temperatures due to the presence of thecentral resonant peak in the density of states. However, when κ exceeds the critical value and the ground-state flux sectorbecomes flux free, C/T only shows a small upturn and is thenstrongly suppressed. The small upturn comes from the in-gapresonant peak at
E > and the suppression is due to the lackof lower energy modes around E = 0 . Note that, if a vacancysite is completely decoupled from the system, the correspond-ing vacancy mode has exactly zero energy and has no contri-bution to the specific heat. However, since we consider thequasivacancy scenario that possibly results from a bond dis-order of Kitaev interactions, those quasivacancy modes havefinite couplings via J (cid:48) and κ and can be hybridized with otherlocalized modes to produce finite-energy eigenmodes, thusleading to the low-temperature upturn in C/T . I P R I P R I P R I P R I P R I P R I P R I P R (a) Bound-flux sector (b) Zero-flux sector FIG. 9. (Color online) Density of states and inverse participation ratio (IPR) for different three-spin couplings κ with fixed J (cid:48) = 0 . and n v = 2% . The IPR results are intensified by a factor of 2 in order to be comparable to the density of in-gap states. (a) Bound-flux sector: threebroad peaks of in-gap states appear when a bulk gap opens for κ > . (b) Zero-flux sector: only two broad peaks appear inside the gap. Thecrucial difference in the number of peaks comes from the hybridization of localized states in the low-energy subspace and, in particular, fromthe presence of Majorana zero modes in the bound-flux sector. The quasivacancy picture and the corresponding
C/T re-sults capture the experimental findings for
C/T in the Kitaevspin-liquid candidate H LiIr O . In Figure 4 of Ref. 8, a pe-culiar scaling law is used for collapsing the C/T data in awide range of magnetic fields h : C/T ∼ h − / T. (12)Since we consider only the leading-order three-spin interac-tion κ emerging from a perturbative treatment of the magneticfield (rather the magnetic field h itself), a simple replacementof h with κ does not provide the correct scaling law for ourdata. Nevertheless, it is possible to assume a general power-law relation h ∼ κ α and test the following scaling behavior: C/T ∼ ( κ α ) − / T. (13)We estimate the optimal α to be around . by collapsing thecurves with various κ below the temperature scale T /κ α ∼ . , as shown in Fig. 10 (b). Qualitatively, our calculation ofthe low-temperature fermionic specific heat is able to capturethe specific heat scaling obtained experimentally in Ref. 8. D. Effect of dangling Majorana fermions
In the previous subsection, we considered the leading-ordereffect of a magnetic field on quasivacancies, where the linkvariables ˆ u ij = i ˆ b αi ˆ b αj emanating from the quasivacancies arewell defined. However, for true vacancies, there are no Ma-jorana fermions on the vacancy sites, and the link variables ˆ u are thus no longer well defined around the vacancies. Conse-quently, the dangling Majorana fermions ˆ b α on the neighbor-ing sites lead to additional terms in the Hamiltonian. Indeed,if we start from the bare Zeeman term on a neighboring site j ,the Majorana-fermion representation immediately gives δ H j = − h α ˆ σ αj = − ih α ˆ b αj ˆ c αj , (14)where ˆ b α is a dangling Majorana fermion if the site j is con-nected to the vacancy site by an α -type bond. The Zeemanterm can then be readily merged into the original tight-bindingHamiltonian as it is equivalent to a hopping term betweenthe dangling Majorana fermion ˆ b α and the matter Majoranafermion ˆ c on the same site. Note that this term is a direct con-sequence of the magnetic field and is not derived from per-turbation theory. As such, it provides the primary effect of amagnetic field in a system with true vacancies.Without loss of generality, this effect is demonstrated in0 (a)(b) FIG. 10. (Color online) Specific heat
C/T calculated from the Ma-jorana fermion spectrum in Fig. 9. (a) Temperature dependence of
C/T for various three-spin couplings κ . The green (blue) dashed(solid) lines are calculated for the bound-flux (zero-flux) sector,which is the ground-state flux sector for small (large) κ . (b) Scalingplot of the curves with large values of κ for which the ground-stateflux sector is the zero-flux sector. Fig. 11 for a magnetic field h z applied in the z direction. En-ergetic considerations show that the zero-flux sector becomesthe ground-state flux sector for h z ≥ . . Again, the clearsuppression in C/T can then be attributed to the formationof an energy gap. Indeed, the dangling Majorana fermions ˆ b z are gapped out through h z , which is similar to the discussionon Fig. 8 (b). Since the field strength h z represents the ac-tual magnetic field instead of the three-spin coupling κ , weare able to see the scaling law C/T ∼ ( h z ) − / T from thedata collapse in the inset of Fig. 11 (b). Interestingly, eventhough we consider the effect of the magnetic field only onthe dangling Majorana fermions (but not on the bulk system),the peculiar scaling behavior can be reproduced in the high-field limit. More detailed results on the contributions of truevacancies to the thermodynamics and dynamical responses ofthe Kitaev model will be reported elsewhere. (b)(a) FIG. 11. (Color online) Effect of dangling Majorana fermions for a2% concentration of true vacancies. The specific heat
C/T is shownfor the bound-flux sector (a) and the zero-flux sector (b). The mag-netic field is applied along the z axis so that it couples the danglingMajorana fermion ˆ b z to the rest of the system. The inset of (b) showsthe field scaling and the collapse of the curves for higher fields. E. Effect of bond disorder
In order to address the peculiar low-energy behavior ofH LiIr O , it was proposed that bond disorder [10, 11] mayplay a major role in generating the upturns found in both C/T and the density of states. The physical origin of bond random-ness can be traced back to the random positions of the protonsbetween the honeycomb layers, leading to a local distortion ofthe oxygen octahedral cage, a change in the Ir-O-Ir bondingangle, and hence the local variation of Kitaev couplings alongthe Ir-Ir links [10, 44, 45]. Here, we consider the combinedeffect of bond randomness and a concentration of vacan-cies. A random coupling term is added to the nearest-neighborcoupling on the normal lattice sites (the first term in Eq. (7)): J (cid:104) ij (cid:105) α −→ J (cid:104) ij (cid:105) α + δJ (cid:104) ij (cid:105) α . (15)In the Gaussian bond disorder model, this replacement is ap-plied to all the bonds except the bonds emanating from the va-1 (b)(a) FIG. 12. (Color online) Effect of bond disorder for a 2% concentra-tion of vacancies. For Gaussian disorder, a normal random distribu-tion of coupling strengths with standard deviation 0.25 is applied toall bonds. For binary disorder, 25% of the bonds obtain additionalcoupling strength δJ = ± . so that the total strength becomes ei-ther . or . . All results are calculated in the bound-flux sector. cancy sites, and δJ is assigned randomly from a Gaussian dis-tribution with mean value 0 and standard deviation σ J = 0 . .An additional constraint of J + δJ ≥ is implemented toprevent couplings with mixed signs which would correspondto random flux insertion [3]. On the other hand, in the bi-nary disorder case, the random coupling term is fixed to be δJ = ± . , and only of the bonds are randomly cho-sen to be disordered. This implementation of binary disorderis consistent with the previous work [11], apart from usingthe random-flux sector. For both true vacancies and quasiva-cancies in our system, the bound-flux sector has lower energythan the zero-flux and random-flux sectors. Thus, the resultspresented in Fig. 12 are all calculated for the bound-flux sec-tor.For true vacancies, the additional bond randomness onlyleads to a small change in the power-law exponent of the spe-cific heat C/T , implying that the primary cause of the low-temperature upturn is the presence of vacancies. However,a non-zero quasivacancy coupling in conjunction with bondrandomness can result in a more noticeable change up to apower law
C/T ∼ T − , indicating that the distinction be-tween true vacancies and quasivacancies may be of great im-portance when studying the Kitaev spin-liquid materials. V. HYBRIDIZATION OF LOW-ENERGY MODES
For a large enough three-spin interaction κ , we can see inFig. 9 that the number of broad peaks inside the energy gap isnot the same in the bound-flux sector as in the zero-flux sec-tor. The spectrum of those in-gap states can be understood byconsidering the hybridization among the low-energy localizedmodes introduced by the quasivacancies.Let us first discuss the zero-flux sector. For κ = 0 , eachquasivacancy induces two quasilocalized modes around its va-cancy site. The first one is the fully localized vacancy mode( v-mode ) which can couple to its neighbors through the weakcoupling J (cid:48) (Fig. 1(a)). The second one is the quasilocalizedmode whose wave function is restricted to the opposite sub-lattice and decays as /r with the distance r from the vacancysite (Fig. 1(b)). For κ (cid:54) = 0 , this quasilocalized mode extendsto both sublattices and becomes properly localized with mostof its wave function distributed over the periphery of the va-cancy plaquette ( p-mode ). For a large enough gap, the in-gapspectrum can then be approximated with the hybridization ofthese localized modes. For example, we can consider a simplemodel with only two vacancies and four localized modes: δ H zero ≈ iµ ( c v, ˜ c p, + c v, ˜ c p, ) + iη (˜ c p, ˜ c p, + c v, c v, + c v, ˜ c p, + c v, ˜ c p, ) + h.c., (16)where c v,i is the v-mode Majorana and ˜ c p,i is the p-mode Ma-jorana corresponding to the i th vacancy ( i = 1 , ). We addthe tilde on the p-mode Majorana to highlight that its wavefunction is not confined to a single site. Two energy scales µ and η are defined to represent the couplings between lo-calized modes around the same vacancy and around differentvacancies, respectively. Importantly, µ increases linearly withboth J (cid:48) and κ , while η decays exponentially with the distancebetween the two vacancies. Thus, for well-separated vacan-cies, the first energy scale is much larger than the second one: µ (cid:29) η . The eigenspectrum of δ H zero then consists of twodoublets with energies ± µ and a small splitting ∼ η withineach doublet (see Fig. 13 (b)). We verified this simple pictureby obtaining the exact energy levels of a single L = 20 sys-tem with only two random vacancies and confirming that theresulting in-gap spectrum (see Fig. 13 (a)) is consistent withthe one obtained from δ H zero (see Fig. 13 (b)). For a finiteconcentration of vacancies, there is further hybridization onthe scale of η which broadens the two peaks at ± µ but leavesthe overall two-peak structure intact (see Fig. 9 (b)).On the other hand, for the bound-flux sector with finite κ ,one additional low-energy Majorana mode is introduced bythe flux on each vacancy plaquette ( f-mode ). When the fluxesare far apart from each other, these modes become Majoranazero modes and can be interpreted as anyons with non-Abelianstatistics [27, 46]. However, if the fluxes are closer and inter-act with each other, these modes hybridise and broaden into amini-band [47, 48]. Thus, in the simple two-vacancy model ofthe low-energy subspace, we must consider the hybridization2 (e)(d) (f) ↓ ↓ μ ημ p v f p v f (b)(a) ↓ μ μη (c) p v p v FIG. 13. (Color online) Hybridization between the low-energy modes around two vacancies in the zero-flux sector (a)-(c) and the bound-fluxsector (d)-(f). The densities of states and the real-space wave functions of the low-energy modes are calculated for L = 20 systems with J (cid:48) = 0 . and κ = 0 . . Note that in the density of states (a) and (d), each peak inside the gap (pointed by a red arrow) contains two stateswith a very small energy-level split. This split corresponds to the energy scale η in the spectra of the simple model Hamiltonians δ H zero (16)and δ H bound (17), which are shown in (b) and (e). In (c) and (f), we present cartoon pictures of the hybridization. The labels p , v , and f referto the peripheral mode ( p -mode), the vacancy mode ( v -mode), and the flux mode ( f -mode). of six localized modes: δ H bound ≈ iµ (cid:88) a (cid:54) = b ( c a, c b, + c a, c b, ) + iη (cid:88) a ( c a, c a, )+ iη (cid:88) a (cid:54) = b ( c a, c b, + c a, c b, ) + h.c., (17)where the summations in a and b are over the three dis-tinct types of modes ( v , p , f ). As in the zero-flux case, thelarger energy scale µ represents the couplings between modesaround the same vacancy, while the smaller energy scale η represents the couplings between modes around different va-cancies. The eigenspectrum of δ H bound has three doublets atenergies and ± µ , which explains the presence of the addi-tional zero-energy peak in the density of states (see Fig. 9 (a)).In general, the number of resonant peaks inside the bulk gapcorresponds to the number of localized modes around eachvacancy.The distinction between the bound-flux and the zero-fluxsectors in the presence or absence of the central peak leadsto very different low-temperature behaviors in the fermionicspecific heat. It also clearly distinguishes the vacancies in theKitaev model from those in graphene. VI. CONCLUSION
In this work, we have demonstrated that introducing a smallconcentration of vacancies in the Kitaev spin liquid leads to apileup of low-energy modes which cause a distinctive power-law divergence in the fermionic DOS. Since the vacancies areknown to bind the fluxes of the emergent Z gauge field, wepropose an algorithm to construct the appropriate bound-fluxsector for each random-vacancy configuration.Dilute vacancies preserve most of the spin-liquid behaviorbut lead to distinct changes in the low-energy physics.First, vacancy-induced Majorana modes are accumulated ina low-energy peak of the density of states. The form of thispeak across a broad window at low energies is well fitted bya power-law DOS of the form N ( E ) ∼ E − ν with ν ≈ . .Consequently, the power characteristic to the ’pure’ Dirac dis-persion is lost, i.e., this smoking gun of a Z Dirac spin liquidis not robust to the inclusion of disorder. We remark that ourresults do not preclude a crossover to yet more intricate behav-ior at – possibly much – lower energies, as has been discussedfor graphene [41].Second, the low-energy modes in question include the qua-sivacancy modes and the quasilocalized modes familiar fromsite-diluted graphene. However, the Kitaev spin liquid has ad-ditional flux degrees of freedom which affect the low-energymodes. Furthermore, the IPR results and the real-space wave3functions indicate that the vacancy-induced modes are local-ized, especially in the presence of an external field breakingtime-reversal symmetry. When a bulk gap is opened by such afield, these localized modes survive in the gap, hybridize witheach other, and become disconnected from the bulk modes. Inparticular, we were able to show that a simple hybridizationpicture can largely account for the in-gap spectrum.Third, a flux-sector transition from the bound-flux sector tothe zero-flux sector is found when increasing the strength κ ofthe field. Unlike the bound-flux sector, the zero-flux sector hasno flux-induced modes that can form a band around E = 0 ,and the DOS is therefore gapped.Our work was mainly motivated by the experimental find-ings in the Kitaev spin liquid candidate H LiIr O [8]. Bydemonstrating a robust vacancy-induced divergence in theDOS, our work provides a basic explanation for the specificheat results in H LiIr O . In particular, the power-law scal-ing N ( E ) ∼ E − ν leads to a C/T ∝ T − ν divergence, and ournumerical results with good fit for ν ≈ . are consistent withthe experiment. This implies that such a functional form arisesrather robustly for the energy window under consideration.This effective power-law exponent ν ≈ . of the site-diluted system only changes weakly over a relatively largeenergy (or temperature) window with respect to the additionof bond or flux randomness. Hence, we argue that vacanciesplay a major role in the low-energy physics of the Kitaev spinliquid, which is somewhat surprising and complementary toprevious theories. In the future, it would be desirable to address thermally ac-tivated fluxes, though the concurrence of thermal and disorderaverages is a potential challenge for numerics. In addition,dynamical probes such as Raman or neutron spectroscopy andmagnetic susceptibility measurements can be used for observ-ing signatures of the vacancy-induced low-energy modes boththeoretically and experimentally. In the limit of extremely di-lute vacancies in a magnetic field, the Majorana zero modesbound to vacancy-induced fluxes are far away from each other,which points to the intriguing possibility to observe and po-tentially even manipulate Majorana zero modes in a magneticmaterial. ACKNOWLEDGEMENTS
We acknowledge discussions with Kedar Damle and JohnChalker. W.-H. Kao and N. B. Perkins acknowledges the sup-port from NSF DMR-1929311. This work was in part sup-ported by the Deutsche Forschungsgemeinschaft under grantsSFB 1143 (project-id 247310070) and the cluster of excel-lence ct.qmat (EXC 2147, project-id 390858490). The workof G. B. H. was supported by the Laboratory Directed Re-search and Development Program of Oak Ridge National Lab-oratory, managed by UT-Battelle, LLC, for the US Depart-ment of Energy. [1] A. J. Willans, J. T. Chalker, and R. Moessner, Phys. Rev. Lett. , 237203 (2010).[2] A. J. Willans, J. T. Chalker, and R. Moessner, Phys. Rev. B ,115146 (2011).[3] J. Knolle, Dynamics of a Quantum Spin Liquid (Springer,2016).[4] F. Zschocke and M. Vojta, Phys. Rev. B , 014403 (2015).[5] G. J. Sreejith, S. Bhattacharjee, and R. Moessner, Phys. Rev. B , 064433 (2016).[6] L. Savary and L. Balents, Phys. Rev. Lett. , 087203 (2017).[7] I. Kimchi, J. P. Sheckelton, T. M. McQueen, and P. A. Lee,Nat. Commun. , 4367 (2018).[8] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano,Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi,Nature , 341 (2018).[9] K. Slagle, W. Choi, L. E. Chern, and Y. B. Kim, Phys. Rev. B , 115159 (2018).[10] Y. Li, S. M. Winter, and R. Valent´ı, Phys. Rev. Lett. ,247202 (2018).[11] J. Knolle, R. Moessner, and N. B. Perkins, Phys. Rev. Lett. ,047202 (2019).[12] S. K. Takahashi, J. Wang, A. Arsenault, T. Imai, M. Abram-chuk, F. Tafti, and P. M. Singer, Phys. Rev. X , 031047 (2019).[13] S.-H. Do, C. H. Lee, T. Kihara, Y. S. Choi, S. Yoon, K. Kim,H. Cheong, W.-T. Chen, F. Chou, H. Nojiri, and K.-Y. Choi,Phys. Rev. Lett. , 047204 (2020).[14] M. G. Yamada, arXiv:2004.06257 (2020).[15] J. Nasu and Y. Motome, arXiv:2004.07569 (2020). [16] H. Yamaguchi, M. Okada, Y. Kono, S. Kittaka, T. Sakakibara,T. Okabe, Y. Iwasaki, and Y. Hosokoshi, Sci. Rep. , 16144(2017).[17] M. A. de Vries, D. Wulferding, P. Lemmens, J. S. Lord, A. Har-rison, P. Bonville, F. Bert, and P. Mendels, Phys. Rev. B ,014422 (2012).[18] K. Mehlawat, G. Sharma, and Y. Singh, Phys. Rev. B ,134412 (2015).[19] J. Paddison, M. Daum, Z. Dun, G. Ehlers, Y. Liu, M. Stone,H. Zhou, and M. Mourigal, Nat. Phys. , 117 (2017).[20] Y. Li, D. Adroja, R. I. Bewley, D. Voneshen, A. A. Tsirlin,P. Gegenwart, and Q. Zhang, Phys. Rev. Lett. , 107202(2017).[21] J.-J. Wen, S. M. Koohpayeh, K. A. Ross, B. A. Trump, T. M.McQueen, K. Kimura, S. Nakatsuji, Y. Qiu, D. M. Pajerowski,J. R. D. Copley, and C. L. Broholm, Phys. Rev. Lett. ,107206 (2017).[22] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Annu. Rev. Condens.Matter Phys. , 195 (2016).[23] S. Trebst, arXiv:1701.07056 (2017).[24] M. Hermanns, I. Kimchi, and J. Knolle, Annu. Rev. Condens.Matter Phys. , 17 (2018).[25] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E.Nagler, Nat. Rev. Phys. , 264 (2019).[26] Y. Motome and J. Nasu, J. Phys. Soc. Jpn , 012002 (2020).[27] A. Kitaev, Annals of Physics , 2 (2006).[28] 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 ,041112 (2014). [29] M. Majumder, M. Schmidt, H. Rosner, A. A. Tsirlin, H. Ya-suoka, and M. Baenitz, Phys. Rev. B , 180401 (2015).[30] R. D. Johnson, S. C. Williams, A. A. Haghighirad, J. Singleton,V. Zapf, P. Manuel, I. I. Mazin, Y. Li, H. O. Jeschke, R. Valent´ı,and R. Coldea, Phys. Rev. B , 235119 (2015).[31] J. A. Sears, M. Songvilay, K. W. Plumb, J. P. Clancy, Y. Qiu,Y. Zhao, D. Parshall, and Y.-J. Kim, Phys. Rev. B , 144420(2015).[32] A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li,M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle,S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant,M. D. G., and S. E. Nagler, Nat. Mater. , 733 (2016).[33] F. Bahrami, W. Lafargue-Dit-Hauret, O. I. Lebedev,R. Movshovich, H.-Y. Yang, D. Broido, X. Rocquefelte,and F. Tafti, Phys. Rev. Lett. , 237203 (2019).[34] J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. Lett. ,197205 (2014).[35] J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. B , 115122(2015).[36] J. Yoshitake, J. Nasu, and Y. Motome, Phys. Rev. Lett. ,157203 (2016).[37] V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. R.Peres, and A. H. Castro Neto, Phys. Rev. Lett. , 036801 (2006).[38] V. M. Pereira, J. M. B. Lopes dos Santos, and A. H. Cas-tro Neto, Phys. Rev. B , 115109 (2008).[39] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Rev. Mod. Phys. , 109 (2009).[40] V. H¨afner, J. Schindler, N. Weik, T. Mayer, S. Balakrishnan,R. Narayanan, S. Bera, and F. Evers, Phys. Rev. Lett. ,186802 (2014).[41] S. Sanyal, K. Damle, and O. I. Motrunich, Phys. Rev. Lett. ,116806 (2016).[42] E. H. Lieb, Phys. Rev. Lett. , 2158 (1994).[43] G. B. Hal´asz, J. T. Chalker, and R. Moessner, Phys. Rev. B ,035145 (2014).[44] R. Yadav, R. Ray, M. S. Eldeeb, S. Nishimoto, L. Hozoi, andJ. van den Brink, Phys. Rev. Lett. , 197203 (2018).[45] K. Geirhos, P. Lunkenheimer, M. Blankenhorn, R. Claus,Y. Matsumoto, K. Kitagawa, T. Takayama, H. Takagi,I. K´ezsm´arki, and A. Loidl, Phys. Rev. B , 184410 (2020).[46] D. A. Ivanov, Phys. Rev. Lett. , 268 (2001).[47] V. Lahtinen and J. K. Pachos, Phys. Rev. B , 245132 (2010).[48] V. Lahtinen, A. W. W. Ludwig, J. K. Pachos, and S. Trebst,Phys. Rev. B86