Vacancy-induced low-energy density of states in the Kitaev spin liquid

The Kitaev honeycomb model has attracted significant attention due to its exactly solvable spin-liquid ground state with fractionalized Majorana excitations and its possible materialization in magnetic Mott insulators with strong spin-orbit couplings. Recently, the 5d-electron compound H$_{3}$LiIr$_{2}$O$_{6}$ has shown to be a strong candidate for Kitaev physics considering the absence of any signs of a long-range ordered magnetic state. In this work, we demonstrate 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 measurements of H$_{3}$LiIr$_{2}$O$_{6}$. Physically, the vacancies can originate from various sources such as missing magnetic moments or the presence of non-magnetic impurities (true vacancies), or from local weak couplings of magnetic moments due to strong but rare bond randomness (quasivacancies). We show numerically that the vacancy effect is readily detectable even at low vacancy concentrations and that it is not very sensitive neither to nature of vacancies nor to different flux backgrounds. We also study the response of the site-diluted Kitaev spin liquid to the three-spin interaction term, which breaks time-reversal symmetry and imitates an external magnetic field. We propose a field-induced flux-sector transition where the ground state becomes flux free for larger fields, resulting in a clear suppression of the low temperature specific heat. Finally, we discuss the effect of dangling Majorana fermions in the case of true vacancies and show that their coupling to an applied magnetic field via the Zeeman interaction can also account for the scaling behavior in the high-field limit observed in H$_{3}$LiIr$_{2}$O$_{6}$.


