Analytic calculation of the vison gap in the Kitaev spin liquid

TL;DR

This paper derives an analytical expression for the vison gap in the Kitaev spin liquid by calculating the phase shift of Majorana fermions scattering off a bond-flip, providing a precise value and potential for broader lattice applications.

Contribution

It introduces an analytical method to compute the vison gap in the Kitaev spin liquid, improving upon previous numerical approaches.

Findings

Analytical expression for the vison gap: Δ = 0.2633J

Phase shift calculation for Majorana fermions

Applicability to complex lattices like the hyper-octagonal lattice

Abstract

Although the ground-state energy of the Kitaev spin liquid can be calculated exactly, the associated vison gap energy has to date only been calculated numerically from finite size diagonalization. Here we show that the phase shift for scattering Majorana fermions off a single bond-flip can be calculated analytically, leading to a closed-form expression for the vison gap energy $\Delta = 0.2633J$. Generalizations of our approach can be applied to Kitaev spin liquids on more complex lattices such as the three dimensional hyper-octagonal lattice.

Figures

Equations

$$H_K=\sum _{<ij>}{J}_{{\alpha }_{ij}}{\sigma }_i^{{\alpha }_{ij}}{\sigma }_j^{{\alpha }_{ij}},$$

$$H_{KSL}=2\sum _{<ij>}{J}_{{\alpha }_{ij}}{\hat{u}}_{ij}(i{c}_i{c}_j),$$

$${\mathcal{W}}_p=\prod _{<i,j>\in p}{u}_{ij}(i\in A,j\in B).$$

$$H_0=H_{KSL}[u_{ij}\to 1].$$

$$H_0=\frac{1}{2}\sum _{\mathbf{k}\in \mathbf{BZ}}{\psi }_{\mathbf{k}}^{†}({\vec{\gamma }}_{\mathbf{k}}\cdot \vec{\tau }){\psi }_{\mathbf{k}},$$

$${\psi_{\bf k}}=\frac{1}{\sqrt{N_{c}}}\sum_{j}\left(\begin{array}[]{cc}c_{j,A}\cr c_{j,B}\end{array}\right)e^{-i{\bf{k}}\cdot{\bf R}_{j}}$$

$$\begin{array}{rl}{\gamma }_{\mathbf{k}}& =2i(J_z+J_x{e}^{i{k}_1}+J_y{e}^{i{k}_2}),\\ \mathbf{k}& =\frac{k_1}{2\pi }{\mathbf{b}}_1+\frac{k_2}{2\pi }{\mathbf{b}}_2,k_1,k_2\in [0,2\pi ].\end{array}$$

$$H_{KSL+2v}=H_0+V,$$

$$\hat{V}=-4J_z(i{c}_{0,A}{c}_{0,B})$$

$$H_{KSL+2v}=\frac{1}{2}\sum _{\mathbf{k}\in \mathbf{BZ}}{\psi }_{\mathbf{k}}^{†}({\vec{\gamma }}_{\mathbf{k}}\cdot \vec{\tau }){\psi }_{\mathbf{k}}+\frac{1}{2}{\mathbf{c}}_0^T(V{\tau }_2){\mathbf{c}}_0,$$

$${\bf c}_{0}=\left(\begin{array}[]{cc}c_{0,A}\cr c_{0,B}\end{array}\right)=\frac{1}{\sqrt{N_{c}}}\sum_{{\bf k}\in{\bf BZ}}{\psi_{{\bf{k}}}}$$

$$G_0(i{\omega }_n,\mathbf{k})=[i{\omega }_n-{\vec{\gamma }}_{\mathbf{k}}\cdot \vec{\tau }{]}^{-1}.$$

$$\beta F=-\frac{1}{2}Tr[\ln (-G^{-1})]=-\frac{1}{2}Tr\ln [-G_0^{-1}+\hat{V}],$$

