Theory of the Kitaev model in a [111] magnetic field

Abstract

Recent numerical studies indicate that the antiferromagnetic Kitaev honeycomb lattice model undergoes a magnetic-field-induced quantum phase transition into a new spin-liquid phase. This intermediate-field phase has been previously characterized as a gapless spin liquid. By implementing a recently developed variational approach based on the exact fractionalized excitations of the zero-field model, we demonstrate that the field-induced spin liquid is gapped and belongs to Kitaev’s 16-fold way. Specifically, the low-field non-Abelian liquid with Chern number C = ±1 transitions into an Abelian liquid with C = ±4. The critical field and the field-dependent behaviors of key physical quantities are in good quantitative agreement with published numerical results. Furthermore, we derive an effective field theory for the field-induced critical point which readily explains the ostensibly gapless nature of the intermediate-field spin liquid.

Introduction

The exactly solvable Kitaev model on the honeycomb lattice¹ has deepened our insight into quantum spin liquids and helped us in identifying strongly spin-orbit-coupled 4d and 5d materials that may host these exotic quantum phases of matter^2,3. Indeed, recent years have seen a flurry of such “Kitaev materials” in which the microscopic spin Hamiltonian is believed to approximately realize the Kitaev honeycomb model^4,5,6,7. The most famous ones include the honeycomb iridates, Na₂IrO₃^{8,9,10,11,12,13}, α-Li₂IrO₃^14,15, and H₃LiIr₂O₆¹⁶, as well as the honeycomb halide α-RuCl₃^{17,18,19,20,21,22,23,24,25}.

While most of these materials are magnetically ordered at the lowest temperatures, the zigzag magnetic order in α-RuCl₃ can be suppressed with an in-plane magnetic field^{26,27,28,29,30,31,32,33,34,35}. Also, there are some experimental indications for an intermediate-field spin-liquid phase between the low-field magnetically ordered phase and the high-field spin-polarized phase. Most importantly, a recent experimental work³⁶ reported a half-integer-quantized thermal Hall conductivity in the intermediate-field regime just beyond the transition out of zigzag order. Though the exact nature of this regime is still an open question, the ongoing experimental efforts reveal the importance of precisely characterizing field-induced spin-liquid phases.

Motivated in large part by the intriguing experimental observations, the behavior of the Kitaev model in a magnetic field has been extensively studied³⁷ by various approaches, including exact diagonalization^{38,39,40,41,42}, density-matrix renormalization group (DMRG)^40,41,42,43, infinite DMRG (iDMRG)⁴⁴, tensor-network methods⁴⁵, continuous-time quantum Monte Carlo techniques⁴⁶, and slave-particle mean-field theories⁴⁷. These approaches all give consistent results. While the ferromagnetic Kitaev model has a single transition into a polarized phase, the antiferromagnetic Kitaev model includes a new intermediate-field spin liquid between the low-field non-Abelian spin liquid¹ and the high-field polarized phase.

In this work, we implement a novel variational approach⁴⁸ to investigate the ground-state phase diagram of the antiferromagnetic Kitaev model in a magnetic field parallel to the [111] direction. This approach is based on the exact fractionalized Majorana-fermion (“spinon”) and gauge-flux (“vison”) excitations of the pure Kitaev model at zero field¹. It accounts for two effects of the magnetic field: the renormalization of the Majorana dispersion through a hybridization with pairs of fluxes (see Fig. 1a) and the finite dispersion acquired by the flux pairs themselves (see Fig. 1b). Remarkably, we find a continuous quantum phase transition, induced by a softening of a hybridized excitation, at a critical field h_c ≃ 0.50, which is very close to the critical field h_c ≃ 0.44 reported by a recent iDMRG study⁴⁴. The critical point signals the transition of the non-Abelian spin liquid¹ with Chern number C = ± 1 into an Abelian spin liquid with C = ± 4. The predicted field dependence of the flux expectation value and the second derivative of the ground-state energy is also in good quantitative agreement with the iDMRG results. Moreover, the effective field theory of the quantum critical point, as derived from the microscopic Hamiltonian, predicts a low-energy ring of gapped excitations in momentum space, which is difficult to be distinguished from a gapless Fermi surface in finite systems. We conjecture that this is the main reason why previous works^{38,40,41,42,43} characterized the phase at h ≳ h_c as a gapless spin liquid.

Fig. 1: Effect of the magnetic field.

a Hybridization between a flux pair and a fermion. b Hopping of a flux pair between two neighboring bonds. Definitions of the plaquettes p, the lattice sites 1−6, the two sublattices A and B, and the nearest-neighbor bond vectors ${\hat{{{{{{\bf{r}}}}}}}}_{x,y,z}$ are also shown.

Model

We consider the antiferromagnetic Kitaev model¹ in an external magnetic field along the [111] direction,

$${{{{{\mathcal{H}}}}}}=\mathop{\sum}\limits_{\alpha }\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in A}{\sigma }_{{{{{{\bf{r}}}}}}}^{\alpha }{\sigma }_{{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha }+h\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}}({\sigma }_{{{{{{\bf{r}}}}}}}^{x}+{\sigma }_{{{{{{\bf{r}}}}}}}^{y}+{\sigma }_{{{{{{\bf{r}}}}}}}^{z}),$$

