Introduction

A quantum spin liquid (QSL) is an exotic state of matter characterized by many-body quantum entanglement1,2,3. In contrast to weakly entangled magnetic states, QSLs host emergent fractionalized quasiparticles described by bosonic/fermionic spinons and gauge fields4,5. The exactly solvable honeycomb model by Kitaev reveals the exact ground and excited states featured with Majorana fermions and \({{\mathbb{Z}}}_{2}\) gauge fluxes, so-called Kitaev quantum spin liquid (KQSL)6. Strong spin-orbit coupled systems with 4d and 5d atoms such as α-RuCl3 are proposed to realize KQSL7,8,9,10,11,12,13,14,15,16, and related spin models have been studied intensively17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48.

Recent advances in experiments have unveiled characteristics of QSLs. For α-RuCl3, signatures of Majorana fermion excitations have been observed in various different experiments of neutron scattering, nuclear magnetic resonance, specific heat, magnetic torque, and thermal conductivity49,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 \({\kappa }_{xy}/T=(\pi /12)({k}_{B}^{2}/\hslash )\) may be interpreted as the hallmark of the presence of Majorana fermions and the non-abelian KQSL60,63. At higher magnetic fields, a significant reduction of κxy/T also suggests a topological phase transition60,66. Thermal Hall measurements are known to be not only highly sensitive to sample qualities64 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 non-abelian KQSL applicable to the candidate material α-RuCl3.

Results

Model and symmetries

We consider a generic spin-1/2 model on the honeycomb lattice with edge-sharing octahedron crystal structure,

$${H}_{KJ\varGamma {\varGamma }^{\prime}}= \, \mathop{\sum}\limits_{{\langle jk\rangle }_{\gamma }}\left[K{S}_{j}^{\gamma }{S}_{k}^{\gamma }+J{{{{{{{{\bf{S}}}}}}}}}_{j}\cdot {{{{{{{{\bf{S}}}}}}}}}_{k}+\varGamma \left({S}_{j}^{\alpha }{S}_{k}^{\beta }+{S}_{j}^{\beta }{S}_{k}^{\alpha }\right)\right.\\ \, +\, \left.{\varGamma }^{\prime}\left({S}_{j}^{\alpha }{S}_{k}^{\gamma }+{S}_{j}^{\gamma }{S}_{k}^{\alpha }+{S}_{j}^{\beta }{S}_{k}^{\gamma }+{S}_{j}^{\gamma }{S}_{k}^{\beta }\right)\right],$$

so-called K-J-Γ-\(\varGamma ^{\prime}\) model10,11,13,16. Nearest neighbor bonds of the model are grouped into x, y, z-bonds depending on the bond direction (Fig. 1a). Spins (Sj,k) at each bond are coupled via the Kitaev (K), Heisenberg (J), and off-diagonal-symmetric (\(\varGamma ,\varGamma ^{\prime}\)) 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

$$H(\theta ,\phi )={H}_{KJ\varGamma {\varGamma }^{\prime}}-{{{{{{{\bf{h}}}}}}}}(\theta ,\phi )\cdot \mathop{\sum}\limits_{j}{{{{{{{{\bf{S}}}}}}}}}_{j}.$$
(1)

We specify the magnetic field direction with the polar and azimuthal angles (θ, ϕ) as defined in Fig. 1b. \({H}_{KJ\varGamma \varGamma ^{\prime} }\) possesses the symmetries of time reversal, spatial inversion, C3 rotation about the normal axis to each hexagon plaquette, and C2 rotation about each bond axis (Fig. 1a). The C3 and C2 rotations form a dihedral group D3. Under each of these symmetries, H(θ, ϕ) is transformed to \(H(\theta ^{\prime} ,\phi ^{\prime} )\) with a rotated magnetic field \({{{{{{{\bf{h}}}}}}}}(\theta ^{\prime} ,\phi ^{\prime} )\); see Supplementary Notes 1 and 2.

Fig. 1: Field angle dependence of the pure Kitaev model.
figure 1

