Dynamical spin structure factors of α-RuCl3

Honeycomb-lattice magnet α-RuCl3 is considered to be a potential candidate of realizing Kitaev spin liquid, although this material undergoes a phase transition to the zigzag magnetically ordered state at TN ∼ 7 K. Quite recently, inelastic neutron-scattering experiments using single crystal α-RuCl3 have unveiled characteristic dynamical properties. We calculate dynamical spin structure factors of three ab-initio models for α-RuCl3 with an exact numerical diagonalization method. We also calculate temperature dependences of the specific heat by employing thermal pure quantum states. We compare our numerical results with the experiments and discuss characteristics obtained by using three ab-initio models.


Introduction
Since Kitaev presented an exactly solvable model for a spin-liquid ground state [1], Kitaev spin liquid (KSL) has attracted much attention because of its fascinating properties. Great efforts have been payed for realizing KSL in real materials. A layered honeycomb-lattice compound α-RuCl 3 is considered to be a promising candidate for applying Kitaev physics [2]. In the specific heat measurements for α-RuCl 3 , a double-peak structure has been observed with decreasing temperature (T ) [3,4,5]. The double-peak structure of the specific heat plays a key role in investigating the Kitaev physics [6,7,8], which has been first pointed out in the pure Kitaev model [6,7]. When temperature is decreased, the spin fractionalization to the itinerant Majorana fermions and localized gauge fluxes is realized, leading to the double-peak structure. At each peak temperature, the entropy ∼ (R/2) log 2 is released, where R is the gas constant [6,7]. The double-peak structure is caused by the independent growth of the short-range and longrange spin correlation functions with decreasing temperature [7,8]. Beyond the pure Kitaev model, it has been shown numerically that the specific heat shows a double-peak structure in the magnetically ordered state, when the system is located close to the KSL phase [8]. Thus, α-RuCl 3 is probably located close to the KSL phase in spite of a phase transition to the zigzag magnetically ordered state at T N ∼ 7 K [9]. Several effective models for α-RuCl 3 have been proposed so far. Although some effective models show qualitatively different features with each other, they have succeeded in explaining experimentally observed thermodynamic quantities and/or low-lying excitations. The adequate effective model is still controversial.
Inelastic neutron scattering (INS) experiments for single crystals provide us with ample information on excitation spectra via dynamical spin structure factors. Quite recently, INS experiments for single crystal α-RuCl 3 were performed [5,10]. The INS experiments have revealed following features, which offer clues to determine the adequate effective model.
• The magnetic Bragg peaks stand at the M and Y points (see Fig. 2 (d)) in the reciprocal space in T < T N , reflecting the zigzag magnetic order [10]. This contributes to the 2 1234567890 ''"" characteristic six-pointed star-shaped intensity centered at the Γ point in constant-energy cuts in low temperatures [5,10]. • Around the Γ point, two peaks of the scattering intensity centered at ω ∼ 3 meV and ∼ 6 meV appear in T < T N [10]. The scattering intensity extends to ω ∼ 15 meV [10]. • In T > T N , the scattering intensity around the Γ point has a broad maximum centered at ω ∼ 6 meV [10]. The Γ-point intensity survives up to T ∼ 120 K [10].
In this study, we calculate the dynamical spin structure factor (DSF) of three effective models [11,12,13] with an exact numerical diagonalization method. These effective models are obtained from ab-initio calculations. We call these three ab-initio models MODELs 1-3. We also calculate the temperature dependence of the specific heat by employing thermal pure quantum states [14,15,16]. By comparing the experimental results with our numerical results, we discuss characteristics of the three ab-initio models.