(1)

where h is the magnetic field (in units of the Kitaev energy) and ${\hat{{{{{{\bf{r}}}}}}}}_{\alpha }$ is the nearest-neighbor vector from an A site to a B site along an α bond (see Fig. 1). For the exactly solvable Kitaev model in the h = 0 limit, the low-energy spectrum comprises gapless matter fermions (i.e., spinons) with a single Dirac cone and gapped dispersionless ${{\mathbb{Z}}}_{2}$ gauge fluxes. These elementary excitations are described in terms of four Majorana fermions ${c}_{{{{{{\bf{r}}}}}}}$ and ${b}_{{{{{{\bf{r}}}}}}}^{\alpha }$ with α = x, y, z at each site r, where ${c}_{{{{{{\bf{r}}}}}}}$ are the matter fermions, and ${b}_{{{{{{\bf{r}}}}}}}^{\alpha }$ are bond fermions associated with the ${{\mathbb{Z}}}_{2}$ gauge field ${u}_{{{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha }\equiv i{b}_{{{{{{\bf{r}}}}}}}^{\alpha }{b}_{{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha }=\pm 1$. The gauge fields are conserved bond variables that commute with each other; their product around any plaquette p (see Fig. 1a) is gauge invariant and expressible in terms of the physical spins:

$${W}_{p}={u}_{12}^{z}{u}_{32}^{x}{u}_{34}^{y}{u}_{54}^{z}{u}_{56}^{x}{u}_{16}^{y}={\sigma }_{1}^{x}{\sigma }_{2}^{y}{\sigma }_{3}^{z}{\sigma }_{4}^{x}{\sigma }_{5}^{y}{\sigma }_{6}^{z}.$$

(2)

Thus, W_p = ±1 can be identified as static ${{\mathbb{Z}}}_{2}$ gauge fluxes. In each flux sector, {W_p = ±1}, represented with an appropriate gauge-field configuration, $\{{u}_{{{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha }=\pm 1\}$, the zero-field model then reduces to a quadratic matter-fermion problem.

While the model in Eq. (1) is not exactly solvable for a finite field, we can derive a low-energy effective model by projecting ${{{{{\mathcal{H}}}}}}$ into the low-energy sector of the pure Kitaev model (corresponding to h = 0) generated by single matter-fermion and/or flux-pair excitations⁴⁸. We focus on flux pairs because, unlike single fluxes, they are coherent fermionic quasiparticles⁴⁸ and can readily hybridize with matter fermions (see Fig. 1a). The fermionic flux-pair excitations can be represented with dressed bond-fermion operators ${({\tilde{\chi }}_{{{{{{\bf{r}}}}}}\in A}^{\alpha })}^{{{\dagger}} }=\frac{1}{2}({\tilde{b}}_{{{{{{\bf{r}}}}}}}^{\alpha }-i{\tilde{b}}_{{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha })$ that have the same projective symmetries as the bare bond-fermion operators ${({\chi }_{{{{{{\bf{r}}}}}}\in A}^{\alpha })}^{{{\dagger}} }=\frac{1}{2}({b}_{{{{{{\bf{r}}}}}}}^{\alpha }-i{b}_{{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha })$. The operator ${({\tilde{\chi }}_{{{{{{\bf{r}}}}}}\in A}^{\alpha })}^{{{\dagger}} }$ turns the ground state of the pure Kitaev model into an excited state with a single-flux pair on the α bond connected to the site r ∈ A by not only creating a bond fermion but also distorting the matter-fermion state: ${({\tilde{\chi }}_{{{{{{\bf{r}}}}}}\in A}^{\alpha })}^{{{\dagger}} }\left|\omega \right\rangle \otimes \left|0\right\rangle =\ \left|{\phi }_{{{{{{\bf{r}}}}}}}^{\alpha }\right\rangle \otimes \left|{\chi }_{{{{{{\bf{r}}}}}}}^{\alpha }\right\rangle$, where $\left|\omega \right\rangle$ and $\left|{\phi }_{{{{{{\bf{r}}}}}}}^{\alpha }\right\rangle$ are the matter-fermion vacua of the gauge-field configurations $\left|0\right\rangle$ and $\left|{\chi }_{{{{{{\bf{r}}}}}}}^{\alpha }\right\rangle$ that correspond to the flux-free sector and the single-flux-pair sector, respectively. (Mathematically, $\left|{\chi }_{{{{{{\bf{r}}}}}}}^{\alpha }\right\rangle ={({\chi }_{{{{{{\bf{r}}}}}}}^{\alpha })}^{{{\dagger}} }\left|0\right\rangle$, while $\left|0\right\rangle$ is the bare-bond-fermion vacuum with ${u}_{{{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha }=-1$ for all bonds.) If we project the pure Kitaev model [i.e., the first term of Eq. (1)] to its low-energy sector containing at most one matter-fermion or flux-pair excitation, the resulting low-energy Hamiltonian reads

$${\tilde{{{{{{\mathcal{H}}}}}}}}_{h = 0}=\mathop{\sum}\limits_{\alpha }\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in A}i{c}_{{{{{{\bf{r}}}}}}}{c}_{{{{{{\bf{r}}}}}}+{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}+{{{\Delta }}}_{\chi }\mathop{\sum}\limits_{\alpha }\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in A}{({\tilde{\chi }}_{{{{{{\bf{r}}}}}}}^{\alpha })}^{{{\dagger}} }({\tilde{\chi }}_{{{{{{\bf{r}}}}}}}^{\alpha }),$$

(3)

where the first term is the quadratic matter-fermion problem within the flux-free sector¹, while the second term accounts for the finite energy (Δ_χ ≃ 0.26) of a flux pair. The Zeeman term [i.e., the second term of Eq. (1)] can then either hybridize a flux pair with a matter fermion (see Fig. 1a) or hop a flux pair to a neighboring bond (see Fig. 1b). By summing ${\tilde{{{{{{\mathcal{H}}}}}}}}_{h = 0}$ and the most general symmetry-allowed Hamiltonians describing these two processes, the effective low-energy Hamiltonian for the full model in Eq. (1) becomes

$$\tilde{{{{{{\mathcal{H}}}}}}}= \, {\tilde{{{{{{\mathcal{H}}}}}}}}_{h = 0}+h\mathop{\sum}\limits_{\alpha }\mathop{\sum}\limits_{{{{{{\bf{R}}}}}}}{p}_{{{{{{\bf{R}}}}}},\alpha }\left[\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in A}i{\tilde{b}}_{{{{{{\bf{r}}}}}}}^{\alpha }{c}_{{{{{{\bf{r}}}}}}+{{{{{\bf{R}}}}}}}+\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in B}i{\tilde{b}}_{{{{{{\bf{r}}}}}}}^{\alpha }{c}_{{{{{{\bf{r}}}}}}-{{{{{\bf{R}}}}}}}\right]\\ \,-ihq\mathop{\sum}\limits_{\alpha ,\beta }{\epsilon }_{\alpha \beta }\left\{\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in A}{({\tilde{\chi }}_{{{{{{\bf{r}}}}}}}^{\alpha })}^{{{\dagger}} }{\tilde{\chi }}_{{{{{{\bf{r}}}}}}}^{\beta }+\mathop{\sum}\limits_{{{{{{\bf{r}}}}}}\in B}{({\tilde{\chi }}_{{{{{{\bf{r}}}}}}-{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}^{\alpha })}^{{{\dagger}} }{\tilde{\chi }}_{{{{{{\bf{r}}}}}}-{\hat{{{{{{\bf{r}}}}}}}}_{\beta }}^{\beta }\right\},$$

(4)

where R is a lattice vector, ϵ_αβ = ∑_γϵ_αβγ is an antisymmetric symbol based on the Levi-Civita symbol ϵ_αβγ, while p_R,α and q are dimensionless parameters to be determined. Notice that some p_R,α are identical due to the threefold rotation symmetry acting simultaneously in real space and spin space.

Since the effective Hamiltonian $\tilde{{{{{{\mathcal{H}}}}}}}$ is quadratic, it can be straightforwardly diagonalized in momentum space:

$$\tilde{{{{{{\mathcal{H}}}}}}}= \mathop{\sum}\limits_{{{{{{\bf{k}}}}}}}\left[i{\lambda }_{{{{{{\bf{k}}}}}}}{C}_{{{{{{\bf{k}}}}}},A}^{{{\dagger}} }{C}_{{{{{{\bf{k}}}}}},B}+{{{{{\rm{H}}}}}}.{{{{{\rm{c}}}}}}.\right]\\ +\mathop{\sum}\limits_{{{{{{\bf{k}}}}}},\alpha ,\beta }\left\{{{{\Delta }}}_{\chi }{\delta }_{\alpha \beta }-ihq{\epsilon }_{\alpha \beta }\left[1+{{{{\rm{e}}}}}^{i{{{{{\bf{k}}}}}}\cdot ({\hat{{{{{{\bf{r}}}}}}}}_{\alpha }-{\hat{{{{{{\bf{r}}}}}}}}_{\beta })}\right]\right\}{({\tilde{X}}_{{{{{{\bf{k}}}}}}}^{\alpha })}^{{{\dagger}} }({\tilde{X}}_{{{{{{\bf{k}}}}}}}^{\beta })\\ +\frac{h}{\sqrt{2}}\mathop{\sum}\limits_{{{{{{\bf{k}}}}}},\alpha }\left\{i{P}_{{{{{{\bf{k}}}}}},\alpha }\left[{\tilde{X}}_{-{{{{{\bf{k}}}}}}}^{\alpha }+{({\tilde{X}}_{{{{{{\bf{k}}}}}}}^{\alpha })}^{{{\dagger}} }\right]{C}_{{{{{{\bf{k}}}}}},A}\right.\\ \left.+\,{P}_{-{{{{{\bf{k}}}}}},\alpha }\ {{{{\rm{e}}}}}^{i{{{{{\bf{k}}}}}}\cdot {\hat{{{{{{\bf{r}}}}}}}}_{\alpha }}\left[{\tilde{X}}_{-{{{{{\bf{k}}}}}}}^{\alpha }-{({\tilde{X}}_{{{{{{\bf{k}}}}}}}^{\alpha })}^{{{\dagger}} }\right]{C}_{{{{{{\bf{k}}}}}},B}+{{{{{\rm{H}}}}}}.{{{{{\rm{c}}}}}}.\right\},$$

(5)

where ${\lambda }_{{{{\bf{k}}}}}={\sum }_{\alpha }{{{{\rm{e}}}}}^{{i{{{{{\bf{k}}}}}}\cdot {\hat{{{{{{\bf{r}}}}}}}_{\alpha} }}}$ and P_k,α = ∑_Rp_R,α e^ik⋅R, while

$$C_{{{{\bf{k}}}},\nu }=\frac{1}{\sqrt{2N}}\mathop{\sum}\limits_{{{{{\bf{r}}}}\in \nu }} c_{{{{\bf{r}}}}} e^{{-i{{{\bf{k}}}}\cdot {{{\bf{r}}}}}},\ \ {\tilde{X}}_{{{{\bf{k}}}}}^{\alpha } = \frac{1}{\sqrt{N}} \mathop{\sum}\limits_{{{{{\bf{r}}}}\in A}} {\tilde{\chi }}_{{{{\bf{r}}}}}^{\alpha } e^{{-i{{{\bf{k}}}}\cdot {{{\bf{r}}}}}}$$

(6)

are momentum-space matter and bond fermions in terms of the sublattice index ν = A, B and the system size N. By considering the matrix elements of the Zeeman term ∝ h in Eq. (1) within the low-energy sector of the pure Kitaev model⁴⁸, we relate the dimensionless parameters in Eq. (5) to matter-fermion matrix elements of this exactly solvable model (Note: see the Supplementary Information for more details on the dimensionless parameters of the effective Hamiltonian, the expectation value of the flux operator, the coefficients of the effective field theory, and the nonanalytic behavior of the ground-state energy):

$$q= \, \langle {\phi }_{{{{{{\bf{0}}}}}}}^{\beta }| (1+i{c}_{{{{{{\bf{0}}}}}}}{c}_{{\hat{{{{{{\bf{r}}}}}}}}_{\alpha }})| {\phi }_{{{{{{\bf{0}}}}}}}^{\gamma }\rangle ,\quad \alpha \ne \beta \ne \gamma ,\\ {P}_{{{{{{\bf{k}}}}}},\alpha }= \, \langle {\phi }_{{{{{{\bf{0}}}}}}}^{\alpha }| \omega \rangle +\frac{1}{2}\mathop{\sum}\limits_{{{{{{{\bf{k}}}}}}}^{\prime}}(1-{{{{\rm{e}}}}}^{-i{{{{{{\bf{k}}}}}}}^{\prime}\cdot {\hat{{{{{{\bf{r}}}}}}}}_{\alpha }+i{\varphi }_{{{{{{{\bf{k}}}}}}}^{\prime}}})\langle {\phi }_{{{{{{\bf{0}}}}}}}^{\alpha }| {\psi }_{{{{{{{\bf{k}}}}}}}^{\prime}}^{{{\dagger}} }{\psi }_{{{{{{\bf{k}}}}}}}^{{{\dagger}} }| \omega \rangle ,$$

(7)

where r = 0 is an A site, while ${\psi }_{{{{{{\bf{k}}}}}}}=({C}_{{{{{{\bf{k}}}}}},A}+i{{{{\rm{e}}}}}^{i{\varphi }_{{{{{{\bf{k}}}}}}}}{C}_{{{{{{\bf{k}}}}}},B})/\sqrt{2}$ in terms of ${{{{\rm{e}}}}}^{i{\varphi }_{{{{{{\bf{k}}}}}}}}={\lambda }_{{{{{{\bf{k}}}}}}}/| {\lambda }_{{{{{{\bf{k}}}}}}}|$ are the matter fermions diagonalizing the flux-free sector of the pure Kitaev model. For a finite honeycomb lattice with N = 121 × 121 unit cells, we numerically find q ≃ 0.0494 and P_0,α ≃ 0.722.

Results

We study the low-energy effective model in Eq. (4) as a function of the magnetic field h. At zero field, the spectrum coincides with that of the pure Kitaev model and contains one dispersive matter-fermion band as well as the three flat bond-fermion bands (see Fig. 2a). For a small field, h ≪ Δ_χ, the hybridization between these four bands gives rise to a finite energy gap, Δ_K(h) ∝ h³, at the K point of the Brillouin zone (BZ). The slow field dependence of Δ_K(h), which is expected from a perturbative argument by Kitaev¹, explains why the global minimum of the band structure remains at the K point up to a large field, h₀ ≃ 0.46. As shown in Fig. 3a, the global minimum switches from the K point to the Γ point at h = h₀, and the corresponding gap, Δ_Γ(h), closes at a slightly larger field, h_c ≃ 0.50 (see Fig. 2b). Since the little group of the Γ point includes the threefold rotation C₃, the fermion eigenmodes at the Γ point can be classified according to their C₃ eigenvalues. The natural bond-fermion modes, corresponding to C₃ eigenvalues 1 and e^{∓2πi/3}, respectively, are then

$$ {\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{0}=\left({\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{x}+{\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{y}+{\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{z}\right)/\sqrt{3},\\ {\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{\pm }=\left({\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{x}+{{{{\rm{e}}}}}^{\pm 2\pi i/3}{\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{y}+{{{{\rm{e}}}}}^{\mp 2\pi i/3}{\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{z}\right)/\sqrt{3}.$$

(8)

Since the matter-fermion mode ψ₀ is invariant under C₃, it can only hybridize with the bond-fermion mode ${\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{0}$. At the critical field, ${h}_{c}=3\sqrt{{{{\Delta }}}_{\chi }/2}\ {({\sum }_{\alpha }{P}_{{{{{{\boldsymbol{0}}}}}},\alpha })}^{−1}\simeq 0.50$, one of the resulting hybridized eigenmodes is gapless. In contrast, there is a higher critical field, $h_c^{\prime} ={{{\Delta }}}_{\chi }/(2\sqrt{3}q)\simeq 1.52$ (not shown in Fig. 3), at which the pure bond-fermion eigenmode ${\tilde{X}}_{{{{{{\boldsymbol{0}}}}}}}^{+}$ has vanishing energy. We note that a complete diagonalization over the full BZ reveals yet another critical point at h_c″ ≃ 1.0 due to the softening of a hybridized mode at the M point. We emphasize, however, that the effective model is no longer expected to be valid when h is significantly larger than h_c.

Fig. 2: Low-energy spectrum of the effective Hamiltonian.

The fermion dispersions correspond to h = 0.05 in (a) and h = h_c ≃ 0.50 in (b). The color scale shows the matter-fermion weight, 0 < Z_ψ < 1, of the given fermion eigenmode; red (blue) color indicates predominantly matter-fermion (bond-fermion) character. The insets show the spectrum over the full energy range. Note that a hybridization decay length, ξ = 25, is used to regularize the K-point behavior (See the “Note” above earlier).

Fig. 3: Field dependence of key physical quantities.

a Overall energy gap. The insets show the dispersion of the low-energy fermion eigenmode on both sides of the phase transition, with the black solid hexagon marking the Brillouin zone. The red (blue) line corresponds to the gap Δ_K (Δ_Γ), while the black line corresponds to the gap at the six wave vectors Q_j, the corners of the blue hexagon in the right-hand-side inset. b Second derivative of the ground-state energy. c Expectation value of the ${{\mathbb{Z}}}_{2}$ gauge flux.

Figure 3a shows the overall energy gap as a function of the magnetic field h. As expected, the gap is proportional to h³ at the smallest fields, h ≪ Δ_χ. Just below h_c, the global minimum of the excitation spectrum switches from the K point to the Γ point, and the gap vanishes at h_c ≃ 0.50^38,40,44. Importantly, the zero-energy mode at h = h_c has dominant bond-fermion character with a large bond-fermion weight 6/(6 + Δ_χ) ≃ 0.96 (see also Fig. 2b), which is consistent with the numerical closing of the vison gap in the specific heat³⁸. In contrast, the gap reopens for h ≳ h_c, which appears to be in contradiction with the same numerical results and the corresponding conjecture of a gapless U(1) spin liquid at intermediate fields. However, our analytic approach can also explain the numerical similarity between the gapped spin liquid at h ≳ h_c and a gapless spin liquid with a circular spinon Fermi surface. Indeed, as we explain below, the phase transition at h = h_c gives rise to a low-energy ring at h ≳ h_c (see the inset of Fig. 3a) which expands from the Γ point and corresponds to a small energy gap $\propto {(h-{h}_{c})}^{3/2}$. This low-energy ring naturally explains the large low-energy density of states found by exact diagonalization^38,40. The emergence of the low-energy ring and the nature of the h ≳ h_c phase are explained in the next section, where we derive an effective field theory to describe the continuous topological phase transition at h = h_c.

Figures 3b and c plot the second derivative of the ground-state energy, ${E}_{G}^{^{\prime\prime} }={{{{\rm{d}}}}}^{2}{E}_{G}/{{{\rm{d}}}}{h}^{2}$, and the expectation value of the ${{\mathbb{Z}}}_{2}$ gauge flux, 〈W_p〉, against the magnetic field. As we explain below, the discontinuity of ${E}_{G}^{^{\prime\prime} }$ at h = h_c is a generic property of the corresponding phase transition. This discontinuity leads to a peak in ${E}_{G}^{^{\prime\prime} }$ at h = h_c, which is qualitatively and quantitatively consistent with the iDMRG results⁴⁴. We note that our result for 〈W_p〉 (See the “Note” above earlier) (see Fig. 3c) is also consistent with iDMRG.

We argue that our effective model in Eq. (4) remains valid up to a field h ≳ h_c just beyond the first phase transition. Indeed, the fractionalized excitations of the pure Kitaev model remain well defined throughout the low-field phase at h < h_c; however, after the first phase transition induced by their softening, these original excitations are superseded by the emergent excitations of the higher-field phase. Therefore, we focus on the first phase transition at h = h_c throughout the rest of this work.

Remarkably, the critical field h_c ≃ 0.50 is only 10% higher than the corresponding iDMRG result, h_c ≃ 0.44⁴⁴. Also, the slight overestimation of h_c is not surprising because the inclusion of higher-energy (E ≃ 2Δ_χ) states with four fluxes and one matter fermion would lead to a reduction of h_c. Finally, at h = h_c, the dynamical spin structure factor from iDMRG indicates that the spin excitation gap closes at the Γ point, which is in agreement with our results. Indeed, since a spin excitation fractionalizes into a pair of fermion excitations, and the fermions at h = h_c are gapless at the Γ point (see Fig. 2b), a pair of gapless fermions has zero total momentum, corresponding to a vanishing spin gap at the Γ point. These similarities between the iDMRG results and those obtained from our effective Hamiltonian $\tilde{{{{{{\mathcal{H}}}}}}}$ indicate that our variational low-energy manifold captures the essence of the phase transition at h = h_c and the new spin-liquid phase at h ≳ h_c.

Field theory of topological phase transition

In the vicinity of the critical field, h ≃ h_c ≃ 0.50, the low-energy fermion eigenmodes belong to the trivial representation of C₃, and the long-wavelength limit of $\tilde{{{{{{\mathcal{H}}}}}}}$, corresponding to the region around the Γ point, can be written as

$${\tilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{eff}}}}}}}=\mathop{\sum}\limits_{{{{{{\boldsymbol{k}}}}}}}{f}_{{{{{{\boldsymbol{k}}}}}}}^{{{\dagger}} }[{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x}{\tau }_{x}+{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{y}{\tau }_{y}+{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}{\tau }_{z}]{f}_{{{{{{\boldsymbol{k}}}}}}},$$

(9)

where τ_x,y,z are the Pauli matrices, and ${f}_{{{{{{\boldsymbol{k}}}}}}}={({f}_{1,{{{{{\boldsymbol{k}}}}}}},{f}_{2,{{{{{\boldsymbol{k}}}}}}})}^{{{{\rm{T}}}}}$ is a two-component fermionic operator corresponding to the two zero-energy modes of $\tilde{{{{{{\mathcal{H}}}}}}}$ at the critical field:

$${f}_{1,{{{{{\boldsymbol{k}}}}}}}=\sqrt{\frac{6}{6+{{{\Delta }}}_{\chi }}}{\tilde{X}}_{{{{{{\boldsymbol{k}}}}}}}^{0}-i\sqrt{\frac{{{{\Delta }}}_{\chi }}{6+{{{\Delta }}}_{\chi }}}{\psi }_{{{{{{\boldsymbol{k}}}}}}},\,\,\,\,{f}_{2,{{{{{\boldsymbol{k}}}}}}}={f}_{1,-{{{{{\boldsymbol{k}}}}}}}^{{{\dagger}} }.$$

(10)

The coefficients ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x,y,z}$ in Eq. (9) must be C₃ invariant real polynomials. Up to cubic order in k, there are only four such polynomials: the trivial polynomial 1, the quadratic polynomial ${k}^{2}={k}_{x}^{2}+{k}_{y}^{2}$, and the cubic polynomials ${g}_{{{{{{\boldsymbol{k}}}}}}}^{x}={k}_{x}(3{k}_{y}^{2}-{k}_{x}^{2})$ and ${g}_{{{{{{\boldsymbol{k}}}}}}}^{y}={k}_{y}(3{k}_{x}^{2}-{k}_{y}^{2})$. Moreover, the particle-hole symmetry of the original Hamiltonian ${{{{{\mathcal{H}}}}}}$ dictates that ${\tilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{eff}}}}}}}$ must remain invariant under ${f}_{{{{{{\boldsymbol{k}}}}}}}\to {\tau }_{x}{({f}_{-{{{{{\boldsymbol{k}}}}}}}^{{{\dagger}} })}^{{{{\rm{T}}}}}$, implying that the polynomials ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{\mu }$ must satisfy the following relationships:

$${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x}=-{\beta }_{-{{{{{\boldsymbol{k}}}}}}}^{x},\,\,{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{y}=-{\beta }_{-{{{{{\boldsymbol{k}}}}}}}^{y},\,\,{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}={\beta }_{-{{{{{\boldsymbol{k}}}}}}}^{z}.$$

(11)

These symmetry considerations then lead to the general forms

$$ {\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}={c}_{0}+{c}_{z}{k}^{2},\\ {\beta }_{{{{{{\boldsymbol{k}}}}}}}^{\eta }=\mathop{\sum}\limits_{\nu =x,y}{c}_{\eta \nu }{g}_{{{{{{\boldsymbol{k}}}}}}}^{\nu },\quad \eta =x,y,$$

(12)

where c₀, c_z, and c_ην are, in general, functions of h. Since the phase transition at h = h_c is driven by a sign change in c₀, we assume that c_z and c_ην are constants, while we write ${c}_{0}=c_{0}^{\prime} (h-{h}_{c})$ with a constant $c_{0}^{\prime}$. Starting from Eqs. (5) and (7), and defining all lengths in units of the lattice vector (i.e., the distance between two neighboring A sites), the constants are derived to be $c_{0}^{\prime} \simeq -1.00$, c_z ≃ 0.0125, c_xx ≃ − 0.00268, c_yy ≃ − 0.00088, and c_xy = c_yx = 0 (See the “Note” above earlier). Then, using Eq. (9), the fermion dispersion is given by

$${\omega }_{{{{{{\boldsymbol{k}}}}}}}=\sqrt{{\left({\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x}\right)}^{2}+{\left({\beta }_{{{{{{\boldsymbol{k}}}}}}}^{y}\right)}^{2}+{\left({\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}\right)}^{2}}$$

(13)

and becomes gapless at k = 0 for h = h_c. For h < h_c, the dispersion is dominated by the function ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}$ and is largely quadratic: ${\omega }_{{{{{{\boldsymbol{k}}}}}}}\simeq | c_{0}^{\prime} | ({h}_{c}-h)+{c}_{z}{k}^{2}$. In contrast, for h > h_c, the function ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z}$ vanishes for $| {{{{{\boldsymbol{k}}}}}}| =\sqrt{| c_{0}^{\prime} | (h-{h}_{c})/{c}_{z}}$. Thus, along this ring of radius ∣k∣, the energy gap is determined by the small cubic contributions from ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x,y}$ and has a slow field dependence: ${{\Delta }}\propto {(h-{h}_{c})}^{3/2}$. The net result is a ring of low-energy fermions around the Γ point (see the inset of Fig. 3a).

The effective field theory in Eq. (9) describes a continuous topological phase transition. The phases on both sides of the transition belong to Kitaev’s 16-fold way¹ and are characterized by the fermion Chern number. The contribution from the low-energy fermions to this Chern number is given by⁴⁹

$$C=\frac{1}{4\pi }\int {{{\rm{d}}}}{{{{{\boldsymbol{k}}}}}}\ {{{{{{\boldsymbol{d}}}}}}}_{{{{{{\boldsymbol{k}}}}}}}\cdot [{\partial }_{{k}_{x}}{{{{{{\boldsymbol{d}}}}}}}_{{{{{{\boldsymbol{k}}}}}}}\times {\partial }_{{k}_{y}}{{{{{{\boldsymbol{d}}}}}}}_{{{{{{\boldsymbol{k}}}}}}}],$$

(14)

where d_k = β_k/∣β_k∣ and ${{{{{{\boldsymbol{\beta }}}}}}}_{{{{{{\boldsymbol{k}}}}}}}=({\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x},{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{y},{\beta }_{{{{{{\boldsymbol{k}}}}}}}^{z})$. Geometrically, C is simply the skyrmion number of the vector field d_k. Figure 4 depicts the vector field d_k around the Γ point on both sides of the phase transition at h = h_c. While the field configuration is topologically trivial for h < h_c, it includes six merons (three skyrmions) for h > h_c. The corresponding change in the Chern number, ΔC = 3, is then a generic property of the phase transition described by ${\tilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{eff}}}}}}}$. To understand the emergence of the six merons around the Γ point, we first note that ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{\eta }\propto \,{{\mbox{Im}}}\,({k}_{+}^{3}{{{{\rm{e}}}}}^{-i{\phi }_{\eta }})$ with k₊ = k_x + ik_y and ${\phi }_{\eta }=\arctan ({c}_{\eta x}/{c}_{\eta y})$. Each function ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{\eta }$ (with η = x, y) possesses three nodal lines corresponding to ${k}_{y}/{k}_{x}=\tan ({\phi }_{\eta }/3+\varphi )$ with φ = 0, π/3, 2π/3. Ignoring the ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{y}$ function, the low-energy spectrum then contains six Dirac nodes Q_j (with j = 1, 2, . . . , 6) at the intersections of the nodal lines of ${\beta }_{{{{{{\boldsymbol{k}}}}}}}^{x}$ and the ring of radius $| {{{{{\boldsymbol{k}}}}}}| =\sqrt{| c_0^{\prime} | (h-{h}_{c})/{c}_{z}}$. The vorticity of the vector field d_k around each Dirac node Q_j is (−1)^j. Assuming ϕ_x ≠ ϕ_y (which is true in our case), the finite value of ${\beta }_{{{{{{{\boldsymbol{Q}}}}}}}_{j}}^{y}\propto {(-1)}^{j}$ generates a mass term for each Dirac node in such a way that the Dirac nodes all give identical contributions (+1/2 each or −1/2 each) to the change in the Chern number. The net change in the Chern number is then

