The field-angle anisotropy of proximate Kitaev systems under an in-plane magnetic field

We have investigated the field-angle behaviors of magnetic excitations under an in-plane magnetic field for proximate Kitaev systems. By employing the exact diagonalization method in conjunction with the linear spin wave theory, we have demonstrated that the magnetic excitation gap in the polarized phase is determined by the magnon excitation at $M$ points and has a strong anisotropy with respect to the field direction in the vicinity of the critical field limit. The specific heat from this magnon excitation bears qualitatively the same anisotropic behaviors as expected one for the non-Abelian spin liquid phase in the Kitaev model and experimentally observed one of the intermediate phases in $\alpha$-RuCl$_3$.

We have investigated the field-angle behaviors of magnetic excitations under an in-plane magnetic field for proximate Kitaev systems. By employing the exact diagonalization method in conjunction with the linear spin wave theory, we have demonstrated that the magnetic excitation gap in the polarized phase is determined by the magnon excitation at M points and has a strong anisotropy with respect to the field direction in the vicinity of the critical field limit. The specific heat from this magnon excitation bears qualitatively the same anisotropic behaviors as expected one for the non-Abelian spin liquid phase in the Kitaev model and experimentally observed one of the intermediate phases in α-RuCl3.

I. INTRODUCTION
Quantum fluctuation in frustrated spin systems can prevent any classical magnetic orders and induce exotic quantum phases such as a quantum spin liquid (QSL). The Kitaev model, the ideal S = 1 2 quantum spin system with bond-directional Ising interactions in a honeycomb lattice (see Fig. 1(a)), is one of intriguing systems which host the QSL phase as a ground state. In this exactly solvable model, spin dynamics can be interpreted as free Majorana fermions in a static Z 2 flux 1 . Majorana fermions with a gapless energy spectrum leads to the ground QSL phase and the fractionalized magnetic excitation.
Majorana fermions in the Kitaev model acquire a mass gap ∆ under the magnetic field. The gapped spectrum stabilizes the topological non-Abelian spin liquid (NASL) phase characterized by the Chern number of C = ±1 and protected chiral edge modes 1,2 . Because the mass gap is proportional to the multiplication of three components of the magnetic field with respect to local spin coordinate axes, the sign of C and ∆ shows strong field-angle dependency. For the in-plane field, the sign change and gap closing occur for the field applied to the bond direction between nearest neighboring (NN) spins as shown in Fig. 1(b) [3][4][5] . The thermal Hall conductivity, specific heat, and magnetotropic coefficient emulate characteristic features of the Majorana energy spectrum. The fieldangle behaviors of such quantities have been known as the key feature for the experimental identification of the NALS phase [3][4][5][6] .
A lot of theoretical and experimental studies have been performed to find realistic Kitaev materials [7][8][9] . Among the candidates, α-RuCl 3 has been turned out to be the best proximate Kitaev system with dominant Kitaev interaction 10,11 . The ground state is not the QSL but the antiferromagnetic zigzag order due to the additional non-Kitaev interactions such as the Heisenberg interactions and two types of off-diagonal exchange interactions, * bomisu@kias.re.kr called as Γ and Γ ′ terms [12][13][14] (see the general Hamiltonian for the proximate Kitaev system in Eq. 1). Observed magnetic continuum excitations and two-step magnetic entropy release have been regarded as the evidence of the fractionlized Majorana fermions [15][16][17][18][19][20] . Moreover, recent studies 2,3,6,21-38 have revealed that the zigzag order is easily destroyed when the external magnetic field is applied and the intermediate phase (IP), putative QSL, emerges between the zigzag-order phase and polarized phase. The NASL phase has been highly anticipated as the IP candidate because recent thermal Hall conductivity experiments detected the half-integer plateau and its sign signature under both in-plane and out-of-plane magnetic fields 2,3 . Further, specific heat measurement also probed the expected field-angle anisotropy under an in-plane field 6 . On the other side, strong sample dependency of the measured quantities has been reported 39 . A very recent thermal conductivity experiment has reported the emergence of different types of QSL phase 40 , The connection of three types of bond-dependent Ising interactions, and coordinate axes of local spins (x, y, z) and global honeycomb lattice (a, b, c) in the Kitaev model. The Ising directions in the x-, y-, and z-type bonds drawn with green, blue, and red lines are parallel to the x, y, and z coordinate axes, respectively. The anticipated field-angle variations of (b) Majorana gap of the Kitaev model and (c) magnon gap in the polarized phase of proximate Kitaev systems under an in-plane field. ϕ is the angle between the inplane field and the a axis. The radius of curves refers to the size of gap. The dark violet (sea green) curve denotes the Chern number C = +1 (−1). When the gap is determined at the M1, M2, and M3 points (see the Brillouin zone in inset), the curve is drawn with green, red, and blue colors in (c), respectively. and the detailed thermodynamics study 41 has purported the possibility of the absence of the IP itself. All in all, the nature of the IP and/or its existence are still under debate.
In this study, we investigated the field-angle dependence of magnetic excitation and magnetic specific heat for proximate Kitaev systems under an in-plane magnetic field. Using the exact diagonalization (ED) method and linear spin wave theory (LSWT), we demonstrated that the low-energy excitation features in the polarized phase for various models relevant to α-RuCl 3 can be interpreted in terms of the field-angle anisotropy of magnon gap, determined at one of three M points depending on the field direction (see Fig. 1(c)). The magnetic specific heat dictated by this magnon dynamics shows qualitatively the same anisotropic behaviors as those in the Kitaev model. Our result suggests that the anisotropic behaviors in thermodynamic quantities alone are not a smoking gun of the intermediate NASL phase in α-RuCl 3 under the magnetic field, but requires further considerations.

