Dynamical and thermal magnetic properties of the Kitaev spin liquid candidate $\alpha$-RuCl$_3$

What is the correct low-energy spin Hamiltonian description of $\alpha$-RuCl$_3$? The material is a promising Kitaev spin liquid candidate, but is also known to order magnetically, the description of which necessitates additional interaction terms. The nature of these interactions, their magnitudes and even signs, remain an open question. In this work we systematically investigate dynamical and thermodynamic magnetic properties of proposed effective Hamiltonians. We calculate zero-temperature inelastic neutron scattering (INS) intensities using exact diagonalization, and magnetic specific heat using a thermal pure quantum states method. We find that no single current model satisfactorily explains all observed phenomena of $\alpha$-RuCl$_3$. In particular, we find that Hamiltonians derived from first principles can capture the experimentally observed high-temperature peak in the magnetic specific heat, while overestimating the magnon energy at the zone center. In contrast, other models reproduce important features of the INS data, but do not adequately describe the magnetic specific heat. To address this discrepancy we propose a modified ab initio model that is consistent with both magnetic specific heat and low-energy features of INS data.


INTRODUCTION
Quantum spin liquids (QSL) are long-sought-after states of matter without magnetic order, but with nontrivial topological and potentially exotic properties [1,2]. Much of the search has been focused on frustrated lattice systems [3,4], but in an important development in 2006 Kitaev [5] introduced a novel exactly solvable paradigmatic QSL with bond-directional Ising terms on the bipartite honeycomb lattice. Importantly, this Kitaev model hosts anyonic excitations [6], which are of interest both for fundamental reasons and for their proposed application in topological quantum computing [7,8]. It was realized that such interaction terms naturally appear [9] -and can be large -in Mott-insulating transition-metal systems with edge-sharing octahedra and strong spin-orbit coupling, such as in A 2 IrO 3 (A=Na,Li) [10,11], α-RuCl 3 [12] and other materials [13,14].
However, the three mentioned materials are all found to order magnetically at low temperatures, and hence cannot be perfect realizations of the Kitaev model. Na 2 IrO 3 [15][16][17] and α-RuCl 3 [18][19][20][21] both develop a zigzag order, as illustrated in Fig. 1 (a), while Li 2 IrO 3 displays an incommensurate spiral order [22]. Despite the zigzag order, α-RuCl 3 has emerged as a particularly promising Kitaev QSL candidate. Initial strong evidence came in the form of a strong and unusually stable scattering continuum at the zone center as observed in inelastic neutron scattering experiments (INS) [23][24][25][26], which has been interpreted as evidence for the presence of fractional Majorana excitations. However, in an alternative picture it has been proposed that the continuum may consist of incoherent excitations due to spontaneous magnon decays [27]. More recently, half-integer quantization of the thermal Hall conductivity was reported [28], also consistent with Majorana excitations. The quantization occurs in a narrow range of in-plane magnetic fields, where the magnetic order is melted, possibly uncovering an intermediate QSL state [24,[29][30][31][32]. Further Altogether, these experiments strongly suggest that α-RuCl 3 can be described by a generalized Kitaev-Heisenberg Hamiltonian [13,14], including off-diagonal and furtherrange interactions. Theoretical work leads to the same picture, whether the model is deduced from ab initio methods [12,[44][45][46][47][48][49][50], or from more phenomenological or ab initio inspired approaches [27,[51][52][53]. Unfortunately, these different works, as well as experimental fits [26,40], have led to a veritable zoo of proposed realistic spin Hamiltonian descriptions for α-RuCl 3 , and it's not currently clear which description is most accurate. Moreover, the proposed models disagree in terms of included spin-spin interaction terms, magnitudes of interaction parameters, and even signs. We illustrate this situation in Fig. 1 (c) by a scatter plot of the values for just two relevant interaction terms, and in Table I.
In this work, we adopt a systematic approach to address this uncertainty. We calculate static spin structure factors (SSFs) S (q) and INS intensities I (q, ω) for all models listed in Table I using Lanczos exact diagonalization (ED) [58] on 24-site clusters. We also use ED to explore the evolution of the INS spectra away from the ferromagnetic Kitaev limit as new perturbations are introduced in a generalized Kitaev-Heisenberg model. We then calculate the magnetic specific heat C mag for the models using the thermal pure quantum state (TPQ) method [59] on the 24-and 32-site clusters shown in Fig. 2. A few of the models considered here have previously been studied using similar methods in Refs. [27,53,60]. For clarity we will restrict the discussion in the main text to six particularly relevant models. These models all have a ferromagnetic Kitaev coupling (K 1 < 0), which is expected in α-RuCl 3 [13,14]. Results for the other models are included in the Supplementary Information.
Our key finding is that none of the studied models manages to fully capture the salient features of both the INS and magnetic specific heat data. The energy scales obtained in first principles approaches appear to be needed to reproduce a high-temperature peak in C mag , but the parameters proposed in the literature push the INS intensity at the Γ point to higher energies. On the other hand, models obtained by fits to INS data run the risk of missing significant off-diagonal interactions, and fail to reproduce the experimentally observed temperature dependence of C mag . By modifying one of the ab initio models, we are able to find results consistent with both C mag and low-energy features of the INS spectrum. Our results thus provide important clues for an accurate and realistic description of of α-RuCl 3 . Table I are special cases of a proposed minimal model [13,47] for α-RuCl 3 ,

Several of the Hamiltonians listed in
where . . . and . . . denote nearest and third-nearest neighbors, respectively. γ =X,Y,Z is the bond index shown in Fig. 1 (b), and α, β are the two other bonds. The Γ 1 term is required to explain the moment direction, and J 3 > 0 helps stabilize the zigzag order. Ab initio and DFT studies also tend to report a sizable symmetric off-diagonal Γ 1 interaction, which originates from trigonal distortion [61,62]. Since the crystal structure of α-RuCl 3 features trigonal compression [13,19,20,25] perturbative calculations [47,62] would predict that Γ 1 < 0, which provides an additional mechanism to stabilize the zigzag order [62,63]. In the absence of other crystal distortions, the most general nearest-neighbor Hamiltonian would thus be the Combining these two proposed minimal models results in the Further proposed extensions include second-nearestneighbor Kitaev and Heisenberg terms, third-nearest-neighbor Kitaev terms, and additional symmetry-allowed anisotropies [44,47]. In particular, α-RuCl 3 does not have a perfect honeycomb lattice, which allows the parameters for the Z bond to deviate from those on the X,Y bonds. In Table I we have bond averaged such anisotropies for the sake of clarity, but we will use the full parameter sets in our calculations when appropriate.
Spin structure factors and neutron scattering intensities Fig. 3 shows predicted zero-temperature neutron scattering intensity spectra, I (q, ω), for the six central models. All models feature sharp low-frequency peaks at the M points, indicating the zigzag order. The intensity at the M points is significantly higher than the intensity at the Γ point, which is inconsistent with experimental observations [24]. However, the M peaks could potentially be suppressed at finite temperatures [51]. The models in Fig. 3 (b), (c), (d) and (e) all show clear signs of the scattering continuum at the Γ point at frequencies comparable to the position of the M peak, whereas the two ab initio models in (a) and (f) display a sizable gap up to any noticeable scattering at the Γ point.
In the top row of Fig. 4 we have plotted the static spin structure factors S (q). As shown, all six models are consistent TABLE I. The spin Hamiltonians for α-RuCl 3 considered in this work. Dashes (-) indicate that the value is unavailable or negligible. The bolded models are considered in the main text, and results for the other models are given in the Supplementary Information. Asterisks in the 'BA' column signify that the full Hamiltonian has different values for the X/Y bonds compared with the Z bonds, and that the parameter values given in the row have been bond averaged.

Reference
Method  [57] Spin wave fit / THz spectroscopy +0.46 −3.50 +2.35 ----a Using the proposed minimal model, which is bond averaged and neglects small Γ 1 = −0.9 meV. Values for the monoclinic (C2/m) crystal structure. b We use the sign convention in Refs. [53,55]. c This work gives values for several values of U. Here we use the U = 3.5eV parameters. d Values for the C2 structure. e These are the parameters from the preprint version in Ref. [56]. They were revised in the published version, Ref. [44]. In Supplementary Note 4 we show that this slight modification does not affect our conclusions. f Case 0, corresponding to P3 structure and weaker Hund's coupling than in Model 15. g Values for P3 structure. with a zigzag ordering with some weight at the zone center. The model shown in Fig. 4 (d) showcases how different interaction strengths for the Z bond results in weakly broken C 3 symmetry in the structure factors. The bottom three rows of Fig. 4 shows integrated INS intensities for three different energy windows. Experimentally, a star-like pattern with strong weight at Γ and arms extending to the M points was observed in the ω ∈ [4.5, 7.5] meV energy window [23]. We dub this pattern the M star. The only model displaying this pattern in the right energy window is the one due Yadav et al. [48], in Fig. 4 (e). In contrast, (b), (c) and (d) have star-like patterns where the arms extend towards the K points -K star shapes. The two ab initio models in (a) and (f) do not capture the weight at Γ at all, instead forming a flower-like shape that we would expect to see for lower frequencies, since the peak at Γ is observed to be higher energy than the M point peak  (2.69 ± 0.11 meV vs. 2.2 ± 0.2 meV from INS data [24]). The high-energy window ω ∈ [10.5, 20.0] meV is expected to be dominated by the continuum at Γ, which is consistent with (b), (c), (d), (e), but not the "lotus root-like" shapes in (a), and (f).
We summarize our computed neutron scattering intensity results in Table II, which also includes results for the models discussed in the Supplementary Information. The ab initioinspired approach of Winter et al. [27] was constructed to reproduce certain features in the INS spectrum, and thus does particularly well. It reproduces the Γ point intensity profile well (see Supplementary Figure 3), and has an M star shape in the [5.5, 8.5] meV window, but not in [4.5, 7.5] meV. It is thus natural to use this model as a starting point for INS datacompatible effective Hamiltonians, as done in the THz spectroscopy fit of Ref. [40], and an analysis of the magnon thermal Hall conductivity in Ref. [52]. The latter work (for results see the Supplementary Note 3) proposes a particularly minor change -only reducing the magnitude of J 3 from 0.5 meV to 0.1125 meV while keeping other parameters fixed -which actually leads to an M star shape in the relevant window, but also significantly alters the intensity profile at the Γ point. We mention this fact explicitly as an example of a more general observation: for these models even small changes to the parameters can result in significantly different spectra, while even significantly different models can produce very similar SSFs and the same magnetic order. This difficulty calls for other methods to constrain the possible effective Hamiltonians, which is why we will later study the magnetic specific heat.

Evolution of INS spectra
Above we have provided results for spin Hamiltonians with multiple interaction parameters. To untangle the roles and effects of different interaction terms, we now focus on a minimal parametrization where When χ = 0 this reduces to the notation used in Ref. [61]. Following the typical hierarchy of interaction strengths in Table I (|K 1 | ≥ |Γ 1 | ≥ |J 1 | ∼ |Γ 1 | ∼ |J 3 |) we begin by assuming a dominant FM Kitaev interaction, and introduce other terms one by one. A representative selection of the resulting spectra is shown in Fig. 5. Additional parameter values are studied in Supplementary Note 5. As Fig. 5 (a) shows, the FM Kitaev limit has a flat spectrum with intensity peaked at the Γ point, consistent with the exact theoretical result [64]. The spectral evolution away from this point can be qualitatively understood using previously ob-tained phase diagrams [61,62]. For Γ 1 /K 1 = −1/2 shown in Fig. 5 (b), a sharp low-energy peak develops at the center point between Γ and K 1 , "K 1 /2", signaling a tendency towards spiral order. The peak at the Γ point remains strong, however, preventing a clear signal of the spiral phase in the static spin structure factor. We also note that the resolution of the spiral phase ordering vector is limited by the finite cluster size, and that some zigzag correlations remain at the M points. As Γ 1 is increased to Γ 1 /K 1 = −1.0 in Fig. 5 (c) the intensities at K 1 /2 and the M points become comparable in strength. As can be seen by comparison with the Kitaev limit, the presence of Γ 1 > 0 in (b) and (c) tends to produce a stronger excitation continuum, that stretches to higher frequencies.
We next introduce the nearest neighbor Heisenberg exchange J 1 . Figs. 5 (d) and (e) show the effect of adding antiferromagnetic and ferromagnetic J 1 , respectively (J 1 /K 1 = ∓0.1). In this case J 1 > 0 produces a stronger peak at K 1 /2 and weaker intensity at the Γ point, signaling a stabilized spiral order phase. This is consistent with the classical phase We focus on i) the positions ω Γ and ω M of the initial spin wave peaks in the INS intensity at the Γ and M points, respectively, ii) the shape of the neutron scattering intensity map (IM) in momentum space integrated over [4.5, 7.5] meV, and iii) the position of the high-temperature peak in the magnetic specific heat. K (M) star denotes a star-like shape pointing towards the K (M) points. The bolded models are considered in the main text. Results for the other models are given in the Supplementary  diagram of Ref. [61] and cluster mean field theory results of Ref. [63]. In contrast, J 1 < 0 sees the peaks at the Γ point move down in frequency, a strengthening of the M peaks and significant reduction in intensity at K 1 /2. These observations are consistent with moving into the regime of ferromagnetic ordering. We note that a strong continuum remains at the zone center. In the case of Γ 1 /K 1 = 0 we would have had a weaker continuum, and a stronger low-energy peak at the zone center. Finally we also introduce Γ 1 < 0 and J 3 > 0 in Figs. 5 (f) and (g). These interactions both stabilize the zigzag order, while pushing the Γ point peaks to higher frequency, and generally weakening the continuum nature of the excitation spectra. This suggests that the unrealistically large gaps at the Γ points in Figs. 3 (a) and (f) may be due to an overestimation of the Γ 1 and J 3 parameters.