Model and methods
We consider a generalized Kitaev-Heisenberg model on a honeycomb lattice to describe MODELs 1-3 [11,12,13]. The Hamiltonian is written by whereĴ µν Γp is a 3 × 3 matrix expressing the exchange coupling between p-th neighboring i and j sites on the bond Γ p . For instance, the matrix elements of the nearest-neighbor pairs on a Z bond (see Fig. 1 where J µ p and K µ p denote coupling constants of the Heisenberg-type and the Kitaev-type interactions, respectively. We summarize the parameters of the three ab-initio models in table 1.
In the three models we calculate the DSF, which is defined as where Z(T ) is a partition function, ϕ n denotes an eigenstate of H with the eigenvalue E n , andŜ eigenstates and eigenvalues to calculate the DFSs at finite temperatures, the system size is limit to a small cluster. In this study, we calculate the DSF at finite temperatures for the N =8 clusters shown in Fig. 1 the energy E 0 . ϕ 0 and E 0 are calculated by the Lanczos method, and then S µν (Q, ω; 0) is obtained by a continued fraction expansion [17]. We calculate the DSF at zero temperature for the N =24 clusters shown in Fig. 1 (c). In a similar manner to the INS experiments [5,10], we evaluate the sum of the diagonal elements and set S(Q, ω; T ) ≡ ∑ µ=x,y,z S µµ (Q, ω; T ). The specific heat is calculated by employing thermal pure quantum states [15,16] for the N = 24 clusters. The construction of thermal pure quantum states is summarized in Ref. [8].  , ω; 0). Solid curves are dispersion curves calculated from linearized spin-wave theory. In the calculations with linearized spin-wave theory, we start from the Ising-like zigzag state. Note that the ground state of the spherical model sometimes happens to yield a different ordered state from that obtained by the numerical exact diagonalization. Horizontal axis in (a)-(c) runs along arrows shown in Fig. 2 (d).
We first show the DSF at T = 0 in Fig. 2. The area of circle represents log S(Q, ω; 0). The solid curves are dispersion curves calculated from the linearized spin-wave theory by assuming  the Ising-like zigzag state. We find that in MODELs 1-3 the linearlized spin-wave theory fails in explaining the DSFs. In MODELs 1 and 2, the largest intensities appear at the M and Y points at ω ∼ 2 meV, which stem from the zigzag magnetic order in agreement with the experiments [5,10]. On the other hand, in MODEL 3, the largest intensities appear neither at the M nor Y points, but appear at the K points at ω ∼ 0.3 meV. This result disagrees with the experiments. In MODEL 3, the Heisenberg-type long-range interactions J 2 and J 3 were introduced so as to reproduce the magnetization curves in the zigzag ordered state [13]. Another fitting of J 2 and J 3 that can induce strong intensities at the M and Y points may be effective. In Fig. 3, we show the DSF at the Γ point, S(0, ω; 0). In MODEL 1, lower two peaks appear at ω ∼ 4 meV and ∼ 6 meV, which seem to describe the experimental results at the Γ point in T < T N . In MODEL 2, only a single peak appears at ω ∼ 7.5 meV in this energy scale, which deviates from the experiments. In MODEL 3, it is difficult to see the notable lower two peaks in S(0, ω; 0). There are multipeaks in 3 meV < ω < 9 meV, which we consider to form an excitation continuum in the thermodynamic limit. We show the temperature dependences of the DSFs at the Γ point for the N = 8 clusters in Fig. 4. In this calculation, it is difficult to get a good resolution with respect to ω because of the small-cluster calculations and the broad half width of used Lorenzians. Thus, in MODEL 1, the lower two peaks does not appear in S(0, ω; T ), but a broad maximum appears centered at ω ∼ 7 meV. In MODEL 2, S(0, ω; T ) has a broad maximum centered at ω ∼ 8 meV, which follows the peak in S(0, ω; 0) shown in Fig. 3 (b). In MODEL 3, the peak intensity of S(0, ω; T ) is smaller by approximately one order than that in MODEL 1. This is because, in MODEL 3, only the excitation continuum appears at the Γ point.  [15,16]. All results are drawn with error bars that are smaller than symbol size.

Temperature dependence of S(0, ω; T )
We show the specific heat C(T ) for the N = 24 clusters obtained by using the thermal pure quantum states in Fig. 5. The double-peak structure appears in MODELs 1-3. We evaluate the temperature ratio of the low-temperature peak (T l ) to the high-temperature one (T h ). We obtain T l /T h ∼ 0.12 in MODEL 1, T l /T h ∼ 0.20 in MODEL 2, and T l /T h ∼ 0.12 in MODEL 3. According to the criterion that has been elucidated in Ref. [8], MODELs 1 and 3 are located closer to the KSL phase than MODEL 2. Note that the ratios obtained here are larger than the experimentally observed ratios of α-RuCl 3 : T l /T h ∼ 0.09 [3,4] and ∼ 0.065 [5].
In MODELs 1-3, T h 's are evaluated to be T h ∼ 30 K, ∼ 50 K, and ∼ 20 K, respectively. We compare these values with the temperature regions where the large intensities of S(0, ω; T ) appear in Fig. 4. We find that the large intensity of S(0, ω; T ) appears in T < T h in each MODEL. It has been shown that the nearest-neighbor spin correlation begins to develop at T ∼ T h , as T is lowered [7,8]. The development of the nearest-neighbor spin correlation is a hallmark of the KSL [7,8,18]. Thus the large intensity of the DSF around the Γ point is associated with the Kitaev interaction, as has been pointed out in Refs. [5,10].
The values of T h in MODELs 1-3 shown above are lower than the experimental results: T h ∼ 85 K [3,4] and ∼ 100 K [5]. From the numerical calculations, we find that T h is mainly controlled by the dominant Kitaev interaction. These findings probably mean that the nearestneighbor ferromagnetic Kitaev interaction is underestimated in these ab-initio calculations for explaining T h observed in the experiments.

Summary
We have calculated the DSF and C(T ) in the three ab-initio models. In Refs. [11,13]  the DSFs and C(T ) by using such parameter sets, and have employed the results that describe the experiments best in each model. According to our calculations, the INS experimental features are qualitatively described by MODEL 1, while the double-peak structure of C(T ) is qualitatively reproduced by MODELs 1-3. By comparing the numerical results for MODELs 1-3, we consider that the growths of the long-range spin correlation in MODELs 1 and 3 are suppressed as compared with that in MODEL 2 because of the absence of and the weak third neighbor Heisenberg interaction, respectively. This suppression makes T l /T h of MODELs 1 and 3 small. In the DSF of MODEL 3, the excitation continuum appears in the low-energy region and it is difficult to expect the appearance of the localized mode observed in the experiments at ω ∼ 3 meV. The difference between MODEL 1 and MODEL 3 comes from the signs of the nearest-neighbor Heisenberg terms and the Γ terms. We consider that the signs of both terms affect the stability of the low-lying excitation mode. Thus, we conclude that MODEL 1 well describes the qualitative experimental features of α-RuCl 3 that we have turned our attention to. We believe that it is effective to construct a phenomenological model for α-RuCl 3 that can reproduce the experimental features of the DSF and C(T ) not only qualitatively but also quantitatively.