a Honeycomb lattice enclosed by edge-sharing octahedra. Red, green, blue lines denote the x,y,z-bonds, and the six numbers indicate the numbering convention for sites in each hexagon plaquette (p). Black arrows depict C3 and C2 rotation axes. b Convention for the angular representation of an applied magnetic field h (yellow). The polar and azimuthal angles (θ, ϕ) are measured from the out-of-plane [111] axis and the bond-perpendicular \([11\bar{2}]\) axis, respectively. c Color map of the mass function M(h) = hxhyhz/K2 on the (θ, ϕ) plane. The dashed lines highlight the topological phase transitions between the ν = ±1 states, i.e., the quantum critical lines of the energy gap Δ(h) ∝ ∣M(h)∣ = 0. The black dots mark the bond directions (θ = 90∘, ϕ = 30∘ + n⋅60∘), where n = 0, 1, ..., 5. d A schematic of the field angle dependence of the heat capacity cv at a fixed temperature, where peak positions determine the critical lines of the non-abelian KQSL.

Full size image

In the pure Kitaev model, parton approach provides the exact wave function of KQSL together with gapped \({{\mathbb{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) = hxhyhz/K2, and the topological invariant (Chern number) of the KQSL is given by the sign of the mass function, \({{{{{{{\rm{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, \({\kappa }_{xy}/T=\nu (\pi /12)({k}_{B}^{2}/\hslash )\), where \({\nu}\) 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 \({\nu({{\bf{h}}})}\) on the plane of the magnetic field angles (θ, ϕ). Our major finding is that critical lines of \({\nu({{\bf{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 \(\nu({{{\bf{h}}}})\rightarrow −\nu({{{\bf{h}}}})\).

  • Time reversal operation also reverses the magnetic field direction: h(θ, ϕ) → −h(θ, ϕ) = h(π − θ, ϕ + π).

  • Topologically distinct regions with {\({+\nu({{\bf{h}}})}\), +h} and {\({−\nu({{\bf{h}}})}\), −h} exist on the (θ, ϕ) plane.

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 D3 symmetry of the system. The topological invariant \({\nu}\) and thermal Hall conductivity κxy are A2 representations of the D3 group, i.e., even under C3 rotations but odd under C2 rotations, which reveals the generic form,

$$\nu ({{{{{{{\bf{h}}}}}}}})={{{{{{{\rm{sgn}}}}}}}}[{{{\Lambda }}}_{1}({h}_{x}+{h}_{y}+{h}_{z})+{{{\Lambda }}}_{3}{h}_{x}{h}_{y}{h}_{z}],$$
(2)

where Λ1,3 are field-independent coefficients. The h-linear term (hx + hy + hz) and h-cubic term (hxhyhz) are the leading order A2 representations of magnetic fields. Conducting third order perturbation theory, we find the coefficients

$${{{\Lambda }}}_{1}=-\frac{4\varGamma ^{\prime} }{{{{\Delta }}}_{{{{{{{{\rm{flux}}}}}}}}}}+\frac{6J\varGamma ^{\prime} }{{{{\Delta }}}_{{{{{{{{\rm{flux}}}}}}}}}^{2}}-\frac{4\varGamma \varGamma ^{\prime} }{{{{\Delta }}}_{{{{{{{{\rm{flux}}}}}}}}}^{2}}+\frac{5\varGamma {^{\prime} }^{2}}{2{{{\Delta }}}_{{{{{{{{\rm{flux}}}}}}}}}^{2}}\,\& \,{{{\Lambda }}}_{3}=\frac{18}{{{{\Delta }}}_{{{{{{{{\rm{flux}}}}}}}}}^{2}},\quad$$
(3)

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, \(\nu ({{{{{{{\bf{h}}}}}}}})={{{{{{{\rm{sgn}}}}}}}}({h}_{x}{h}_{y}{h}_{z})\). The dashed lines highlight the critical lines representing the topological phase transitions between the phases with \({\nu({{\bf{h}}})}\) = ±1, where the energy gap of Majorana fermion excitations is closed: Δ(h) = 0.

Exploiting the symmetry analysis, we stress two universal properties of the KQSL with the D3 symmetry.

  1. 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. 2.

    Cubic dependence: for in-plane fields, the h-cubic term governs low field behaviors of the KQSL, e.g., \(\nu ({{{{{{{\bf{h}}}}}}}}) \sim {{{{{{{\rm{sgn}}}}}}}}({h}_{x}{h}_{y}{h}_{z})\) & Δ(h) ~ ∣hxhyhz∣ for θ = 90∘.

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

$${\hat{\chi }}_{p}={S}_{2}^{x}{S}_{1}^{z}{S}_{6}^{y}+{S}_{5}^{x}{S}_{4}^{z}{S}_{3}^{y}+{C}_{3}\,{{{{{{{\rm{rotated}}}}}}}}\,{{{{{{{\rm{terms}}}}}}}}$$
(4)

at each hexagon plaquette p and investigate the expectation value of the chirality operator, shortly the chirality,

$$\chi ({{{{{{{\bf{h}}}}}}}})\equiv \frac{1}{N}\mathop{\sum}\limits_{p}\left\langle {{{\Psi }}}_{{{{{{{{\rm{KQSL}}}}}}}}}({{{{{{{\bf{h}}}}}}}})\right|{\hat{\chi }}_{p}\left|{{{\Psi }}}_{{{{{{{{\rm{KQSL}}}}}}}}}({{{{{{{\bf{h}}}}}}}})\right\rangle ,$$
(5)

and its sign,

$$\bar{\nu }({{{{{{{\bf{h}}}}}}}})\equiv {{{{{{{\rm{sgn}}}}}}}}[\chi ({{{{{{{\bf{h}}}}}}}})],$$
(6)

where \(\left|{{{\Psi }}}_{{{{{{{{\rm{KQSL}}}}}}}}}({{{{{{{\bf{h}}}}}}}})\right\rangle\) is the ground state of the full Hamiltonian H(θ, ϕ) in the KQSL phase (N is the number of unit cells). The chirality operator \({\hat{\chi }}_{p}\) produces the mass term of Majorana fermions and determines the topological invariant in the pure Kitaev limit6. More precisely, how magnetic fields couple to the chirality operator determines the topological invariant and the Majorana energy gap. The chirality χ and its sign \(\bar{\nu }\) are in the A2 representation of the D3 group as of the topological invariant \({\nu}\). 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. Figures 2 and 3 display our major results, the ED calculations of the chirality χED(h) for the KQSL. The used parameter sets are listed in Table 1. The zero lines [χED(h) = 0; dashed lines in the figures] exist in all the cases.

Fig. 2: Chirality of the non-abelian KQSL.
figure 2

Upper: color maps of the chirality χED(h) on the plane of the field angles (θ, ϕ), where the magnetic field strength is fixed by h = 0.01 (horizontal axis: ϕ[°], vertical axis: θ[°]). The dashed lines highlight the zero lines χED(h) = 0, and the black dots mark the bond directions. Lower: χED(h) as a function of h3 for the in-plane fields (θ = 90∘, ϕ = 0∘, 10∘, 20∘, 30∘), illustrating the universality of the h3 behavior in the KQSL. In the case #3, the bending at h3 > 0.5 × 10−6 is an effect of higher order contributions (h5, h7, ...). The parameter sets used in the four cases (#1~4) are listed in Table 1.

Full size image
Table 1 Parameter sets for exact diagonalization and spin wave calculations.
Full size table

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, hx + hy + hz, vanishes for in-plane fields, and the cubic term, hxhyhz, 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 \(\bar{\nu }\) 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: \(\bar{\nu }({{{{{{{\bf{h}}}}}}}})=\nu ({{{{{{{\bf{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 model. This result is completely consistent with the perturbative parton analysis [Eq. (3) and Supplementary Fig. 2b], indicating the validity of our strategy.

Figure 2c, d presents consequences of the other non-Kitaev interactions, \(\varGamma ^{\prime}\) and Γ. The zero lines tend to be flatten around the equator θ = 90∘, which can be attributed to the h-linear term induced by the non-Kitaev interactions: \({h}_{x}+{h}_{y}+{h}_{z}=h\sqrt{3}\cos \theta\). In other words, the zero lines substantially deviate from those of the pure Kitaev model by the non-Kitaev couplings, \(\varGamma ^{\prime}\) 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, \(\varGamma ^{\prime}\) and Γ, the shape of the critical lines is governed by the cubic term hxhyhz, as shown in Fig. 3a, b. In presence of the \(\varGamma ^{\prime}\) or Γ coupling, the h-cubic term competes with the h-linear term as illustrated in Fig. 3c, d. Namely, the linear term dominates over 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.

Fig. 3: Field evolutions of the zero lines.
figure 3

Color maps of the chirality χED(h) on the plane of the field angles (θ, ϕ) with increasing magnetic field h = 0.004, 0.008 (horizontal axis: ϕ[°], vertical axis: θ[°]). The dashed lines highlight the zero lines χED(h) = 0, and the black dots mark the bond directions. The parameter sets used in the four cases (#1~4) are listed in Table 1.

Full size image

Similarities and differences between the topological invariant, \({\nu({{\bf{h}}})}\), and the sign of the chirality, \(\bar{\nu }({{{{{{{\bf{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 D3 group, so \(\bar{\nu }({{{{{{{\bf{h}}}}}}}})\) and \({\nu({{\bf{h}}})}\) have in common the symmetric zeroes. Third, differences between the two quantities may be understood by considering other possible A2 representation spin operators that 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 \({\nu}\) 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 \({\nu}\) (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 D3 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 D3 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 D3 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.

Fig. 4: Comparison of the KQSL with magnetically ordered phases.
figure 4

a KQSL: Majorana energy gap Δparton as a function of h3 for the in-plane fields, (θ = 90∘, ϕ = 0∘, 10∘, 20∘, 30∘), obtained by a parton theory. b–f Ferromagnetic (FM), stripy, vortex, Neel, and zigzag phases: Magnon gap ΔSW as a function of h for the in-plane fields, (θ = 90∘, ϕ = 0∘, 30∘), obtained with a spin wave theory for the parameter sets #5~9 in Table 1.

Full size image

The universal properties of the KQSL can be observed by heat capacity experiments. Figure 5a illustrates the calculated specific heat cv(ϕ) 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 α-RuCl3 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 non-abelian KQSL and zigzag states. Remarkably, such behaviors were observed in the recent heat capacity experiments with in-plane magnetic fields65. Covering the polar angle (θ) in the heat capacity measurements will provide more detailed information on the critical lines and spin interactions in α-RuCl3 (see Fig. 1d).

Fig. 5: Magnetic field angle dependence of specific heat in the KQSL and zigzag states.
figure 5

a KQSL state: Majorana gap Δparton and specific heat cv as functions of in-plane field angle ϕ, obtained by a parton theory. b Zigzag state: magnon gap ΔSW and specific heat cv as functions of in-plane field angle ϕ, obtained by a spin wave theory with \((K,\varGamma ,\varGamma ^{\prime} ,h)=(-1,0.8,-0.05,0.03)\,\& \,{k}_{B}T=0.01\). Magnon gaps for other values of Γ are shown together to highlight the generality of the field angle dependence.

Full size image

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 non-abelian 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 α-RuCl3. 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 quantities61,65,68,69,70,71.

Methods

Exact diagonalization

The KQSL and other magnetic phases of H(θ, ϕ) are mapped out by using the flux operator \({\hat{W}}_{p}={2}^{6}{S}_{1}^{z}{S}_{2}^{y}{S}_{3}^{x}{S}_{4}^{z}{S}_{5}^{y}{S}_{6}^{x}\), the second derivative of the ground state energy \({\partial }^{2}{E}_{{{{{{{{\rm{gs}}}}}}}}}/\partial {\xi }^{2}\,(\xi =J,\varGamma ,\varGamma ^{\prime} )\), and the spin structure factor \(S({{{{{{{\bf{q}}}}}}}})=\frac{1}{N}{\sum }_{i,j}\langle {{{{{{{{\bf{S}}}}}}}}}_{i}\cdot {{{{{{{{\bf{S}}}}}}}}}_{j}\rangle {e}^{i{{{{{{{\bf{q}}}}}}}}\cdot ({{{{{{{{\bf{r}}}}}}}}}_{i}-{{{{{{{{\bf{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 \({{{{{{{\mathcal{S}}}}}}}}\) matrix6,38,72.

Figures 6 and 7 display the phase diagrams of \({H}_{KJ\varGamma {\varGamma }^{\prime}}\). 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 \({\varGamma}\) 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 states10,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 \(\langle {\hat{W}}_{p}\rangle\), 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 \(\varGamma ^{\prime}\) is changed to negative. A zigzag antiferromagnetic order instead sets in under negative sign of \(\varGamma ^{\prime}\) (Supplementary Fig. 6).

Fig. 6: Ferromagnetic KQSL and nearby magnetic states.
figure 6

a Phase diagram of \({H}_{KJ\varGamma {\varGamma }^{\prime}}\) at K = − 1 and \(\varGamma ^{\prime} =0.05\). The color encodes the flux operator expectation value \(\langle {\hat{W}}_{p}\rangle\), and the dashed lines denote phase boundaries determined by the ground state energy second derivatives − ∂2Egs/∂ξ2(ξ = J, Γ). b–e Color maps of the spin structure factor S(q) for the KQSL, QSL, ferromagnetic (FM), and stripy states. The inner and outer hexagons denote the first and second Brillouin zones of the honeycomb lattice.

Full size image

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 \({\varGamma}\) coupling, and surrounded by magnetically ordered states such as the vortex, zigzag, and Neel states10,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 \(\langle {\hat{W}}_{p}\rangle\) 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 \(\varGamma ^{\prime}\) are provided in Supplementary Fig. 6.

Fig. 7: Antiferromagnetic KQSL and adjacent magnetic states.
figure 7

a Phase diagram of \({H}_{KJ\varGamma {\varGamma }^{\prime}}\) at K = 1 and \(\varGamma ^{\prime} =-0.05\). The color encodes the flux operator expectation value \(\langle {\hat{W}}_{p}\rangle\), and the dashed lines denote phase boundaries determined by the ground state energy second derivatives −∂2Egs/∂ξ2(ξ = J, Γ). b–e Color maps of the spin structure factor S(q) for the KQSL, vortex, zigzag, and Neel states. The inner and outer hexagons denote the first and second Brillouin zones of the honeycomb lattice.

Full size image

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 degeneracy6,73 and modular \({{{{{{{\mathcal{S}}}}}}}}\) matrix38,72,74. As an example, the \({{{{{{{\mathcal{S}}}}}}}}\) matrix

$${{{{{{{{\mathcal{S}}}}}}}}}_{{{{{{{{\rm{ED}}}}}}}}}= \, \left\langle {{{\Psi }}}^{{{{{{{{\rm{MES-I}}}}}}}}}| {{{\Psi }}}^{{{{{{{{\rm{MES-II}}}}}}}}}\right\rangle \\ = \, \left[\begin{array}{ccc}0.45{e}^{-i0.08}&0.53{e}^{-i0.03}&0.70{e}^{i0.04}\\ 0.53{e}^{-i0.03}&0.50&-0.71{e}^{-i0.01}\\ 0.70{e}^{i0.04}&-0.71{e}^{-i0.01}&0.02{e}^{-i2.09}\end{array}\right]\\ \approx \, \left[\begin{array}{ccc}1/2&1/2&1/\sqrt{2}\\ 1/2&1/2&-1/\sqrt{2}\\ 1/\sqrt{2}&-1/\sqrt{2}&0\end{array}\right]$$
(7)

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.