Anderson–Kitaev spin liquid

The bond-disordered Kitaev model attracts much attention due to the experimental relevance in α-RuCl3 and A3LiIr2O6 (A = H, D, Ag, etc.). Applying a magnetic field to break the time-reversal symmetry leads to a strong modulation in mass terms for Dirac cones. Because of the smallness of the flux gap of the Kitaev model, a small bond disorder can have large influence on itinerant Majorana fermions. The quantization of the thermal Hall conductivity κxy/T disappears by a quantum Hall transition induced by a small disorder, and κxy/T shows a rapid crossover into a state with a negligible Hall current. We call this immobile liquid state Anderson–Kitaev spin liquid (AKSL). Especially, the critical disorder strength δJc1 ~ 0.05 in the unit of the Kitaev interaction would have many implications for the stability of Kitaev spin liquids.


INTRODUCTION
The Kitaev model 1 is one of the greatest examples of twodimensional (2D) solvable models of quantum spin liquids (QSLs) [2][3][4] , especially in the perspective of spin-orbital-entangled physics 5,6 . This model has a bond-dependent anisotropic interaction, which brings about exchange frustration and realizes gapped and gapless spin liquid states depending on its parameters. Amazingly, this interaction can be furnished in materials with a strong spin-orbit coupling 7 . Iridates and α-RuCl 3 are prominent examples of candidate materials for the Kitaev model [8][9][10] , and thus are called Kitaev materials. As a hallmark, the half-integer quantization of the thermal Hall conductance has been observed in α-RuCl 3 with a magnetic field 11 , suggesting the fractionalization of spin degrees of freedom predicted from the Kitaev model. However, it is also known that these honeycomb materials cannot fully be understood by the original (pure) Kitaev model 12,13 . While other diagonal or offdiagonal interactions might be important in real materials 14 , the importance of disorder has been ignored in these materials until recently [15][16][17] . Indeed, experiments in A 3 LiIr 2 O 6 (A = H, Ag, etc.) show a universal scaling in the field dependence of the heat capacity 5,18 , which strongly suggests the existence of disorder 19,20 . The candidate ground states must be disordered QSLs, and the absence of long-range order can be attributed to the critical role of disorder.
In fact, the role of disorder in QSLs itself is a long-standing problem because of the absence of a solvable model, except for limited cases 21 . We propose a disordered Kitaev model as a "numerically" solvable model for the disordered QSL, where we can treat the magnetic field effect within the perturbation theory. Thus, this study is not only a model investigation for the disordered Kitaev materials like A 3 LiIr 2 O 6 (A = H, D, Ag, etc.) 5,18,22 , but also a systematic examination of a numerically solvable disordered QSL, which would be an attempt towards the universal understanding of various disordered QSLs. Especially, since most QSLs are unsolvable, an unbiased study of disordered QSLs was impossible in the previous method.
Moreover, the bond-disordered Kitaev model with an applied magnetic field is not an "analytically" solvable problem. While it seems that the problem is reduced to the Anderson localization of symmetry class D with a short-range correlated disorder 23 , we discovered that this is not the end of the story. The nontrivial part comes from the third-order perturbation in an applied field because the computation of the "disorder" in the time-reversalbreaking term requires diagonalization 1 . This nonlocal operation makes the disorder long-range correlated. Such a long-range correlated random mass disorder could be a relevant perturbation 24 , and a simple renormalization group analysis may fail. Treating the long-range correlated disorder essentially requires a large-scale calculation and the extrapolation to the thermodynamic limit, which were numerically difficult in the finite-size Monte Carlo methods 25,26 . For these reasons, we invented a powerful numerical method based on the kernel polynomial method (KPM) 27 to overcome the difficulty.
From the materials side, the dispersion of low-energy excitations of H 3 LiIr 2 O 6 is a mystery. Although the behavior of the heat capacity of H 3 LiIr 2 O 6 without a magnetic field 5 is partially explained by the bond-disordered Kitaev model 17 , the T 2dependence of the heat capacity of H 3 LiIr 2 O 6 with a magnetic field 5 , which suggests the linear dispersion of the low-energy density of states (DOS) 5 , has never been described by an unbiased calculation. Such heat capacity behaviors are phenomenologically explained by the random-singlet theory with a Dzyaloshinskii-Moriya (DM) interaction 19 . However, it is not clear whether the DM interaction is important in H 3 LiIr 2 O 6 , and symmetric Kitaev-Γ interactions may be a dominant consequence of the spin-orbit coupling. We rather try to explain the observed linear dispersion with an applied field from the bond-disordered Kitaev model. Indeed, an unbiased numerical simulation of the bond disorder is not only theoretically nontrivial but also experimentally important.
The Anderson transition (or quantum Hall transition) is visible in a Kitaev spin liquid (KSL) by applying a small magnetic field. The quantized thermal Hall conductivity serves as an order parameter, which must go to zero in the large-disorder limit. This limit without a linear slope of the thermal Hall conductivity is named Anderson-Kitaev spin liquid (AKSL). The transition between KSL and AKSL is nothing but a transition between different Hall plateaus. Although we defined the transition in 2D, the transition may be observable even in three-dimensional (3D) generalizations 28,29 (mobility edge).
In this article, we simulate the bond-disordered Kitaev model to see a crossover between KSL and AKSL, especially from the topological transition in the thermal Hall effect 11,25,30 . We discovered that the quantized thermal Hall effect is not as stable as expected from the Anderson transition with a short-range correlated disorder.

