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

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[1] 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-7]. The most famous ones include the honeycomb iridates, Na2IrO3[8-13], α-Li2IrO3[14,15], and H3LiIr2O6[16], as well as the honeycomb halide α-RuCl3[17-25].

While most of these materials are magnetically ordered at the lowest temperatures, the zigzag magnetic order in α-RuCl3 can be suppressed with an in-plane magnetic field[26-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[36] 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[37] by various approaches, including exact diagonalization[38-42], density-matrix renormalization group (DMRG)[40-43], infinite DMRG (iDMRG)[44], tensor-network methods[45], continuous-time quantum Monte Carlo techniques[46], and slave-particle mean-field theories[47]. 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[1] and the high-field polarized phase.

In this work, we implement a novel variational approach[48] 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[1]. 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 ≃ 0.50, which is very close to the critical field h ≃ 0.44 reported by a recent iDMRG study[44]. The critical point signals the transition of the non-Abelian spin liquid[1] 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-43] characterized the phase at h ≳ h 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 are also shown.

Model

We consider the antiferromagnetic Kitaev model[1] in an external magnetic field along the [111] direction,where h is the magnetic field (in units of the Kitaev energy) and 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 gauge fluxes. These elementary excitations are described in terms of four Majorana fermions and with α = x, y, z at each site r, where are the matter fermions, and are bond fermions associated with the gauge field . 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:Thus, W = ±1 can be identified as static gauge fluxes. In each flux sector, {W = ±1}, represented with an appropriate gauge-field configuration, , 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 into the low-energy sector of the pure Kitaev model (corresponding to h = 0) generated by single matter-fermion and/or flux-pair excitations[48]. We focus on flux pairs because, unlike single fluxes, they are coherent fermionic quasiparticles[48] and can readily hybridize with matter fermions (see Fig. 1a). The fermionic flux-pair excitations can be represented with dressed bond-fermion operators that have the same projective symmetries as the bare bond-fermion operators . The operator 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: , where and are the matter-fermion vacua of the gauge-field configurations and that correspond to the flux-free sector and the single-flux-pair sector, respectively. (Mathematically, , while is the bare-bond-fermion vacuum with 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 readswhere the first term is the quadratic matter-fermion problem within the flux-free sector[1], 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 and the most general symmetry-allowed Hamiltonians describing these two processes, the effective low-energy Hamiltonian for the full model in Eq. (1) becomeswhere R is a lattice vector, ϵ = ∑ϵ is an antisymmetric symbol based on the Levi-Civita symbol ϵ, while p and q are dimensionless parameters to be determined. Notice that some p are identical due to the threefold rotation symmetry acting simultaneously in real space and spin space.

Since the effective Hamiltonian is quadratic, it can be straightforwardly diagonalized in momentum space:where and P = ∑p e, whileare 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[48], 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):where r = 0 is an A site, while in terms of 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.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, Δ(h) ∝ h3, at the K point of the Brillouin zone (BZ). The slow field dependence of Δ(h), which is expected from a perturbative argument by Kitaev[1], explains why the global minimum of the band structure remains at the K point up to a large field, h0 ≃ 0.46. As shown in Fig. 3a, the global minimum switches from the K point to the Γ point at h = h0, and the corresponding gap, ΔΓ(h), closes at a slightly larger field, h ≃ 0.50 (see Fig. 2b). Since the little group of the Γ point includes the threefold rotation C3, the fermion eigenmodes at the Γ point can be classified according to their C3 eigenvalues. The natural bond-fermion modes, corresponding to C3 eigenvalues 1 and e∓2, respectively, are thenSince the matter-fermion mode ψ is invariant under C3, it can only hybridize with the bond-fermion mode . At the critical field, , one of the resulting hybridized eigenmodes is gapless. In contrast, there is a higher critical field, (not shown in Fig. 3), at which the pure bond-fermion eigenmode has vanishing energy. We note that a complete diagonalization over the full BZ reveals yet another critical point at h″ ≃ 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.

Fig. 2

Low-energy spectrum of the effective Hamiltonian.

The fermion dispersions correspond to h = 0.05 in (a) and h = h ≃ 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 Δ (ΔΓ), while the black line corresponds to the gap at the six wave vectors , 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 gauge flux.

Figure 3a shows the overall energy gap as a function of the magnetic field h. As expected, the gap is proportional to h3 at the smallest fields, h ≪ Δ. Just below h, the global minimum of the excitation spectrum switches from the K point to the Γ point, and the gap vanishes at h ≃ 0.50[38,40,44]. Importantly, the zero-energy mode at h = h 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[38]. In contrast, the gap reopens for h ≳ h, 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 and a gapless spin liquid with a circular spinon Fermi surface. Indeed, as we explain below, the phase transition at h = h gives rise to a low-energy ring at h ≳ h (see the inset of Fig. 3a) which expands from the Γ point and corresponds to a small energy gap . 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 phase are explained in the next section, where we derive an effective field theory to describe the continuous topological phase transition at h = h.