I. INTRODUCTION
Various types of disorder in quantum spin liquids (QSLs) have recently attracted a lot of attention from both experimental and theoretical points of view [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. There are three main reasons for this interest. First, some level of disorder in various forms of dislocations, vacancies, impurities, and bond disorder is inevitable in real materials. Second, disorder can significantly affect the low-energy properties of these systems. In particular, if the system is close to a QSL state, quenched disorder on top of the quantum disordered strongly correlated spin state of a QSL can give rise to diverse and often puzzling behaviors [8,12,13,16]. Namely, sometimes disorder is detrimental for the QSL state since it either localizes the resonating spin singlets or induces competing glassy states instead of entangled ones [17][18][19][20]. However, in some other cases, e.g., in classical spin ice materials, certain forms of disorder can instead enhance the quantum dynamics of spins throughout the system and generate a QSL with long-range entanglement [6,16,21]. Third, given that the properties of QSLs are difficult to detect directly because such states lack any local order parameter, much additional information can be obtained by studying the distinctive responses to local perturbations, such as static defects, dislocations, and magnetic or non-magnetic impurities. In particular, these perturbations may nucleate exotic excitations characteristic to the spin liquid under consideration [1,2].
Of specific interest is the role of disorder in the materials * nperkins@umn.edu that have been suggested to be potential candidates [22][23][24][25][26] to realize the Kitaev QSL [27]. In a flurry of recent experiments on the honeycomb ruthenium chloride α-RuCl 3 , it was shown that both bond disorder and stacking disorder are not negligible [28][29][30][31][32]. Perhaps, disorder also plays a crucial role for a potential proximity of Ag 3 LiIr 2 O 6 to a Kitaev QSL state [33]. However, perhaps the most remarkable and intriguing consequences of disorder have been observed in a presumptive quantum spin liquid state of the hydrogen intercalated iridate H 3 LiIr 2 O 6 [8]: (i) the specific heat displays a lowtemperature divergence of C/T ∝ T −1/2 ; (ii) only a small fraction of the total magnetic entropy is released at these low temperatures; and (iii) there is a non-vanishing contribution down to the lowest temperature in the NMR rate 1/T 1 and an almost flat Knight shift. All of these observations signal the presence of abundant low-energy density of states (DOS) related to magnetic excitations. However, despite the presence of dominant Kitaev exchange, this phenomenology is at odds with the thermodynamics of the pure Kitaev QSL [34][35][36], which has a vanishing specific heat C/T ∝ T and a significant release of half of its total entropy at low T . Motivated by these experimental findings, some of us have recently considered a minimal model of a bond disordered Kitaev QSL [11] that can account for these salient experimental observations in H 3 LiIr 2 O 6 [8]. However, in order to recover the low temperature scaling of the specific heat, the Kitaevlike model of Ref. 11 assumed a somewhat ad-hoc form of binary bond disorder and invoked a random flux background even at very low temperatures.
In this work, we show that a finite density of vacancies in the Kitaev model induces a pileup of the low-energy DOS, N (E), and a low temperature divergence of the specific heat, C/T . These results are consistent with an algebraic divergence with an exponent around ν = 1/2 over a broad range of low energies/temperatures. As a finite density of static, randomly-located vacancies is always present in the two dimensional limit of layered materials, i.e., similar to the case of graphene [37][38][39][40][41], our simple model provides a natural explanation for the experimental observations in H 3 LiIr 2 O 6 [8].
We also carefully treat the flux background and show that the energy of a random-vacancy configuration is minimized when a flux is bound to each vacancy. Hence, we resolve the problem of a random flux background since now the flux configuration is determined by the vacancy distribution. Finally, we show numerically that the vacancy effect for the lowenergy DOS is robust with respect to the addition of bondrandomness or a random flux background.

II. THE MODEL
The extended Kitaev honeycomb model is defined in terms of localized spin 1/2 degrees of freedom that are coupled in a bond-anisotropic manner on the honeycomb lattice [27]. The spin Hamiltonian reads whereσ α i denotes Pauli spin operators with α = x, y, z and ij α labels the nearest-neighbor sites i and j along an αtype bond. The second term is the three-spin interaction with strength κ ∼ hxhyhz J 2 that imitates an external magnetic field and breaks time-reversal symmetry while preserving the exact solvability [27]. By rewriting each spin operator in terms of four Majorana fermions,σ α i = ib α iĉ i , and defining the link operatorsû ij = ib α ib α j , the Hamiltonian takes the form The solvability of the Kitaev model relies on the extensive number of conserved fluxes defined on each hexagonal plaquette,Ŵ p =σ x 1σ y 2σ z 3σ x 4σ y 5σ z 6 = ij ∈pû ij α , which can block-diagonalize the Hamiltonian (1) into flux sectors since fluxes commute with each other, [Ŵ p ,Ŵ p ] = 0, and with the Hamiltonian, [Ŵ p , H] = 0. Both the flux operatorsŴ p and the link operatorsû ij α have eigenvalues ±1. Once the link variable is specified for each bond, the physically relevant flux sector is determined, and the Hamiltonian can be solved exactly as a tight-binding model of Majorana fermions. For the pure Kitaev model (κ = 0), it has been proven that the groundstate 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 contains two sublattice sites labeled as A and B. Thus, for a system with N unit cells, there are 2N lattice sites and N hexagonal plaquettes. Under periodic boundary conditions (PBC), the fluxes can be excited only in pairs since flipping one link variable in the zero-flux sector results in W p = −1 on both sides of that link. As a result, there exists a global constraint for the fluxes, p W p = 1, such that the number of independent fluxes is reduced by one. Since two additional flux degrees of freedom are introduced by the toric topology, the total number of different flux sectors is then 2 N +1 .
With this decomposition, the spin degrees of freedom in the original Hamiltonian are now fractionalized into itinerant Majorana fermions and static Z 2 gauge fluxes. In a given flux sector, the gauge can be fixed and all link variables can be specified, leading to a Hamiltonian of non-interacting Majorana matter fermions, where the hopping amplitudes between sites on sublattice A and sublattice B are M ij = J αû ij α , while the entries in the diagonal blocks F and −D contain hopping amplitudes between sites on the same sublattice. Two adjacent Majorana fermions from the same unit cell can be combined into a complex matter fermion, such that the Hamiltonian can be written in a Bogoliubov−de-Gennes form and diagonalized in the standard way, . Therefore, for a given flux configuration, the fermionic ground-state energy reads E 0 ({u ij }) = − 1 2 n n , and the global density of states is given by

A. Vacancies, quasivacancies, and fluxes
A vacancy is usually a simple absence of an atom at a given site. However, in the present work, we will use this term more generally since it can also correspond to a non-magnetic impurity or a magnetic moment that is weakly connected to its neighbors due to extremely strong but relatively rare bond randomness. To distinguish it from the simple absence of an atom (which we call a true vacancy), we will refer to the latter type of defect as a quasivacancy.
In order to introduce randomly distributed vacancies into the Kitaev honeycomb model (2), we first consider the timereversal symmetric case at κ = 0. In this case, the second term in Eq. (2) is absent, while the first term can be written as where P denotes the subset of normal lattice sites and V denotes the subset of vacancy sites. We consider a compensated case with equal numbers of vacancies on the two sublattices of the honeycomb lattice. By taking the limit of J α J α , sites belonging to V behave as quasivacancies. For such a quasivacancy, a Majorana fermionĉ remains on the vacancy site, but its nearest-neighbor hopping amplitudes are removed in the limit of J α → 0. Therefore, the way we diagonalize the pure Kitaev Hamiltonian remains valid, even though the number of flux degrees of freedom is effectively reduced.
For each vacancy with J α → 0, there are three hexagonal plaquettes around the vacancy site, resulting in 2 3 = 8 distinct flux sectors labeled by W 1,2,3 = ±1. However, the Majorana Hamiltonian in each flux sector is completely determined by , which corresponds to a large vacancy plaquette merged from the three individual plaquettes around the vacancy. Since the remaining two degrees of freedom do not affect the Majorana DOS and can even be partially "gauged away" in the case of true vacancies [43], we ignore them in the rest of this work and characterize the vacancy with a single vacancy fluxŴ v =Ŵ 1Ŵ2Ŵ3 . Consequently, for a system with N v isolated vacancies, the number of flux degrees of freedom is effectively reduced by 2N v . With periodic boundaries, the global constraint for the fluxes then becomes p W h,p q W v,q = 1, where W h correspond to bulk hexagonal plaquettes. Since there are two physically distinct flux operatorsŴ h andŴ v , introducing a flux pair falls into one of three situations: both fluxes on hexagonal plaquettes, both fluxes on vacancy plaquettes, and one flux on each kind of plaquette. We will show that, in order to minimize the total energy, fluxes must be bound to the vacancy plaquettes.

B. Quasilocalized eigenmodes and flux binding
When considering vacancies in the Kitaev model, many results are completely analogous to those in graphene [37][38][39][40][41]. For example, it was found for both systems that introducing a vacancy leads to a zero-energy eigenmode with a quasilocalized wavefunction on the other sublattice around the vacancy site [1,2,37,38].
Indeed, if we consider only the nearest-neighbor couplings and ignore the possible flux excitations in the Kitaev model, then there is a one-to-one mapping between the Majorana hopping in the Kitaev model and the fermionic hopping in graphene (see Fig. 1 (b)).
It was shown by Willans et al. through analytical calculations that a single vacancy binds a flux in the gapped phase of the Kitaev model (J x , J y J z ) [1,2]. In the gapless phase (including the isotropic point J x = J y = J z ), the flux-binding effect can be verified numerically. Practically, we consider two vacancies, one on an A sublattice site and the other on a B sublattice site (see Fig. 1 (a)), which are separated by a distance ∼ L/2, where L is the linear dimension of the system, and then calculate the energy difference between the boundflux ( Fig. 1 (c)) and the zero-flux sectors ( Fig. 1 (a)): In Fig. 2 (a), different system sizes up to L = 120 are considered and, by extrapolation, we show that the flux-binding energy converges to E bind = −0.0268J, which is consistent with the previous result [1]. In the same way, we also calculate the binding energies for non-zero J and show (see Fig. 2 (b)) that the flux-binding effect remains for J /J < 0.0544, which can be compared to another previous result with a slightly different setup [4]. The flux-binding effect for non-zero J indicates that we can extend the quasivacancy picture to a bond-disordered model where the vacancy sites are not truly removed or replaced by non-magnetic ions but the coupling strengths are strongly suppressed by structural disorder.
Note that, in this setting, the boundary conditions play a crucial role in the flux binding effect. For open boundary conditions, it is possible to create a single flux on each vacancy plaquette so that the energy is minimized. However, for periodic boundary conditions, the fluxes must be created in pairs. For a system with only one vacancy, there exists only one vacancy plaquette to bind the flux and thus the other flux must be bound to a hexagonal plaquette. This arrangement results in a higher total energy since the flux excitation energy on a hexagonal plaquette is 0.1536J [27] which can not be fully compensated by E bind . Therefore, the ground-state sector is still flux free. In contrast, when the system contains an even number of isolated vacancies, it is possible to bind fluxes to all vacancy plaquettes, as we describe in the following.

C. Bound-flux sector
In order to minimize the total energy of a random-vacancy configuration, fluxes must be introduced and bound to each vacancy plaquette. The most unbiased approach is to apply a Markov chain Monte Carlo simulation that samples the flux configurations at low temperatures [34,35]. However, this approach is not easily realized for large systems with disorder since the tight-binding Hamiltonian (2L 2 × 2L 2 matrix) must be diagonalized for each update. Thus, here we discuss how to generate the low-energy bound-flux sector in which each vacancy has a flux attached to it. Note that when two or more vacancies are connected to each other and thus the discussion in Sec. III A is no longer valid, our approach may not provide the actual ground-state flux sector. However, in the dilute limit, most vacancies are isolated and the ground state is well approximated with the bound-flux sector.
In all our calculations, we use periodic boundary conditions because open boundaries lead to additional zero-energy eigenmodes and make it difficult to examine the direct consequences of adding vacancies into a finite-size system. As discussed in the previous subsection, the consequence of PBC is 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 vacancy plaquette, one flux emerges on this vacancy plaquette and the other one emerges on a neighboring hexagonal plaquette. The former changes the energy by ∆E v = E bind < 0, while the latter changes the energy by ∆E f > 0.
Numerical studies of the single-vacancy effect show that the energy decrease of a flux binding to a vacancy cannot compensate the energy increase by the other flux [2]. Indeed, at the isotropic point, ∆E v converges to −0.0268J, while ∆E f converges to 0.1536J. In contrast, for a system with two isolated vacancies, it is possible to generate a flux pair on hexagonal plaquettes, and then propagate the two single fluxes individually until they bind to the vacancies, lowering the total energy by approximately 2|∆E v | compared to the zero-flux case (see Fig. 3 However, finding the ground-state flux configuration for more than two vacancies in terms of link variables is a nontrivial task because fluxes may create or annihilate during the flipping of link variables. Instead, we propose an algorithm that generates a vacancy configuration along with the bound fluxes. First, we randomly choose a position on the lattice and place a vacancy pair. The pair contains one vacancy on an A-sublattice site and the other one on a nearby B-sublattice site, as shown in Fig. 3(a). One shared link of the two vacancy plaquettes is flipped such that each plaquette binds a flux. Second, we randomly move one vacancy in the directions of the two primitive vectors. When a vacancy migrates, the link variables are flipped along the same path, keeping the flux bound to the vacancy plaquette ( Fig. 3 (b)). Thus, instead of propagating fluxes to minimize the energy of a given vacancy configuration, as is shown in Fig. 1 (c), we propagate composites of one vacancy and one flux to generate a random configuration of vacancies with bound fluxes.
Note that this method does not give the ground-state flux sector for all vacancy configurations. When two or more vacancies are connected after migration, the binding fluxes may annihilate each other such that no flux binds to the merged vacancy plaquette. Nevertheless, this method works relatively well for a low density of vacancies since connected vacancies are then rare. In the following sections, we denote the finiteflux sector generated by this method as the bound-flux sector.
D. Time-reversal symmetry broken case (κ = 0) The effective three-spin interaction term in Eq. (1) breaks time-reversal symmetry and changes the energetics of the model by simultaneously gapping out the fermionic spectrum and introducing localized zero-energy Majorana modes in the presence of isolated fluxes [27]. Here we test the disorderaveraged total energy of the system in the bound-flux and zero-flux sectors for various κ. By comparing the energies of the two flux sectors, a ground-state transition from the bound-flux to the zero-flux sector is observed. This transition is shown in Fig. 4, where we plot the difference between the two energies, each obtained by averaging over 4000 disor-  dered samples with 2% quasivacancies at J = 0.01. While the bound-flux sector is lower in energy for 0 ≤ κ ≤ 0.1, the zero-flux sector becomes energetically favorable for κ ≥ 0.1.

IV. DENSITY OF STATES AND SPECIFIC HEAT
In this section, we discuss how the presence of vacancies results in a low-temperature divergence in the specific heat, C/T , which might be related to the recent experimental observation by Kitagawa [8]. We consider both the case of true vacancies with J = 0 and the case of quasivacancies with J > 0.
At finite temperatures, both itinerant Majorana fermions and fluxes contribute to the specific heat and 1 2 ln 2 in the thermal entropy [34,35]. However, in H 3 LiIr 2 O 6 , it is reported that the low-energy excitations release only 1%-2% of ln 2 entropy at 5K, implying that the flux degrees of freedom might be frozen. Therefore, we assume that the temperature dependence of the specific heat is solely due to the thermal occupation of itinerant Majorana fermions in static flux sectors, where n F (E, T ) = (e E/T + 1) −1 is the Fermi function, and we use the definition of the density of states in Eq. (6) to reach the final expression.
As before, we first consider the time-reversal symmetric case with κ = 0 in Secs. IV A and IV B, and then discuss the effect of a finite κ in Sec. IV C.

A. True vacancies
We introduce a certain amount of vacancies into the Majorana problem, as described by Eq. (7), at the isotropic point (J x = J y = J z ≡ J = 1) and see how the fermionic specific heat and density of states are affected. The vacancy concentration n v is defined as the number of vacancies (N v ) divided by the number of lattice sites (2L 2 ). Half of the vacancies are on the A sublattice and the other half are on the B sublattice. For example, in a system with L = 20 and n v = 2%, we randomly put 8 vacancies on A sites and 8 vacancies on B sites. When generating the random vacancy configurations, we avoid removing the same site twice or more. Thus, for a given concentration, each disorder realization contains the same amount of vacancies. In the remainder of the paper, the results are presented for systems with linear dimension between L = 20 and L = 40 on a torus, and all the data are averaged over 10 3 to 10 4 disorder realizations.
First, we demonstrate the pileup of low-energy states in a system with 2% true vacancies. The density of states in the bound-flux sector, averaged over 4000 realizations (L = 40), is shown in Fig. 5. Since the low concentration of vacancies acts like a weak disorder on top of the Kitaev spin liquid, the overall behavior of the density of states is similar to the analytical result for the pure Kitaev model in the zero-flux sector, except for the low-energy region. Note that the states at exactly zero energy are removed from the density of states, so that the pileup in the low-energy region is exclusively from states with small but non-zero energies. The upturn is clearly seen in the log-log plot (inset of Fig. 5), and it fits well to a power-law form N (E) ∼ E −ν with ν ≈ 1/2.
Next, the fermionic specific heat (9) is calculated from the density of states and shown in Fig. 6. As expected, the pileup of low-energy states approximated by the power law The system-size dependence of C/T . The low-temperature upturn becomes more robust for larger systems due to the reduction of the finite-size effect. The curves are averaged over 4000-10000 disorder realizations, depending on the system size.
N (E) ∼ E −ν brings about a similar behavior in the specific heat, C/T ∝ T −ν with an exponent close to 0.5. This powerlaw behavior is consistent with what has been found in the Kitaev spin-liquid candidate H 3 LiIr 2 O 6 [8], where the authors mentioned that an approximately 2% density of magnetic impurities or vacancies in the material may be responsible for the magnetization and specific heat results. Thus, even without the presence of bond randomness [11], the low-energy fermionic states produced by the vacancies in the ground-state flux sector can give rise to a similar power-law behavior.
In Fig. 6 (b), the size dependence of this power-law upturn is presented. Because of the gapless nature of the Kitaev spin liquid at the isotropic point, the finite-size effects are considerable at low temperatures. We found that L = 40 is a reasonable choice in practice such that the finite energy of the vacancy states can be extended to the scale of 10 −3 . In the rest of this paper except Sec. V, systems with L = 40 will be In addition to the bound-flux sector, we also considered the zero-flux and random-flux sectors for the same set of randomvacancy configurations. The corresponding densities of states and fermionic specific heats are presented in Fig. 7 for various vacancy concentrations. The pileup of vacancy-induced states (Fig. 7, upper row) and the corresponding upturn in the specific heat ( Fig. 6 (a) and Fig. 7, lower row) appears in all the flux sectors, indicating that the effect of vacancies plays a major role in the low-energy region. Besides, we also see that the upturn power is slightly dependent on the vacancy concentration. In the bound-flux and zero-flux sectors, densities of states with different n v start splitting below a characteristic energy scale, which is similar to the tight-binding model of graphene with compensated vacancies [40].

B. Quasivacancies
In order to study the effect of quasivacancies introduced in the Kitaev model, we computed the density of states for the model (7) with different coupling strengths J α = J up to 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 the energetic analysis shown in Fig. 2 (b), the bound-flux sector is lower in energy than the zero-flux and the random-flux sectors for all vacancy couplings J ≤ 0.05.
In the limit of J → 0, the resonant peak around zero energy is present in all three cases, indicating that the lowenergy physics is governed by vacancies rather than fluxes. In the bound-flux sector, a finite value of J leads to a coupling of the quasivacancy mode to its surroundings and results in a larger width of the zero-energy peak. In the zero-flux sector, however, a tiny energy gap opens when increasing the magnitude of J . This phenomenon comes from the hybridization of the two different zero-energy modes corresponding to the same quasivacancy. When switching on the coupling J , the quasivacancy mode ( Fig. 1(a)) begins to hybridize with the zero-energy quasilocalized mode ( Fig. 1(b)), leading to a splitting of the corresponding energy levels. The formation of this pseudo-gap is also reported in site-diluted graphene [38]. The hybridization picture of low-energy modes will become even more clear when we open a larger bulk gap by adding the time-reversal symmetry breaking term to the system.

C. Three-spin interaction
In the previous sections we have demonstrated how the vacancy-induced low-energy states give rise to upturns in both the fermionic DOS and the specific heat C/T , similar to the experimental findings in H 3 LiIr 2 O 6 [8]. Now we turn our attention to the remaining question from the experiment: how does this low-energy upturn get suppressed in the presence of an external magnetic field? As discussed above, the threespin interaction with strength κ ∼ hxhyhz J 2 in Eq. (1) represents the leading-order perturbation effect of the Zeeman term [27]. This interaction breaks time-reversal symmetry and introduces zero-energy Majorana modes in the presence of fluxes. Therefore, the presence of vacancies in conjunction with the flux-binding effect provides a natural scenario for creating Majorana zero modes at low temperatures.
By introducing the three-spin term into the pure Kitaev honeycomb model, the gapless spin liquid becomes gapped due to the next-nearest-neighbor hopping of Majorana fermions. This effect is clearly seen in Fig. 9, where the fermonic den- sity of states is shown for different three-spin couplings κ in both the bound-flux and the zero-flux sectors for a 2% concentration of vacancies (recall that the bound-flux sector is the ground state flux sector only for 0 ≤ κ ≤ 0.1). Fig. 9 also clearly shows that there is a pronounced difference between the two flux sectors in the number of resonant peaks inside the bulk energy gap. In the zero-flux sector, two broad peaks appear inside the energy gap and move away from E = 0 with increasing κ. In contrast, one additional resonant peak is present around E = 0 in the bound-flux sector whose position is independent of κ. Note that, due to its finite width, this central peak contains a number of eigenmodes with small but nonzero energy, resulting in measurable signatures in thermodynamic quantities such as the specific heat at low temperatures. The presence or absence of this peak along with the flux-sector transition shown in Fig. 4 plays a crucial role in understanding the C/T results. A cartoon picture of the eigen-mode hybridization leading to different numbers of peaks in the two flux sectors will be discussed in Sec. V.
The localized nature of these vacancy-induced eigenmodes can be illustrated by the inverse participation ratio (IPR). This quantity is defined as 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 each eigenmode is shown by red dots. For a delocalized mode, the IPR scales roughly as ∼ 1/N in a system with N sites since the wavefunction is spread out uniformly over the entire lattice. This behavior is precisely what we see for the fermionic bulk modes. However, for the in-gap modes introduced by the vacancies, the IPR is significantly larger since the wave function is confined to a small portion of the lattice. Similarly to graphene, when κ = 0, each vacancy leads to a zero-energy eigenmode with a quasilocalized wave function on the other sublattice around the vacancy site. For a single vacancy, this wave function can be written in an analytical form [37,38]: The wavevectors K and K denote the two different Dirac points on the corner of the first Brillouin zone. While the analytical form of the quasilocalized wave function is no longer available for a finite density of vacancies, we can still relate the enhanced values of the IPR for the in-gap states in our numerical calculations to the 1/r decay of the quasilocalized wavefunctions Ψ(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 presented for various values of κ. According to the flux-sector transition shown in Fig. 4, the ground-state flux sector is the bound-flux sector for κ < 0.10 and the zero-flux sector for κ ≥ 0.10. As discussed in the previous sections, the upturn of the curve for κ = 0 can be extended to very low temperatures for large system sizes, and it is comparable to the experimental result for H 3 LiIr 2 O 6 without magnetic field. For κ = 0.08, a small number of localized modes appear in the DOS (see Fig. 9 (a)), giving rise to a steeper upturn in C/T . Still, there is no suppression at the lowest temperatures due to the presence of the central resonant peak in the density of states. However, when κ exceeds the critical value and the ground-state flux sector becomes flux free, C/T only shows a small upturn and is then strongly suppressed. The small upturn comes from the in-gap resonant peak at E > 0 and the suppression is due to the lack of lower energy modes around E = 0. Note that, if a vacancy site is completely decoupled from the system, the corresponding vacancy mode has exactly zero energy and has no contribution to the specific heat. However, since we consider the quasivacancy scenario that possibly results from a bond disorder of Kitaev interactions, those quasivacancy modes have finite couplings via J and κ and can be hybridized with other localized modes to produce finite-energy eigenmodes, thus leading to the low-temperature upturn in C/T . The quasivacancy picture and the corresponding C/T results capture the experimental findings for C/T in the Kitaev spin-liquid candidate H 3 LiIr 2 O 6 . In Figure 4 of Ref. 8, a peculiar scaling law is used for collapsing the C/T data in a wide range of magnetic fields h: Since we consider only the leading-order three-spin interaction κ emerging from a perturbative treatment of the magnetic field (rather the magnetic field h itself), a simple replacement of h with κ does not provide the correct scaling law for our data. Nevertheless, it is possible to assume a general powerlaw relation h ∼ κ α and test the following scaling behavior: We estimate the optimal α to be around 1.5 by collapsing the curves with various κ below the temperature scale T /κ α ∼ 1.2, as shown in Fig. 10 (b). Qualitatively, our calculation of the low-temperature fermionic specific heat is able to capture the specific heat scaling obtained experimentally in Ref. 8.

D. Effect of dangling Majorana fermions
In the previous subsection, we considered the leading-order effect of a magnetic field on quasivacancies, where the link variablesû ij = ib α ib α j emanating from the quasivacancies are well defined. However, for true vacancies, there are no Majorana fermions on the vacancy sites, and the link variablesû are thus no longer well defined around the vacancies. Consequently, the dangling Majorana fermionsb α on the neighboring 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 whereb α is a dangling Majorana fermion if the site j is connected to the vacancy site by an α-type bond. The Zeeman term can then be readily merged into the original tight-binding Hamiltonian as it is equivalent to a hopping term between the dangling Majorana fermionb α and the matter Majorana fermionĉ on the same site. Note that this term is a direct consequence of the magnetic field and is not derived from perturbation theory. As such, it provides the primary effect of a magnetic field in a system with true vacancies. Without loss of generality, this effect is demonstrated in  Fig. 11 for a magnetic field h z applied in the z direction. Energetic considerations show that the zero-flux sector becomes the ground-state flux sector for h z ≥ 0.4. Again, the clear suppression in C/T can then be attributed to the formation of an energy gap. Indeed, the dangling Majorana fermionsb z are gapped out through h z , which is similar to the discussion on Fig. 8 (b). Since the field strength h z represents the actual magnetic field instead of the three-spin coupling κ, we are able to see the scaling law C/T ∼ (h z ) −3/2 T from the data collapse in the inset of Fig. 11 (b). Interestingly, even though we consider the effect of the magnetic field only on the dangling Majorana fermions (but not on the bulk system), the peculiar scaling behavior can be reproduced in the highfield limit. More detailed results on the contributions of true vacancies to the thermodynamics and dynamical responses of the Kitaev model will be reported elsewhere. . The magnetic field is applied along the z axis so that it couples the dangling Majorana fermionb z to the rest of the system. The inset of (b) shows the 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 of H 3 LiIr 2 O 6 , it was proposed that bond disorder [10,11] may play a major role in generating the upturns found in both C/T and the density of states. The physical origin of bond randomness can be traced back to the random positions of the protons between the honeycomb layers, leading to a local distortion of the oxygen octahedral cage, a change in the Ir-O-Ir bonding angle, and hence the local variation of Kitaev couplings along the Ir-Ir links [10,44,45]. Here, we consider the combined effect of bond randomness and a 2% concentration of vacancies. A random coupling term is added to the nearest-neighbor coupling on the normal lattice sites (the first term in Eq. (7)): In the Gaussian bond disorder model, this replacement is applied to all the bonds except the bonds emanating from the va- cancy sites, and δJ is assigned randomly from a Gaussian distribution with mean value 0 and standard deviation σ J = 0.25. An additional constraint of J + δJ ≥ 0 is implemented to prevent couplings with mixed signs which would correspond to random flux insertion [3]. On the other hand, in the binary disorder case, the random coupling term is fixed to be δJ = ±0.8, and only 25% of the bonds are randomly chosen to be disordered. This implementation of binary disorder is consistent with the previous work [11], apart from using the random-flux sector. For both true vacancies and quasivacancies in our system, the bound-flux sector has lower energy than the zero-flux and random-flux sectors. Thus, the results presented in Fig. 12 are all calculated for the bound-flux sector.
For true vacancies, the additional bond randomness only leads to a small change in the power-law exponent of the specific heat C/T , implying that the primary cause of the lowtemperature upturn is the presence of vacancies. However, a non-zero quasivacancy coupling in conjunction with bond randomness can result in a more noticeable change up to a power law C/T ∼ T −1 , indicating that the distinction between true vacancies and quasivacancies may be of great importance when studying the Kitaev spin-liquid materials.

V. HYBRIDIZATION OF LOW-ENERGY MODES
For a large enough three-spin interaction κ, we can see in Fig. 9 that the number of broad peaks inside the energy gap is not the same in the bound-flux sector as in the zero-flux sector. The spectrum of those in-gap states can be understood by considering the hybridization among the low-energy localized modes introduced by the quasivacancies.
Let us first discuss the zero-flux sector. For κ = 0, each quasivacancy induces two quasilocalized modes around its vacancy site. The first one is the fully localized vacancy mode (v-mode) which can couple to its neighbors through the weak coupling J (Fig. 1(a)). The second one is the quasilocalized mode whose wave function is restricted to the opposite sublattice and decays as 1/r with the distance r from the vacancy site ( Fig. 1(b)). For κ = 0, this quasilocalized mode extends to both sublattices and becomes properly localized with most of its wave function distributed over the periphery of the vacancy plaquette (p-mode). For a large enough gap, the in-gap spectrum can then be approximated with the hybridization of these localized modes. For example, we can consider a simple model with only two vacancies and four localized modes: where c v,i is the v-mode Majorana andc p,i is the p-mode Majorana corresponding to the ith vacancy (i = 1, 2). We add the tilde on the p-mode Majorana to highlight that its wave function is not confined to a single site. Two energy scales µ and η are defined to represent the couplings between localized modes around the same vacancy and around different vacancies, respectively. Importantly, µ increases linearly with both J and κ , while η decays exponentially with the distance between the two vacancies. Thus, for well-separated vacancies, the first energy scale is much larger than the second one: µ η. The eigenspectrum of δH zero then consists of two doublets with energies ±µ and a small splitting ∼ η within each doublet (see Fig. 13 (b)). We verified this simple picture by obtaining the exact energy levels of a single L = 20 system with only two random vacancies and confirming that the resulting in-gap spectrum (see Fig. 13 (a)) is consistent with the one obtained from δH zero (see Fig. 13 (b)). For a finite concentration of vacancies, there is further hybridization on the scale of η which broadens the two peaks at ±µ but leaves the 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 by the flux on each vacancy plaquette (f-mode). When the fluxes are far apart from each other, these modes become Majorana zero modes and can be interpreted as anyons with non-Abelian statistics [27,46]. However, if the fluxes are closer and interact with each other, these modes hybridise and broaden into a mini-band [47,48]. Thus, in the simple two-vacancy model of the low-energy subspace, we must consider the hybridization (e) (d) where the summations in a and b are over the three distinct types of modes (v, p, f ). As in the zero-flux case, the larger energy scale µ represents the couplings between modes around the same vacancy, while the smaller energy scale η represents the couplings between modes around different vacancies. The eigenspectrum of δH bound has three doublets at energies 0 and ±µ, which explains the presence of the additional zero-energy peak in the density of states (see Fig. 9 (a)). In general, the number of resonant peaks inside the bulk gap corresponds to the number of localized modes around each vacancy.
The distinction between the bound-flux and the zero-flux sectors in the presence or absence of the central peak leads to very different low-temperature behaviors in the fermionic specific heat. It also clearly distinguishes the vacancies in the Kitaev model from those in graphene.

VI. CONCLUSION
In this work, we have demonstrated that introducing a small concentration of vacancies in the Kitaev spin liquid leads to a pileup of low-energy modes which cause a distinctive powerlaw divergence in the fermionic DOS. Since the vacancies are known to bind the fluxes of the emergent Z 2 gauge field, we propose an algorithm to construct the appropriate bound-flux sector for each random-vacancy configuration.
Dilute vacancies preserve most of the spin-liquid behavior but lead to distinct changes in the low-energy physics.
First, vacancy-induced Majorana modes are accumulated in a low-energy peak of the density of states. The form of this peak across a broad window at low energies is well fitted by a power-law DOS of the form N (E) ∼ E −ν with ν ≈ 0.5. Consequently, the power characteristic to the 'pure' Dirac dispersion is lost, i.e., this smoking gun of a Z 2 Dirac spin liquid is not robust to the inclusion of disorder. We remark that our results do not preclude a crossover to yet more intricate behavior at -possibly much -lower energies, as has been discussed for graphene [41].
Second, the low-energy modes in question include the quasivacancy modes and the quasilocalized modes familiar from site-diluted graphene. However, the Kitaev spin liquid has additional flux degrees of freedom which affect the low-energy modes. Furthermore, the IPR results and the real-space wave functions indicate that the vacancy-induced modes are localized, especially in the presence of an external field breaking time-reversal symmetry. When a bulk gap is opened by such a field, these localized modes survive in the gap, hybridize with each other, and become disconnected from the bulk modes. In particular, we were able to show that a simple hybridization picture can largely account for the in-gap spectrum.
Third, a flux-sector transition from the bound-flux sector to the zero-flux sector is found when increasing the strength κ of the field. Unlike the bound-flux sector, the zero-flux sector has no 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 findings in the Kitaev spin liquid candidate H 3 LiIr 2 O 6 [8]. By demonstrating a robust vacancy-induced divergence in the DOS, our work provides a basic explanation for the specific heat results in H 3 LiIr 2 O 6 . In particular, the power-law scaling N (E) ∼ E −ν leads to a C/T ∝ T −ν divergence, and our numerical results with good fit for ν ≈ 0.5 are consistent with the experiment. This implies that such a functional form arises rather robustly for the energy window under consideration.
This effective power-law exponent ν ≈ 0.5 of the sitediluted system only changes weakly over a relatively large energy (or temperature) window with respect to the addition of bond or flux randomness. Hence, we argue that vacancies play a major role in the low-energy physics of the Kitaev spin liquid, which is somewhat surprising and complementary to previous theories.
In the future, it would be desirable to address thermally activated fluxes, though the concurrence of thermal and disorder averages is a potential challenge for numerics. In addition, dynamical probes such as Raman or neutron spectroscopy and magnetic susceptibility measurements can be used for observing signatures of the vacancy-induced low-energy modes both theoretically and experimentally. In the limit of extremely dilute vacancies in a magnetic field, the Majorana zero modes bound to vacancy-induced fluxes are far away from each other, which points to the intriguing possibility to observe and potentially even manipulate Majorana zero modes in a magnetic material.