RESULTS
Magnetic field effect The pure Kitaev model on the honeycomb lattice can be defined from Fig. 1a. The bonds parallel to the red, green, and blue ones are x-, y-, and z-labeled bonds, respectively. The Hamiltonian is where 〈jk〉 means a nearest-neighbor (NN) bond, and J > 0. γ = x, y, or z is determined by a bond label. This model is known to be solvable by representing σ γ j by Majorana fermions ib γ j c j . This representation still works even if we introduce bond disorder as follows: where J jk = J ± δJ is a bond-dependent hopping, and δJ > 0 is the strength of bond disorder. This model is still numerically solvable if we can assume that the ground state is 0-flux when δJ is in the perturbative regime. Under this assumption, all the states with a vison, which is defined as a pair of neighboring π-flux vortices excited by flipping the sign of a single bond, are assumed to be the "first" excited states from the ground state flux sector. This is how the perturbation theory works for this Kitaev model. We employ Kitaev's trick to solve these Hamiltonians with an applied magnetic field 1 . Let us first review the perturbation for the pure Kitaev model with a magnetic field where h ¼ ðh x ; h y ; h z Þ t is an applied magnetic field. It is well known that V can be treated by the third-order perturbation 1 . The result after introducing itinerant Majorana fermions c j is X hhklii c k c l þ ðfour À fermion termsÞ: where α 0 = 0.262433 in the thermodynamic limit for the 0-flux state, and α 0 J is a vison gap in the uniform case. The determination of the prefactor follows a mean-field solution 31 .
The direction of the next-nearest-neighbor (NNN) bond 〈〈kl〉〉 is defined clockwise as shown in Fig. 1a around the site j. A site connected by the γ-bond from j is called γ[j] for γ = x, y, and z, as shown in Fig. 1a.
Kitaev's trick Next, let us include binary disorder as H = H bond + V. Following Kitaev 1 , we can always do perturbation from any random H bond by a formula where Π 0 is a projection onto the ground state flux sector, G 0 0 ðEÞ is an unperturbed Green function constructed from H bond with the ground state flux sector excluded from the Hilbert space, and E 0 is an initial energy. Since H bond is solvable by Majorana fermions, it is in principle possible to calculate G 0 0 ðEÞ numerically to exhaust every term appearing in the third order. For example, a Green function for excited states is efficiently obtained by KPM 27 numerically. However, this strategy is surely overkill for our problem.
A much simpler solution is to use a trick introduced by Kitaev. Though we still need an O(N 4 ) calculation cost to decide all terms by usual matrix diagonalization, where N is the number of sites, there is no need for matrix exponentiation or integration. Kitaev's trick is done by replacing G 0 0 ðE 0 Þ by − (1 − Π 0 )/Δ vison , assuming that the virtual state energy is constant determined just by a vison gap Δ vison . This is a bold approximation to simplify the problem drastically, but as we will see essential features, such as the modulation of the mass term, can be captured even within Kitaev's approximation.
In this way, a typical third-order term is as follows: whereκ kl depends on the intermediate site j in the third-order perturbation process. Seen from j,κ can be calculated by replacing where Δ γ is the energy of the vison excited by flipping the sign of the bond between j and γ[j].
We note that three bonds have the same value ofκ kl around j. Thus, the disorder simply modulate the mass term of Dirac cones via random NNN hoppings, and the problem is still solvable numerically.
In this case, four-fermion terms are short-ranged and we have just ignored them as we are only interested in the Hall conductivity in the h → 0 limit. Though we will assume the ground state of H bond to be 0-flux in the following discussions, the perturbation can be done from any flux configuration. We note that a second-order perturbation in h is ignored because it just renormalizes bond-dependent hoppings J jk and does not break the time-reversal symmetry. Indeed, there is a priori no way to determine the ratio of the coefficients of the second-and thirdorder perturbations, although we can always use a mean-field solution of the pure Kitaev model to estimate it 31 . As stressed in the Introduction, the computation by Kitaev's trick is an essential part of this study, producing long-range correlation in the mass disorder.
Thermal conductivity We only consider zero temperature and ignore thermal flux fluctuations above the 0-flux sector. Lieb's theorem 32 no longer applies, but we can expect it to be applicable on average. Anyway, the calculation is relevant only in the regime where the flux gap is not closed by thermal fluctuation or bond disorder (δJ < δJ c2 in Fig. 1b, c). We employed Kitaev's trick to calculate a Majorana spectrum with an external magnetic field for each quenched bond disorder. From this, we can compute an inplane thermal Hall conductivity κ xy (T), especially behavior of κ xy /T at T → 0. Here xy does not coincide with the Cartesian axis but means a transverse component of the thermal conductivity.
As described in Methods, the thermal Hall conductivity can be computed directly by the Kubo formula. However, we can alternatively use the so-called noncommutative Chern number (NCCN) 33 , which is defined by a spectral projector for occupied free fermions. This formula is advantageous because it is proven to become integer after disorder average with some conditions, whereas it only makes sense in the thermodynamic limit.
where P F is a spectral projector, i.e. the projection operator onto the occupied states, and L is a linear scale of the system. This formula must agree with the Kubo formula in the thermodynamic limit by a well-known relation κ xy =T ¼ πk 2 B Ch=ð12_Þ for Majorana fermions. We note that there are other ways to detect the topological nontriviality [34][35][36] .
After taking an average of Ch ∝ κ xy /T over a number of disorder configurations, we plot a physical thermal Hall conductivity as a function of δJ. The error bar is estimated from a statistical deviation. From now on we set h = k B = J = 1.