II. SPIN HAMILTONIAN
Let S A i and S B i be two base spins at the i-th unit cell in a honeycomb lattice. The general spin Hamiltonian of proximate Kitaev systems with first, second, and third NN interactions and external magnetic field are given as following: where γ 1 , γ 2 , and γ 3 represent the bond types of first, second, and third NN spins, respectively. i and i γn are the unit cell indices of two spins in the γ n bond (see Fig. 2(a)). γ (= x, y, z) can be characterized by three coordinate axes of spins.γ 2 refers to the bond with the opposite direction of the γ 2 bond. g is the g factor of spins and µ B is the Bohr magneton.J γn is the superexchange dyadic tensor of the γ-type n-th NN bond defined as where α, β, and γ are cyclically ordered coordinate axes of local spins andγ is the unit vector along the γ axis. J n , K n , Γ n , and Γ ′ n are the exchange parameters of the Heisenberg interaction, Kitaev interaction, and two types of off-diagonal interactions between n-th NN spins, respectively. The global coordinate axes a, b, and c can be determined so that the a (b) axis is parallel (perpendicular) to the displacement between two spins in the z 1 -type  [45] bond and the c axis is normal to the honeycomb lattice (See Fig. 1(a)). Unit vectorsâ,b, andĉ are given aŝ in terms of local spin coordinate axes. In the text, the exchange parameters between first NN spins are simply denoted as J, K, Γ, and Γ ′ omitting neighbor indices.
To investigate the field-angle anisotropy of proximate Kitaev model under an in-plane magnetic field, we first adopted the simple K-Γ-J 3 model with J 3 /|K| = Γ/|K| = 0.1 and K < 0, which also accounts for the magnetic phase of the system well. Then we extended our study to more realistic models presented in Table I.

