Anderson-Kitaev spin liquid

The bond-disordered Kitaev model attracts much attention due to the experimental relevance in $\alpha$-RuCl$_3$ and $A_3$LiIr$_2$O$_6$ ($A=$ H, D, Ag, $\textit{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, and Majorana fermions will be in the Anderson localization state immediately. We call this immobile liquid state Anderson-Kitaev liquid state with two localized Majorana fermions, one frozen by gauge fluctuations and the other localized by disordered mass terms. The quantization of the thermal Hall conductivity $\kappa/T$ disappears by a quantum Hall transition induced by a small disorder, and $\kappa/T$ shows a rapid crossover into the Anderson-Kitaev liquid with a negligible Hall current. Especially, the critical disorder strength $\delta J_{c1} \sim 0.05$ in the unit of the Kitaev interaction would have many implications for the stability of Kitaev spin liquids.

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, and Majorana fermions will be in the Anderson localization state immediately. We call this immobile liquid state Anderson-Kitaev liquid state with two localized Majorana fermions, one frozen by gauge fluctuations and the other localized by disordered mass terms. The quantization of the thermal Hall conductivity κ/T disappears by a quantum Hall transition induced by a small disorder, and κ/T shows a rapid crossover into the Anderson-Kitaev liquid with a negligible Hall current. 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 two-dimensional (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], but it is also known that these honeycomb materials cannot fully be understood by the original (pure) Kitaev model [11,12]. While other diagonal or offdiagonal interactions might be important in real materials [13], the importance of disorder has been ignored in these materials until recently [14][15][16]. 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,17], which strongly suggests the existence of disorder [18,19]. The candidate ground state 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 longstanding problem because of the absence of a solvable model, except for limited cases [20]. 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,17,21], 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. We invented a powerful numerical method based on kernel polynomial method (KPM) [22] to do a large-scale investigation (O(10000) sites) for QSL.
Specifically, a Kitaev spin liquid (KSL) [1] is characterized by the fractionalization of the spin into two types of Majorana fermions. As such, there is a possibility that an itinerant part of Majorana fermions will be localized by the Anderson transition after introducing a quenched disorder. This effect is strongest in 2D, but may be observable even in three-dimensional (3D) generalizations [23,24] (mobility edge). These states with Majorana fermions in an Anderson (weak) localization is named Anderson-Kitaev (AK) spin liquid, or AK liquid in short. We try to investigate the crossover between KSL and AK liquid by the bond-disordered Kitaev model. The pure Kitaev model is described by the following Hamiltonian: 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 pair of π-flux vortices (vison) is 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].
In this Letter, we simulate the bond-disordered Kitaev model to see a crossover between KSL and AK liquid, especially from the topological transition in the thermal Hall effect [25][26][27]. We discovered that quantized thermal Hall effect is not as stable as expected, and Majorana fermions are very easily localized by disorder. Utilizing an approximation trick introduced by Kitaev, a largescale calculation up to O(10000) sites is possible. Important information for the Anderson transition like density of states (DOS) has been calculated.
Magnetic field effect. -The Kitaev model on the honeycomb lattice can be defined from Fig. 1(a). The bonds parallel to the red, green, and blue ones are x-, y-, and z-labeled bonds. We first consider the pure Kitaev model with a magnetic field as follows.
where h = (h x , h y , h z ) t is an applied magnetic field. We define a position operator r α for the n α -direction for α = 1, 2.
It is well-known that V 0 can be treated by the thirdorder perturbation [1]. The result after introducing itinerant Majorana fermions c j is jk c j c k + (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 [28]. The direction of the bond kl is defined clockwise as shown in Fig. 1(a) around the site j. A site connected by the γ-bond from j is called γ[j] for γ = x, y, and z, as shown in Fig. 1(a). We defineh = h x h y h z /48. Kitaev's trick. -Next, let's 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 (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 (E) numerically to exhaust every term appearing in the third order. For example, a Green function for excited states is efficiently obtained by the KPM [22] 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 (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 like the following: whereκ kl depends on the intermediate site j in the thirdorder perturbation process. From j,κ can be calculated by replacing 3/(α 0 J) 2 by 1/(∆ x ∆ y ) + 1/(∆ y ∆ z ) + 1/(∆ z ∆ x ), where ∆ γ is a vison gap for 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 irrelevant, so 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 [29].
Thermal conductivity. -We only consider zero temperature and ignore thermal flux fluctuations above the 0-flux sector. Lieb's theorem [30] 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. 1 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 in-plain thermal Hall conductivity κ xy (T ), especially a 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. A Kubo formula for κ xy at zero temperature is reduced to the generalized Thouless-Kohmoto-Nightingale-den Nijs (TKNN) formula [31] for noninteracting Majorana Hamiltonians [32]: where m and n label eigenvalues of H, ε m and ε n , corresponding to eigenstates |m and |n , respectively [33,34]. ϑ(x) is a Heaviside theta and v α = i[H, r α ]/ is a velocity operator along the α-direction. This Kubo-TKNN formula [35] is nothing but a real-space formulation of the Chern number calculation. We can alternatively use the so-called noncommutative Chern number (NCCN) [36], 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 = n ϑ(−ε n ) |n n| is a spectral projector. These two formulae must agree in the thermodynamic limit by a well-known relation κ xy /T = πk 2 B Ch/(12 ) for Majoranas. The finite-size effect is suppressed exponentially by an artificial k-space quantization of a size L × L and by replacing the commutator [36] as follows: where ∆ = 2π/L, c 0 = 0 and c q = −c −q are determined to hold x − L/2 q=−L/2 c q e iq∆x = O(∆ L ), and Q ≤ L/2. When Q = L/2, this formula exponentially converges to the thermodynamic limit with a self-converging property.
Thus, we can expect that these two methods may agree with a large L, while the Hall conductivity and the Chern number are a priori different quantities. We note that there are other ways to detect the topological nontriviality [37][38][39].
After taking an average of κ 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 = 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 [40], 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. 2(a)]. We note that Haar-random vectors used in this calculation show large errorbars and are not used in the following as described in Supplemental Material (SM) [41]. From here we will prefer the NCCN because we can use the KPM to approximate the spectral projector P F to avoid the diagonalization [42]. We fixed Q = 15 for L > 30 because otherwise the calculation cost becomes O(N 3 ). KPM can reproduce the vison gap approximately and at most reduce the computational cost to O(N ) with a truncation [43][44][45]. However, later we found that the truncation cause a problem in our simulation, and thus we used the O(N 2 ) algorithm [22,46,47].
Next, we would move on to a large-scale calculation by Kitaev's trick. From now on, κ xy is always calculated through NCCN. 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 where ∆ min is the minimum vison gap as a vison gap has spatial dependence on each bond, for simplicity [48]. In order to reduce the finite-size effect, we adopt Kitaev's torus basis where the finite-size effects cancel out, which is defined from a torus basis (L n 1 , L n 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 SM [41]. 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 [49]. The mean and minimum value of vison gaps are plotted for each δJ in Fig. 2(b). When δJ > δJ c2 ∼ 0.3, the vison gap approaches 0 for some plaquette, and the 0-flux ground state is destabilized. From here, the perturbation from the 0-flux sector cannot be justified. Moreover, after the gap closing, some flux sectors get almost degenerate and the first-order 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. 2(c) 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. 2(d)  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. 3(a)] the LDOS becomes gapless as the disorder strength increases. In the gapless region, DOS behaves linearly around ε = 0 [see Fig. 3 The localization in Fig. 3(b)-(d) is clear from the discrepancy between the arithmetic and geometric averages of LDOS. Details are included in SM [41].
Discussions. -Though we only did a finite-size calculation, the transition between KSL and AK liquid was well-observed and the schematic phase diagram in Fig. 1(c) was confirmed. From the extrapolation, δJ c1 is very small and δJ c1 /J ∼ 0.05. This fragility may be related to the long-range correlation in the mass term disorder [50], and reflects the nonlocality of the definition of Majorana fermions. We note that the vortex disorder is known to be relevant, so the introduction of random vortices may change the universality [51]. 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].
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 [52]. 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 [53]. The sensitivity also resembles unconventional superconductors [54]. It might be universal in strongly correlated systems. Thus, we need to reconsider the importance of cleanness for the topological order in general. Last but not least, we fixed h/∆ min = 1.0 for simplicity, so it is necessary to check another parameter region for comparison.
We thank Y. Akagi  We would like to introduce the kernel polynomial method (KPM) [1][2][3]. For this approximation, We wrote a program in the Julia 1.3.1 language.
First, let's 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 skewsymmetric 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. (2) 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. From the connectivity of the Hamiltonian we can set s = 6(J + δJ) without a magnetic field, but it is more convenient to use Arpack.jl or ArnoldiMethod.jl to compute the maximum absolute value E max of eigenvalues, and set s = E max + 0.1.

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 ∆, 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) [1]. Fortunately, Julia has a wrapper for FFTW [4] called FFTW.jl. Using this, the integral is reduced to a discrete weighted sum of the Fermi(-Majorana) function evaluated at some specific points. FFT (type-III DCT for the real diagonal part) decreases the computational cost drastically for evaluating the Chebyshev polynomials. The discretization sizeM was setM = 2M for simplicity. We tried M = 64, 128, 256, 512, 1024, and 2048. We found that M = 1024 has the best performance for our purpose, where the error is always about 0.01J.
As for the estimation of a noncommutative Chern number (NCCN) [5], it is better to expand the Fermi function directly instead of an FFT scheme. The spectral projector P F can be written as: where C is a contour which encloses every negative eigenvalue of H. The integrant is nothing but a Green function, so this 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 is a cutoff of the expansion for NCCN [6]. We used the Jackson kernel for the Chern number: where m = 0, . . . , M − 1. We here used M = 512, instead. Thus, We will again use M = 1024 for local density of states (LDOS) later.

O(N ) APPROXIMATION OF THE CHERN NUMBER
If the method above is directly applied, the estimation of P F still requires O(N 2 ). We propose a stochastic method to evaluate the Chern number based on the randomized algorithm to estimate a trace [1]. By picking up R normalized random vector |r , the trace is approximated as follows.
where B is some N × N matrix. We can use a Haarrandom vector, a Z 2 -random vector, etc. In Fig. S1, R = 100 Haar-random vectors are used. Fig. S1 shows the convergence of the Kubo formula and NCCN better than in the main text. An unphysically large magnetic field of h/∆ min = 3 √ 2 is applied. However, it is actually better to use a site basis |i for a random site i to take a trace. Actually, in the thermodynamic limit NCCN is expected to be independent of i [5]. R = 24 is enough for our purpose, and the result gets accurate as the system size increases.
In this O(N ) method, q = 1, . . . , Q terms are all expanded. In the large scale we fixed Q = 15, so we need to evaluate Q 2 = 225 terms at the same time. We note that terms with a negative q can easily be computed from the positive ones.
Instead of calculating every element of P F , we can now estimate P F |α for any (sparse or dense) vector |α with an O(N M ) cost. Since r 1 and r 2 are both diagonal in the original basis, the trace can be computed by sequential estimations of P F |α after decomposing the summation over q. Eventually, the calculation costs becomes O(N M QR + N Q 2 R), and it scales as O(N ) as long as M, Q, and R are kept constant. We note that the preparation of the look-up tables has been ignored.

RANDOM-WALK O(N ) APPROXIMATION FOR THE CHEBYSHEV SERIES
Though in the main text we only used a conventional O(N 2 ) method for calculating the Chebyshev series, there exists a stochastic algorithm which can approximate the series with a computational cost of O(N ) without any truncation. This algorithm can work as a Markov chain simulation of the matrix product process and the probability that random walkers go back to the original site gives each (diagonal) element of the moment of the Hamiltonian. Thus, we can use this weighted random walk to obtain every (diagonal or offdiagonal) element of the moments, and we will quickly demonstrate the algorithm for diagonal cases, and show that all the diagonal parts can be obtained in O(N ) time. This is nothing but an importance sampling among every con-tour which goes from the site i to the same site i, and the weight stays constant, i.e. the importance sampling goes well, when the disorder strength is small enough. Examples are shown in Fig. S2(a).
The algorithm is very simple. The number of walkers N walker should be kept constant. for j ∈ {1, . . . , N walker } do 11: ν(pj) ← nearest neighbor sites of pj 12: ω(pj) ← hopping amplitudes from pj 13: wj ← wj×(the sum of ω(pj))/a 14: sample pj from ν(pj) with weights ω(pj) 15: if xj = i then 16: hm ← hm + wj In the weak-disorder regime, the growth of weight can be kept almost unity, so the histogram itself can be approximated by a Poisson distribution. Thus, we put an errorbar for the number of eventsÑ as Ñ . The calculation is compared with the exact results in Fig. S2(b).
One of the advantages of this method is that by replacing H by H + sI, where I is an identity matrix, the transition matrix can always be made positive-definite. The convergence is guaranteed by the Perron-Frobenius theorem, and the relaxation time is determined from the difference between the first and second largest eigenvalues. Thus, after the relaxation, all the moments can be replaced by the equilibrium value and we can reconstruct the Chebyshev series until M → ∞.
However, the prefactor of the cost of this algorithm is very large, so we rather prefer the O(N 2 ) deterministic algorithm, instead, in this work. This O(N ) algorithm is advantageous only in a very large scale N ∼ 10 6 and we would pursue its usefulness in the future. Especially, in the strongly disordered case, the weight w j should be balanced by a node-weighted algorithm.
In the flux sector except for the 0-flux one, this algorithm strongly suffers from a "sign problem". Contours enclosing a π flux will be affected by this sign problem.
Dealing with a general flux based on this algorithm is also an important future problem. We note that this method can be combined with a truncation by excluding walkers which go too far from the origin.

KITAEV'S PERIODIC BOUNDARY CONDITION
In our method, the Chern number calculation assumes infinite supercells of the periodic boundary condition, and also the unitarity of the discrete Fourier transformation. Thus, if we employ Kitaev's periodic boundary condition (KPBC) [7], we must do a linear transformation to retain the periodicity of plane waves. The coordinate transformation is necessary from the usual crystallographic coordinate (r 1 , r 2 ) to (r 1 − r 2 /L, r 2 ). In this new basis, the computation of NCCN is possible. Apparently, this basis change does not affect the thermodynamic limit.

EVALUATION OF LOCAL DENSITY OF STATES
LDOS at site i is defined as follows [1].
where |k is an eigenstate of H corresponding to an eigenvalue E k . Apparently, the average over every i would be the density of states (DOS). As long as the system is uniform, LDOS coincides with DOS. LDOS can be calculated through KPM. The expansion is like We note M = 1024 is used for LDOS. Usually, i is randomly chosen in the same way as the Chern number calculation, and the arithmetic mean of ρ i (E) is denoted by ρ ave . Its geometric mean ρ typ over a number of samples is also important.
where the disorder average is taken in the expression. The reason why the difference in the arithmetic and geometric averages of LDOS captures the localization is as follows. If the localization happens, the LDOS strongly depends on the site, and LDOS at site i deviates from the true DOS (averaged LDOS). Thus, the typical LDOS (ρ typ ) deviates from the averaged LDOS (ρ ave ), and it signals the localization transition. Though we have plotted LDOS for one-body Majorana Hamiltonians, which is not necessarily the same as the many-body (L)DOS, the results can always be transformed into the many-body language, which will be done in the future work.
As for another signal, the localization length is also calculable using the real-space Kubo-Greenwood formalism [8], but we would leave it as future work because LDOS already worked as the proof of the localization.

FUTURE DIRECTION
In reality, a problem of solving the disordered Kitaev model is many-body localization (MBL) [9] because there is a four-fermion term [7]. The effect of a finite interaction should be studied in the future to find out some MBL criterion. The problem itself resembles physics of the Sachdev-Ye-Kitaev model and its MBL [10,11]. Especially, the localization of the protected edge states is not discussed in this Letter.
In our calculation, there exists an intermediate gapless phase between the non-Abelian topologically ordered phase and the trivial phase. This corresponds to the observed crossover between the Kitaev spin liquid and the Anderson-Kitaev liquid, and the thermal Hall conductivity takes a nonquantized value in the crossover regime. In order to characterize this phase, the longitudinal component κ xx is also important [12], which would be an important future study. We note that the use of a binary disorder makes the transition steep, so the calculation by Anderson disorder is also necessary.
As for the computational complexity, a complete O(N ) method is possible without a truncation by a probabilistic approximation of the moments of the Hamiltonian using multiple random walks. This new approach was tested and we confirmed its accuracy. Applying this stochastic O(N ) method to our problem setting is also interesting future work, where essentially O(10 5-6 ) sites are achievable. However, unfortunately, this approach in the original form only works for the 0-flux sector and suffers from the "sign problem" in other flux sectors, although the problem at hand is completely classical. This means that a full O(N ) method is not universal and solving this classical sign problem is also a future problem.
In the 0-flux calculation the thermal fluctuation to destroy the quantization is underestimated, so we need an unbiased Monte Carlo simulation [13,14] to go to finite temperature correctly. From the previous study, the thermal Hall conductivity is quantized only in the regime T /J < O(0.01) [14], and this suggests that the quantization is very weak under bond disorder as well as with thermal fluctuation. The combined effect is not sought in this Letter, and thus Monte Carlo simulations are worth doing. However, in three dimensions non-Abelian topological phase could potentially survive under both bond and flux disorders because the flux fluctuation is frozen at finite T c by the second-order transition [13], although we do not know a good lattice where the thermal conductivity can be quantized in three dimensions.