Identification of a Kitaev quantum spin liquid by magnetic field angle dependence

Quantum spin liquids realize massive entanglement and fractional quasiparticles from localized spins, proposed as an avenue for quantum science and technology. In particular, topological quantum computations are suggested in the non-abelian phase of Kitaev quantum spin liquid with Majorana fermions, and detection of Majorana fermions is one of the most outstanding problems in modern condensed matter physics. Here, we propose a concrete way to identify the non-abelian Kitaev quantum spin liquid by magnetic field angle dependence. Topologically protected critical lines exist on a plane of magnetic field angles, and their shapes are determined by microscopic spin interactions. A chirality operator plays a key role in demonstrating microscopic dependences of the critical lines. We also show that the chirality operator can be used to evaluate topological properties of the non-abelian Kitaev quantum spin liquid without relying on Majorana fermion descriptions. Experimental criteria for the non-abelian spin liquid state are provided for future experiments.

Recent advances in experiments have unveiled characteristics of QSLs. For α-RuCl 3 , signatures of Majorana fermion excitations have been observed in various different experiments of neutron scattering, nuclear magnetic resonance, specific heat, magnetic torque, and thermal conductivity [49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65] . Among them, the half quantization of thermal Hall conductivity κ xy =T ¼ ðπ=12Þðk 2 B =_Þ may be interpreted as the hallmark of the presence of Majorana fermions and the non-abelian KQSL 60,63 . At higher magnetic fields, a significant reduction of κ xy /T also suggests a topological phase transition 60,66 . Thermal Hall measurements are known to be not only highly sensitive to sample qualities 64 but also very challenging due to the required precision control of heat and magnetic torque from strong magnetic fields. This strongly motivates an independent way to detect the Majorana fermions and non-abelian KQSL.
In this work, we propose that the non-abelian KQSL may be identified by the angle dependent response of the system under applied magnetic fields. As a smoking gun signature of the KQSL, quantum critical lines are demonstrated to occur on a plane of magnetic field directions whose existence is protected by topological properties of the KQSL. The critical lines vary depending on the microscopic spin Hamiltonian, which we show by investigating a chirality operator via exact diagonalization. We further propose that the critical lines can be detected by heat capacity measurements and provide experimental criteria for the nonabelian KQSL applicable to the candidate material α-RuCl 3 .