III. ED CALCULATION
Employing the ED method with the periodic 24-site cluster as shown in Fig. 2(a), we investigated the magnetic phase transition of the K-Γ-J 3 model with J 3 /|K| = Γ/|K| = 0.1 and K < 0 under an in-plane field. With the help of the thick-restarted Lanczos method 49 , we obtained the ground state |Ψ g and its energy E g . Figure 2(b) and (c) present the static spin structure factors (SSSFs) Ψ g | S −q · S q |Ψ g at q = Γ, M 1 , and M 2 points (see the Brillouin zone in Fig. 1) for the a-and b-axis fields. S q is defined as S q = 1 √ N N j=1 e −iq·rj S j , where r j is the position vector of the j-th spin S j in the honeycomb lattice and N is the total number of spin sites. Note that the SSSFs at Γ and three M points characterize the polarized phase and three types of the zigzag order, respectively. Results indicate that the magnetic phase transition happens from the zigzag-order phase to polarized phase at around gµ B h/|K| = 0.1 without hosting any IP under both a-and b-axis fields. In the ED calculation, the IP is only feasible for an out-of-plane  field. For the c-axis field, the IP appears in the range of 0.396 ≤ gµ B h/|K| ≤ 0.416 (not shown here).
Further, we numerically explored the magnetic excitation features as a function of strength and direction of an in-plane magnetic field by calculating the dynamical spin structure factor (DSSF) χ q (ω) as following: where S q,ν is the ν (= x, y, z) component of S q and δ (= 0.01|K|) is the broadening parameter. When a magnetic field is applied along the a axis (ϕ = 0 • , where ϕ is the angle between the in-plane field and the a axis), the minimum excitation energies at both the Γ and M 1 points decrease in a weak field limit but they start to increase at around the critical field corresponding to the zigzag order to the polarized phase transition. The excitation gap is determined at the M 1 point regardless of the field strength. In contrast, the minimum excitation energy at the M 2 point, originally degenerate with those at other M points without the field, monotonically increases with losing its spectral weight when the field strength increases (see Fig. 2(d), (e), and (f)).
The 24-site cluster has the dihedral D 3 symmetry which includes both C 2 rotation along three NN bond directions and C 3 rotation along the c axis. Due to the C 2 rotation symmetry, the excitation spectra at three M points have the 180 • periodicity for the field angle ϕ. In addition, the excitation spectra at three M points have a cyclic relation under the C 3 rotation. These symmetric features are well captured in the polarized phase region of gµ B h/|K| = 0.12 as shown in Fig. 2(g) and (h). The DSSFs χ M1 , χ M2 , and χ M3 have the 180 • periodicity, and the excitation gap is determined from χ M1 , χ M2 , and χ M3 when the field angle ϕ is located at 0 • ≤ ϕ ≤ 60 • , 60 • ≤ ϕ ≤ 120 • , and 120 • ≤ ϕ ≤ 180 • , respectively. The minima of the excitation gap appear whenever a magnetic field is parallel to one out of three NN bond directions (ϕ = (2n + 1) × 30 • where n is the integer number).
As shown in Fig. 2(g) and (h), the excitation gap at ϕ = 0 • is higher than that at ϕ = 90 • , which means that the magnetic entropy under the b-axis field (ϕ = 90 • ) is easier to be thermally populated in a low-temperature limit than the a-axis-field case (ϕ = 0 • ). Hence, the magnetic specific heat C m has a lower gap for the field along the b axis than the a axis, which is well captured in Fig. 3(a). Here we employed the K-Γ-J 3 model with gµ B h/|K| = 0.12 (see Appendix A for the calculation detail). Also, we found the six-fold symmetricity of C m upon the field angle ϕ for the low-temperature case (k B T /|K| 0.01) as shown in Fig. 3(b). C m periodi- cally varies with minimum values at ϕ = n × 60 • and maximum values at ϕ = (2n + 1) × 30 • . The anisotropic behaviors are suppressed when the field strength is increased from gµ B h/|K| = 0.12 to 0.2, which is consistent with the experimental observations of α-RuCl 3 6 .

IV. SPIN WAVE THEORY
To get further insight on the field-angle anisotropy under an in-plane magnetic field, we examined the excitation features in terms of the LSWT (see the detail in Appendix B). We performed the LSWT calculation within the polarized phase in which all spins are ferromagnetically ordered along the field direction. By reducing the field strength from the infinity, we traced the variation of magnon dispersions in the K-Γ-J 3 model.
In the classical magnetic phase diagram for J 3 /|K| = Γ/|K| = 0.1, the polarized phase is stabilized when gµ B h/|K| is larger than about 0.295 (0.305) for the a-axis (b-axis) field as shown in the Fig. 4 (see the calculation detail in Appendix C). The calculated magnon bands are always gapped in this field limit. The gap is monotonically diminished as the field strength is reduced. Eventually, the magnon bands condense at the M 2,3 points (M 1 point) under the a-axis (b-axis) field with the critical field strength of gµ B h/|K| = 0.2930 (0.3048) (see Fig. 4 and Fig. 5(a)). For the field-angle dependency, as analogous to the previous DSSFs, the magnon bands have the D 3 character in the critical field limit: the lowlying excitation gap at 0 • ≤ ϕ ≤ 60 • , 60 • ≤ ϕ ≤ 120 • , and 120 • ≤ ϕ ≤ 180 • , is determined at M 1 , M 2 , and M 3 points, respectively. The gap minima are located at ϕ = (2n + 1) × 30 • as shown in Fig. 5(b).
Based on the magnon dispersions, we calculated the magnon specific heat C m which is attributed to onemagnon excitation. Figure 5 as a function of temperature T under both a-and b-axis fields at the critical field strength of gµ B h/|K| = 0.3048. For the b-axis field, the C m is gapless since the magnon is gappless at the critical strength, while it is still gapped for the a-axis field. The six-fold symmetric behaviors of C m , at low T , can be also observed such that minimum and maximum values at ϕ = n × 60 • and (2n + 1) × 30 • as shown in Fig. 5(d). The anisotropic behaviors are progressively suppressed as the field strength is increased beyond the critical value. All behaviors are consistent with both results of ED calculation and experimental observations of α-RuCl 3 6 .