$$\displaystyle\Delta F=-\frac{1}{2\beta}{\hbox{Tr}}[\ln(1-\hat{V}G_{0})]=\frac{1}{2\beta}\sum_{r=1}^{\infty}\frac{1}{r}{\rm Tr}\bigl{[}(\hat{V}G_{0})^{r}\bigr{]}.$$

$$\Delta F=\frac{1}{2\beta}\sum_{i\omega_{n}}\sum_{r=1}^{\infty}\frac{1}{r}{\rm tr}\biggl{[}\left(\frac{V\tau_{2}}{N_{c}}\sum_{{\bf{k}}}G_{0}(i\omega_{n},{\bf{k}})\right)^{r}\biggr{]},$$

$$g(z)=\frac{1}{N_c}\sum _{\mathbf{k}\in \mathbf{BZ}}{G}_0(z,\mathbf{k}),$$

$$\Delta F$$

$$=$$

$${\gamma }_c(\mathbf{k})$$

$${\gamma }_s(\mathbf{k})$$

$$g(z)=\frac{1}{N_c}\sum _{\mathbf{k}\in \mathbf{BZ}}\frac{z-({\gamma }_c(\mathbf{k}){\tau }_2+{\gamma }_s(\mathbf{k}){\tau }_1)}{z^2-|{\gamma }_{\mathbf{k}}{|}^2}.$$

$$\begin{array}{rl}1-\hat{V}g(z)& =1-\frac{2V}{N_c}\sum _{\mathbf{k}\in \mathbf{BZ}}\frac{z{\tau }_2-{\gamma }_c(\mathbf{k})}{z^2-|{\gamma }_{\mathbf{k}}{|}^2}\\ & =1-V({\tau }_2{g}_0(z)-{\mathbb{I}}_2{g}_2(z))\end{array}$$

$$g_0(z)$$

$$g_2(z)$$

$$\Delta F$$