Magnetic specific heat
As shown in Fig. 6 (a), the magnetic specific heat of the pure Kitaev model features two characteristic, well-separated peaks at T l and T h , where the low-T one is due to thermal fluctuations of localized Majorana fermions, and the high-T peak is related to itinerant Majoranas [65,66]. This two-peak structure appears to be stable to small perturbations away from the Kitaev point [6,67,68]. Note that the presence of two peaks is not itself a unique signature of Kitaev physics [67,69] and occurs also for e.g. the Γ model [68,70], see Fig. 6 (a).
A similar two-peak structure has been found in α-RuCl 3 experiments, using both RhCl 3 [34] and ScCl 3 [25,71,72] as nonmagnetic analogue compounds. In clean samples, a sharp low-T peak representing the magnetic ordering occurs at T l ≈ 6.5 K [30,34], and then a broader peak occurs at a higher temperature T h , followed by a (non-magnetic) structural transition of α-RuCl 3 near 165 − 170 K [25,34] [72] find a broad maximum around 80 − 100 K), but it appears to be an order of magnitude larger than T l . Whether or not the T h peak can be attributed to fractionalized excitations due to a proximate Kitaev QSL, the feature appears to be real and ought to be captured by a realistic spin Hamiltonian. Two additional comments are in order. First, accurately determining the magnetic specific heat at higher temperatures is challenging, and sensitive to details of the analysis. This may partly explain the range of T h values mentioned above. Second, an optical spectroscopy study [73] extracted a crossover temperature for magnetic correlations, T ∼ 35 K, which, in an analysis relying on the pure Kitaev model, was equated to T h . In the (a) FM Kitaev limit, Evolution of the INS intensity I (q, ω) away from the ferromagnetic Kitaev limit. The limit is shown in (a). In (b), (c) a Γ 1 > 0 term is added. The cases of AFM J 1 and FM J 1 are considered in (d) and (e), respectively. Finally Γ 1 < 0 and J 3 > 0 are introduced in (f) and (g). These results were obtained using the 24-site cluster.
following we will rely mainly on data from Widmann et al. [34], but there is clearly some uncertainty to the value of T h . In Figs. 6 (b)-(c) we plot the magnetic specific heat, C mag (T ), for the six considered models on 24-site clusters, along with the excess heat capacity determined in Ref. [34]. We note that the finite-size clusters are far from the thermodynamic limit, so we cannot expect to numerically observe a sharp magnetic transition, but the location of the peaks can provide useful information. We see that the models plotted in (b) are clearly inconsistent with the experimental data. However, the two ab initio models in (c) (which did not capture the INS data) actually have peak positions that are consistent with the data. In fact, the model fully determined from firstprinciples in Eichstaedt et al. [44] has a peak at T h ≈ 66 K, while the experimental data is centered around 70 K, with a peak at 68 K.
In Fig. 7 we provide 32-site cluster TPQ results for a subset of the models, and show the finite size scaling tendencies. The two cluster sizes have different symmetry properties, which could explain part of the differences. Unfortunately, going to even larger cluster sizes (for better scaling or to preserve symmetries) using the TPQ method becomes computationally prohibitive. We find that the position of the high-temperature peak changes only marginally (see Supplementary Note 2 for details), while the low-temperature behavior is much less well-converged. We thus conclude that the two ab initio models describe the magnetic specific heat bet-ter than the other models.