$${{\Delta }}C=3\ {{{{{\rm{sgn}}}}}}\left[\det \hat{{{{{{\mathcal{C}}}}}}}\right],\quad \hat{{{{{{\mathcal{C}}}}}}}=\left(\begin{array}{ccc}{c}_{xx}&{c}_{xy}&0\\ {c}_{yx}&{c}_{yy}&0\\ 0&0&{c}_{z}\end{array}\right).$$

(15)

Using the constants c_z and c_ην given above, we obtain ΔC = 3 at the critical field h = h_c. Since the low-field phase at h < h_c is well known¹ to have Chern number 1, we conclude that the higher-field phase at h ≳ h_c has Chern number 4.

Fig. 4: Topological phase transition.

Configuration of the unit-vector field d_k on the two sides of the phase transition, a h < h_c and b h > h_c. The color scale shows the component ${d}_{{{{{{\boldsymbol{k}}}}}}}^{y}$, while the black arrows represent the components $({d}_{{{{{{\boldsymbol{k}}}}}}}^{z},{d}_{{{{{{\boldsymbol{k}}}}}}}^{x})$. The green circle marks the low-energy ring.

We next consider the second derivative of the ground-state energy ${E}_{G}^{^{\prime\prime} }$ with respect to the magnetic field h. The universal critical behavior at h = h_c is determined by the low-energy modes ∣k∣≤Λ, where the cutoff Λ can be made arbitrarily small (corresponding to an infrared singularity). While the contribution of these modes to ${E}_{G}^{^{\prime\prime} }$ is ∝ Λ² for $h\to {h}_{c}^{-}$, it is an ${{{{{\mathcal{O}}}}}}(1)$ constant for $h\to {h}_{c}^{+}$. In particular, there is a contribution from the neighborhood of the low-energy ring at h ≳ h_c which is independent of the cutoff Λ. Therefore, we obtain a discontinuity in ${E}_{G}^{^{\prime\prime} }$ at the critical field (See the “Note” above earlier):

