Identification of magnetic interactions and high-field quantum spin liquid in α-RuCl3

The frustrated magnet α-RuCl3 constitutes a fascinating quantum material platform that harbors the intriguing Kitaev physics. However, a consensus on its intricate spin interactions and field-induced quantum phases has not been reached yet. Here we exploit multiple state-of-the-art many-body methods and determine the microscopic spin model that quantitatively explains major observations in α-RuCl3, including the zigzag order, double-peak specific heat, magnetic anisotropy, and the characteristic M-star dynamical spin structure, etc. According to our model simulations, the in-plane field drives the system into the polarized phase at about 7 T and a thermal fractionalization occurs at finite temperature, reconciling observations in different experiments. Under out-of-plane fields, the zigzag order is suppressed at 35 T, above which, and below a polarization field of 100 T level, there emerges a field-induced quantum spin liquid. The fractional entropy and algebraic low-temperature specific heat unveil the nature of a gapless spin liquid, which can be explored in high-field measurements on α-RuCl3.


Supplementary Note 1. Determination of the Effective α-RuCl 3 Hamiltonian
Fitting the model parameters from thermodynamic measurements. In this section we show the workflow of the model parameter fittings. By performing exponential tensor renormalization (XTRG) calculations on a YC4×4×2 lattice, we scan the parameter space spanned by the couplings [K, Γ, Γ , J], and fit the specific heat C m (T ), in-plane (χ ab ) and out-of-plane susceptibilities (χ c * ) measured under a small magnetic field µ 0 H 1 T. From the susceptibility simulations, we also determine the corresponding g-factors, g ab and g c * , along with other couplings. Given the determined parameters, we can compute the static spin structure factors, and confirm the appearance of zigzag magnetic order at low temperature (see   [1][2][3], which are shown in the background. The energy scale |K| is tuned adaptively in different cases to fit the C m curves. To be concrete, we show in Supplementary Figs. 1, 2 part of our simulated data in the thermodynamic properties. In Supplementary Fig. 1(a,b), we start with scanning over various Γ values with ferromagnetic K < 0, by setting Γ = J = 0 at first. As shown in Supplementary Fig. 1(a), the C m curves are sensitive to Γ for Γ < 1, while the curves do not change much for Γ ≥ 1 in Supplementary Fig. 1(b), given that the energy scale |K| is properly tuned. Therefore, to uniquely pinpoint the parameter Γ/|K|, we need to include more thermodynamic measurements like the magnetic susceptibility.
As shown in Supplementary Fig. 2(a,b), as the Γ interaction increases, the height of computed susceptibility peak decreases accordingly and deviates the experimental in-plane susceptibility χ ab curves. After a thorough scanning, we find Γ = 0.3 constitutes an overall optimal parameter in our model fittings of the specific heat and susceptibility. With the given Γ, we find in Supplementary Fig. 2(c) the susceptibility curves turn out to be not very sensitive to the small Γ interactions, while, on the other hand, the low-T peak of C m moves towards higher temperatures as |Γ | increases in Supplementary Fig. 1(c). This can be understood as the Γ term is crucial for stabilizing the zigzag order when K < 0 and thus has strong influences on the low-T peak. Therefore, we fix Γ = −0.02 from the specific heat fittings.
Lastly, we determine the nearest-neighboring Heisenberg term J. In Supplementary Fig. 2(d), we find the height of the susceptibility χ ab peak enhances and thus approaches the experimental curves as J changes from −0.06 to −0.1, while the specific heat C m is not sensitive to J, as shown in Supplemen- neighboring Heisenberg interaction J 3 has also been suggested to stabilize the zigzag order for K < 0 [4][5][6], which can play a similar role to the off-diagonal Γ interaction. In order to explore its effects, we introduce an additional J 3 term to our α-RuCl 3 model and show the computed results in Supplementary   Fig. 3.
In the spin structure factors in Supplementary Fig. 3(a,b), we find the M-point intensity clearly enhances when J 3 increases from 0.01 to 0.05. At the same time, even a small J 3 can have a considerable impact on the low-temperature thermodynamics and makes our model fittings deviate from the experimental data. In Supplementary Fig. 3(c,d), we find a J 3 ≥ 0.05 clearly spoils the fittings to C m and χ: The low-T specific heat peak moves towards higher temperature, and the height of the susceptibility peak gets reduced when increasing J 3 . In addition, from Supplementary Fig. 3(b), J 3 = 0.05 interaction evidently cripples the spin intensity at the Γ point of the Brillouin zone (BZ). Overall, we find J 3 indeed plays a very similar role as Γ in stabilizing the zigzag order, and, such a term, if exists, should have a small coupling strength. We therefore leave out J 3 in the main text as well as the discussion below, and mostly focus on the minimal Exact diagonalization results of dynamical spin structure factors. In the dynamical simulations, the 24-site cluster with C 3 symmetry (denoted as 24a henceforth) has been widely adopted in ED calculations of the Kitaev model (see, e.g., Refs. [5,[7][8][9]). Amongst other geometries that are accessible by ED, the 24a cluster is unique as it contains all the high-symmetry points in the BZ [c.f., inset in Supplementary   Fig. 4(c)], which is important for dynamical property calculations. Therefore, we choose the 24a cluster in presenting our dynamical data in the main text. Nevertheless, in Supplementary Fig. 4 we also perform dynamical ED simulations on a different 24-site geometry (24b) and compare the results to those of the 24a cluster. It can be seen in Supplementary Fig. 4 that the positions ω M and ω Γ of the intensity peaks as well as the double-peak structure of Γ-intensity curves are qualitatively consistent. Note the 24b cluster shown in Supplementary Fig. 4(b,d) are not C 3 symmetric and thus does not possess all high symmetry moment points.
Magnetic anisotropy and the off-diagonal Γ term. The off-diagonal Γ interaction has been suggested to be responsible for the strong magnetic anisotropy between in-and out-of-plane directions [10,11]. In a pure Kitaev model without the Γ term, it has been shown that the in-plane (χ ab ) and out-of-plane (χ c * ) susceptibilities are of similar magnitudes [12]. On the other hand, as shown in Supplementary Fig. 2, when the Γ term is introduced, we find the two susceptibility curves are clearly separated, suggesting that the offdiagonal Γ interaction indeed constitutes a resource of strong easy-plane anisotropy observed in α-RuCl 3 .
Automatic searching of the Hamiltonian parameters. Through the fitting process described above, we have manually determined an accurate set of model parameters of α-RuCl 3 . As a double-check, below we exploit the recently developed automatic parameter searching technique based on the Bayesian approach [13] to conduct a global optimization of the model parameters in Supplementary Fig. 5. The target is to minimize the fitting loss of the simulated results to measured thermodynamic quantities, including the specific heat C m as well as the in-(χ ab ) and out-of-plane (χ c * ) susceptibilities. Since the specific heat data vary greatly in various measurements and the lower-T peak diverges likely due to 3D effects, we take the mean values of several measurements as the specific heat data C exp m (T ) used in the fittings. In practice, the loss function L is designed as follows: where λ C = 1/max(C exp m ) 2 and λ χ = 1/max(χ exp ab ) 2 controls the weights of specific heat and magnetic susceptibility in defining the fitting loss function. N χ ab , N χ c * , and N Cm are the data point numbers in three experimental curves, respectively. To focus on the fittings of the locations of two specific heat peaks at 100 K and 8 K, respectively, we further introduce a penalty factor P defined as K} are the position of the peaks, with {C 1 , C 2 } the specific heat value at corresponding higher and lower T . As our many-body calculations are performed on a finite-size system, we need to introduce the T cut (T cut ) as the lowest temperature involved in the fittings, and in practice it is set as T cut = 25 K for magnetic susceptibility and T cut = 3 K for the specific heat data. The factor P emphasizes the double-peak structure as well as the locations and heights of each peaks in the simulated C m , where we set an empirical factor a = 2 to enlarge the loss function as a penalty for the simulated specific heat curves without double-peak structure. Figure 5. The estimated fitting loss L shown in (a) K-Γ, (b) K-Γ , and (c) K-J planes. The results are obtained by Bayesian optimization through 120 iterations of XTRG calculations on the YC4×4×2 cylinder. The asterisk represents the optimal parameter set found to also reproduce other major experimental observations including the magnetization curve and dynamical spin structures, which is located within the bright regime with small fitting loss of thermodynamics. The two Landé factors g ab and g c * at each parameter point are tuned to their optimal values (between 1.5 and 3) to minimize the loss L. Note, when plotting the cross sectors, like the K-Γ plane in (a), the other model parameters are fixed at their optimal values.
In Supplementary Fig. 5, we show the estimated loss L in the K-Γ, K-Γ , and K-J planes, respectively.
From the color map of L, we see clearly optimal (bright) regimes, and the determined optimal parameter set in the main text is indicated by the asterisk, i.e., K = −25 meV, Γ = 0.3 |K|, Γ = −0.02 |K|, and J = −0.1 |K|, with the in-and out-of-plane Landé factors found to be g ab = 2.5 and g c * = 2.3, respectively, which can fit excellently both the thermodynamic and dynamic measurements in Fig. 2  Refs.
‡ Check if the model exhibits a low-T zigzag order. * Check if the dynamical M-star structure exists (taken from Ref. [9], except for our model in the last row).
Our α-RuCl 3 model proposed in the main text.
Supplementary Table 1. Checklist of various candidate models of α-RuCl 3 on typical experimental features.
The "! " symbol indicates a good fitting between model calculations and experimental measurements, while "% " represents a disagreement between the two.
Specific heat and susceptibility curves. In Supplementary Fig. 6 we show the simulated magnetic specific heat C m and susceptibility χ of various candidate models, and compare them to the experimental measurements. From Supplementary Fig. 6(a-c), we find only the models Suzuki2019 [18] and Ran2017 [19] exhibit a double-peaked C m curve, each located at the characteristic temperature precisely as in experiments. For the rest of revisited candidate models, however, we did not observe the desired double-peak feature in the right temperature window.
Overall, as concluded in Ref. [9], we also find no single candidate model revisited here can satisfactorily explain all observed phenomena of experiments on α-RuCl 3 . For example, both Cookmeyer2018 [16] and Ran2017 [19] have M-star shape in the intermediate-energy dynamical spin structure, however, the former does not fit the specific heat and susceptibility measurements well and the latter does not host a zigzag oder at low temperature. Nevertheless, we note that the model Suzuki2019 [18] can reproduce the prominent thermodynamic features such as the double-peaked specific heat, highly anisotropic susceptibility, and zigzag static spin structures (although the dynamical M star was not found in this candidate model according to Ref. [9]). Our model, on the other hand, accurately describes the spin interactions and explain thus major experimental findings in the compound α-RuCl 3 .

Supplementary Note 3. α-RuCl 3 Model Simulations under In-plane and Tilted Magnetic Fields
In this Note, we show various thermodynamic properties of α-RuCl 3 under in-plane and tilted c fields (see Fig. 1b in the main text), and compare the simulated results to experimental measurements.
Magnetization curves along the a-and c -axis. In Supplementary Fig. 8(a,b)  derivative dM/dH, where the zigzag order becomes suppressed and the system enters the polarized phase.
As the field further increases, the uniform magnetization M gradually approaches the saturation value.
Specific heat under magetic fields and the Zeeman energy scale T l . In Supplementary Fig. 8(c,d), we show the low-T part of the magnetic specific heat curves on YC4×6×2 systems calculated by XTRG.
For µ 0 H [112] < 7 T, we find the height of the C m peaks is suppressed by fields, and the low-temperature scale T l associated with the zigzag order, as indicated by the red arrow in Supplementary Fig. 8(c), decreases to zero when the field approaches the critical value, in agreement with experiments [1]. Above that, a new low-temperature scale T l builds up and moves towards higher temperatures almost linearly vs H [112] , as indicated by the black arrows in Supplementary Fig. 8(c). We find T l is intimately related to the uniform magnetization, as observed from the temperature variation of spin structure intensity S(Γ). In Supplementary Fig. 8(e), we show the derivative dS(Γ)/dT exhibits clear peaks, whose locations are in agreement with the temperature scale T l determined from the low-temperature peak of C m in Supplementary Fig. 8(c). The peaks moving towards higher temperature as the magnetic field increases, and we thus relate T l to the Zeeman energy scale. Beside in-plane fields, as depicted in Supplementary Fig. 8(d,f) Energy spectra under in-plane fields. In Supplementary Fig. 9 we show ED results of the energy spectra on the 24-site cluster (24a), under two in-plane fields H [112] a and H [110] b axis [c.f., Supplementary   Fig. 9(a)]. The energy spectra results under various fields are plotted in Supplementary Fig. 9(b)  , and the spectra appear differently along a and b directions, consistent with the observation in recent experiments [10,25]. 8 T, respectively, as indicated by the thick vertical lines.

Supplementary Note 4. Quantum Spin Liquid under Out-of-plane Magnetic Fields
Finite-temperature Spin structures under out-of-plane fields. In Supplementary Fig. 10, we plot the spin structure factors at two typical momentum points in the BZ, i.e., k =M and Γ. From the S(M) curves in Supplementary Fig. 10(a), we find the zigzag order [cf. 13 T and 26 T lines in Supplementary   Fig. 10(a)] gets suppressed in the quantum spin liquid (QSL) phase under strong out-of-plane magnetic fields [see, e.g., 45.5 T, 78 T, and 91 T lines in Supplementary Fig. 10(a)]. In the QSL phase, there exists enhanced M-point spin correlation near the emergent lower temperature scale T l , as indicated by the black arrows in Supplementary Fig. 10(a). In Supplementary Fig. 10(b), we show the temperature derivative −dS(Γ)/dT , whose peak position signals the low-temperature scale T l , below which the field-induced uniform magnetization gets established, i.e., T l is a crossover temperature scale to the polarized state. As shown in Supplementary Fig. 10(b), T l moves to higher temperatures as the field increases.
In the main text, we have discussed the static spin-structure factorsS zz (k) under an out-of-plane field of 45.5 T (cf. the insets of Fig. 5b of the main text). In Supplementary Fig. 10(c,d), we supplement with 10 100  Kitaev interaction (∼ −25 meV), which becomes more and more clear as the system cools down near the low-T scale T l , around which the brightness at Γ and M points are slightly enhanced due to the strong spin fluctuations. When T < T l the Γ peak in the spin structures becomes flattened, yet the stripy background remains distinct, which is in excellent agreement with the DMRG results, supporting strongly the existence of QSL state under high out-of-plane fields.
Ground-state static spin structures under a field of 78 T. In Supplementary Fig. 11 and Supplementary Fig. 12, we show the ground-state static structure factorS αβ (k) = j∈bulk, j =i 0 e ik·(r j −r i 0 ) ( S α i 0 S β j − S α i 0 S β j ) of our α-RuCl 3 model on cylinders of different sizes in the high-field QSL phase (under a magnetic fields of 78 T), obtained by DMRG simulations with bond dimension up to D = 2048. Here we fix a central site i 0 and compute the correlations with respect to this reference site, instead of computing the all-to-all correlations, to speed up the computations. From Supplementary Fig. 11, we find strong stripy background inS xx ,S yy , andS zz , divided by N bulk , where N bulk is the number of bulk sites (i.e., with the left-and right-most columns of the cylinder skipped), representing intrinsic Kitaev spin liquid characteristics. There also exists a bright Γ point in the BZ, and we depict the profiles with fixed k y = 0 in Supplementary Fig. 12, to check whether the structure factor diverges at any k point. In Supplementary Fig. 12 we find both the diagonal partS diag (k) = α,βS αβ (k)δ αβ and the off-diagonal part S v . Considering a bipartition of the system into two parts A and B, the reduced density matrix of A can be obtained by tracing out the degrees of freedom of B, i.e., ρ A = Tr B ρ, where ρ = |ψ 0 ψ 0 | is density matrix of the normalized ground state |ψ 0 of the whole system. The von Neumann entropy is defined as Here, we consider the cut along circumference direction and measure S v for each cut at l. By virtue of the conformal mapping: l →l = (L/π) sin(πl/L), a central charge c on the open cylindrical geometry can be formally extracted using originally proposed in the (1+1)-D conformal field theory. The prominent "dome"-like S v curves in Supplementary Fig. 13(a) together with the linear fitting with effective central charge c in Supplementary   Fig. 13(b), indicates a gapless intermediate phase. This is consistent with the gapless feature observed in the energy spectra as shown in Fig. 5d and algebraic magnetic specific heat in Fig. 6d of the main text.
Furthermore, the effective central charge c increases with the system width W , i.e. the allowedŷ momenta in the Brillouin zone, suggesting the increased gapless modes.

Supplementary Note 5. Variational Monte Carlo Results
Spin-liquid ansatz based on D 3d × Z T 2 symmetry. We start with the extended Kitaev honeycomb model containing K, J, Γ and Γ interactions. The most general expression of the mean-field Hamiltonian ansatz [26][27][28] with nearest-neighbor couplings reads Owing to the subtle SU (2) gauge symmetry [29] in the fermionic spinon representation, the particle number constraintN i = 1 is extended into the SU (2) gauge invariant three-component form Tr(ψ ψ ψ i τ τ τ ψ † i ) = 0 which are ensured by the three Lagrangian multipliers λ x,y,z , where τ x,y,z are the generators of the SU (2) gauge group. The matrices U (0,1,2,3) ji can be expanded with the identity matrix and τ x,y,z , where the expanding coefficients form a subset of R R R. For the present model, the SU (2) gauge symmetry breaks down to its subgroup Z 2 (which is called the invariant gauge group, IGG) and the resultant quantum spin liquids are Z 2 QSLs.
A QSL ground state preserves the whole space group symmetry whose point group is D 3d × Z T 2 . However, the symmetry group of a spin liquid mean-field Hamiltonian is the projective symmetry group (PSG) [30,31] whose group elements are space group operations followed by SU (2) gauge transformations. Although there are more than 100 classes of PSGs for Z 2 QSLs, here we adopt the PSG that describes the symmetry of the mean-field Hamiltonian of the pure Kitaev model (called Kitaev PSG). Under this constraint, the allowed mean-field Hamiltonian contains the following parameters.
Firstly, recalling that in KSL the Kitaev interactions are decoupled in a way that the c Majorana fermions are not mixed with the b γ Majorana fermions, we have Similarly, the Γ interactions are decoupled as and the Γ interactions are decoupled as ji , in which j and i specify γ.
Secondly, we consider a more general mean-field decoupling which means c Majorana fermions can mix with b γ Majorana fermions. In addition, the most general coefficients preserving the C 3 rotation symmetry (in the PSG sense) also contain multiples of the uniform (I) and τ x + τ y + τ z gauge components,Ũ ji = η 6 +iη 7 (τ x +τ y +τ z ). If the full symmetry group, G = D 3d ×Z T 2 , is preserved, then only three parameters η 0 , η 3 , and η 5 are allowed. Thus a spin-liquid ansatz that preserves the full Kitaev PSG contains the variables with seven real parameters, ρ a , ρ c , ρ d , ρ f , η 0 , η 3 and η 5 .
Magnetically ordered states. Because fermions do not condense alone to form a magnetic order, we introduce the classical order under single-Q Q Q approximation [32] to describe the magnetic order of the spinsymmetry-breaking phases of the K-J-Γ-Γ model. where Q Q Q is the ordering momentum,ê e e x,y,z are the local spin axes (not to be confused with the global spin axes), and φ is the canting angle. π/2 − φ describes the angle by which the spins deviate from the plane spanned byê e e x andê e e y . The classical ground state is obtained by minimizing the energy of the trial states.
In our VMC calculations, the static order is treated as a background field coupling to the spins as sitedependent Zeeman field, hence the complete mean-field Hamiltonian for the K-J-Γ-Γ model reads Finally, since the magnetic order essentially breaks the C 3 symmetry, the parameter η 7 and an additional parameter θ describing the ratio of certain parameters in the z-bond and x-(y-) bonds are allowed (similar to the model discussed in Ref. [33]).
Field-induced chiral spin liquid. By comparing various ansatzes, including the spin-liquid and magnetically ordered states, our VMC calculations reveals an exotic phase under the out-of-plane magnetic field B B B ≡ µ 0 H [111] . As the zigzag order being suppressed by a large field, an intermediate chiral spin liquid (CSL) phase with Chern number ν = 2 emerges before the system enters the polarized phase. The phase transition between the zigzag ordered state and the the CSL phase is of first order. According to Kitaev's seminal work [34], the CSL is an Abelian phase whose quasiparticle excitations include a,ā, ε, 1, where 1 denotes the vacuum, ε is the fermion, a andā are two types of vortices with topological spin e iπ/4 . The fusion rules are ε × ε = 1, a ×ā = 1, a × ε =ā,ā × ε = a,ā ×ā = ε. The edge of the CSL state is gapless and contains 2 branches of chiral Majorana excitations, each branch carries a chiral central charge of 1/2.
Therefore, the thermal Hall conductance is integer quantized, with κ xy /T = πk 2 B /6h at low temperature.
In the following, we discuss whether the field-induced Z 2 CSL is topologically nontrivial or not. The confinement or deconfinement of the Z 2 gauge field is reflected in the ground state degeneracy (GSD) of the Gutzwiller projected state when placed on a torus. Therefore, we calculate the density matrix of the projected states from the wave-function overlap ρ αβ = P G ψ α |P G ψ β = ρ * βα , where α, β ∈ {++, +−, −+, −−} are the boundary conditions (+ stands for the periodic boundary condition and − for antiperiodic boundary condition) along a a a 1 and a a a 2 direction, respectively. If ρ has only one significant eigenvalue, with the others vanishingly small, then the GSD is 1, indicating that the Z 2 gauge field is confined. Otherwise, for ρ with more than one (near-degenerate) nonzero eigenvalues, then the GSD is nontrivial and hence the Z 2 gauge fluctuations are deconfined. In the VMC calculations, we find the field-induced CSL (with Chern number ν = 2) is deconfined, with GSD equals 4 (the eigenvalues of the overlap matrix ρ are given by 0.5915, 0.8066, 1.0419, 1.5600 under g c * µ B B/|K| = 0.52 for a system on an 8×8×2 torus). Actually, we have performed a finite-size scaling calculation (not shown), which indeed indicates that the GSD of this state is 4 in the large-size limit. Indeed, such a GSD matches the number of topologically distinct quasiparticle types.