Modified ab initio model
Having established that the two ab initio models are consistent with the experimental specific heat, we now ask whether they can be modified to better describe the INS data. As discussed in the section on evolution of INS spectra, the large gaps at the zone centers in Figs. 3 (a) and (f) may be due to overestimated Γ 1 or J 3 values. Since Γ 1 is sensitive to the degree of trigonal distortion, it can be expected to vary between crystal samples. It is thus likely the parameter with the highest degree of uncertainty. For these reasons, we consider the effect of reducing |Γ 1 | in the model of Ref. [44], while leaving other parameters (J 3 included) unchanged. We use bondaveraged interaction parameters.
We again find that the spin wave gap at the Γ point depends strongly on the value of Γ 1 , and that the low-energy features of the INS spectrum can be well explained when |Γ 1 | is significantly reduced but still finite. Specifically, we take Γ 1 → 0.05Γ 1 . (The full parameter set is given in Supplementary Table 2.) The spectrum for this case is shown in Fig. 8 (a). In this case, we find ω Γ = 2.5 meV, ω M = 2 meV, close to the values obtained in Ref. [24]: ω Γ = 2.69 ± 0.11 meV, ω M = 2.2 ± 0.2 meV. We note that our value for ω Γ is consistent with the the low-energy magnon energy of 2.5 meV shows C mag for the six chosen models for α-RuCl 3 . For comparison, the experimentally determined excess heat capacity from Ref. [34] is plotted using black dots. The peak in the experimental data near 6.5 K signals the magnetic ordering, and the strong peak at 170 K is a structural transition, unrelated to the magnetic specific heat. Finally, the peak near 70 K may correspond to itinerant Majorana quasiparticles [34,65,66]. The peak position is inconsistent with the models plotted in (b), but consistent with the ab initio models plotted in (c), with higher interaction strengths. . The data is averaged over 15 initial TPQ vectors, and the shaded regions show the standard deviations. We generically observe a two-peak structure, the high-T peak of which appears to be well converged. The lower temperature peak corresponds to magnetic ordering, and may be quite sensitive to finite size effects, or to the difference in symmetry between the 24-and 32-site clusters.
observed using THz spectroscopy [40,41]. As we noted earlier, the relative intensity of the Γ point peak may be enhanced at finite temperature [51]. Fig. 8 (b) shows the SSF, which remains consistent with zigzag order, and I(q, ω) integrated over the [4.5, 7.5] meV range. From this integrated intensity and Fig. 8 (a) we see that there is a lack of intensity at the zone center within the chosen energy range, unlike the notable star shape in Ref. [23]. The source of this discrepancy is not clear, and calls for further study and parameter refinement. Finally, in (c) we show the magnetic specific heat calculated for this modified model. We do find T h ≈ 83 K, which is higher than the ≈ 70 K reported in Ref. [34], yet consistent with the broader range of values proposed (70 − 100 K). DISCUSSION We have found that there is a considerable qualitative difference between proposed spin Hamiltonians that describe the INS data well, and realistic models derived using ab initio methods, which are consistent with the reported magnetic specific heat observations. This difference is accompanied by a significant discrepancy in overall energy scales. The specific heat measurements probe the energy density of states, and should represent a good guide to the energy scale, provided the phonon background is handled adequately. In light of our results we thus expect the Kitaev and off-diagonal couplings strengths to be larger, and that α-RuCl 3 may be closer to the QSL regime than previously believed. (We note that recent shows I(q, ω) integrated over [4.5, 7.5] meV, with less intensity at the zone center than is experimentally observed. (c) shows the magnetic specific heat calculated using TPQ and 15 random initial vectors compared with the experimental data from Ref. [34]. The shaded region shows the standard deviation. The calculations were all done using the 24-site cluster.
anisotropic susceptibility [74] and THz spectroscopy experiments [41] are also consistent with higher Kitaev strengths than in Model 2.) In contrast, the calculated dynamical spin structure factors and INS intensities are much more sensitive to the relative strengths of different interaction terms. They are particularly useful probes for models with fewer degrees of freedom, such as the At the same time, static properties such as magnetic order or SSFs, are clearly insufficient to fully constrain the α-RuCl 3 spin Hamiltonians. In this respect, properties in the presence of magnetic fields, such as phase transitions and the magnon thermal Hall effect [52], present a particularly promising direction for both theory and experiment. By using one of the ab initio models (Eichstaedt et al. [44]) as a starting point and reducing the magnitude of Γ 1 , we are able to identify a set of parameters (full values are given in Supplementary Table 2) that partly resolve the discrepancy between the two classes of models mentioned above. We find low-energy peaks in the INS spectrum and a high-temperature peak in the magnetic specific heat that are consistent with experiment. However, this should not yet be considered a fully accurate model, as there is an unexplained lack of intensity at the zone center at intermediate frequencies. Instead, we consider it a new starting point.
With these results in mind, we now return to the question we posed in the abstract, about the nature of the correct spin Hamiltonian for α-RuCl 3 . From a variety of ab initio and DFT calculations, we expect a minimal model to include ferromagnetic nearest-neighbor Kitaev and Heisenberg couplings, Γ 1 > 0, a Γ 1 < 0 term, and a small J 3 > 0. Our results for the modified ab initio model further suggest that Γ 1 should be small, but finite. Since both Γ 1 and J 3 act to stabilize the zigzag order, small values are consistent with the fact that a relatively weak in-plane magnetic field can take α-RuCl 3 out of the ordered phase. Alternatively, α-RuCl 3 might be close to a quantum critical point [32,33,51], which would be a very exciting scenario. Anisotropic susceptibility measurements [74] point towards significant off-diagonal Γ 1 and Γ 1 terms, which may also help stabilize the purported spin liquid phase at finite magnetic fields [68,75,76]. At this point it is not clear whether anisotropies between bonds or the interlayer coupling play a qualitative role, but they are also expected in a full model.
We hope that our results can help guide further theory development and interpretation of experimental results going forward, both for α-RuCl 3 and other Kitaev spin liquid candidate materials. Since we found the INS predictions to be particularly sensitive to small parameter changes, it would be very useful to consider additional modeling techniques and additional observables. For example, machine learning methods may be a promising way to efficiently handle the highdimensional parameter space. In addition, further experiments in applied magnetic fields can help constrain the Hamiltonian by suppressing fluctuations.