Results
Model and symmetries. We consider a generic spin-1/2 model on the honeycomb lattice with edge-sharing octahedron crystal structure, 10,11,13,16 . Nearest neighbor bonds of the model are grouped into x, y, z-bonds depending on the bond direction (Fig. 1a). Spins (S j,k ) at each bond are coupled via the Kitaev (K), Heisenberg (J), and off-diagonal-symmetric (Γ; Γ 0 ) interactions. The index γ ∈ {x, y, z} denotes the type of bond, and the other two α, β are the remaining components in {x, y, z} other than γ. Under an applied magnetic field (h), the Hamiltonian becomes We specify the magnetic field direction with the polar and azimuthal angles (θ, ϕ) as defined in Fig. 1b. H KJΓΓ 0 possesses the symmetries of time reversal, spatial inversion, C 3 rotation about the normal axis to each hexagon plaquette, and C 2 rotation about each bond axis (Fig. 1a). The C 3 and C 2 rotations form a dihedral group D 3 . Under each of these symmetries, H(θ, ϕ) is transformed to Hðθ 0 ; ϕ 0 Þ with a rotated magnetic field hðθ 0 ; ϕ 0 Þ; see Supplementary Notes 1 and 2.
In the pure Kitaev model, parton approach provides the exact wave function of KQSL together with gapped Z 2 flux and gapless Majorana fermion excitations. Application of magnetic fields drives the KQSL into the non-abelian phase by opening an energy gap in Majorana fermion excitations. The gap size is proportional to the mass function, M(h) = h x h y h z /K 2 , and the topological invariant (Chern number) of the KQSL is given by the sign of the mass function, sgnðh x h y h z Þ 6 .
Topologically protected critical lines. Topological invariant of non-abelian phases with Majorana fermions can be defined from the quantized thermal Hall conductivity, κ xy =T ¼ νðπ=12Þðk 2 B =_Þ, where ν is the topological invariant representing the total number of chiral Majorana edge modes (T: temperature) 6 . While the topological invariant in the pure Kitaev model is exactly calculated by the Chern number of Majorana fermions, it is a nontrivial task to analyze the topological invariant for the generic model H(θ, ϕ).
Our strategy to overcome the difficulty is to exploit symmetry properties of the topological invariant and find characteristic features of the non-abelian KQSL. Concretely, we focus on the landscape of νðhÞ on the plane of the magnetic field angles (θ, ϕ). Our major finding is that critical lines of νðhÞ must arise as an intrinsic topological property of the non-abelian KQSL.
We first consider time reversal symmetry and note the following three facts: • Time reversal operation reverses the topological invariant as νðhÞ ! − νðhÞ.
These properties enforce the two regions to meet by hosting critical lines where Majorana fermion excitations become gapless. In other words, topological phase transitions must occur as the field direction changes. We propose that the very existence of critical lines can be used in experiments as an identifier for the KQSL.
We further utilize the D 3 symmetry of the system. The topological invariant ν and thermal Hall conductivity κ xy are A 2 representations of the D 3 group, i.e., even under C 3 rotations but odd under C 2 rotations, which reveals the generic form, where Λ 1,3 are field-independent coefficients. The h-linear term (h x + h y + h z ) and h-cubic term (h x h y h z ) are the leading order A 2 representations of magnetic fields. Conducting third order perturbation theory, we find the coefficients where Δ flux = 0.065|K| means the flux gap in the Kitaev limit. See Supplementary Notes 3-5 and ref. 67 for more details of the perturbation theory.
Notice that the h-linear term is completely absent in the pure Kitaev model (Λ 1 = 0). Figure 1c visualizes the topological invariant, νðhÞ ¼ sgnðh x h y h z Þ. The dashed lines highlight the critical lines representing the topological phase transitions between the phases with νðhÞ = ±1, where the energy gap of Majorana fermion excitations is closed: Exploiting the symmetry analysis, we stress two universal properties of the KQSL with the D 3 symmetry.
1. Symmetric zeroes: for bond direction fields, topological transition/gap closing is guaranteed to occur by the symmetry, i.e., Δ(h) = 0 for (θ = 90 ∘ , ϕ = 30 ∘ + n⋅60 ∘ ). 2. Cubic dependence: for in-plane fields, the h-cubic term governs low field behaviors of the KQSL, e.g., The universal properties and critical lines of the KQSL are numerically investigated for the generic Hamiltonian H(θ, ϕ) in the rest of the paper.
Chirality operator. We introduce the chirality operator at each hexagon plaquette p and investigate the expectation value of the chirality operator, shortly the chirality, and its sign, where Ψ KQSL ðhÞ is the ground state of the full Hamiltonian H(θ, ϕ) in the KQSL phase (N is the number of unit cells). The chirality operatorχ p produces the mass term of Majorana fermions and determines the topological invariant in the pure Kitaev limit 6 . More precisely, how magnetic fields couple to the chirality operator determines the topological invariant and the Majorana energy gap. The chirality χ and its sign ν are in the A 2 representation of the D 3 group as of the topological invariant ν. We note that the relation between the chirality and the Majorana energy gap, χ~Δ, holds near the symmetric zeroes in a generic KQSL beyond the pure Kitaev model due to the symmetry properties (Supplementary Note 6).
Exact diagonalization. The Hamiltonian H(θ, ϕ) is solved via exact diagonalization (ED) on a 24-site cluster with sixfold rotation symmetry and a periodic boundary condition (Fig. 1a). Resulting phase diagrams are provided in the section of Methods. The two universal features of the KQSL are well captured by the chirality (Figs. 2 and 3). Firstly, the zero lines of the chirality χ ED (h) always pass through the bond directions (marked by black dots), i.e., the symmetric zeroes. Secondly, the chirality shows the cubic dependence for in-plane fields (θ = 90 ∘ ). The linear term, h x + h y + h z , vanishes for in-plane fields, and the cubic term, h x h y h z , determines the chirality at low fields, which is confirmed in our ED calculations (lower panels of Fig. 2). Below, we show how non-Kitaev interactions affect topological properties of the non-abelian KQSL.
Most of all, we find that ν becomes identical to ν for the pure Kitaev model in Fig. 2a. It is remarkable that the two different methods, ED calculations of the chirality sign and the parton analysis, show the complete agreement: νðhÞ ¼ νðhÞ. The agreement indicates that the topological phase transitions can be identified by using the chirality operator, which becomes a sanity check of our strategy to employ the chirality operator. Figure 2b illustrates effects of the Heisenberg interaction (J) on the chirality. The shape of the critical lines is unaffected by the Heisenberg interaction, remaining the same as in the pure Kitaev  (3) and Supplementary Fig. 2b], indicating the validity of our strategy. Figure 2c, d presents consequences of the other non-Kitaev interactions, Γ 0 and Γ. The zero lines tend to be flatten around the equator θ = 90 ∘ , which can be attributed to the hlinear term induced by the non-Kitaev interactions: cos θ. In other words, the zero lines substantially deviate from those of the pure Kitaev model by the non-Kitaev couplings, Γ 0 and Γ. We point out that the signs of the chirality are opposite to the Chern numbers of the third order perturbation parton analysis ( Supplementary Fig. 2c, d). The opposite signs indicate that the two methods have their own valid conditions, calling for improved analysis (Supplementary Note 8).
Impacts of the non-Kitaev interactions also manifest in the field evolution of the zero lines (Fig. 3). Without the non-Kitaev couplings, Γ 0 and Γ, the shape of the critical lines is governed by the cubic term h x h y h z , as shown in Fig. 3a, b. In presence of the Γ 0 or Γ coupling, the h-cubic term competes with the h-linear term as illustrated in Fig. 3c, d. Namely, the linear term dominates over  Table 1. the cubic term at low fields while the dominance gets reversed at high fields (see Supplementary Note 11 and Supplementary  Table 4 for more results). The competing nature may be used to quantitatively characterize the non-Kitaev interactions.
Similarities and differences between the topological invariant, νðhÞ, and the sign of the chirality, νðhÞ, are emphasized. First, the two quantities are identical in the Kitaev limit while they can be generally different by non-Kitaev interactions. Second, the two quantities are in the same representation of the D 3 group, so νðhÞ and νðhÞ have in common the symmetric zeroes. Third, differences between the two quantities may be understood by considering other possible A 2 representation spin operators that  Table 1. may contribute to the topological invariant. For example, linear and higher-order spin operators exist in addition to the chiral operator. Since the topological invariant ν is related with the thermal Hall conductivity κ xy , the associated energy current operator directly informs us of relevant spin operators to the topological invariant. We find that linear spin operator is irrelevant to κ xy and ν (Supplementary Note 7). We also evaluate the expectation values of higher-order spin operators for the KQSL and confirm that their sizes are substantially small compared to the chirality (Supplementary Note 8). Therefore, we argue that the critical lines of the non-abelian KQSL are mainly determined by the zero lines of the chirality.