$${{\Delta }}{E}_{G}^{^{\prime\prime} }=\mathop{\lim }\limits_{h\to {h}_{c}^{-}}{E}_{G}^{^{\prime\prime} }-\mathop{\lim }\limits_{h\to {h}_{c}^{+}}{E}_{G}^{^{\prime\prime} }=\frac{\sqrt{3}{(c_{0}^{\prime} )}^{2}}{8\pi {c}_{z}}.$$

(16)

Remarkably, this discontinuity in ${E}_{G}^{^{\prime\prime} }$, as shown in Fig. 3b, is entirely determined by two coefficients of the effective field theory. From the constants $c_{0}^{\prime}$ and c_z given above, it is found to be ${{\Delta }}{E}_{G}^{^{\prime\prime} }\simeq 5.5$, which is consistent with the corresponding result for a finite lattice (see Fig. 3b). The quantitative agreement between this value and the one obtained from iDMRG⁴⁴ indicates that the effective field theory at h = h_c is both qualitatively correct and quantitatively accurate.

Discussion

Our simple and accurate variational approach to extended Kitaev models⁴⁸ indicates that the antiferromagnetic (AFM) Kitaev model undergoes a continuous quantum phase transition driven by a magnetic field parallel to the [111] direction. According to this approach, the new phase, which has been reported in previous numerical works^{38,39,40,41,42,43,44}, is a gapped chiral spin liquid with a ring of low-energy excitations. Due to its large low-energy density of states, it is difficult for numerical simulations to distinguish this low-energy ring from a gapless Fermi surface. In particular, while DMRG may, in principle, detect gapless modes via a finite value of the central charge^41,42,43,44, different studies find conflicting values⁴³ or even unphysical non-integer values⁴⁴, thereby indicating that the currently available system sizes cannot be used to determine whether the new phase is gapped or gapless⁴⁴.

