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.

where k = k /| k | and k = k·r , the hybridization term in Eq. (4) of the main text can then be written as The two terms in Eq. (S2) can be matched with the microscopic model by considering the matrix elements of the bare Zeeman term, H ℎ = ℎ r r r , between appropriate states of the pure Kitaev model: the ground state, |Ω = | ⊗ |0 , the matterfermion eigenstate, | k = ( † k | ) ⊗ |0 , the flux-pair eigenstate, |˜ r = | r ⊗ [( r ) † |0 ], and the approximate eigenstate containing both a matter fermion and a flux pair, |˜ r −k = ( † −k | r ) ⊗ [( r ) † |0 ], where r∈ = 1 2 ( r + r+r ) are the bond fermions of the pure Kitaev model, |0 is the vacuum of these bond fermions, | is the matter-fermion ground state of the flux-free sector, and | r is the matter-fermion ground state of the flux sector with a single flux pair around the bond connecting the sites r ∈ and (r +r ) ∈ . According to Eq. (S2), these matrix elements must be equal to and, therefore, the hybridization parameters k, are found to be where { , } ≡ + . Setting r = 0 without loss of generality, these hybridization parameters then become Importantly, each matter fermion k with a given energy k belongs to a degenerate set of matter fermions whose momenta k are related by the various symmetries of the Kitaev model. From the perspective of these matter fermions, the presence of a flux pair is a local perturbation which only affects the two sites 0 ∈ andr ∈ connected by the corresponding bond. Therefore, we can form appropriate linear combinations of the degenerate matter fermions such that only two linear combinations couple to the perturbation while the remaining ones have vanishing wave functions at both sites 0 andr . Exploiting the residual inversion symmetry around the flux pair (i.e., the corresponding bond), the natural choice for these two linear combinations is where {k} is the set of all momenta k satisfying k = . Since the matter fermions ,± have eigenvalues ± under the residual inversion symmetry (which acts projectively on the matter fermions), the two-fermion matrix element in Eq. (S5) then becomes the momentum-space hybridization parameters in Eq. (S5) finally take the form It can be shown numerically (and argued analytically) that, at the lowest energies, → 0 and k → 0, the leading-order behavior of the two-fermion matrix element in Eq. (S9) is given by ( ) ∼ is the low-energy matter-fermion density of states. Consequently, while ( k ) is finite at k = 0, it is not analytic because its derivative diverges. Moreover, the factor − k has a nontrivial phase winding around the K point. This nonanalytic behavior of the hybridization function k, reflects the gapless Dirac cone of the matter fermions. However, the matter fermions are known to be gapped out by an infinitesimally small magnetic field. Thus, a finite field should remove the nonanalytic behavior by generating an exponential decay for the real-space hybridization function R, . On a phenomenological level, we can account for this exponential decay by multiplying R, with exp(−|R|/ ), which is equivalent to taking a convolution between k, and a Lorentzian function of width ∼ 1/ in momentum space. We emphasize that, while this regularization procedure is important for producing the correct field dependence of the gap opening at the K point, it has negligible effects on all the other results of this work. In practice, we take = 25 for the hybridization decay length.
To determine the hopping parameter , we match the flux-pair hopping term in Eq.
(3) of the main text with the microscopic model by considering the matrix element of the bare Zeeman term, H ℎ = ℎ r r , where r = r r = − 2 , r r , between the flux-pair eigenstates |˜ 0 = | 0 ⊗ [( 0 ) † |0 ] and |˜ 0 = | 0 ⊗ [( 0 ) † |0 ] of the pure Kitaev model. The hopping parameter is then found to be Also, we can straightforwardly determine the signs of 0 and by considering appropriate products of the corresponding matrix elements ˜ 0 | 0 |˜ 0 = − 0 and ˜ 0 | r |˜ 0 = − (along with their cyclic permutations in , , ): Since the individual matrix elements are expected to be O (1) due to the absence of an orthogonality catastrophe, it is a reasonable approximation to neglect the projectors to the intermediate states. In this approximation, the products in Eq. (S11) become

II. EXPECTATION VALUE OF THE FLUX OPERATOR
Here we describe how the expectation value of the Z 2 gauge-flux operator in Eq.
(2) of the main text can be computed for the effective HamiltonianH in Eq. (3) in the main text. We first recognize that the dressed bond-fermion operators˜ r and the bare bond-fermion operators r have exactly the same effect on the flux degrees of freedom as they only differ in an appropriate distortion of the matter-fermion state. Therefore, in terms of the dressed bond variables˜ = ˜ ˜ and the ground state |Ω of the quadratic fermion HamiltonianH , the expectation value of the flux operator becomes where the subscript = 1, 2, ..., 6 labels the six sites around the hexagon (see Fig. 1 of the main text). In turn, this expectation value can be computed by means of Wick's theorem, which reduces the 12-fermion expectation value to products of two-fermion expectation values, ˜ ˜ = Ω |˜ ˜ |Ω . If we then write the bond-fermion operators˜ in terms of the fermion eigenmodes , (with = 1, 2, ..., 8) of the effective Hamiltonian, where , , are obtained from a straightforward diagonalization ofH , each two-fermion expectation value takes the form While the ground state |Ω contains no bond fermions for ℎ = 0, corresponding to = 1, the hybridization between the bond fermions and the matter fermions leads to a finite density of bond fermions for ℎ > 0, which corresponds to < 1.

III. COEFFICIENTS OF THE EFFECTIVE FIELD THEORY
Here we provide the coefficients of the effective field theory in Eq. (11) of the main text. These coefficients can be computed by projecting the effective Hamiltonian in Eq. (4) of the main text to the two low-energy fermion bands corresponding to Eq. (9) of the main text. For concreteness, the momentum ≡ ( , ) is described in Cartesian coordinates defined by the unit vectors e r −r and e r [see Fig. 1 of the main text for definitions ofr ], while the length unit is taken as the lattice constant (i.e., the distance between two neighboring sites).

IV. NONANALYTIC BEHAVIOR OF THE GROUND-STATE ENERGY
Here we analyze the nonanalytic behavior of the ground-state energy (ℎ) at the critical field ℎ = ℎ . Specifically, we show that the second derivative, = 2 / ℎ 2 , is discontinuous at ℎ = ℎ , and provide an expression for its discontinuity, Δ , in terms of the effective field theory [see Eqs. (8) and (11) of the main text]. Given the infrared nature of the singularity, it is useful to write the ground-state energy as a sum of two contributions, =˜ (Λ) + (Λ), which correspond to the long-wavelength modes with momentum ≡ |k| < Λ and the remaining modes with momentum > Λ, respectively: where Ω = 8 2 / √ 3 is the area of the Brillouin zone (in units of −2 ), Λ 1 is an arbitrary cutoff, and k is the energy of the mode at momentum k [see Eq. (12) of the main text]. The second contribution (Λ) is analytic at ℎ = ℎ because all of its derivatives are well defined. Therefore, In other words, the discontinuity in at the critical field can be completely extracted from the first contribution in the Λ → 0 limit (i.e., the effective field theory) because it does not depend on the cutoff Λ.