Numerical results
We first note that, since we only include the third-order perturbation, the results here are not simply comparable with experiments. However, it was proposed that the contribution from h x h y h z can be picked up by applying an inplane magnetic field 37 , so we only take an odd component under every sign change (h x ↦ − h x , h y ↦ − h y , and h z ↦ − h z ) of the three components of h from total κ xy . From now on we denote κ xy as an odd component under every sign change and ignore other components.
The approximate correspondence between the Kubo formula and NCCN is confirmed for the pure Kitaev model (see Fig. 2a). We note that Haar-random vectors are used in this calculation, but are not used in the following as described in Supplementary Methods and Supplementary Fig. 1. From here we will prefer NCCN to the Kubo formula because we can use KPM to approximate the spectral projector P F to avoid diagonalization. We note that the application of KPM to the Kubo formula requires efforts 38 .
Next, we would move on to a large-scale calculation by Kitaev's trick. We only take (Kitaev's) L × L periodic boundary condition (for spins) from L = 10, where the vison gap gets close to the thermodynamic limit. As long as we are interested in the topological property, the h → 0 limit does not have to be taken. We set h x ¼ h y ¼ h z h ¼ Δ min , where Δ min is the minimum vison gap for each disorder configuration. In order to reduce the finitesize effect, we adopt Kitaev's torus basis where the finite-size effects cancel out, which is defined from a torus basis (Ln 1 , Ln 2 + n 1 ) 1 . We call it Kitaev's periodic boundary condition (KPBC) for simplicity. The NCCN formula for KPBC has to be modified as described in Supplementary Methods. This arbitrary choice of boundary conditions does not matter in the thermodynamic limit. The averaged 〈κ xy 〉/T for T → 0 is shown as a function of δJ, and drops rapidly to 0 from the quantized value as the disorder strength δJ grows. From here 〈κ xy 〉/T is plotted in the unit of a quantum π/12. We used R = 24 vectors to approximate the trace 39 .
The mean and minimum values of vison gaps are plotted for each δJ in Fig. 2b. When δJ > δJ c2~0 .3, the vison gap approaches 0 for some plaquette, and the 0-flux ground state is destabilized. We tentatively define δJ c2 as the point where Δ min ¼ 0:05. From here, the perturbation from the 0-flux sector cannot be justified. Moreover, some flux sectors get almost degenerate and the firstorder perturbation in h now becomes relevant. Beyond this point, a quantized thermal Hall current is no longer a well-defined notion. Flux excitations and (itinerant) Majorana fermions are not separable, and the discussion based only on free Majorana fermions breaks down.
When δJ ≪ δJ c2 , the calculation by Kitaev's trick can be justified. Fig. 2c shows NCCN calculated by diagonalization (line plot) and KPM (scatter plot). These two methods agree well. From the data of KPM we extrapolated the thermodynamic limit. The finite-size data are fit by exponential functions, and extracted the converged value for L → ∞. The extrapolation is plotted in Fig. 2d  Local density of states When δJ ≪ δJ c2 , free Majorana fermions are only relevant lowenergy excitations, and we can use many tools of free fermions to discuss properties of the transition, such as DOS and a localization length. DOS around the ground state can be measured from the information of the 0-flux sector. As is often the case, we only calculated local density of states (LDOS), instead. LDOS for a site i is defined as follows 27 : where i j i is a site basis, and k j i is an eigenstate of the Hamiltonian with an energy ε = ε k .
The nonlocality of Majorana fermions does not matter as averaged LDOS approximates DOS well enough. Both of the quantities are easily computed using KPM, and LDOS is enough for our purpose. The ratio of the arithmetic and geometric means of LDOS also works as the order parameter of an Anderson transition instead of the localization length. From the gapped Dirac spectrum (see Fig. 3a) the LDOS becomes gapless as the disorder strength increases. In the gapless region, DOS behaves linearly around ε = 0 (see Fig. 3b-d). The localization in Fig. 3b-d is clear from the discrepancy between the arithmetic and geometric averages of LDOS. Details are included in Supplementary Methods.