Discussion
Intrinsic topological properties of the non-abelian KQSL including the critical lines, the symmetric zeroes, and the cubic dependence are highlighted in this work by exploiting the symmetries of time reversal and D 3 point group. The chirality operator is used to evaluate the topological properties of the KQSL via the ED calculations. Now we discuss how the properties are affected by lattice symmetry breaking such as stacking faults in real materials. First, the existence of the critical lines relies on time reversal, thus it is not destroyed by lattice symmetry breaking. The symmetric zeroes appearing at the bond directions, however, are a consequence of the D 3 symmetry. The locations of the zeroes get shifted upon breaking the symmetry, which is confirmed in ED calculations of the chirality.
The cubic dependence for in-plane fields also originates from the D 3 symmetry and topology in the KQSL. The characteristic nonlinear response is not expected in magnetically ordered phases, which we check by performing spin wave calculations. Figure 4 contrasts the KQSL with magnetically ordered phases in terms of excitation energy gap (Majorana gap vs. magnon gap). The magnetic phases show completely different behaviors from the h-cubic dependence. Hence, the characteristic cubic dependence under in-plane magnetic fields may serve as an experimentally measurable signature of the KQSL.
The universal properties of the KQSL can be observed by heat capacity experiments. Figure 5a illustrates the calculated specific heat c v (ϕ) for the KQSL as a function of in-plane field angle ϕ (where magnetic field is rotated within the honeycomb plane). The specific heat is maximized by gapless continuum of excitations when the magnetic field is aligned to the bond directions. For comparison, the zigzag state, observed in α-RuCl 3 at zero field, is investigated by using a spin wave theory. The magnon spectrum is gapped due to completely broken spin rotation symmetry, so there is no critical line on the (θ, ϕ) plane (Fig. 5b). Compared to the KQSL, the zigzag state exhibits reverted patterns  of ϕ dependence in the excitation energy gap and specific heat. The energy gap is maximized and the specific heat is minimized at the bond directions. This behavior is closely related with the structure of spin configuration: all spin moments are aligned perpendicular to a certain bond direction selected by magnetic field direction (Supplementary Note 9). The distinct patterns of ϕ dependence in Fig. 5 characterize differences between the nonabelian KQSL and zigzag states. Remarkably, such behaviors were observed in the recent heat capacity experiments with in-plane magnetic fields 65 . Covering the polar angle (θ) in the heat capacity measurements will provide more detailed information on the critical lines and spin interactions in α-RuCl 3 (see Fig. 1d).
Lastly, we have examined the chirality and critical behaviors of excitation energy gap for magnetically ordered phases of H(θ, ϕ). It is found that the associated magnon gap does not have any critical lines, and there is no resemblance/correlation between the magnon gap and the chirality (Supplementary Note 9).
To summarize, we have uncovered characteristics of the nonabelian Kitaev quantum spin liquid, including the topologically protected critical lines, the symmetric zeroes, and the h-cubic dependence for in-plane fields, by using ED calculations with the chirality operator. Furthermore, we characterize the topological fingerprints of the KQSL in heat capacity. We expect our findings to be useful guides for identifying the KQSL in candidate materials such as α-RuCl 3 . Investigation of the universal properties in field angle dependence of thermodynamic quantities such as spin susceptibility is highly desired, and it would be also useful to apply our results to the recently studied field angle dependence of thermodynamic quantities 61,65,[68][69][70][71] .