V. APPLICATION TO α-RUCL3
We extended our study to a few representative models (models 1 to 7), which have been proposed for the magnetic and thermal properties of α-RuCl 3 , as shown in Table I.
We investigated their magnetic phase transition under an in-plane magnetic field. As shown in Fig. 6, the second derivative of their ground state energy with respect to the field strength (−d 2 E g /dh 2 ) shows single peak in the ED calculation. All models seem to exhibit the phase transition from the zigzag order to the polarized phase without hosting any IP under both a-and b-axis fields like the K-Γ-J 3 model. However, the ED calculation with small-size cluster is limited due to the finite size effect. The absence of IP under an in-plane field is still questionable.
We examined the behaviors of DSSFs, magnon dispersions, and magnon specific heat under both a-and b-axis fields for models 1-7. Results are presented in Fig. 7 and 8. Overall, we found all considered models bear the essentially same field-angle anisotropy. More specifically, the DSSF, magnon condensation, and magnon specific heat behaviors under an in-plane magnetic field are the same between models 1-4 and the K-Γ-J 3 model (see   Table I under (a) a-and (b) b-axis fields (ϕ = 0 • and ϕ = 90 • ). The ground state energy is obtained with the ED calculation of the periodic 24-site cluster shown in Fig. 2(a).   Fig. 7). Small differences in DSSFs and/or magnon condensation are found for models 5-7: (i) DSSFs χ M1 and χ M2 do not determine the excitation gap for a-and baxis fields for model 7. (ii) The magnon condensation point is slightly shifted from the M 2 to Γ points under the b-field for models 5 and 6. Still, the quantitative features remain very similar (see Fig. 8). Therefore, our results suggest the potential role of magnon dynamics in explaining the field-angle anisotropy of low-temperature thermal properties for α-RuCl 3 .
Previously, we proposed the K-J 3 model with K < 0 and J 3 > 0 as the simplest model to account for the IP for the a-axis field 38 , where, we found the magnon simultaneously condenses at both the M 1 and M 2 points with the suppression of the anisotropy in the magnon specific heat C m (T ) (see Appendix D). Hence, we conjecture that the Γ-term plays a significant role in the inplane anisotropic feature in the proximate Kitaev model with ferromagnetic K (K < 0).

VI. DISCUSSION
In contrast with the ED calculations of various models to show no IP, plausible evidence of intermediate NASL phase such as the half plateau thermal Hall conductivity and field-angle anisotropy of specific heat has been reported in α-RuCl 3 in the intermediate range of the a-axis field (7 ∼ 10 T) 3,6 . The lower bound field is coincident with the critical field at which the zigzag order totally disappears. However, the upper phase boundary is not evident. The thermodynamic quantities such as magnetic susceptibility, specific heat, and magnetic Grüneisen parameter show no clear anomalies at around the upper field 32,35,41 . Moreover, the recent electron spin resonance (ESR) experiment unveiled that the single-magnon excitation is present across the upper field 36 . The magnon dynamics of polarized phase is certainly anticipated to emerge even for the intermediate region of α-RuCl 3 .
An important remark on the magnon dynamics of the polarized phase is that the excitation gap is determined at not the Γ point but M points in contrast with the NASL phase in which the magnetic excitation gap appears at the Γ point. The momentum of the excitation gap could be good measure the origin of the fieldangle anisotropy of specific heat. The full magnetic excitation features, which can be measured by the inelastic neutron scattering (INS) and resonant inelastic x-ray scattering experiments, are crucial for determining the nature/existence of IP of α-RuCl 3 . According to recent INS, Raman spectroscopy, and ESR experiments 18,36,37,53 , the minimum excitation energy at the Γ point evidently decreases before the phase boundary and increases again after the boundary. In the experimental resolution, however, it is hardly resolved whether the gap is genuinely closed or not. It was also verified that the excitation energy at the M 2 point in the INS is almost constant in the zigzag ordering limit and that it increases by losing its spectral weight after the phase boundary. This behavior is quite similar to our ED calculation. However, only magnetic excitation along the Γ-M 2 line under the a-axis field has been reported in the hitherto INS measurement. The momentum resolved excitations except for the Γ point are limited in Raman spectroscopy and ESR experiments. Experimental observations with different momentum lines and various field directions are highly required.
The recent theoretical study has pointed out that the magnon topology in the polarized phase of the proximate Kitaev system can give rise to the field-angle variation in the thermal Hall conductivity which has the same sign structure as that in the NASL phase 52 . Our result supports novel significant feature of this magnon dynamics. Thus, the field-angle anisotropy of both specific heat and thermal Hall conductivity cannot be taken as a key evidence for the NASL phase. For the ultimate identification, one should test more comprehensively a few characteristic features which are inherent to Majorana fermion dynamics, e.g., the continuum excitation at the Γ point, the half-integer plateau of thermal Hall conductivity under the a-axis field, and T 2 -behavior of the specific heat at low temperature under the b-axis field.