DISCUSSION
Even in a finite-size calculation, the transition between KSL and AKSL was well-observed and the schematic phase diagram in Fig.  1c was confirmed. From the extrapolation, δJ c1 is very small and δJ c1 /J~0.05. δJ c1 /J ≪ α 0~0 .26 suggests the fragility of KSL. The fragility potentially explains the observed sample dependence of the half-integer quantized thermal Hall effect in α-RuCl 3 40 . There is a possibility that δJ c1 /J → 0 as h=Δ min ! 0, which means that the bond disorder in the Kiteav model is a relevant perturbation 24 .
After the transition, the V-shaped behavior of DOS completely agrees with an observed linear low-energy DOS for H 3 LiIr 2 O 6 with an applied magnetic field 5 . It was proposed that the spin-polarized scanning tunneling microscopy (SP-STM) could be a local probe for edge states of the Kitaev model 41 . During the crossover to AKSL, edge states should disappear gradually and thus SP-STM should be a good probe for the topological transition.
The h=Δ min ! 0 limit must be taken carefully. A small h=Δ min suffers from the finite-size effect, so we need to increase N according to the h=Δ min ! 0 limit, as described in Supplementary Methods and Supplementary Fig. 2. Whether or not δJ c1 /J → 0 in this limit is an important future problem. We note that the vortex disorder is known to be relevant, so the introduction of random vortices may change the universality 42 .
The fragility of the quantization has many implications to experiments. Disorder always exists in real materials, especially in any 2D layered system, and even in clean samples of α-RuCl 3 stacking faults must exist 43 . Thus, the situation is quite similar to that of the fractional quantum Hall effect (FQHE). The observation of FQHE requires a really clean sample, and the recently observed quantized thermal Hall current of FQHE is more sensitive to disorder 44 . The sensitivity also resembles unconventional superconductors 45 . It might be universal in strongly correlated systems. Thus, we need to reconsider the importance of cleanness for the topological order in general. Further discussions can be found in Supplementary Discussion.

Note added
After the early version of the present paper on arXiv, a finitetemperature calculation of the disordered Kitaev models has been reported 46 .

METHODS
We would like to introduce the KPM 27,47,48 . KPM is used twofold in this work. The first is for the vison gap calculation and the second is for the Chern number calculation. The calculation cost for the former is O(N 2 ) and for the latter O(N).