METHODS
We use the HΦ [77] library for numerical calculations on finite-size systems. We employ a 24-site cluster with C 3 symmetry, and a rhombic 32-site cluster, see Fig. 2. The momenta compatible with the finite-size clusters are shown in Supplementary Note 1. Finite-temperature specific heat is computed using the microcanonical thermal pure quantum state (TPQ) method [59,67,78], and averaged over ≥ 15 random initial vectors. The key idea behind the TPQ method is that a quantum system at thermal equilibrium can be reliably described by a single, iteratively constructed state. Utilizing this fact allows for a significant reduction in computational cost compared with finite-temperature exact diagonalization methods.
Zero temperature properties are calculated using the Lanczos exact diagonalization method, and the continued fraction expansion (CFE) [58] is used to compute the dynamical quantities. A total of 500 Lanczos steps are used to calculate the CFE. We take 1 meV as a representative value for the experimental energy resolution at full-width/half-maximum [23,24], and emulate it in the exact diagonalization calculations by using a Lorentzian broadening of 0.5 meV. For the INS spectral evolution calculations we used a Lorentzian broadening of 0.05 in the fixed energy scale. The neutron scattering intensity I (q, ω) is defined [27,51] where f (q) is the magnetic form factor, q a is the projection of the momentum vector onto the spin components in the local cubic coordinate system also used for the spin Hamiltonian, and S µν (q, ω) is the dynamical spin structure factor at momentum q and frequency ω, Note that the off-diagonal elements of S µν (q, ω) contribute significantly for most models studied here, due to the presence of Γ 1 and Γ 1 interactions. The static spin structure factor, S (q) = S (q, ω) dω, is evaluated separately. We note that neutron scattering experiments probe the magnetization M, while the Hamiltonians in Eqs. (1) and (2) are expressed in terms of the pseudospin S. Hence the form of Eq. (9) amounts to an assumption that M and S are approximately parallel. However, trigonal distortion can induce an angle between the two vectors, and a resulting g-factor anisotropy [79]. While there have been conflicting reports about the degree of anisotropy [48,72,80], a more recent X-ray absorption spectroscopy study [81] found α-RuCl 3 to have only weak trigonal distortion and a nearly isotropic g-factor. In light of this result and the relatively weak Γ 1 interactions in Table I we assume that M and S are indeed approximately parallel in α-RuCl 3 . f (q) is assumed to be isotropic, which is justified for small scattering wave numbers [23]. The magnetic form factor for Ru 3+ was calculated using DFT in the Supplementary Material of Ref. [25]. By fitting their data to a Gaussian we have obtained the analytical approximation We integrate over the momentum direction perpendicular to the honeycomb plane, following the experiment [23]. Since the ED calculation is necessarily two-dimensional we assume that S µν (q, ω) is constant along the perpendicular direction during the integration step. We expect this to be a reasonable approximation due to the strong two-dimensionality of α-RuCl 3 [82], and the relatively small interlayer coupling [32].

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. ACKNOWLEDGMENTS We thank C. Balz, A. Banerjee, T. Berlijn, S. E. Nagler, A. M. Samarakoon, and D. A. Tennant for helpful discussions, and A. Loidl for providing the magnetic specific data. We thank Y. Yamaji both for useful discussions and assistance with HΦ. We are also grateful to an anonymous referee for suggesting a systematic approach to study the spectral evolution away from the Kitaev limit, which lead us to the improved interaction parameters proposed in this work.
The In this Supplement, we present results for the Hamiltonians not covered extensively in the main text, as well as more details for the models covered in the main text. In particular, we include comparisons of the predicted inelastic neutron scattering (INS) profiles at the Γ and M points with experimental data from Ref. [1] for all models, as well as additional results for the magnetic specific heat and structure factors. All results reported in the Supplement, except where otherwise explicitly noted, were obtained using the N = 24 cluster.