VII. CONCLUSION
Based on the numerical ED calculation and LSWT analysis, we have explored the field-angle anisotropy of proximate Kitaev systems under an in-plane magnetic field. We have found that the low-lying excitation gap, interpreted as the magnon excitation in the polarized phase, is determined at not the Γ but M points in the vicinity of the phase boundary between the zigzag order and polarized phase. Also, the excitation gap has the 60 • periodicity for the field angle with its minimum when the field is along the NN bond direction. The anisotropy of the low-energy magnon gap can reproduce the field-angle anisotropy of the specific heat, which was considered a hallmark of the NASL phase. Our results provide a novel We performed the finite-temperature Lanczos method (FTLM) 50,51 to calculate the magnetic specific heat of proximate Kitaev systems with the 24-site cluster. The specific heat at temperature T is given as where N is the total number of spin sites, H is the spin Hamiltonian, and k B is the Boltzmann constant, respectively. The expectation values of H and H 2 can be approximately estimated as following: where N st , N sc , and N L are the size of Hilbert space, the number of initial random states, and the number of Lanczos iteration steps, respectively. ε r,m refers to the m-th eigenvalue calculated by the Lanczos method with the r-th initial random state. Z is the partition function which is also approximately calculated as where |φ r is the normalized r-th initial random state and |ψ r,m is the m-th eigenstate calculated by the Lanczos method with the initial state |φ r . In the calculation, we set N sc = 250 and N L = 200. In the polarized phase, all spins are ferromagnetically ordered along the magnetic field direction. We performed the liner spin wave theory (LSWT) assuming that an inplane magnetic field is applied away from the a axis with the field angle ϕ. According to the Holstein-Primakoff transformation, the spin operators can be written in terms of two bosonic operators as following: where a † i (a i ) and b † i (b i ) are the bosonic creation (annihilation) operators of magnons at the i-th A and B sublattices, respectively, and x ′ , y ′ , and z ′ are coordinate axes defined asx where a, b, and c are global coordinate axes of the lattice. The magnetic interactions between first, second, and third NN spins in Eq. 1 can be approximately expressed in terms of the bosonic operators as following: and where HereJ αβ γn =α ·J γn ·β. We define the spinor operators can be given as following: where where r γn = −rγ n = r iγ n − r i . According to the Bogoliubov transformation, Equation B8 is extended as following: H ≈ −N u S [SE(h, ϕ) + ǫ(h, ϕ)] We obtained the dynamic equation of motion of spinor operators as follows: where ∆(ϕ, k) = 1 2 [d(ϕ, k) + d(ϕ, −k) ⊺ ]. By solving the general eigenvalue problem of Eq. B12, we calculated magnon dispersions. Let ǫ n (k) be the n-th magnon band at a given momentum k (n = 1, 2). The average energy per unit cell can be given as To determine the classical phase diagram of proximate Kitaev systems under an in-plane field, we performed the classical Monte Carlo (MC) calculation with the standard Metropolis algorithm. By considering a periodic 2 × 36 × 36 cluster, we ran the 40000 MC steps after 20000 MC steps for thermalization at k B T /|K| = 0.02 and calculated the expectation value of order parameters for both the zigzag order and the polarized phase as a function of the field strength.
Appendix D: K-J3 model The K-J 3 model with K < 0 and J 3 > 0 has been proposed as the simplest proximate Kitaev model which shows the genuine intermediate phase under both the aand c-axis fields 38 . We have explored the behaviors of DSSFs, magnon dispersions, and magnon specific heat of the K-J 3 model under in-plane magnetic fields. As shown in Fig. 9, the field-angle anisotropy in the K-J 3 model is less evident than those in other models. The DSSFs χ M1 , χ M2 , and χ M3 are almost constant in the range of 0 • ≤ ϕ ≤ 60 • , 60 • ≤ ϕ ≤ 120 • , and 120 • ≤ ϕ ≤ 180 • , respectively. The magnon bands in the polarized phase condense almost simultaneously at the three M points and the field-angle anisotropy of magnon specific heat is almost suppressed.