$$\Delta F={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\left(\frac{1}{2}-f(\omega )\right){\delta }_v(\omega ),$$

$$\delta_{v}(\omega)={\rm Im}\ \ln\biggl{[}(1+2Kg_{2}(z))^{2}-(2Kg_{0}(z))^{2}\biggr{]}_{z=\omega-i\delta}$$

$${\Delta }_v$$

$$\begin{array}{r}{g}_0(z)={\int }_0^{2\pi }\frac{d{k}_1}{2\pi }{\int }_0^{2\pi }\frac{d{k}_2}{2\pi }\frac{z}{z^2-|{\gamma }_{\mathbf{k}}{|}^2},\\ g_2(z)={\int }_0^{2\pi }\frac{d{k}_1}{2\pi }{\int }_0^{2\pi }\frac{d{k}_2}{2\pi }\frac{{\gamma }_c(\mathbf{k})}{z^2-|{\gamma }_{\mathbf{k}}{|}^2},\end{array}$$

$$g_0(z)={\int }_0^{2\pi }\frac{dk}{2\pi }\frac{z}{(z^2-(3+2c))\sqrt{1-\frac{8(c+1)}{(z^2-(3+2c){)}^2}}},$$

$$g_2(z)={\int }_0^{2\pi }\frac{dk}{2\pi }\frac{2c+1}{(z^2-(3+2c))\sqrt{1-\frac{8(c+1)}{(z^2-(3+2c){)}^2}}},$$

$$\Delta \rho (\omega )=\frac{1}{2\pi }\frac{d{\delta }_v}{d\omega }$$

Full text

I Introduction

Kitaev spin liquids (KSL) are a class of exactly solvable quantum spin liquid that exhibit spin fractionalization, anyonic excitations and long-range entanglementKitaev (2006); Hermanns et al. (2018); O’Brien et al. (2016); Trebst (2017); Takagi et al. (2019). The fractionalization of spins into Majorana fermions is accompanied by the formation of emergent $\mathbb{Z}_{2}$ gauge fields, giving rise to $\mathbb{Z}_{2}$ vortex excitations or “visons”. These excitations are gapped, and the energy cost associated with creating two visons on adjacent plaquettes is called the vison gap $\Delta_{v}$ (Fig[1]). Proposals for the practical realization of Kitaev spin liquids in quantum materials, including $\alpha$-RuCl3 Banerjee et al. (2016); Do et al. (2017); Jackeli and Khaliullin (2009); Kasahara et al. (2018); Liu et al. (2020); Takagi et al. (2019); Winter et al. (2017); Wolter et al. (2017); Wulferding et al. (2020); Yamada (2020) and Iridates Chaloupka et al. (2010); Kitagawa et al. (2018); Winter et al. (2017) have renewed interest in the thermodynamics of Kitaev spin liquid Eschmann et al. (2019); Feng et al. (2020); Joy and Rosch (2022); Kato et al. (2017); Li et al. (2020); Nasu et al. (2015); Motome and Nasu (2020); Udagawa and Moessner (2019). The extension of these ideas to Yao-Lee spin liquid Yao and Kivelson (2007); Yao and Lee (2011) and its application to Kondo models, Coleman et al. (2022); Tsvelik and Coleman (2022) motivate the development of an analytical approach to calculate the vison gap $\Delta_{v}$.

The vison gap in KSLs has to date, been determined by numerical diagonalization of finite size systems O’Brien et al. (2016); Kitaev (2006). Here we present a Green’s function approach for the analytical computation of the vison gap $\Delta_{v}$ from the scattering phase shift associated with a a $\mathbb{Z}_{2}$ bond-flip. Our work builds on theoretical developments in the field of Kitaev spin liquids which relate to the interplay between Majorana fermions and visonsKitaev (2006); Baskaran et al. (2007); Knolle et al. (2014); Lahtinen (2011); Lahtinen et al. (2014, 2008, 2012); Théveniaut and Vojta (2017); Kitaev (2006); Joy and Rosch (2022); Kao et al. (2021); Knolle et al. (2014); Lahtinen (2011); Lahtinen et al. (2014, 2008, 2012); Théveniaut and Vojta (2017); Zhang et al. (2019). Using exact calculations, we find the vison gap energy of $\Delta_{v}=0.263313(6)J$ for the Kitaev spin liquid on honeycomb lattice in the gapless phase, extending the accuracy of previous calculations Kitaev (2006); O’Brien et al. (2016). Our calculations reveal the formation of Majorana resonances in the density of states which accompany the formation of two adjacent visons. Our approach can be simply generalized to more complex lattices and are immediately generalizable to Yao-Lee spin liquids.

II Vison Gap in the Kitaev Honeycomb Model

The Kitaev honeycomb lattice model Kitaev (2006) is described by the Hamiltonian

[TABLE]

where the Heisenberg spins $\vec{\sigma}_{i}=(\sigma^{x}_{i},\sigma^{y}_{i},\sigma^{z}_{i})$ at site $i$ interact with their nearest neigbors via an Ising coupling between the $\alpha_{ij}=x,y,z$ spin components, along the coresponding bond directions $\langle ij\rangle$, with strength $J_{\alpha_{ij}}$, as shown in Fig[1]. An exact solution of Kitaev Model Kitaev (2006) is found by representing the spins as products of Majorana fermions, $\sigma^{\alpha}_{j}=2c_{j}b^{\alpha}_{j}$ which satisfy canonical anticommutation algebras, $\{c_{i},b^{\alpha}_{j}\}=0$, $\{b^{\alpha}_{i},b^{\beta}_{j}\}=\delta_{ij}\delta^{\alpha,\beta}$, (taking the convention that $c_{j}^{2}=(b_{j}^{\alpha})^{2}=1/2$). The system is projected into the physical subspace by selecting $\mathcal{D}_{j}\equiv-4ic_{j}b^{x}_{j}b^{y}_{j}b^{z}_{j}=1$ at each site, allowing the Hamiltonian (1) to be rewritten as $\mathbb{Z}_{2}$ gauge theory

[TABLE]

where the gauge fields $\hat{u}_{ij}=2ib^{\alpha_{ij}}_{i}b^{\alpha_{ij}}_{j}=\pm 1$ on bond $ij$ commute with the Hamiltonian, $[\hat{u}_{ij},H_{K}]=0$. The plaquette operators $\mathcal{W}_{p}$

[TABLE]

formed from the product of gauge fields $\hat{u}_{ij}$ around the hexagonal loop $p$ ( plaquette), are gauge invariant and also commute with the Hamiltonian $[\mathcal{W}_{p},H_{K}]=0$ and constraint operators $[\mathcal{W}_{p},\mathcal{D}_{j}]=0$, giving rise to a set of static constants of motion which take values $\mathcal{W}_{p}=\pm 1$. Each eigenstate is characterized by the configurations of $\{\mathcal{W}_{p}\}$; Lieb’s theorem Lieb (1994) specifies that the ground state configuration is flux-free, i.e. $\mathcal{W}_{p}=1$ for all hexagons $p$. In what follows we will choose the gauge $\hat{u}_{ij}=1$ when $i\in A$ and $j\in B$ sublattice, assigning

[TABLE]

Rewriting $H_{0}$ in momentum space, we obtain

[TABLE]

where

[TABLE]

creates a Majorana in momentum space, where $N_{c}$ is the number of unit cells and ${\bf R}_{j}$ is the location of the unit cell and $\vec{\gamma}_{{\bf k}}=(Re({\gamma_{\bf k}}),-Im(\gamma_{\bf k}))$ is expressed in terms of the form factor

[TABLE]

Here we have employed a reciprocal lattice basis ${\bf b_{1},b_{2}}$ to span the momentum ${\bf k}\in{\bf BZ}$, which transforms to rhombus shaped Brillouin zone in the reciprocal lattice (see Fig. 2 ). The Majorana excitation spectrum of the Kitaev spin liquid is given by the eigenvalues of $H_{0}$, $\epsilon_{\bf k}=\pm|\gamma_{\bf k}|$.

We create two adjacent visons by flipping the gauge field in the unit cell at origin to $\hat{u}_{(0,A)(0,B)}=-1$ as shown in Fig [1], resulting in the following Hamiltonian:

[TABLE]

where

[TABLE]

acts as a scattering term for majoranas in the bulk. In this way, the vison gap calculation is formulated as a scattering problem.

For this case, the Hamiltonian is given by

[TABLE]

[TABLE]

creates a Majorana fermion at the origin and $V=4J_{z}$.

We now set up the scattering problem in terms of Green’s functions. The Green’s function of the unscattered majoranas is $G_{0}=G_{0}(i\omega_{n},{\bf k})\delta_{{\bf k},{\bf k^{\prime}}}$, where

[TABLE]

In the presence of the bond-flip at the origin, the Green’s function of the scattered majoranas is given by $G=(G_{0}^{-1}-\hat{V})^{-1}$, where $\hat{V}_{{\bf k},{\bf k^{\prime}}}=(V\tau_{2})/N_{c}$ is the scattering matrix. The total free energy of the non-interacting ground-state in the presence of the scattering is given by the standard formula

[TABLE]

where ${\rm Tr}$ denotes the full trace over Matsubara frequencies, momenta and sublattice degrees of freedom. The change in free energy is then given by

[TABLE]

We now carry out the trace over the Matsubara frequencies and momenta, so that

[TABLE]

where {\rm tr}\bigl{[}\ \ \bigr{]} denotes the residual trace over sublattice degrees of freedom. Now, we can incorporate the summations over momentum by introducing the local Green’s function

[TABLE]

so that

[TABLE]

where we have re-assembled the Taylor series as a logarithm.

We shall illustrate our method for the isotropic case $J_{x}=J_{y}=J_{z}=J$, setting $K=2J$ and $V=4J$. In this case, $\gamma_{\bf k}=iK(1+e^{ik_{1}}+e^{ik_{2}})$. If we divide $\gamma_{{\bf{k}}}=i(\gamma_{c}({\bf{k}})+i\gamma_{s}({\bf{k}}))$ into its even and odd components

[TABLE]

then $G_{0}(i\omega_{n})$ can be rewritten as

[TABLE]

The odd component $\gamma_{s}(\bf k)$ vanishes under momentum summation so that

[TABLE]

where

[TABLE]

Carrying out the trace in the free energy we then obtain

[TABLE]

The Matsubara summation can then be carried out as an anti-clockwise contour integral around the imaginary axis weighted by Fermi function, $f(z)=[e^{\beta z}+1]^{-1}$. Deforming the contour to run clockwise around the real axis we obtain

[TABLE]

where

[TABLE]

is identified as the scattering phase shift. Note that $\delta_{v}(\omega)=-\delta_{v}(-\omega)$ is an antisymmetric function of frequency. At zero temperature the vison gap is then

[TABLE]

where we have rescaled the frequency in units of $K$, setting $z=\omega/K$. In the reciprocal basis

[TABLE]

where we have set $K=1$ in $\gamma({\bf{k}})$, i.e $\gamma_{{\bf{k}}}=1+e^{ik_{1}}+e^{ik_{2}}$ and $\gamma_{c}=\cos(k_{1})+\cos(k_{2})$. The interior integral over $k_{2}$ can be carried out as a complex contour integral over $w=e^{ik_{2}}$ around the unit circle, (Appendix A), giving

[TABLE]

[TABLE]

where $c\equiv\cos(k)$. These integrals were evaluated numerically, to obtain the phase shift $\delta_{v}(\omega)$ (Fig[4]). The phase shift was interpolated over over a discrete set of $N$ points and the integral (28) was carried out numerically on the interpolated phase shift. By extrapolating the limit $1/N\rightarrow 0$, we find the vison gap energy to be $\Delta_{v}=0.1311656(3)K=0.263313(6)J$ for the isotropic case $J_{z}=J_{y}=J_{z}=J$.

This analytically-based calculation improves on the earlier result obtained via numerical diagonalization of finite size systems Kitaev (2006) i.e. $\Delta_{v}\approx 0.267J$. Its main virtue however, is that the method can be easily generalized, and we gain insights from the calculated scattering phase shifts.

From the calculated phase shift, we can calculate the change in density of states (DOS)

[TABLE]

(Fig. 4 c.) associated with a Bond flip, which is seen to contain a resonance centered around $\epsilon_{0}\approx\pm 0.07K$. This resonance can be examined in detail by expanding $g_{0}(z)$ and $g_{2}(z)$ for small $z$:

[TABLE]

Which can be used to evaluate scattering phase shift $\delta_{v}(\omega)$ (27), and the resonant DOS change $\Delta\rho(\omega)$ (32) analytically. The position of the resonance is determined by the integration over the entire band but its width is determined by the density of states at low energies. Since the DOS vanishes inside the spectral gap, the resonance may become sharp in the topological state. The sharp peak in the gapped state signifies the binding of Majorana fermions to the visons formed by $\mathbb{Z}_{2}$ bond flip at origin.

III Discussion

In this work we have presented an analytical method for determination of the vison gap by treating the flipping of the $\mathbb{Z}_{2}$ gauge field as a scattering potential for the Majorana fermions. In this way, we have been able to analytically extend the numerical treatment by Kitaev for the isotropic model on honeycomb lattice Kitaev (2006) to obtain an analytic result for the vison gap energy $\Delta_{v}$.

A key part of our approach is the calculation of the Majorana phase shift for scattering off the bond-flipped configuration. One of the interesting observations is that the scattering contains a Majorana bound-state resonance, located at an energy $\epsilon_{0}\approx\pm 0.07K$. Since this bound-state is formed from scattering throughout the entire Brillouin zone, its location is expected to be quite robust. Thus in those cases where the excitation spectrum acquires a gap, eg through time-reversal symmetry breaking Kitaev (2006); Haldane (1988), we expect this resonance to transform into a sharp in-gap excitation.

While it is possible to extend our method to analytically calculate the energy associated with anyons by flipping $x-x$ bonds along ${\bf{a_{1}}}$ direction, a much simpler derivation of the anyon energy in the KSL can be made by taking making two copies of the KSL, forming a complex fermion Hamiltonian $H_{c}=H_{KSL}+H_{KSL}$. The line of reverse bonds around the torus can then be absorbed by a unitary transformation that redistributes the odd boundary condition into an effective vector potential that shifts all the momenta ${\bf k}=(k_{1},k_{2})\rightarrow(k_{1}+\frac{\pi}{L},k_{2})$, equivalent to introducing a half magnetic flux with vector potential $A_{x}=\frac{\pi}{L}$. Treating the response to the vector potential in an analogous fashion to a superconductor, the putative the energy cost of an anyon would be

[TABLE]

where $\rho_{s}$ is the superfluid stiffness associated with the ground-state, $A=\pi/L$ is the vector potential and the factor of $4$ derives from halving the energy of the complex fermion system. However, since the complex fermion Hamiltonian $H_{c}$ preserves the global $U(1)$ symmetry, its superfluid stiffness $\rho_{s}$ vanishes so it costs no energy to create anyons in the gapless state. From this line of reasoning, we can see that the ground state of the Kitaev spin liquid has a four-fold degeneracy and is topologically ordered.

Finally, we note that our method also admits various generalizations. For example, it can be extended to anisotropic couplings i.e. $J_{x}\neq J_{y}\neq J_{z}$ as well as to higher dimensions, such as the three-dimensional hyperoctagonal lattice. Moreover, our method can be applied to study the impact of spinor order formation as a consequence of hybridization between conduction electrons and Majorana spinons in the CPT model for a Kondo lattice coupled to a Yao-Lee spin liquid Tsvelik and Coleman (2022); Coleman et al. (2022). This allows us to study the stability of Yao-Lee spin liquid against spinor order formation, the subject of a forthcoming article by the authors.

Acknowledgements.

This work was supported by Office of Basic Energy Sciences, Material Sciences and Engineering Division, U.S. Department of Energy (DOE) under Contracts No. DE-SC0012704 (AMT) and DE-FG02-99ER45790 (AP and PC ). All authors contributed equally to this work.

Appendix A Analytic Calculation of Green’s Function in Honeycomb Lattice

Here we show how to simplify the integrals

[TABLE]

where $\gamma_{c}({{\bf{k}}})=1+\cos(k_{1})+\cos(k_{2})$, using a contour integral. We begin by noting that the integrals over $k_{1}$ and $k_{2}$ can be carried out in either order, allowing us to pull the cosines in $\gamma_{c}(k)$ out of the inner integral, so that

[TABLE]

where

[TABLE]

Writing $s=e^{ik_{1}}$ and $w=e^{ik_{2}}$, we can rewrite $I_{0}$ as a counter-clockwise integral around the unit circle $|w|=1$,

[TABLE]

Rewriting the denominator as a quadratic function of $w$,

[TABLE]

where

[TABLE]

We can thus write the integral in the form

[TABLE]

where

[TABLE]

are the poles of the integrand.

Now since $w_{+}w_{-}=s=e^{ik_{1}}$, it follows that $|w_{+}w_{-}|=1$, so that only one of these poles lies inside the contour. (In general, this may depend on the way we treat the branch cuts inside the square root of (43). However, we don’t actually need to know which pole it is, as this we will fix the sign and the branch-cuts in the final expression by demanding that the asymptotic behavior of $I_{0}\sim 1/z^{2}$ is analytic at large $z$.) Lets assume that the pole closest to the origin is at $w=w_{-}$, then we obtain

[TABLE]

Now expanding the denominator, we have

[TABLE]

where we have factorized the final expression, to guarantee that at large $z$, $I_{0}(z,s)\sim 1/z^{2}$ is analytic. Combining the above results, gives us the following expressions for $g_{0}(z)$ and $g_{2}(z)$

[TABLE]

Where $c\equiv\cos(k)$, which are the expressions given in (30) and (31).