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

TL;DR

This study evaluates various effective Hamiltonians for $ ext{RuCl}_3$, revealing none fully explain experimental magnetic properties, and proposes a modified model aligning with both specific heat and neutron scattering data.

Contribution

The paper systematically compares existing models and introduces a modified ab initio Hamiltonian that better captures experimental magnetic behaviors of $ ext{RuCl}_3$.

Findings

No single current model explains all experimental data.

First-principles models match high-temperature specific heat peaks but overestimate magnon energies.

Modified model aligns with both specific heat and inelastic neutron scattering data.

Abstract

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.

Tables

Table 1. Table 1: The spin Hamiltonians for α 𝛼 \alpha -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	$J_1$	$K_1$	${\Gamma }_1$	${\Gamma }_1^{\prime }$	$J_2$	$K_2$	$J_3$	$K_3$	BA
1 Winter et al. PRB PhysRevB.93.214431 111Using the proposed minimal model, which is bond averaged and neglects small ${\Gamma }_1^{\prime }=-0.9$ meV. Values for the monoclinic ($C2/m$) crystal structure.	Ab initio (DFT + exact diag.)	$-1.7$	$-6.7$	$+6.6$	$-0.9$	–	–	$+2.7$	–	⋆
2 Winter et al. NC 10.1038/s41467-017-01177-0	Ab initio-inspired (INS fit)	$-0.5$	$-5.0$	$+2.5$	–	–	–	$+0.5$	–
3 Wu et al. PhysRevB.98.094425	THz spectroscopy fit	$-0.35$	$-2.8$	$+2.4$	–	–	–	$+0.34$	–
4 Cookmeyer and Moore PhysRevB.98.060412	Magnon thermal Hall (sign)	$-0.5$	$-5.0$	$+2.5$	–	–	–	$+0.1125$	–
5 Kim and Kee PhysRevB.93.155143	DFT + $t/U$ expansion	$-1.53$	$-6.55$	$+5.25$	$-0.95$	–	–	–	–	⋆
6 Suzuki and Suga PhysRevB.97.134424 ; Suzuki2019	Magnetic specific heat	$-1.53$	$-24.4$	$+5.25$	$-0.95$	–	–	–	–	⋆
7 Yadav et al. Yadav2016 222We use the sign convention in Refs. PhysRevB.97.134424 ; PhysRevB.96.064430 .	Quantum chemistry (MRCI)	$+1.2$	$-5.6$	$+1.2$	$-0.7$	$+0.25$	–	$+0.25$	–
8 Ran et al. PhysRevLett.118.107203	Spin wave fit to INS gap	–	$-6.8$	$+9.5$	–	–	–	–	–
9 Hou et al. PhysRevB.96.054410 333This work gives values for several values of $U$. Here we use the $U=3.5$eV parameters.	Constrained DFT +$U$	$-1.87$	$-10.7$	$+3.8$	–	–	–	$+1.27$	$+0.63$	⋆
10 Wang et al. PhysRevB.96.115103 444Values for the C2 structure.	DFT + $t/U$ expansion	$-0.3$	$-10.9$	$+6.1$	–	–	–	$+0.03$	–
11 Eichstaedt et al. PhysRevB.100.075110 ; Eichstaedt2019 555These are the parameters from the preprint version in Ref. Eichstaedt2019 . They were revised in the published version, Ref. PhysRevB.100.075110 . In Supplementary Note 4 we show that this slight modification does not affect our conclusions.	Fully ab initio (DFT + cRPA + $t/U$)	$-1.4$	$-14.3$	$+9.8$	$-2.23$	–	$-0.63$	$+1.0$	$+0.03$	⋆
12 Eichstaedt et al. PhysRevB.100.075110 ; Eichstaedt2019 55footnotemark: 5	Neglecting non-local Coulomb	$-0.2$	$-4.5$	$+3.0$	$-0.73$	–	$-0.33$	$+0.7$	$+0.1$	⋆
13 Eichstaedt et al. PhysRevB.100.075110 ; Eichstaedt2019 55footnotemark: 5	Neglecting non-local SOC	$-1.3$	$-13.3$	$9.4$	$-2.3$	–	$-0.67$	$+1.0$	$+0.1$	⋆
14 Banerjee et al. Banerjee2016	Spin wave fit	$-4.6$	$+7.0$	–	–	–	–	–	–
15 Kim et al. PhysRevB.91.241110 ; PhysRevB.96.064430	DFT + $t/U$ expansion	$-12$	$+17$	$+12$	–	–	–	–	–
16 Kim and KeePhysRevB.93.155143 666Case 0, corresponding to P3 structure and weaker Hund’s coupling than in Model 15.	DFT + $t/U$ expansion	$-3.5$	$+4.6$	$+6.42$	$-0.04$	–	–	–	–
17 Winter et al. PRB PhysRevB.93.214431 777Values for P3 structure.	Ab initio (DFT + exact diag.)	$-5.5$	$+7.6$	$+8.4$	$+0.2$	–	–	$+2.3$	–
18 Ozel et al. PRB PhysRevB.100.085108	Spin wave fit / THz spectroscopy	$-0.95$	$+1.15$	$+3.8$	–	–	–	–	–
19 Ozel et al. PRB PhysRevB.100.085108	Spin wave fit / THz spectroscopy	$+0.46$	$-3.50$	$+2.35$	–	–	–	–	–