Vison gap calculation
First, let us consider a Majorana Hamiltonian with the following form.
where H is a Hermitian matrix. For Majorana fermions c j , H has a form H ¼ iA, where A is a real skew-symmetric matrix. From now on, we assume H to be the ones considered in the main text, either with or without a magnetic field. The eigenvalues of the N × N matrix H is denoted by E k with k = 1, …, N.
A Green function can be expanded by a Chebyshev polynomial T m (x) as follows: where λ = 4.0 was used in the Lorentz kernel g m . M is the expansion order and m = 0, …, M − 1. ε is a small parameter which goes to 0 when M → ∞.
The scaling factor s is necessary so that the spectrum of the Hamiltonian H=s falls within the domain of the Chebyshev polynomials [−1, 1]. We note that this expression is for diagonal components, but almost the same is true for offdiagonal components. We simply set s ¼ E max þ 0:1, where E max is the maximal eigenvalue. Elements of Chebyshev moments T m ðH=sÞ can be computed recursively by using cost is required to compute all the necessary elements.
From the expanded Green function, we can compute the energy change by the local modification of the Hamiltonian H ! H þ δH. We define By extending this function to a complex number, the energy difference, i.e. vison gap Δ vison , can be computed as follows: where F(E) is a Fermi-Majorana function at zero temperature.
The evaluation of the integral in the Green function requires fast Fourier transformation (FFT) or discrete cosine transformation (DCT) 27 . Using this, the integral is reduced to a discrete weighted sum of the Fermi(-Majorana) function evaluated at some specific points. The discretization sizeM for the summation was setM ¼ 2M for simplicity. We compared M = 64, 128, 256, 512, 1024, and 2048. Among these, we found that M = 1024 has the best performance for our purpose, where the error is always about 0.01J.
KPM can reproduce the vison gap with an enough precision and reduce the computational cost to O(N) with a truncation algorithm [49][50][51] . However, later we found that the truncation causes a problem in our simulation, and thus we used the abovementioned O(N 2 ) algorithm without a truncation 27,47,48 . Chern number calculation A Kubo formula for κ xy at zero temperature is reduced to the generalized Thouless-Kohmoto-Nightingale-den Nijs (TKNN) formula 52 for noninteracting Majorana Hamiltonians: where m and n label eigenvalues of H, ε m and ε n , corresponding to eigenstates m j i and n j i, respectively 53,54 . ϑ(x) is a Heaviside theta and v α ¼ i½H; r α =_ is a velocity operator along the α-direction. We define a position operator r α for the n α -direction for α = 1, 2. A gravitomagnetic term should be added to derive this formula. We note the generalized formula without a translation symmetry was originally discussed by Kitaev 1 using a flow of unitary matrices and extended by Kapustin and Spodyneiko 55 . This Kubo-TKNN formula is nothing but a real-space formulation of the Chern number calculation.
However, as mentioned in the main text, it is better to use a noncommutative Chern number (NCCN) Eq. (10) 33 to compute the same quantity in the thermodynamic limit. Though Eq. (10) makes sense only in the thermodynamic limit, we can make use of a discretized version, which exponentially converges to the thermodyanamic limit by an artificial kspace quantization of a size L × L. This version can be derived by replacing the commutator 33 as follows: Ài½r α ; P F 7 ! X Q q¼ÀQ c q e ÀiqΔrα P F e iqΔrα ; where Δ = 2π/L, c 0 = 0 and c q = −c −q are determined to hold x À P L=2 q¼ÀL=2 c q e iqΔx ¼ OðΔ L Þ, and Q ≤ L/2. Thus, we can expect that these two methods to compute κ xy may agree with a large L, while the Hall conductivity and the Chern number are a priori different quantities. We fixed Q = 15 for L > 30 because otherwise the calculation cost becomes O (N 3 ).
The spectral projector P F can be rewritten as: where C is a contour which encloses every negative eigenvalue of H. The integrant is nothing but a Green function, and can be expanded by KPM. It is better to use the following Fermi function instead of approximating a spectral projector directly.
Due to the particle-hole symmetry, this deformation also gives a correct Chern number. This expression can again be expanded by Chebyshev polynomials. Especially, the Fermi function has been expanded as where M 0 is a cutoff of the expansion for NCCN 39 . We used the Jackson kernel for the Chern number: where m ¼ 0; ; M 0 À 1. We here used M 0 ¼ 512, instead. Thus, g 0 m f m T m ðHÞ: This can indeed be evaluated with a calculation cost of O(N). Detailed information for the O(N) approximation is included in Supplementary Methods.