Methods
Exact diagonalization. The KQSL and other magnetic phases of H(θ, ϕ) are mapped out by using the flux operatorŴ p ¼ 2 6 S z 1 S y 2 S x 3 S z 4 S y 5 S x 6 , the second derivative of the ground state energy ∂ 2 E gs =∂ξ 2 ðξ ¼ J; Γ; Γ 0 Þ, and the spin structure factor SðqÞ ¼ 1 N ∑ i;j hS i Á S j ie iqÁðr i Àr j Þ . We find that the KQSL differently responds to non-Kitaev interactions depending on the sign of the Kitaev interaction. Furthermore, the non-abelian phase of the KQSL is ensured by checking topological degeneracy and modular S matrix 6,38,72 . Figures 6 and 7 display the phase diagrams of H KJΓΓ 0 . A different structure of phase diagram is found depending on the sign of the Kitaev interaction. With the ferromagnetic Kitaev coupling (K < 0 as in Fig. 6), the KQSL takes an elongated region along the J axis but substantially narrowed along the Γ axis, showing the sensitivity to the Γ coupling of the ferromagnetic KQSL. Crossover-type continuous transitions are mostly observed among the KQSL and nearby magnetically ordered states such as ferromagnetic, stripy and vortex states 10,20,43 . Nature of the phase Y is unclarified within the finite size calculation. Unlike the aforementioned magnetically ordered states (Fig. 6d, e), the phase Y does not exhibit sharp peaks and periodicity in the spin structure factor, from which the phase is speculated to have an incommensurate spiral order or no magnetic order. It is remarkable that another quantum spin liquid phase, characterized by negative hŴ p i, exists in a broad region of the phase diagram (blue region of Fig. 6a) 39 . The QSL and KQSL show similarity in the spin structure factor (Fig. 6b, c). Nonetheless, the QSL as well as the phase Y get suppressed when the sign of Γ 0 is changed to negative. A zigzag antiferromagnetic order instead sets in under negative sign of Γ 0 ( Supplementary  Fig. 6).
In case of the antiferromagnetic Kitaev coupling (K > 0 as in Fig. 7), the KQSL is found to be more sensitive to the Heisenberg coupling rather than the Γ coupling, and surrounded by magnetically ordered states such as the vortex, zigzag, and Neel states 10,20,43 . In contrast to the ferromagnetic KQSL case, phase transitions between the antiferromagnetic KQSL and adjacent ordered states are all discontinuous as shown by hŴ p i in Fig. 7a. We also find that the antiferromagnetic KQSL and ferromagnetic KQSL are distinguished by different patterns of spin structure factor (Figs. 6b and 7b). Further phase diagrams for other values of Γ 0 are provided in Supplementary Fig. 6.
We also examine the phase diagrams at weak magnetic fields, and confirm that the overall structures remain the same as the zero-field results. We find that the chirality is useful for the identification of distinct phase boundaries. In some cases, the chirality performs better than the conventionally used flux (Supplementary Note 12).
We ensure the non-abelian KQSL phase by checking the Ising anyon topological order via threefold topological degeneracy 6,73 and modular S matrix 38,72,74 is obtained for the parameter set #4 in Table 1 with the magnetic field fixed along the [111] direction (θ = 0 ∘ ). See Supplementary Note 10 for the topological degeneracy and modular matrix computation.

Data availability
The data that support the findings of this study are available from the corresponding author on reasonable request.

Code availability
The code used to generate the data in this study is available from the corresponding author upon reasonable request.