Heat transport in the anisotropic Kitaev spin liquid
Abstract
We present a study of longitudinal thermal transport in the Kitaev spin model on the honeycomb lattice, focusing on the role of anisotropic exchange to cover both, gapless and gapped phases. Employing a complementary combination of exact diagonalization on small systems and an average gauge configuration approach for up to spinful sites, we report results for the dynamical energy current autocorrelation function as well as the dc thermal conductivity over a wide range of temperatures and exchange anisotropies. Despite a pseudogap in the current correlation spectra, induced by emergent thermal gauge disorder on any finite system, we find that in the thermodynamic limit, gapless and gapful phases both feature normal dissipative transport, with a temperature dependence crossing over from power law to exponentially activated behavior upon gap opening.
I Introduction
Quantum spin liquids (QSL) represent a rare form of magnetic matter, with fluctuations strong enough to suppress the formation of local order parameters even at zero temperature. QSLs typically result from frustration of magnetic exchange and tend to exhibit many peculiar properties, such as massive entanglement, quantum orders and fractionalized excitations Balents2010 (); Balents2016 (). Among the many models proposed for putative QSLs, Kitaev’s compass exchange Hamiltonian on the honeycomb lattice stands out as one of the few, in which a QSL can exactly be shown to exist Kitaev2006 (). In fact, pairs of spins in this model fractionalize in terms of mobile Majorana fermions coupled to static gauge fields Kitaev2006 (); Feng2007 (); Chen2008 (); Nussinov2009 (); Mandal2012 (). Solid state realizations of Kitaev’s model are based on Mottinsulators with strong spinorbit coupling (SOC) Khaliullin2005 (); Jackeli2009 (); Chaloupka2010 (); Nussinov2015 (), although residual nonKitaev exchange interactions remain an issue at low temperatures and low energies, where most of the present systems still tend to order magnetically Trebst2017 ().
Mobile Majorana matter has been suggested to leave several fingerprints in spectroscopic measurements, like inelastic neutron Banerjee2016 (); Banerjee2016a (); Banerjee2018 (), Raman scattering Knolle2014 (), as well as in local resonance probes Baek2017 (); Zheng2017 (). Apart from spectroscopy, thermal transport is yet another powerful tool, able to discriminate microscopic models of elementary excitations and their scattering mechanisms in quantum magnets  see Ref. Hess2018a () for a review. Among the candidates for Kitaev QSLs,  Plumb2014 () has been under intense scrutiny regarding measurements of the longitudinal thermal conductivity Hirobe2017 (); Leahy2017 (); Hentrich2018 (); Yu2018 (). At first, heat transport by itinerant spin excitations had been inferred from an anomaly in the inplane longitudinal heat conductivity at around 100 K Hirobe2017 (), and from a magnetic fieldinduced lowtemperature enhancement of for fields T parallel to the material’s honeycomb planes Leahy2017 (). These interpretations however, have recently been ruled out by results for the outofplane heat conductivity where both types of anomalies are present as well Hentrich2018 (). The emergent picture for rationalizing the longitudinal heat transport in  material is thus, that its most salient features can be explained by phononic heat transport which is likely to be renormalized by scattering from Majorana matter Hentrich2018 (); Yu2018 (). However, it remains to be settled whether also a subleading direct magnetic contribution to the heat transport is present, masked by the phononic transport. Very recently also transverse thermal transport in magnetic fields, i.e. thermal Hall conductivity , and its potential quantization has been observed in  Kasahara2018 (). This may be a rather strong evidence for Kitaev physics and chiral Majorana edge modes in this material.
Theoretically, thermal transport studies in pure Kitaev QSLs have been performed using quantum Monte Carlo in 2D Motome2017 () or applying ED in 1D and 2D Steinigeweg2016 (); Metavitsiadis2016 (); Metavitsiadis2017 (). Thermal transport in KitaevHeisenberg ladders has been considered recently, employing ED and quantum typicality Metavitsiadis2018a (). Magnetically ordered phases of a KitaevHeisenberg model have also been investigated for transport using spin wave calculations Stamokostas2017 (). Remarkably, in the 2D QSL case, all of these studies have been confined to the point of Csymmetrically sized, isotropic compass exchange, i.e. the gapless phase of the Kitaev model. The main purpose of this work is to step beyond this limitation and investigate the longitudinal thermal conductivity, covering also anisotropic Kitaev exchange, ranging from the gapless to the gapped case. Exploring both phases, we find that regardless of the exchange coupling regime, the system supports a finite dc transport coefficient at all nonzero temperatures investigated, however with the emergence of thermal activation behavior in the gapped phase.
The outline of this paper is as follows: in Sec. II we briefly recall the properties of the Kitaev model, along with basic ingredients of thermal transport calculations, as used in this work. In Sec. III we detail our results for the thermal conductivity versus the anisotropy of the exchange coupling, as obtained using the average gauge configurations (AGC). The method is explained in Sec. III.1, which is then applied to study the dynamical heat current autocorrelation function in Sec. III.2, and the dc heat conductivity in Sec. III.3. We conclude in Sec. IV. Appendix A relates our findings to transportphenomenology.
Ii Model
We study the longitudinal dynamical thermal conductivity of the Kitaev model in the presence of anisotropy, introduced by tuning the relative exchange couplings and considering both, gapless and gapped regimes. The Hamiltonian is
(1) 
where indicates nearestneighbors on a 2D honeycomb lattice (HL) of 2 sites, with and the exchange couplings and the Pauli matrices, respectively, for coordinates , as shown in Fig. 1(a). Model (1) can be mapped onto a sum of Hamiltonians of noninteracting itinerant Majorana fermions, each of which is classified by a fixed, distinct distribution of values of static gauge fields , residing on, e.g., the blue bonds in Fig. 1(a). The Majorana fermions of each gauge sector can be rewritten in terms of spinless noninteracting fermions with occupation numbers 0,1. This implies a total of states of gauges and spinless fermions, consistent with the dimension of the original spin Hilbert space. Each of the fermionic subHamiltonians is a QSL with spin correlations not exceeding nearest neighbor distance. Several routes have been documented to achieve the aforementioned mapping and we refer to Refs. Kitaev2006 (); Feng2007 (); Chen2008 (); Nussinov2009 (); Mandal2012 () for details. In this language (1) reads
(2)  
Here, for convenience, the bonds, which feature only an dependent onbond potential for the fermions are shrunk onto just a point, transforming the HL into a square lattice (SL), with primitive vectors as depicted in Fig. 1(b). We will refer the SL geometry throughout this paper. To quantify the anisotropy of the Kitaev exchange, we employ a parameter , such that , and . Tuning allows to keep the fermionic energy scale constant, while accessing both, the gapless and the gapped regime of the zero temperature spectrum of (2), for and , respectively Kitaev2006 ().
Our discussion of the thermal conductivity is based on calculations within linear response theory, assuming the long wave length limit . For this, thermal transport coefficients are computed from the energy current autocorrelation functions , where is the energy current and refers to the thermal trace FHM2007 (). The real part of thermal conductivity, i.e. , follows from the Fourier transform
(3)  
(4) 
Generically, . The Drude weight (DW) refers to a ballistic channel, and if finite marks a perfect conductor. encodes the dissipative parts of the heat flow. If is finite, the system is a normal heat conductor with a finite DC conductivity. Interacting quantum systems with are rare, see e.g. Ref. FHM2007 (), and for has also been argued for the Kitaev model Metavitsiadis2016 (); Metavitsiadis2017 (). In turn, a prime goal of this work will be the analysis of and its limit for in order to determine the DC conductivity.
To derive the heat current operator in (3) one has to realize that most of the gauge sectors in (2) are not translationally invariant, due to the real space distribution of . In turn, we employ the polarization operator Metavitsiadis2017 (), to obtain the current operator in real space
(5) 
It is worth noting, that in a translationally invariant system from (5) is exactly identical to that obtained from the continuity equation in the limit of . On the SL we get
(6) 
with , and for .
For the remainder of this work we resort to numerical methods to treat (3). For this we note, that not only the Hamiltonian, but also the current operator does not mediate transitions between gauge sectors. In turn, any thermal trace can be decomposed into a classical trace over gauges and a trace over the fermions, separately. For the latter, and for each fixed gauge sector , we define a component spinor of the fermions, with indices referring to the sites on the SL. With this the Hamiltonian and the current for each gauge sector read , and . Using a Bogoliubov transformation onto quasiparticle fermions via , the Hamiltonian is diagonalized as
(7) 
with quasiparticle energies . The contribution of each gauge sector to (3) can then be obtained straightforwardly by a Wickdecomposition of in the quasiparticle basis Metavitsiadis2016 ().
Iii Results
In this section we will present our results for the dynamical energy current autocorrelation function and the dclimit of the heat conductivity versus anisotropy at finite temperatures.
iii.1 Average gauge configuration method
Tracing over the gauge sectors in (3) can be done in various ways and
to various degrees of precision. Among the methods are Markovian sampling of the
gauges by QMC Motome2015b (); Motome2017 (), exact summation over all gauges,
equivalent to ED of the spin model Metavitsiadis2016 (); Metavitsiadis2017 (),
and summing only over a dominant subset of gauge configurations, i.e. the
average gauge configuration (AGC) method Metavitsiadis2016 (); Metavitsiadis2017 (), which will be used here. The following sketches the main ideas
of the AGC method, while for more details we refer the reader to Refs. Metavitsiadis2016 (); Metavitsiadis2017 (). The notion of an average
gauge configuration is trivial in the low and the hightemperature limit. In
the former, the gauges assume a fluxfree configuration Kitaev2006 (), and
in the latter a typical gauge configuration is fully random. These limiting
cases are separated by a temperature scale at which
fluxproliferation occurs upon increasing temperature. Here refers to
the anisotropy parameter. scales with the energy to excite a
single flux, i.e. the flux gap Kitaev2006 (). Because
of the latter inequality, and the crossover regime can be read
off from the specific heat, i.e. from the entropy release of the gauge fields
Motome2015b (), since this is well separated from the fermionic energy
scales . Our strategy will be to first determine and then
to confine all evaluations of to temperatures above this scale, such
that averaging over only completely random gauge states is sufficiently
justified for the gauge trace. Consequently, the remaining effective fermionic
model is a binary disorder problem. It turns out that the temperature range
accessible by this approach is large enough for our purposes.
To approximate we now calculate the specific heat noteSpec (). This is done exactly on finite systems, tracing over both, fermions and gauges. The system sizes we can reach for this, namely , i.e. 72 spins, are larger than those for exact thermodynamics using the many body spin basis, however, still much smaller than those which will be used within the AGC subsequently. Fig. 2 shows the temperature dependence of the specific heat for three different values of anisotropy. As expected, two peaks develop at two different temperature scales, corresponding to the release of entropy from the gauge degrees of freedom and the mobile Majorana fermions, at and , respectively. Obviously, the global fermionic energy scale is only weakly affected by the gauge disorder, and therefore the high peak remains rather insensitive to . However the low peak at is strongly shifted to smaller temperatures, by roughly two orders of magnitude, upon decreasing from 1 to 0.5. This finding is consistent with results from QMC Motome2015b (). While in principle sets an individual temperature scale for each , above which averaging over fully random gauge configurations provides a sufficient approximation for the gauge trace, we will use the maximum of these, i.e. as the lowest temperature to apply the AGC for all subsequent calculations, independent of .
iii.2 Dynamical current correlations
Now we turn to the dynamical current correlation function . Results on up to 6060sites systems are displayed in Fig. 3 for two representative anisotropy values . These place the system into the gapless(ful) phase. Two temperatures, i.e. are shown, which allow to highlight the difference between low and high matter fermion densities. We begin with Fig. 3(a). Here the fermion occupation number at all momenta is small and the intensity is essentially due to twofermion pairbreaking contributions to . For , the oneparticle dispersion displays a gap , at Motome2015b (), with a corresponding twoparticle excitation gap in at . Moving to Fig. 3(b) the matterfermion density is large and direct transport by thermally populated quasiparticle states contributes. In a clean system this transport would show up as a Drude peak strictly at . Due to the gauge disorder however, the Drude peak is smeared over a substantial range of finite frequencies . Interestingly, the randomness provided for by the gauge field is not sufficient to fill in the gap completely, as is obvious from panel (a) and (b). Related observations have been made in Refs. Motome2015b (); Motome2016 (). The ’smeared Drude peak’ in Fig. 3(b) features a very narrow zerofrequency dip, which requires careful finite size analysis. Examples of such analysis are shown in Fig. 3(c,d) for the case of . The main point can be read off best from panel (d), where it is amply clear, that as the zerofrequency dip closes in onto the yaxis. This observation renders the system a normal, dissipative heat conductor in the thermodynamic limit and extends identical findings from Ref. Metavitsiadis2017 () to the case of . Note, that there is a difference of several orders of magnitude between in Figs. 3(c) and 3(d). This is related to the thermal activation behavior in the gapful case, to be discussed in the next section.
iii.3 DClimit of heat conductivity
Based on a series of calculations as in Fig. 3, we have extrapolated the dclimit over a wide range of and , sufficient to acquire a consistent picture of the thermal conductivity both, in the gapless and the gapful regime. This result is depicted in Fig. 4. To arrive at this, the brown dashed lines in Fig. 3(c,d) exemplarily indicate the extrapolation procedure performed, for each of the largest systems, i.e. for . It consists of leastsquare fitting of the low range of the correlation function for each pair of by second order polynomials, from which the dc value is extracted. We have chosen a range of to determine the coefficients of the polynomial in all cases. To obtain some form of finitesize scaling, we have compiled the results of the fitting procedure for two system sizes, i.e. and , into Fig. 4, showing dashed (solid) curves with empty (solid) circles, for . Although the linear dimension of the systems shown differs by almost a factor of 4, displays only minimal differences for most of the temperatures. Only around , a small, systematic correction towards lower values with increasing system size is observed for every . This implies visible, but satisfyingly little finitesize effects, presumably putting the largest system sizes rather close to the thermodynamic limit.
Fig. 4 is a main result of this paper. It extends previous analysis of the thermal conductivity of the Kitaev model Motome2017 (); Metavitsiadis2017 () into the anisotropic regime. It demonstrates a qualitative change of the temperature dependence of , as the system crosses over from the gapless to the gapful spectrum. In the former the thermal conductivity roughly scales with a power law, for , while in the latter an exponential behavior with provides a reasonable description. This is in line with the notion of quasiparticle transport, since the activation gap is consistent with the oneparticle gap, e.g. for as discussed in Sec. III.2.
We would like to emphasize, that while the trend of the curves in Fig. 4 might suggest a strictly vanishing , e.g. for and , this is not so, based on our results. Rather, the yaxis scale in Fig. 3(c) is very small, albeit finite. Such small orders of magnitude are consistent with the rapid low temperature suppression by the exponential activation. In turn, our findings suggest, that the Kitaev model is a normal dissipative heat conductor at all finite temperatures, at least above , for all .
Iv Conclusion
In summary, we have investigated the longitudinal dynamical thermal conductivity of the Kiteav model on the honeycomb lattice with anisotropic exchange couplings. Using an average gauge configuration method, combined with exact diagonalization of the matter fermion sector, systems of up 7200 spinful sites have been analyzed over a wide range of temperatures and anisotropies. We have confirmed both, the applicability and range of validity of our method by complementary exact analysis of thermodynamic properties. Our main finding is that for all temperatures and anisotropies investigated, the systems feature normal dissipative transport, consistent with scattering of the matter fermions from the static fluxes, and also consistent with a temperature dependence set by the matter fermion mass gap. From an experimental point of view, transport in putative Kitaev quantum spin liquids could be interesting under strain in order to tune exchange couplings, while monitoring the temperature dependence of thermal transport.
V Acknowledgments
We thank M. Vojta and C. Hess for comments and discussion. Work of A.P. and W.B. has been supported in part by the DFG through SFB 1143, project A02. Work of W.B. has been supported in part by QUANOMET and the NTHSchool CiN. W.B. also acknowledges kind hospitality of the PSM, Dresden. This research was supported in part by the National Science Foundation under Grant No. NSF PHY1748958.
Appendix A Kinetic model of thermal conductivity
This appendix is a digression into phenomenology, which we add out of curiosity. It is customary to perform kinetic modeling of thermal conductivities, expressing in terms of quasiparticle properties Hess2018a (); Hess2001 ()
(A.1) 
where , and refer to the specific heat, velocity, and mean free path at momentum of quasiparticles of type ’’. In the present case these types are fluxes, , and matter fermions, . Disentangling the contributions to from various types ’’ is only an option, if their spectral supports are well separated. As one can read off directly from Fig. (2), exactly this is possible for strong anisotropy, e.g. , where the flux peak in at is well separated from the dominant contributions to the specific heat of the fermions. Moreover, realizing that the flux velocity , and adopting the usual approximation to drop the momentum summation in (A.1) in favor of momentum integrated quantities, one gets for .
From this it is tempting to analyze the ratio of , as obtained from our calculations, e.g. for . This is shown in Fig. 5. It is this variation of with , which this appendix is meant to speculate on as follows. First, at all , the flux lattice is completely random, setting a mean free path for the fermions, independent of . Second, for sufficiently large , fermions far up in the Dirac cone set the Fermi velocity, which is also independent of . Third, upon lowering , fermions close to the gap set . This velocity approaches zero as . Remarkably, these points are roughly captured in the temperature variation of , seen in Fig. 5.
References
 L. Balents, Nature 464, 199208 (2010).
 L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2017).
 A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
 X.Y. Feng, G.M. Zhang, and T. Xiang, Phys. Rev. Lett. 98, 087204 (2007).
 H.D. Chen and Z. Nussinov, J. Phys. A: Math. Theor. 41, 075001 (2008).
 Z. Nussinov and G. Ortiz, Phys. Rev. B 79, 214440 (2009).
 S. Mandal, R. Shankar and G. Baskaran, J. Phys. A: Math. Theor. 45, 335304 (2012).
 G. Khaliullin, Prog. Theor. Phys. Suppl. 160, 155 (2005).
 G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
 J. Chaloupka, G. Jackeli, and G. Khaliullin Phys. Rev. Lett. 105 027204 (2010).
 Z. Nussinov and J. van den Brink, Rev. Mod. Phys. 87, 1 (2015).
 S. Trebst, arXiv:1701.07056.
 A. Banerjee, C. A. Bridges, J. Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Nat. Mater. 15, 733 (2016).
 A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 6342 (2017).
 A. Banerjee, P. LampenKelley, J. Knolle, C. Balz, A. A. Aczel, B. Winn, Y. Liu, D. Pajerowski, J. Yan, C. A. Bridges, A. T. Savici, B. C. Chakoumakos, M. D. Lumsden, D. A. Tennant, R. Moessner, D. G. Mandrus, and S. E. Nagler, Nat. Part. J. Quantum Mater. 3, 8 (2018).
 J. Knolle, G.W. Chern, D.L. Kovrizhin, R. Moessner, and N.B. Perkins, Phys. Rev. Lett. 113, 187201 (2014).
 S.H. Baek, S.H. Do, K. Y. Choi, Y.S. Kwon, A.U.B. Wolter, S. Nishimoto, J. van den Brink, and B. Büchner, Phys. Rev. Lett. 119, 037201 (2017).
 J. Zheng, K. Ran, T. Li, J. Wang, P. Wang, B. Liu, Z.X. Liu, B. Normand, J. Wen, and W. Yu, Phys. Rev. Lett. 119, 227208 (2017).
 C. Hess, arXiv:1805.01746.
 K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.Y. Kee, and Y.J. Kim, Phys. Rev. B 90, 041112 (2014).
 D. Hirobe, M. Sato, Y. Shiomi, H. Tanaka, E. Saitoh, Phys. Rev. B 95, 241112 (2017).
 I.A. Leahy, C.A. Pocs, P.E. Siegfried, D. Graf, S.H. Do, K.Y. Choi, B. Normand, and M. Lee, Phys. Rev. Lett. 118, 187203 (2017).
 R. Hentrich, A. U. B. Wolter, X. Zotos, W. Brenig, D. Nowak, A. Isaeva, T. Doert, A. Banerjee, P. LampenKelley, D. G. Mandrus, S. E. Nagler, J. Sears, Y.J. Kim, B. Büchner, C. Hess, Phys. Rev. Lett. 120, 117204 (2018).
 Y. J. Yu, Y. Xu, K. J. Ran, J. M. Ni, Y. Y. Huang, J. H. Wang, J. S. Wen, and A. Y. Li, Phys. Rev. Lett. 120, 067202 (2018).
 Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Nature 559, 227231 (2018).
 J. Nasu, J. Yoshitake, and Y. Motome, Phys. Rev. Lett. 119, 127204 (2017).
 R. Steinigeweg and W. Brenig, Phys. Rev. B 93, 214425 (2016).
 A. Metavitsiadis and W. Brenig, Phys. Rev. B 96, 041115(R) (2017).
 A. Metavitsiadis, A. Pidatella, and W. Brenig, Phys. Rev. B 96, 205121 (2017).
 A. Metavitsiadis, C. Psaroudaki, and W. Brenig, arXiv:1806.02344.
 G. L. Stamokostas, P. E. Lapas, and G. A. Fiete, Phys. Rev. B 95, 064410 (2017).
 F. HeidrichMeisner, A. Honecker, and W. Brenig, Eur. Phys. J. Special Topics 151, 135 (2007).
 J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. B 92, 115122 (2015).

The specific heat is evaluated numerically via temperature derivative of entropy
(A.2)  J. Nasu, and Y. Motome, Phys. Conf. Ser. 683, 012037 (2016).
 C. Hess, C. Baumann, U. Ammerahl, B. Büchner, F. HeidrichMeisner, W. Brenig, and A. Revcolevschi, Phys. Rev. B 64, 184305 (2001).