Figures 3b and c plot the second derivative of the ground-state energy, , and the expectation value of the gauge flux, 〈W〉, against the magnetic field. As we explain below, the discontinuity of at h = h is a generic property of the corresponding phase transition. This discontinuity leads to a peak in at h = h, which is qualitatively and quantitatively consistent with the iDMRG results[44]. We note that our result for 〈W〉 (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 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; 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 throughout the rest of this work.

Remarkably, the critical field h ≃ 0.50 is only 10% higher than the corresponding iDMRG result, h ≃ 0.44[44]. Also, the slight overestimation of h 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. Finally, at h = h, 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 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 indicate that our variational low-energy manifold captures the essence of the phase transition at h = h and the new spin-liquid phase at h ≳ h.

Field theory of topological phase transition

In the vicinity of the critical field, h ≃ h ≃ 0.50, the low-energy fermion eigenmodes belong to the trivial representation of C3, and the long-wavelength limit of , corresponding to the region around the Γ point, can be written aswhere τ are the Pauli matrices, and is a two-component fermionic operator corresponding to the two zero-energy modes of at the critical field:The coefficients in Eq. (9) must be C3 invariant real polynomials. Up to cubic order in , there are only four such polynomials: the trivial polynomial 1, the quadratic polynomial , and the cubic polynomials and . Moreover, the particle-hole symmetry of the original Hamiltonian dictates that must remain invariant under , implying that the polynomials must satisfy the following relationships:These symmetry considerations then lead to the general formswhere c0, c, and c are, in general, functions of h. Since the phase transition at h = h is driven by a sign change in c0, we assume that c and c are constants, while we write with a constant . 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.0125, c ≃ − 0.00268, c ≃ − 0.00088, and c = c = 0 (See the “Note” above earlier). Then, using Eq. (9), the fermion dispersion is given byand becomes gapless at  = 0 for h = h. For h < h, the dispersion is dominated by the function and is largely quadratic: . In contrast, for h > h, the function vanishes for . Thus, along this ring of radius ∣∣, the energy gap is determined by the small cubic contributions from and has a slow field dependence: . 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[1] and are characterized by the fermion Chern number. The contribution from the low-energy fermions to this Chern number is given by[49]where  = /∣∣ and . Geometrically, C is simply the skyrmion number of the vector field . Figure 4 depicts the vector field around the Γ point on both sides of the phase transition at h = h. While the field configuration is topologically trivial for h < h, it includes six merons (three skyrmions) for h > h. The corresponding change in the Chern number, ΔC = 3, is then a generic property of the phase transition described by . To understand the emergence of the six merons around the Γ point, we first note that with k+ = k + ik and . Each function (with η = x, y) possesses three nodal lines corresponding to with φ = 0, π/3, 2π/3. Ignoring the function, the low-energy spectrum then contains six Dirac nodes (with j = 1, 2, . . . , 6) at the intersections of the nodal lines of and the ring of radius . The vorticity of the vector field around each Dirac node is (−1). Assuming ϕ ≠ ϕ (which is true in our case), the finite value of 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 thenUsing the constants c and c given above, we obtain ΔC = 3 at the critical field h = h. Since the low-field phase at h < h is well known[1] to have Chern number 1, we conclude that the higher-field phase at h ≳ h has Chern number 4.

Fig. 4

Topological phase transition.

Configuration of the unit-vector field on the two sides of the phase transition, a h < hc and b h > hc. The color scale shows the component , while the black arrows represent the components . The green circle marks the low-energy ring.

We next consider the second derivative of the ground-state energy with respect to the magnetic field h. The universal critical behavior at h = h is determined by the low-energy modes ∣∣≤Λ, where the cutoff Λ can be made arbitrarily small (corresponding to an infrared singularity). While the contribution of these modes to is ∝ Λ2 for , it is an constant for . In particular, there is a contribution from the neighborhood of the low-energy ring at h ≳ h which is independent of the cutoff Λ. Therefore, we obtain a discontinuity in at the critical field (See the “Note” above earlier):Remarkably, this discontinuity in , as shown in Fig. 3b, is entirely determined by two coefficients of the effective field theory. From the constants and c given above, it is found to be , 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[44] indicates that the effective field theory at h = h is both qualitatively correct and quantitatively accurate.

Discussion

Our simple and accurate variational approach to extended Kitaev models[48] 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-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-44], different studies find conflicting values[43] or even unphysical non-integer values[44], thereby indicating that the currently available system sizes cannot be used to determine whether the new phase is gapped or gapless[44].

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[50] or the topological entanglement entropy for a bipartition of an infinite cylinder[51-53], readily available in iDMRG[44]. 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[38], 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 α-RuCl3, has been reported[42]. According to our theory, the key experimental signature of the higher-field spin liquid is a specific quantized value of the thermal Hall conductivity, , 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[48]. 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, . 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[48]. We note that this asymmetry is not apparent in the simplified perturbative analysis of ref. [1] 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 ≃ 0.028[44]. At this first-order phase transition, corresponding to , 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[54] has also found a field-induced chiral spin liquid phase with Chern number C = 4 through a completely different approach.