In contrast to the non-Abelian low-field phase, the new phase at higher fields possesses Abelian topological order with four distinct types of anyons: 1 (vacuum), ε (fermion), as well as e and m (vortices). The two phases can then be distinguished numerically by computing the entanglement spectrum⁵⁰ or the topological entanglement entropy for a bipartition of an infinite cylinder^51,52,53, readily available in iDMRG⁴⁴. However, due to the challenges mentioned above, such a numerical confirmation of our predictions may require the addition of irrelevant Hamiltonian terms that increase the gap in the higher-field phase without generating new phase transitions.

From an experimental perspective, it is important to note that the higher-field spin liquid is known to be stable against both Heisenberg and Gamma interactions³⁸, making it more likely to emerge in real materials. Also, in the presence of ferromagnetic Heisenberg terms, a field-induced transition between the higher-field spin liquid and a lower-field zigzag order, potentially relevant for α-RuCl₃, has been reported⁴². According to our theory, the key experimental signature of the higher-field spin liquid is a specific quantized value of the thermal Hall conductivity, ${\kappa }_{xy}=\pi {k}_{B}^{2}T/(3\hslash )$, which is four times larger than for the low-field non-Abelian spin liquid.

We next remark that our variational approach is still approximately valid in the presence of both a matter-fermion and a flux-pair excitation and that, in the presence of non-Kitaev interactions, it can also be used to describe bound states between these two types of excitations⁴⁸. Since such a bound state corresponds to a spin excitation, its softening leads to a divergent magnetic susceptibility for some wave vector and thus signals the onset of magnetic ordering.