Table 2. Table 2: Summary of results, highlighting important features in the INS and magnetic specific heat predictions on the 24 24 24 -site cluster. We focus on i) the positions ω Γ subscript 𝜔 Γ \omega_{\Gamma} and ω M subscript 𝜔 M {\omega_{\mathrm{M}}} of the initial spin wave peaks in the INS intensity at the Γ Γ \Gamma and M points, respectively, ii) the shape of the neutron scattering intensity map (IM) in momentum space integrated over [ 4.5 , 7.5 ] 4.5 7.5 [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 Information.

Reference	${\omega }_{\Gamma }$ [meV]	${\omega }_M$ [meV]	IM shape	$T_h$ [K]
1 Winter et al. PRB PhysRevB.93.214431	$14.2$	$2.55$	Flower	54
2 Winter et al. NC 10.1038/s41467-017-01177-0	$2.8$	$0.75$	K star	22
3 Wu et al. PhysRevB.98.094425	$1.18$	$0.54$	Dominated by $\Gamma$	17
4 Cookmeyer and Moore PhysRevB.98.060412	$1.67$111There is also a clear, distinct peak at lower frequency, which would be hidden by the elastic scattering continuum.	$0.57$	M star	24
5 Kim and Kee PhysRevB.93.155143	$4.12$	$1.49$	K star-like	34
6 Suzuki and Suga PhysRevB.97.134424 ; Suzuki2019	$2.05$	$1.98$	K star	94
7 Yadav et al. Yadav2016	$2.58$	$0.58$	M star	13
8 Ran et al. PhysRevLett.118.107203	$4.8$	$1.37$	M star	57
9 Hou et al. PhysRevB.96.054410	$5.5$	$1.87$	K star	33
10 Wang et al. PhysRevB.96.115103	$3.13$11footnotemark: 1	$1.03$	K star	53
11 Eichstaedt et al. PhysRevB.100.075110	$11.2$	$2.79$	Flower / M3 222Flower shape when bond-averaged, otherwise dominated by M3. Peak positions are given for the non-bond averaged case.	66
12 Eichstaedt et al. PhysRevB.100.075110	$5.1$	$0.98$	Lotus root	22
13 Eichstaedt et al. PhysRevB.100.075110	$11.9$	$2.43$	Dominated by M3	63
14 Banerjee et al. Banerjee2016	$6.38$	$1.21$	Lotus root	21
15 Kim et al. PhysRevB.91.241110 ; PhysRevB.96.064430	$6.52$	$2.65$	Flower	81
16 Kim and Kee PhysRevB.93.155143	$4.39$	$2.25$	Dominated by $\Gamma$ and K	35
17 Winter et al. PRB PhysRevB.93.214431	$7.35$	$0.43$	Ring	41
18 Ozel et al. PRB PhysRevB.100.085108	$4.39$	$0.31$	Dominated by $\Gamma$	16
19 Ozel et al. PRB PhysRevB.100.085108	$3.18$	$3.92$	Dominated by $\Gamma$	18

Equations

$$H_{J_1-K_1-{\Gamma }_1-J_3}$$

$$+J_3\sum _{⟨⟨⟨i,j⟩⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$H_{{\Gamma }_1^{\prime }}$$

$${\Gamma }_1$$

$$J_1$$

$$K_1$$

$${\Gamma }_1^{\prime }$$

$$J_3$$

$$I\left(\mathbf{q},\omega \right)$$

$$S^{\mu \nu }\left(\mathbf{q},\omega \right)$$

$$f(q)$$

Full text

Introduction

Quantum spin liquids (QSL) are long-sought-after states of matter without magnetic order, but with nontrivial topological and potentially exotic properties Savary2017 ; RevModPhys.89.025003 . Much of the search has been focused on frustrated lattice systems Balents2010 ; RevModPhys.88.041002 , but in an important development in 2006 Kitaev Kitaev2006 introduced a novel exactly solvable paradigmatic QSL with bond-directional Ising terms on the bipartite honeycomb lattice. Importantly, this Kitaev model hosts anyonic excitations Hermanns2018 , which are of interest both for fundamental reasons and for their proposed application in topological quantum computing Kitaev2003 ; RevModPhys.80.1083 . It was realized that such interaction terms naturally appear Khaliullin2005 — and can be large — in Mott-insulating transition-metal systems with edge-sharing octahedra and strong spin-orbit coupling, such as in A2IrO3 (A=Na,Li) PhysRevLett.102.017205 ; PhysRevLett.105.027204 , $\alpha$-RuCl3 PhysRevB.90.041112 and other materials Winter2017 ; Takagi2019 .

However, the three mentioned materials are all found to order magnetically at low temperatures, and hence cannot be perfect realizations of the Kitaev model. Na2IrO3 PhysRevB.83.220403 ; PhysRevLett.108.127204 ; PhysRevB.85.180403 and $\alpha$-RuCl3 PhysRevB.91.144420 ; PhysRevB.92.235119 ; PhysRevB.93.134423 ; Banerjee2016 both develop a zigzag order, as illustrated in Fig. 1 (a), while Li2IrO3 displays an incommensurate spiral order PhysRevB.93.195158 . Despite the zigzag order, $\alpha$-RuCl3 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) Banerjee2017 ; Banerjee2018 ; Do2017 ; PhysRevLett.118.107203 , 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 10.1038/s41467-017-01177-0 . More recently, half-integer quantization of the thermal Hall conductivity was reported 10.1038/s41586-018-0274-0 , 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 PhysRevLett.119.227208 ; PhysRevLett.119.037201 ; PhysRevB.95.180411 ; Banerjee2018 ; Balz2019 . Further evidence for Kitaev physics has been found in experiments reporting magnetic specific heat PhysRevB.96.041405 ; Do2017 ; PhysRevB.99.094415 , NMR PhysRevLett.119.037201 ; Jansa2018 , microwave absorption PhysRevB.98.184408 , Raman scattering PhysRevLett.114.147201 ; Du2018 ; Wang2018 , and THz spectroscopy results PhysRevB.98.094425 ; Reschke2019 ; PhysRevLett.119.227201 ; PhysRevLett.119.227202 .

Altogether, these experiments strongly suggest that $\alpha$-RuCl3 can be described by a generalized Kitaev-Heisenberg Hamiltonian Winter2017 ; Takagi2019 , including off-diagonal and further-range interactions. Theoretical work leads to the same picture, whether the model is deduced from ab initio methods PhysRevB.100.075110 ; PhysRevB.90.041112 ; PhysRevB.91.241110 ; PhysRevB.93.155143 ; PhysRevB.93.214431 ; Yadav2016 ; PhysRevB.96.054410 ; PhysRevB.96.115103 , or from more phenomenological or ab initio inspired approaches 10.1038/s41467-017-01177-0 ; PhysRevLett.120.077203 ; PhysRevB.98.060412 ; PhysRevB.97.134424 . Unfortunately, these different works, as well as experimental fits PhysRevLett.118.107203 ; PhysRevB.98.094425 , have led to a veritable zoo of proposed realistic spin Hamiltonian descriptions for $\alpha$-RuCl3, 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 1.

In this work, we adopt a systematic approach to address this uncertainty. We calculate static spin structure factors (SSFs) $S\left(\mathbf{q}\right)$ and INS intensities $I\left(\mathbf{q},\omega\right)$ for all models listed in Table 1 using Lanczos exact diagonalization (ED) RevModPhys.66.763 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_{\mathrm{mag}}$ for the models using the thermal pure quantum state (TPQ) method PhysRevLett.108.240401 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. 10.1038/s41467-017-01177-0 ; PhysRevB.97.134424 ; Suzuki2018b . 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 $\alpha$-RuCl3 Winter2017 ; Takagi2019 . 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_{\mathrm{mag}}$, but the parameters proposed in the literature push the INS intensity at the $\Gamma$ 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_{\mathrm{mag}}$. By modifying one of the ab initio models, we are able to find results consistent with both $C_{\mathrm{mag}}$ and low-energy features of the INS spectrum. Our results thus provide important clues for an accurate and realistic description of of $\alpha$-RuCl3.

Results

Several of the Hamiltonians listed in Table 1 are special cases of a proposed minimal model PhysRevB.93.214431 ; Winter2017 for $\alpha$-RuCl3,

[TABLE]

where $\langle\dots\rangle$ and $\langle\langle\langle\dots\rangle\rangle\rangle$ denote nearest and third-nearest neighbors, respectively. $\gamma=$X,Y,Z is the bond index shown in Fig. 1 (b), and $\alpha,\beta$ are the two other bonds. The $\Gamma_{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 $\Gamma_{1}^{\prime}$ interaction,

[TABLE]

which originates from trigonal distortion PhysRevLett.112.077204 ; Rau2014 . Since the crystal structure of $\alpha$-RuCl3 features trigonal compression Winter2017 ; PhysRevB.92.235119 ; PhysRevB.93.134423 ; Do2017 perturbative calculations Rau2014 ; PhysRevB.93.214431 would predict that $\Gamma_{1}^{\prime}<0$, which provides an additional mechanism to stabilize the zigzag order Rau2014 ; PhysRevB.99.064425 . In the absence of other crystal distortions, the most general nearest-neighbor Hamiltonian would thus be the $J_{1}-K_{1}-\Gamma_{1}-\Gamma_{1}^{\prime}$ model. Combining these two proposed minimal models results in the $J_{1}-K_{1}-\Gamma_{1}-\Gamma_{1}^{\prime}-J_{3}$ model.

Further proposed extensions include second-nearest-neighbor Kitaev and Heisenberg terms, third-nearest-neighbor Kitaev terms, and additional symmetry-allowed anisotropies PhysRevB.93.214431 ; PhysRevB.100.075110 . In particular, $\alpha$-RuCl3 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 1 we have bond averaged such anisotropies for the sake of clarity, but we will use the full parameter sets in our calculations when appropriate.

.1 Spin structure factors and neutron scattering intensities

Fig. 3 shows predicted zero-temperature neutron scattering intensity spectra, $I\left(\mathbf{q},\omega\right)$, 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 $\Gamma$ point, which is inconsistent with experimental observations Banerjee2018 . However, the M peaks could potentially be suppressed at finite temperatures PhysRevLett.120.077203 . The models in Fig. 3 (b), (c), (d) and (e) all show clear signs of the scattering continuum at the $\Gamma$ 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 $\Gamma$ point.

In the top row of Fig. 4 we have plotted the static spin structure factors $S\left(\mathbf{q}\right)$. As shown, all six models are consistent 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 $\Gamma$ and arms extending to the M points was observed in the $\omega\in[4.5,7.5]$ meV energy window Banerjee2017 . 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. Yadav2016 , 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 $\Gamma$ at all, instead forming a flower-like shape that we would expect to see for lower frequencies, since the peak at $\Gamma$ is observed to be higher energy than the M point peak ($2.69\pm 0.11$ meV vs. $2.2\pm 0.2$ meV from INS data Banerjee2018 ). The high-energy window $\omega\in[10.5,20.0]$ meV is expected to be dominated by the continuum at $\Gamma$, 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 2, which also includes results for the models discussed in the Supplementary Information. The ab initio-inspired approach of Winter et al. 10.1038/s41467-017-01177-0 was constructed to reproduce certain features in the INS spectrum, and thus does particularly well. It reproduces the $\Gamma$ 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 data-compatible effective Hamiltonians, as done in the THz spectroscopy fit of Ref. PhysRevB.98.094425 , and an analysis of the magnon thermal Hall conductivity in Ref. PhysRevB.98.060412 . 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 $\Gamma$ 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.

.2 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 $J_{1}-K_{1}-\Gamma_{1}-\Gamma_{1}^{\prime}-J_{3}$ model. We fix the energy scale $1=\sqrt{J_{1}^{2}+K_{1}^{2}+\Gamma_{1}^{2}+\left(\Gamma_{1}^{\prime}\right)^{2}+J_{3}^{2}}$, and use a hyperspherical parametrization where

[TABLE]

When $\chi=0$ this reduces to the notation used in Ref. PhysRevLett.112.077204 . Following the typical hierarchy of interaction strengths in Table 1 ($|K_{1}|\geq|\Gamma_{1}|\geq|J_{1}|\sim|\Gamma_{1}^{\prime}|\sim|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 $\Gamma$ point, consistent with the exact theoretical result PhysRevLett.112.207203 . The spectral evolution away from this point can be qualitatively understood using previously obtained phase diagrams PhysRevLett.112.077204 ; Rau2014 . For $\Gamma_{1}/K_{1}=-1/2$ shown in Fig. 5 (b), a sharp low-energy peak develops at the center point between $\Gamma$ and K1, “K${}_{1}/2$”, signaling a tendency towards spiral order. The peak at the $\Gamma$ 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 $\Gamma_{1}$ is increased to $\Gamma_{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 $\Gamma_{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}=\mp 0.1$). In this case $J_{1}>0$ produces a stronger peak at K${}_{1}/2$ and weaker intensity at the $\Gamma$ point, signaling a stabilized spiral order phase. This is consistent with the classical phase diagram of Ref. PhysRevLett.112.077204 and cluster mean field theory results of Ref. PhysRevB.99.064425 . In contrast, $J_{1}<0$ sees the peaks at the $\Gamma$ 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 $\Gamma_{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 $\Gamma_{1}^{\prime}<0$ and $J_{3}>0$ in Figs. 5 (f) and (g). These interactions both stabilize the zigzag order, while pushing the $\Gamma$ point peaks to higher frequency, and generally weakening the continuum nature of the excitation spectra. This suggests that the unrealistically large gaps at the $\Gamma$ points in Figs. 3 (a) and (f) may be due to an overestimation of the $\Gamma_{1}^{\prime}$ and $J_{3}$ parameters.

.3 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_{{\mathrm{l}}}$ and $T_{{\mathrm{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 PhysRevLett.113.197205 ; PhysRevB.92.115122 . This two-peak structure appears to be stable to small perturbations away from the Kitaev point PhysRevB.93.174425 ; Hermanns2018 ; Catuneanu2018 . Note that the presence of two peaks is not itself a unique signature of Kitaev physics PhysRevB.68.014424 ; PhysRevB.93.174425 and occurs also for e.g. the $\Gamma$ model Catuneanu2018 ; PhysRevB.98.045121 , see Fig. 6 (a).

A similar two-peak structure has been found in $\alpha$-RuCl3 experiments, using both RhCl3 PhysRevB.99.094415 and ScCl3 Do2017 ; PhysRevB.95.241112 ; PhysRevB.91.094422 as nonmagnetic analogue compounds. In clean samples, a sharp low-$T$ peak representing the magnetic ordering occurs at $T_{{\mathrm{l}}}\approx 6.5$ K PhysRevLett.119.037201 ; PhysRevB.99.094415 , and then a broader peak occurs at a higher temperature $T_{{\mathrm{h}}}$, followed by a (non-magnetic) structural transition of $\alpha$-RuCl3 near $165-170$ K Do2017 ; PhysRevB.99.094415 . So far, there is no clear consensus for the precise value of $T_{{\mathrm{h}}}$ (Widmann et al. report $T_{{\mathrm{h}}}\approx 70$ K PhysRevB.99.094415 , Do et al. find $T_{{\mathrm{h}}}\simeq 100$ K Do2017 , while Hirobe et al. PhysRevB.95.241112 and Kubota et al. PhysRevB.91.094422 find a broad maximum around $80-100$ K), but it appears to be an order of magnitude larger than $T_{{\mathrm{l}}}$. Whether or not the $T_{{\mathrm{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_{{\mathrm{h}}}$ values mentioned above. Second, an optical spectroscopy study PhysRevB.94.195156 extracted a crossover temperature for magnetic correlations, $T^{\star}\sim 35$ K, which, in an analysis relying on the pure Kitaev model, was equated to $T_{{\mathrm{h}}}$. In the following we will rely mainly on data from Widmann et al. PhysRevB.99.094415 , but there is clearly some uncertainty to the value of $T_{{\mathrm{h}}}$.

In Figs. 6 (b)–(c) we plot the magnetic specific heat, $C_{\mathrm{mag}}(T)$, for the six considered models on $24$-site clusters, along with the excess heat capacity determined in Ref. PhysRevB.99.094415 . 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 first-principles in Eichstaedt et al. PhysRevB.100.075110 has a peak at $T_{{\mathrm{h}}}\approx 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 better than the other models.

.4 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 $\Gamma_{1}^{\prime}$ or $J_{3}$ values. Since $\Gamma_{1}^{\prime}$ 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 $|\Gamma_{1}^{\prime}|$ in the model of Ref. PhysRevB.100.075110 , while leaving other parameters ($J_{3}$ included) unchanged. We use bond-averaged interaction parameters.

We again find that the spin wave gap at the $\Gamma$ point depends strongly on the value of $\Gamma_{1}^{\prime}$, and that the low-energy features of the INS spectrum can be well explained when $|\Gamma_{1}^{\prime}|$ is significantly reduced but still finite. Specifically, we take $\Gamma_{1}^{\prime}\rightarrow 0.05\Gamma_{1}^{\prime}$. (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 $\omega_{\Gamma}=2.5$ meV, ${{\omega_{\mathrm{M}}}}=2$ meV, close to the values obtained in Ref. Banerjee2018 : $\omega_{\Gamma}=2.69\pm 0.11$ meV, ${\omega_{\mathrm{M}}}=2.2\pm 0.2$ meV. We note that our value for $\omega_{\Gamma}$ is consistent with the the low-energy magnon energy of $2.5$ meV observed using THz spectroscopy PhysRevB.98.094425 ; Reschke2019 . As we noted earlier, the relative intensity of the $\Gamma$ point peak may be enhanced at finite temperature PhysRevLett.120.077203 .

Fig. 8 (b) shows the SSF, which remains consistent with zigzag order, and $I(\mathbf{q},\omega)$ 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. Banerjee2017 . 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_{{\mathrm{h}}}\approx 83$ K, which is higher than the $\approx 70$ K reported in Ref. PhysRevB.99.094415 , 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 $\alpha$-RuCl3 may be closer to the QSL regime than previously believed. (We note that recent anisotropic susceptibility PhysRevB.98.100403 and THz spectroscopy experiments Reschke2019 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 $J_{1}-K_{1}-\Gamma_{1}-\Gamma_{1}^{\prime}-J_{3}$ model we study. At the same time, static properties such as magnetic order or SSFs, are clearly insufficient to fully constrain the $\alpha$-RuCl3 spin Hamiltonians. In this respect, properties in the presence of magnetic fields, such as phase transitions and the magnon thermal Hall effect PhysRevB.98.060412 , present a particularly promising direction for both theory and experiment.

By using one of the ab initio models (Eichstaedt et al. PhysRevB.100.075110 ) as a starting point and reducing the magnitude of $\Gamma_{1}^{\prime}$, 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 $\alpha$-RuCl3. From a variety of ab initio and DFT calculations, we expect a minimal model to include ferromagnetic nearest-neighbor Kitaev and Heisenberg couplings, $\Gamma_{1}>0$, a $\Gamma_{1}^{\prime}<0$ term, and a small $J_{3}>0$. Our results for the modified ab initio model further suggest that $\Gamma_{1}^{\prime}$ should be small, but finite. Since both $\Gamma_{1}^{\prime}$ 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 $\alpha$-RuCl3 out of the ordered phase. Alternatively, $\alpha$-RuCl3 might be close to a quantum critical point PhysRevB.96.041405 ; PhysRevLett.120.077203 ; Balz2019 , which would be a very exciting scenario. Anisotropic susceptibility measurements PhysRevB.98.100403 point towards significant off-diagonal $\Gamma_{1}$ and $\Gamma_{1}^{\prime}$ terms, which may also help stabilize the purported spin liquid phase at finite magnetic fields Catuneanu2018 ; Gordon2019 ; PhysRevB.99.224409 . 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 $\alpha$-RuCl3 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 high-dimensional parameter space. In addition, further experiments in applied magnetic fields can help constrain the Hamiltonian by suppressing fluctuations.

Methods

We use the $\mathcal{H}\Phi$ Kawamura2017 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 PhysRevLett.108.240401 ; PhysRevB.93.174425 ; Sugiura2017 , and averaged over $\geq{}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) RevModPhys.66.763 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 Banerjee2017 ; Banerjee2018 , 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\left(\mathbf{q},\omega\right)$ is defined 10.1038/s41467-017-01177-0 ; PhysRevLett.120.077203

[TABLE]

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^{\mu\nu}\left(\mathbf{q},\omega\right)$ is the dynamical spin structure factor at momentum $\mathbf{q}$ and frequency $\omega$,

[TABLE]

Note that the off-diagonal elements of $S^{\mu\nu}\left(\mathbf{q},\omega\right)$ contribute significantly for most models studied here, due to the presence of $\Gamma_{1}$ and $\Gamma_{1}^{\prime}$ interactions. The static spin structure factor, $S\left(\mathbf{q}\right)=\int S\left(\mathbf{q},\omega\right)\,\mathrm{d}\omega$, is evaluated separately. We note that neutron scattering experiments probe the magnetization $\mathbf{M}$, while the Hamiltonians in Eqs. (1) and (2) are expressed in terms of the pseudospin $\mathbf{S}$. Hence the form of Eq. (9) amounts to an assumption that $\mathbf{M}$ and $\mathbf{S}$ are approximately parallel. However, trigonal distortion can induce an angle between the two vectors, and a resulting $g$-factor anisotropy PhysRevB.94.064435 . While there have been conflicting reports about the degree of anisotropy PhysRevB.91.180401 ; PhysRevB.91.094422 ; Yadav2016 , a more recent X-ray absorption spectroscopy study PhysRevB.96.161107 found $\alpha$-RuCl3 to have only weak trigonal distortion and a nearly isotropic $g$-factor. In light of this result and the relatively weak $\Gamma_{1}^{\prime}$ interactions in Table 1 we assume that $\mathbf{M}$ and $\mathbf{S}$ are indeed approximately parallel in $\alpha$-RuCl3.

$f(q)$ is assumed to be isotropic, which is justified for small scattering wave numbers Banerjee2017 . The magnetic form factor for Ru3+ was calculated using DFT in the Supplementary Material of Ref. Do2017 . By fitting their data to a Gaussian we have obtained the analytical approximation

[TABLE]

We integrate over the momentum direction perpendicular to the honeycomb plane, following the experiment Banerjee2017 . Since the ED calculation is necessarily two-dimensional we assume that $S^{\mu\nu}\left(\mathbf{q},\omega\right)$ 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 $\alpha$-RuCl3 PhysRevB.98.205110 , and the relatively small interlayer coupling Balz2019 .

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 $\mathcal{H}\Phi$. 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 research by P.L. and S.O. was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by the US Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725, and of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is managed by UT-Battelle and supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. An award of computer time was provided by the INCITE program. A portion of the work was conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility.

Competing interests

The authors declare that there are no competing interests.

Author contributions

P.L. performed the ED and TPQ calculations. S.O. supervised the study. Both authors contributed to the writing of the manuscript.