SUPPLEMENTARY NOTE 1: FINITE-SIZE CLUSTER MOMENTA
The finite-size clusters are compatible with the momenta shown in Supplementary Fig. 1. To plot the static spin structure factors (SSFs), and energy slices of the INS intensity, I (q, ω), over reciprocal space, we use cubic interpolation over a hexagonal grid centered in the small hexagons.

SUPPLEMENTARY NOTE 2: ADDITIONAL MAGNETIC SPECIFIC HEAT RESULTS
Supplementary Fig. 2 shows the magnetic specific heat for the models not considered in the main text. Two of the panels also include comparisons to a few models that were discussed in the main text. This is the case for (a), where the Winter et al. Nat. Comms. model [2] is shown along with two related models, due to Cookmeyer and Moore [3], and Wu et al. [4], respectively. In (c) the different cases considered in Eichstaedt et al. [5] are shown. Both these panels clearly indicate that, for related models, the high-temperature peak T h and the C(T ) profile vary mostly smoothly with the overall energy scale of the model. In (a) we also note that the Cookmeyer and Moore parameters produce a more prominent low-temperature peak than the Winter et al. Nat. Comms. model. We interpret this as an effect of the weaker third-nearest Heisenberg exchange, J 3 , which results in a less stable zigzag ordering and puts the system closer to the Kitaev limit, for which a two-peak structure is expected.
In Supplementary Fig. 2 (b), the models of Ran et al. [9] and Wang et al. [7] both have high-temperature peaks near ≈ 55 K, which is relatively close to the experimental value (≈ 70 K) obtained in Ref. [15]. Meanwhile, the Suzuki and Suga model [6,16] has its peak at ≈ 94 K, consistent with the peak position of 100 K obtained in Ref. [17]. In (d) results for four Hamiltonians with K 1 > 0 are shown. We note that the Kim et al. model [12,18] has its peak at a reasonable position (T h = 81 K), but a rather flat temperature dependence. Supplementary Fig. 2 (e) shows the magnetic specific heat for the two spin wave fits to THz spectroscopy data obtained by Ozel et al. [14]. The fit with K 1 > 0 has T h ≈ 16 K, and the fit with K 1 < 0 has T h = 17.5 K. Both models are inconsistent with the magnetic specific heat data, but we note that these T h values are very consistent with the similarly obtained model of Wu et al. [4], which has T h ≈ 17 K.
In Fig. 7 of the main text we contrasted the 24-and 32site cluster results for four models, arguing that the high-T peak appears to be fairly stable. To quantify this, Supplementary Table I lists the numerical positions for the high-T peaks. Overall, the relative difference ∆T h between the two clusters considered is on the order of a few percent, despite the cluster changing both size and shape. This relatively small difference suggests that the TPQ method is a useful numerical tool, particularly at high temperatures, despite being limited to small clusters. This is in line with previous results for nearest-neighbor honeycomb Kitaev-Heisenberg [19] and kagome Heisenberg [20] models. Supplementary Figure 2. Magnetic specific heat calculated using the TPQ method for various proposed α-RuCl 3 Hamiltonians, compared with experimentally determined excess heat capacity from Ref. [15]. The model numbers match those in Table I of the main text. The solid lines show the calculated average value over 15 initial vectors, and the shaded areas show the standard deviation.