We also emphasize that our approach straightforwardly generalizes to the ferromagnetic (FM) Kitaev model. In this case, the first term in Eq. (3) has a negative sign, and the flux-pair-hopping parameter in Eq. (7) is found to be q ≃ 1.35, i.e., about 30 times larger than for the AFM Kitaev model. Therefore, the lowest-field phase transition is driven by a softening of a pure flux-pair mode and happens at a much smaller critical field, $h^{\prime} ={{{\Delta }}}_{\chi }/(2\sqrt{3}q)\simeq 0.056$. The strong asymmetry between the FM and AFM Kitaev models is due to opposite (constructive and destructive) interference effects between the two processes contributing to flux-pair hopping⁴⁸. We note that this asymmetry is not apparent in the simplified perturbative analysis of ref. ¹ because it neglects the energy dispersions of the intermediate states. We further remark that our results for the FM Kitaev model are also consistent with numerical studies that report a single first-order transition into a trivial polarized phase at a critical field h_p ≃ 0.028⁴⁴. At this first-order phase transition, corresponding to ${h}_{p}\lesssim h^{\prime}$, the fluxes suddenly proliferate and confine all fractionalized excitations.

Finally, going back to the AFM Kitaev model, it is interesting to note that a recent work⁵⁴ has also found a field-induced chiral spin liquid phase with Chern number C = 4 through a completely different approach.