SUPPLEMENTARY NOTE 3: ADDITIONAL INS RESULTS
The T = 0 INS intensity profiles at the Γ and M 1 points for the six models considered in the main text are plotted in Supplementary Fig. 3 as a function of energy. The ED calculation is compared with experimental data from Ref. [1], in which the positions of the first spin-wave peaks were estimated to be 2.69 ± 0.11meV at the Γ point, and 2.2 ± 0.2meV at the M point.
The ab initio inspired Winter et al. Nat. Comms. parameters reproduce the intensity profile at the Γ point particularly well. This is no surprise, given that the parameters were chosen to reproduce broad features of the INS spectrum, especially near the Γ point. In contrast, from panels (a) and (f), we see that the two models that most accurately predicted the high-temperature peak in the magnetic specific heat have the Supplementary Table I. Cluster dependence on the location of the high-temperature peak in the magnetic specific heat calculated using the TPQ method. Both cluster size and symmetry properties may influence the results. The larger difference for model 11 may be due to the fact that the full ab initio Hamiltonian, without C 3 symmetry was used, including bond anisotropies and second nearest neighbor interactions. The relative difference is defined intensity at the Γ point shifted to much higher frequencies.
We next study the fully ab initio-determined parameters from Eichstaedt et al. [5] more closely. These parameters were derived using density functional theory, constrained RPA, and perturbation theory (t/U expansion), and the calculations included nonlocal Coulomb interaction and spin-orbit coupling terms. This fact makes it interesting to consider approximate versions of the full Hamiltonian, to better understand the role of different effects. The I (q, ω) spectra are shown in Supplementary Fig. 4, while the static spin structure factor S (q) and energy-integrated slices of I (q, ω) are shown in Supplementary Fig. 5. The full Hamiltonian is bondanisotropic, which results in inequivalent M peaks. Bond averaging makes the model C 3 symmetric, and the M points equivalent. As can be seen in Table I of the main text, the non-local SOC only has a small effect on the spin interaction parameters, which leads to minor changes to I (q, ω). Finally, the INS intensity profiles for the four cases are plotted in Supplementary Fig. 6.
We turn next to the nine of the remaining Hamiltonians. Their respective I (q, ω) spectra are shown in Supplementary  Fig. 7. Five (four) of these models have ferromagnetic (antiferromagnetic) nearest-neighbor Kitaev interactions. Further results for models with K 1 < 0 are shown in Supplementary  Figs. 8 and 9. See Supplementary Figs. 10 and 11 for the models with K 1 > 0.
Among the first five models, we note that the Cookmeyer and Moore parameters [3] produce the expected SSF and intensity map in the [4.5, 7.5] meV window. In this sense it improves on the Winter et al. Nat. Comms. parameters [2]. However, the latter model yields intensity profiles at the Γ and M points that better fit the experiment. The SSF for Suzuki and Suga's model [6,16] is dominated by the Γ point peak, due to a pronounced scattering continuum extending to large energies. However, note that Supplementary Fig. 7 (b) shows that the M point peak in I (q, ω) occurs at lower energy than the Γ point peak, and is stronger in intensity, which is consistent with the zigzag order.
The Kitaev-Γ model proposed by Ran et al. [9], model 8, displays strong low-energy peaks at both the M points, and the center between Γ and K points. This spectral structure is similar to that of the pure "antiferromagnetic" Γ model. Model 9 from Hou et al. [8] produces a Γ point peak at too high ω, similar to the ab initio parameters for models 1 and 11. Model 10 from Wang et al. [7] is an interesting case. While the strongest I (q, ω) peak occurs at the M point, it also has a low-lying peak at the Γ point. Such a peak could potentially be hidden by the elastic continuum in an INS experiment, but we note that THz spectroscopy experiments down to 0.3 − 0.4meV find no signs of such a peak [4,21]. The SSF for this model has higher intensity at the Γ point than the M points. However, we note that the trace over diagonal components of the SSF, S µµ (q), has the opposite pattern, while the off-diagonal components are weaker. This suggests that the pattern in momentum space for this model may be quite sensitive to the radial dependence of the magnetic form factor.
We next turn to the four models with antiferromagnetic NN Kitaev coupling. Models 14 and 16 [10,22] both have I (q, ω) dominated by low-energy M point peaks. The Γ point peaks are too high in energy, and we do not observe extended scattering minima. Models 15 and 16 [11,12] are found to have strong weight at the K points at low frequencies, as well as in the SSF. This shows that they order in the 120 • configuration, which is consistent with the phase diagrams obtained in Ref. [23].
Finally we turn to two models obtained by Ozel et al. [14] by linear spin wave fits to THz spectroscopy data. Their spectra and intensity profiles are shown in Supplementary Fig. 12, and the static spin structure factors and energy slices of their spectra are shown in Supplementary Fig. 13. Since THz spectroscopy is a probe of the physics at the Γ point, we may expect good agreement there. We do indeed find peaks that reasonably resemble the shape of the experimental INS data. The frequency shifts in the peak position may be explained by quantum renormalization. The M point intensities are, however, not captured. Despite prominent low-energy peaks at the M and K 1 /2 points, respectively, the static spin structure factors both suggest ferromagnetic phases.

SUPPLEMENTARY NOTE 4: COMPARING THE UPDATED AND ORIGINAL EICHSTAEDT ET AL. PARAMETERS
As mentioned in Table I of the main text, the Eichstaedt et al. interaction parameters were revised in the published version [24]. Above and in the main text (except for the case of the modified ab initio model) we have used the preliminary values from the preprint version [5]. As Supplementary  Table II shows  (e) 7. Yadav et al.   We also find visually near-identical INS spectra. We plot the intensity profiles at the Γ and M 1 points in Supplementary  Fig. 15, from which it is clear that the strong low-energy features in the spectrum are nearly unchanged. We thus conclude that there is no qualitative difference between the revised and original interaction parameters.   [5] considered in the main text (in bold), compared with the revised, published parameters from Ref. [24]. The differences are only minor. The last row shows the parameters for the modified ab initio model, which is based on the revised, fully ab initio model. Dashes (-) indicate that the value is unavailable or negligible. Asterisks in the 'BA' column signify that the full Hamiltonian has different values for the X/Y bonds compared to the Z bonds, and that the parameter values given in the row have been bond averaged.