Finite-momentum energy dynamics in a Kitaev magnet
FFinite-momentum energy dynamics in a Kitaev magnet
Alexandros Metavitsiadis ∗ and Wolfram Brenig † Institute for Theoretical Physics, Technical University Braunschweig, D-38106 Braunschweig, Germany
We study the energy-density dynamics at finite momentum of the two-dimensional Kitaev spin-model on the honeycomb lattice. Due to fractionalization of magnetic moments, the energy re-laxation occurs through mobile Majorana matter, coupled to a static Z gauge field. At finitetemperatures, the Z flux excitations act as an emergent disorder, which strongly affects the energydynamics. We show that sufficiently far above the flux proliferation temperature, but not yet inthe classical regime, gauge disorder modifies the coherent low-temperature energy-density dynamicsinto a form which is almost diffusive, with hydrodynamic momentum scaling of a diffusion-kernel,which however remains retarded, primarily due to the presence of two distinct relaxation channelsof particle-hole and particle-particle nature. Relations to thermal conductivity are clarified. Ouranalysis is based on complementary calculations in the low-temperature homogeneous gauge and amean-field treatment of thermal gauge fluctuations, valid at intermediate and high temperatures. I. INTRODUCTION
Ever since the discovery of large magnetic heattransport in quasi one-dimensional (1D) local-momentsystems , the dynamics of energy in quantummagnets , has been a topic of great interest . Un-fortunately, rigorous theoretical progress has essentiallyremained confined to 1D . Above 1D, understandingenergy dynamics in quantum magnets remains an openissue at large. If magnetic long range order (LRO) ispresent and magnons form a reliable quasi particle ba-sis, various insights have been gained for antiferromag-nets and cuprates . If LRO is absent, and in particu-lar in quantum spin liquids (QSL) , energy transporthas recently come into focus as a probe of potentiallyexotic elementary excitations. In fact, experiments inseveral quantum disordered, frustrated spin systems in D ≥ suggest unconventional magnetic energy dynam-ics. For bulk transport, this pertains, e.g., to quasi 2Dtriangular organic salts , or to 3D quantum spin icematerials . For boundary transport, i.e., the magneticthermal Hall effect, recent examples include Kagomémagnets and spin ice . A microscopic descriptionof such observations is mostly lacking.In this context, Kitaev’s compass exchange Hamilto-nian on the honeycomb lattice is of particular interest, asit is one of the few models, in which a Z QSL can exactlybe shown to exist . The spin degrees of freedom of thismodel fractionalize in terms of mobile Majorana fermionscoupled to a static Z gauge field . From a material’sperspective, α -RuCl may be a promising candidate forthis model . Free mobile Majorana fermions have beeninvoked to interpret ubiquitous unconventional continuain spectroscopies on various Kitaev materials, like inelas-tic neutron and Raman scattering , as well aslocal resonance probes .Clearly, Majorana fermions should also play a rolein the energy dynamics in putative Kitaev materials.However, bulk thermal conductivity, i.e., κ xx in α -RuCl , seems to be governed primarily by phonons,with phonon-Majorana scattering as a potential dissipa- tion mechanism . Yet, Majorana fermions may havebeen observed in the transverse energy conductivity κ xy in magnetic fields, i.e., in the thermal Hall effect, and itsalleged quantization . In view of this, energy trans-port in Kitaev spin liquids has been the subject of severaltheoretical studies . Previous investigations however,have focused on energy- current correlation functions andat zero momentum only.In this work, we take a different perspective and con-sider the energy- density correlation function directly andat finite momentum. In particular we will be interestedin the role, which the emergent disorder introduced bythermally excited Z gauge field plays in the long wave-length regime. Therefore, we map out the energy diffu-sion kernel of the Kitaev model and its momentum andenergy dependence, ranging from low, up to intermedi-ate temperatures and we contrast this with expectationfor simple diffusion in random systems. The paper is or-ganized as follows. In Section II, we briefly recapitulatethe Kitaev model. In Section III, details of our calcula-tions are provided, for the homogeneous and the randomgauge state in Subsections III A and III B, respectively,and the extraction of the diffusion kernel is described inSubsection III C. We discuss our results in Section IVand summarize in Section V. An Appendix, App. A, ongeneralized Einstein relations is included. II. MODEL
We consider the Kitaev spin-model on the two dimen-sional honeycomb lattice H = (cid:88) l ,α J α S α l S α l + r α , (1)where l = n R + n R runs over the sites of the trian-gular lattice with R = (1 , , [( , √ )] for lattice con-stant a ≡ , and r α = x,y,z = ( , √ ) , ( − , √ ) , (0 , − √ ) refer to the basis sites α = x, y, z , tricoordinated to eachlattice site of the honeycomb lattice. As extensive lit-erature, rooted in Ref. , has clarified, Eq. (1) can be a r X i v : . [ c ond - m a t . s t r- e l ] F e b mapped onto a bilinear form of Majorana fermions inthe presence of a static Z gauge η l = ± , residing on,e.g., the α = z bonds H = − i (cid:88) l ,α J α η l ,α a l c l + r α . (2)Here we introduce η l ,α to unify the notation, with η l ,x ( y ) = 1 and η l ,z = η l . There are two types of Majo-rana particles, corresponding to the two basis sites. Wechose to normalize them as { a l , a l (cid:48) } = δ l , l (cid:48) , { c m , c m (cid:48) } = δ m , m (cid:48) , and { a l , c m } = 0 . For each gauge sector { η l } ,Eq. (2) represents a spin liquid.For the purpose of this work a local energy density h c on some repeating “unit” cluster c has to be cho-sen. Obviously H = (cid:80) c h c . As for any local den-sity, the latter does not fix a unique expression for h c .Different shapes of the real-space clusters c support-ing h c will typically lead to differing high-frequency andshort wave-length spectra for its autocorrelation func-tion. However the low-frequency, long-wave-length dy-namics is governed by energy conservation and will notdepend on a particular choice of h c qualitatively . Forthe remainder of this work we therefore set c = l , fo-cusing on h l = − ( i/ (cid:80) α J α η l ,α a l c l + r α , i.e. the energydensity formed by the tricoordinated bonds around eachsite on the triangular lattice. Its Fourier transform is h q = (cid:80) l e i q · l h l with h † q = h − q and h = H III. ENERGY SUSCEPTIBILTY
In this section, we present our evaluation of the dy-namical energy susceptibility. We focus on two temper-ature regimes, namely T (cid:46) ( (cid:38) ) T (cid:63) . Here T (cid:63) is the socalled flux proliferation temperature. In the vicinity ofthis temperature the gauge field and therefore fluxes getthermally excited. Previous analysis has shown,that the temperature range over which a complete prolif-eration of fluxes occurs is confined to a rather narrow re-gion, less than a decade centered around T (cid:63) ≈ . J forisotropic exchange, J z,y,x = J , used in this work, and de-crease rapidly with anisotropy . Our strategy there-fore is to consider a homogeneous ground state gauge,i.e., η l = 1 for T (cid:46) T (cid:63) and a completely random-gaugestates for T (cid:38) T (cid:63) . This approach has proven to workvery well on a quantitative level in several studies of thethermal conductivity of Kitaev models . Note thatfor α -RuCl , the case of T (cid:46) T (cid:63) is more of conceptual,than of practical interest, since it refers to very low tem-peratures, assuming a generally accepted | J | ∼ K . A. Homogeneous gauge for T (cid:46) T (cid:63) For η l = 1 , the Hamiltonian (1) can be diagonalized analytically in terms of complex Dirac fermions. Map-ping from the real Majorana fermions to the latter can be achieved in various ways, all of which require some typeof linear combination of real fermions in order to formcomplex ones. Here we do the latter by using Fouriertransformed Majorana particles, a k = (cid:80) l e − i k · l a l / √ N with momentum k and analogously for c k . The momen-tum space quantization is chosen explicitly to comprise ± k for each | k | . Other approaches, involving reshapedlattice structures , may pose issues regarding the dis-crete rotational symmetry of the susceptibility.The fermions introduced in momentum space are com-plex with a † k = a − k , i.e., with only half of the momentumstates being independent. This encodes, that for eachDirac fermion, there are two Majorana particles. Stan-dard anticommutation relations apply, { a k , a † k (cid:48) } = δ k , k (cid:48) , { c k , c † k (cid:48) } = δ k , k (cid:48) , and { a ( † ) k c ( † ) k (cid:48) } = 0 . From this, the di-agonal form of H reads H = ∼ (cid:88) k ,γ =1 , sg γ (cid:15) k d † k ,γ d k ,γ , (3)where the ˜ (cid:80) sums over a reduced “positive” half ofmomentum space and sg γ =1(-1) for γ =1(2). Thequasiparticle energy is (cid:15) k = J [3 + 2 cos( k x ) +4 cos( k x /
2) cos( √ k y / / / . In terms of reciprocal lat-tice coordinates x, y ∈ [0 , π ] , this reads (cid:15) k = J [3 +2 cos( x ) + 2 cos( y ) + 2 cos( x − y )] / / with k = x G + y G , where G = (1 , − √ ) , [(0 , √ )] . The quasiparti-cles are given by (cid:20) c k a k (cid:21) = (cid:20) u ( k ) u ( k ) u ( k ) u ( k ) (cid:21) (cid:20) d k d k (cid:21) , (4) u ( k ) = − u ( k ) = i (cid:80) α e − i k · r α / (cid:15) k ,u ( k ) = u ( k ) = 1 √ . From the sign change of the quasiparticle energy betweenbands γ =1,2 in Eq. (3) it is clear that the relations a † k = a − k and c † k = c − k for reversing momenta of the originalMajorana fermions has to change into d † k = d − k ,switching also the bands. Indeed this is also born out ofthe transformation (4). Inserting the latter into h q , theenergy density in the quasiparticle basis reads h q = 12 ˜ (cid:88) k (cid:110)(cid:104) d † k + q d † k + q (cid:105) × (cid:20) (cid:15) k + q + (cid:15) k (cid:15) k + q − (cid:15) k (cid:15) k − (cid:15) k + q − (cid:15) k + q − (cid:15) k (cid:21) (cid:20) d k d k (cid:21)(cid:27) . (5)As to be expected, h q = = H from (3) and the off-diagonal interband transitions vanish in that limit.The energy density susceptibility χ ( q , z ) is ob-tained from Fourier transformation of the imag-inary time density Green’s function χ ( q , z ) = (cid:82) β dτ (cid:104) T τ ( h q ( τ ) h − q ) (cid:105) e iω n τ /T by analytic continuation ofthe Bose Matsubara frequency iω n = i πnT → z ∈ C and eventually z → ω + i + . In order to ease geomet-rical complexity, we refrain from confining the complexfermions to only a reduced “positive” region of momen-tum space. This comes at the expense of additionalanomalous anticommutators like, e.g., { d k , d k (cid:48) } = δ − k , k (cid:48) and their corresponding contractions. Simple al-gebra yields χ ( q , z ) = χ ph ( q , z ) + χ pp ( q , z ) (6) χ ph ( q , z ) = 1 N (cid:88) k ( (cid:15) k + q + (cid:15) k ) f k + q ( T ) − f k ( T ) z − (cid:15) k + q + (cid:15) k χ pp ( q , z ) = 12 N (cid:88) k ( (cid:15) k + q − (cid:15) k ) { [ f k + q ( T )+ f k ( T ) − × (cid:18) z − (cid:15) k + q − (cid:15) k − z + (cid:15) k + q + (cid:15) k (cid:19)(cid:27) , where the superscripts ph(pp) indicate particle-hole(particle-particle) or, synonymous intra(inter)band typeof intermediate states of the fermions, f k ( T ) =1 / ( e (cid:15) k /T + 1) is the Fermi function. This concludes theformal details for T (cid:46) T (cid:63) . B. Random gauge for T (cid:38) T (cid:63) In a random gauge configuration, translational invari-ance of the Majorana system is lost, and we resort toa numerical approach in real space. This has been de-tailed extensively for 1D and 2D models andis only briefly reiterated here for completeness sake. Firsta spinor A † σ = ( a . . . a l . . . a N , c . . . c l + r x . . . c N ) , com-prising the Majoranas on the N sites of the lattice isdefined. Using this, the energy density h q and the Hamil-tonian (2), i.e., h , are rewritten as h q = A † g q A / . Boldfaced symbols refer to vectors and matrices, i.e., g q is a N × N array. Next a spinor D † σ = ( d † . . . d † N , d . . . d N ) of N complex fermions is defined by D = FA using theunitary (Fourier) transform F . The latter is built fromtwo disjoint N × N blocks I i =1 , σρ = e − i k σ · R iρ / √ N , with R iρ = l and l + r x , for a - and c -Majorana lattice sites,respectively. k is chosen such, that for each k , there ex-ists one − k , with k (cid:54) = − k . Finally, for convenience, F is rearranged such as to associate the d † . . . d † N with the N/
2) = N “positive” k -vectors. With this h q = D † ˜g q D / , (7)where ˜ o = FoF † . We emphasize, that (i) F does not diagonalize h q and (ii) that in general, the N × N matrices of Fourier transformed operators ˜ o will containparticle number non-conserving entries of D fermions.As for the case of the homogeneous gauge in Sec. III A,the energy density susceptibility χ ( q , z ) for a particular gauge sector { η l } is obtained by analytic continuation from the imaginary time density Green’s function χ ( q , τ ) = (cid:104) T τ ( h q ( τ ) h − q ) (cid:105) { η l } = 14 (cid:104) T τ [( D † ˜g q D )( τ )( D † ˜g q D ) † ] (cid:105) { η l } . (8)This is evaluated using Wick’s theorem for quasiparticles T = UD , referring to a N × N Bogoliubov transforma-tion U , determined numerically for a given distribution { η l } , such as to diagonalize ˜g , i.e., ( U˜g U † ) ρσ = δ ρσ (cid:15) ρ ,with (cid:15) ρ = ( (cid:15) . . . (cid:15) N , − (cid:15) · · · − (cid:15) N ) . We get χ ( q , z ) = (cid:88) ρσ Π σρ ( z ) w σρ, q [ w (cid:63) ¯ ρ ¯ σ, q − w (cid:63)σρ, q ] , Π σρ ( z ) = f σ ( T ) − f ρ ( T ) z − (cid:15) σ + (cid:15) ρ , (9) w ρσ, q = ( 12 U ˜ g q U † ) ρσ , where f σ ( T ) = 1 / ( e (cid:15) σ /T + 1) , and overbars refer to swap-ping the upper and lower half of the range of 2 N indices,e.g., ¯ ρ = ρ ∓ N for ρ ≷ N .As a final step, χ ( q , z ) from Eq. (9) is averaged overa sufficiently large number of random distributions { η l } .This concludes the formal details of the evaluation of theenergy density susceptibility for T (cid:38) T (cid:63) . C. Diffusion Kernel
We will relate χ ( q , z ) to a diffusion kernel D ( q , z ) bythe following phenomenological Ansatz , rooted in hydro-dynamic theory and memory function approaches χ ( q , z ) = χ q iq D ( q , z ) z + iq D ( q , z ) . (10)This should be viewed as a definition of D ( q , z ) and astatic energy-density susceptibility χ q . Neither does thistake into account fine details concerning differences be-tween static, adiabatic, or isolated susceptibilities nordoes it formulate the momentum scaling in terms of har-monics of the honeycomb lattice instead of the simplerfactor of q . The latter implies, that the momentum de-pendence of D ( q , z ) is adapted best to the small q regime.Since by construction of (10), D ( q , z ) has a properspectral representation, χ q results from the sum rule χ q = (cid:90) ∞−∞ dωπω Im[ χ ( q , ω + )] , (11)with χ ( q , ω + ) = χ ( q , ω + i + ) . Introducing a normalizedsusceptibility ¯ χ ( q , ω + ) = χ ( q , ω + ) /χ q , we will extractthe diffusion kernel from D ( q , ω + ) = 1 iq ω ¯ χ ( q , ω + )1 − ¯ χ ( q , ω + ) . (12) ω/ J z . . . I m χ ( q , ω + ) h m g / J z J x , y = J z hmg , L = 300 q i = [0, 20] T = 0.01 J z i + ≡ i . . i + ≡ i (cid:15) q phpp Figure 1. Spectrum, Im χ ( q , ω + ) , of energy density suscepti-bility versus ω at fixed small q , for low T = 0 . J z (cid:46) T (cid:63) , us-ing Eq. (6) for homogeneous ground state gauge. System size L = 300 × . Momentum q = 2 π/L (cid:80) j =1 , q ij G j . Inset:Blow up of low- ω region with reduced imaginary broadening i + . Dashed black: upper ph-continuum bound at T = 0 . IV. RESULTS
We now discuss the density dynamics as obtained fromthe previous sections. First we consider the low- T be-havior, T (cid:46) T (cid:63) , using the homogeneous gauge groundstate. As from Eq. (6), χ ( q , z ) sums two channels:(i) particle-hole (ph) and (ii) particle-particle (pp) ex-citations. Their spectral support is < | ω | <(cid:15) ˜ q for (ph)and (cid:15) ˜ q < | ω | (cid:46) max(2 (cid:15) k )=3 J at | q |(cid:28) for (pp), where ˜ q = q + k D refers to the wave vector with respect to the lo-cation of the Dirac cone. A typical spectrum Im χ ( q , ω + ) is shown in Fig. 1 at small, albeit finite q . The insetdepicts a blowup of the low- ω region dissecting the spec-trum into its ph and pp contributions. While the imagi-nary broadening in the inset is already reduced such thatfinite size oscillations start to show, the pp-channel stillexhibits some weight below its cut-off at (cid:15) ˜ q (cid:39) . J z .This will vanish as Im z → . For the ph-channel how-ever, the spectral weight in this energy range is not a fi-nite size effect. Due to the Dirac cone, the Fermi volumeshrinks to zero in the Kitaev model as T → , i.e., occu-pied states only stem from a small patch with (cid:15) k (cid:46) T around the Dirac cone. Therefore, the weight of theph-channel decreases rapidly to zero as T → . In thisregime and for small q , because of the linear fermiondispersion close to the cones, only a narrow strip of or-der ω ∈ [max(0 , (cid:15) ˜ q − T ) , (cid:15) ˜ q ] from the spectral supportdominates the ph-continuum. At the upper edge of itssupport the ph DOS is singular. The inset of Fig. 1is consistent with this, considering the finite system sizeand imaginary broadening used.Regarding the pp-channel, the complete two-particlecontinuum is unoccupied and available for excited statesas T → . This leads to the broad spectral hump seenin Fig. 1, which extends out to max(2 (cid:15) k ) = 3 J z , at J x,y = J z and is two orders of magnitude larger than theph-process at this temperature. . . . I m χ ( q , ω + ) h m g / J z T = 0.1 J z ( a ) hmg , i + ≡ i hmg , i + ≡ i . . . . T = 0.5 J z J x , y = J z , L = 60 q i = [0,4] ( b ) hmg , i + ≡ i hmg , i + ≡ i ω/ J z . . . I m χ ( q , ω + ) r nd / J z T = 0.1 J z ( c ) rnd , comb , i + ≡ i rnd , ph , 200 realzrnd , ppmax ( (cid:15) ph ) ω/ J z . . . . T = 0.5 J z ( d ) rnd , comb , i + ≡ i rnd , ph realzrnd , ppmax ( (cid:15) ph ) Figure 2. Spectrum, Im χ ( q , ω + ) , of energy density suscep-tibility versus ω at fixed small q = 2 π/L (cid:80) j =1 , q ij G j , con-trasting homogeneous (a,b) versus random (c,d) gauge on allidentical system sizes L = 60 × for two identical tempera-tures T = 0 . J z (a,c) and T = 0 . J z (b,d). Solid black, green,and red lines in panels (c,d): total, ph-, and pp-spectrum. Ex-tended vertical gray ticks in panels (c,d): upper ph-continuumbound at q in homogeneous gauge. Imaginary smoothing i + for homogeneous case, both, identical (solid black) and signif-icantly larger than random gauge case (solid magenta). Notethe different y-axis scales. A fingerprint of potentially diffusive relaxation of den-sity modes at finite momentum q is the near-linear behav-ior of Im χ ( q , ω + ) ∼ χ q ω/ ( Dq ) at small ω . Definitely,this should neither be expected, nor is it observed in Fig.1, since for T (cid:46) T (cid:63) , the density dynamics is set by co-herent two-particle excitations of the Dirac fermions inthe homogeneous gauge state.For the remainder of this work, we therefore now turnto temperatures above the flux proliferation, i.e., T (cid:38) T (cid:63) ,using random gauge states. To begin, and in Fig. 2,we first describe the impact of the random gauges, bycontrasting the dynamic density susceptibilities againsteach other with, and without random gauges, for other-wise identical system parameters, and for two differenttemperatures, T = 0 . and 0.5, in Figs. 2(a,c) and (b,d),respectively. Note, that while the linear system size issmaller by a factor of five with respect to Fig. 1, thewave vector has also been rescaled accordingly. Thereforethese two figures can be compared directly. Obviously,in the homogeneous gauge, significant degeneracies, evenon 60 ×
60 lattices, lead to visible discretization effects.Therefore, in Figs. 2(a,b) we include spectra with animaginary broadening, increased relative to Figs. 2(c,d),for a better comparison with the latter.Several features can be observed. First, while the ph-channel in the uniform gauge clearly displays the singularbehavior at ω = (cid:15) ˜ q , mentioned earlier and visible becauseof the elevated temperatures, in the random gauge case,it displays a smooth peak. Second, as can be read of fromthe y -axis, the weight of the ph-channel strongly increaseswith increasing T . For the temperatures depicted, the . . . I m χ ( q , ω + ) / χ q ( a ) T = 0.1 J z ω/ J z . . . ( b ) T = 0.5 J z J x , y = J z L = 60200 realz i + ≡ i . . . ( c ) T = 0.9 J z q = [0, 2] q = [0, 4] q = [0, 6] q = [2, 4] q = [2, 6] q q q q q Figure 4. Spectrum of the normalized density susceptibility Im χ ( q , ω + ) /χ q in the random gauge state, with 200 realizations,versus ω for various temperatures T /J z = 0 . , q =2 π/L (cid:80) j =1 , q ij G j within the small- q region, inset panel (b), of the irreducible wedge of the BZ for L = 60 × sites and forimaginary broadening i + . pp-channel is much less T -dependent. Third, the ph-,versus the pp-contributions to Im χ ( q , ω + ) cannot only bedissected in the uniform gauge by virtue of Eqs. (6), butalso for the random gauge, Eq. (9) can be decomposedinto addends with (cid:15) σ (cid:15) ρ ≷ . This evidences, that in thelatter case, the ph-spectrum changes completely. As isobvious from Figs. 2(c,d), the ph-channel spreads intoa broad feature, extending over roughly the entire one-particle energy range. The pp-channel on the other handseems less affected by the gauge disorder, with a shapequalitatively similar to that in the gauge ground state,as can be read off by comparing Figs. 2(a,c).Most remarkably, for intermediate temperatures, as inFig. 2(d) at T /J z = 0 . , the overall shape of the spec-trum is very reminiscent of a diffusion-pole behavior atfixed momentum, i.e., Im χ ( q , ω ) ∝ ω Γ / ( ω + Γ ) , withsome relaxation rate Γ . To clarify this in more detail,we therefore proceed and analyze χ ( q , z ) in terms of thehydrodynamic expression Eq. (10).To this end, we first extract the static susceptibility χ q ,performing the sum rule of Eq. (11) via numerical inte- Figure 3. Static susceptibility χ q versus momentum in therandom gauge state, with 200 realizations. Small solid dots:all of irreducible wedge of the BZ at L = 30 × for threetemperatures T /J z = 0 . , 0.5, and 0.9 (black, orange, andblue). Big solid red dots, finite size effects: χ q at low − q for L = 60 × and simultaneously q ∈ × BZ. gration, using Im χ ( q , z ) from the corresponding randomgauge states. The results are shown in Fig. 3, spanningan irreducible q -wedge of the BZ. Since energy conser-vation renders the dynamic density response singular at q = , this momentum will be excluded hereafter. Obvi-ously χ q is a smooth and featureless function. The figurealso contrasts 30 ×
30 against 60 ×
60 systems at selectedmomenta. The finite size effects are small.Next, and in Fig. 4, we consider the global variationwith momentum, of the normalized spectrum of the dy-namical energy density susceptibility versus ω . Since ourfocus is on the hydrodynamic regime, we remain withsmall momenta. These momenta are indicated on a frac-tion of an irreducible wedge of the BZ in Fig. 4(b). Aspacing of Mod(2) of the momenta has been chosen toallow for later analysis of finite size effects in compar-ison to systems with a linear dimension smaller by afactor of 2. The spectra show significant changes withmomentum. First, the low- ω contributions, which stemprimarily from the ph-channel show a broadening of theirsupport with increasing | q | . Second, the spectral rangeof pp-excitations displays a global increase of weight with | q | . Finally, at intermediate and elevated T in Fig. 4(b)and (c), the former effect dominates the latter regardingthe global shape of the spectrum.Now we extract the diffusion kernel D ( q , ω + ) as inEq. (12). Its real part is depicted in Fig. 5 versus ω for identical momenta and temperatures as in Fig.4. Clearly, at low T , Fig. 5(a), the diffusion kerneldisplays significant q -dependence, implying that the en-ergy currents are not proportional to the gradient ∇ l h l of the energy density and therefore a simply hydrody-namic picture is not applicable. In contrast to that, atintermediate and high T , Fig. 5(b) and (c), the dif-fusion kernel is approximately momentum independent Re D ( | q | (cid:28) , ω + ) (cid:39) Re D ( ω + ) . This is consistent withFick’s law regarding q -scaling. However, the diffusionprocess remains retarded, although not very strong, dis-playing a shoulder in the pp-range and a peak at low ω in the ph-range. R e D ( q , ω + ) / ( a J z ) ( a ) T = 0.1 ω/ J z . . . ( b ) T = 0.5 J x , y = J z L = 60200 realz i + ≡ i ω min = 0.005 J z . . . ( c ) T = 0.9 q = [0, 2] q = [0, 4] q = [0, 6] q = [2, 4] q = [2, 6] q q q q q Figure 5. Real part of the diffusion kernel Re D ( q , ω + ) in the random gauge state, with 200 realizations, versus ω ≥ . J z forvarious temperatures T /J z = 0 . , q = 2 π/L (cid:80) j =1 , q ij G j within the small- q region, inset panel (b), of the irreducible wedge of the BZ for L = 60 × sites and for imaginary broadening i + . As ω/J z → , the numerical accuracy of the trans-form Eq. (12), comprising small numbers in the numer-ator and denominator, is an issue with respect to sys-tem size and imaginary broadening and we have to re-main with ω ≥ ω min = 0 . J z for the parameters used.See also the discussion of Fig. 7. In view of the steepslope in this regime, it may be that on any finite sys-tem Re D ( | q | (cid:28) , ω = 0) vanishes singularly below theunphysically small energy scale ω min , while in the ther-modynamic limit Re D ( | q | (cid:28) , ω = 0) is of the order ofthe low-energy peak height. This behavior is very rem-iniscent of similar findings for the thermal conductivity,see Ref. and App. A. For ω/J z (cid:38) . , the spectralsupport terminates and the diffusion kernel turns purelyimaginary ∝ ω − . In conclusion, at not too low tem-peratures and not too short time scales, gauge disorder T / J z . . . . . R e D ( q , ω + , T ) / ( a J z ) J x , y = J z L = 30200 realz i + ≡ i b ) q i = [1,0] 0 1 2 T / J z χ q / J z ( d ) . . . . . R e D ( q , ω + , T ) / ( a J z ) ( a ) q i = [6,0] ω = 0.1 J z ω = 0.3 ω = 0.5 ω = 1.0 0 1 2 T / J z χ q / J z ( c ) Figure 6. Real part of the diffusion kernel Re D ( q , ω + , T ) inthe random gauge state, with 200 realizations, versus T /J z at q = 2 π/L (cid:80) j =1 , q ij G j at (a) q i = [6 , and (b) [1 , , andvarious fixed energies ω , for L = 30 × sites and imagi-nary broadening i + . Insets (c,d): static susceptibility χ q ( T ) versus T /J z at identical L and q . in the Kitaev magnet leads to an energy density dynam-ics, very similar to conventional diffusion, regarding itsmomentum scaling, with, however, some retardation re-maining. This should be contrasted with the underlyingspin model being a translationally invariant system.Turning to the temperature dependence, we considertwo representative momenta q and several energies. Thecorresponding diffusion kernel Re D ( q , ω + ) and the staticsusceptibility χ q are shown versus T in Figs. 6(a), (b)and (c), (d), respectively. The temperature range hasbeen truncated deliberately to T /J z (cid:38) . , since belowsuch temperatures, and from the q -scaling in Fig. 5 thedensity dynamics is far from hydrodynamic. The figureclearly demonstrates, that for T /J z (cid:38) . , the diffusionkernel rapidly settles onto some almost constant value,set by energy and momentum. This is consistent withFig. 5, which displays only weak overall change betweenthe diffusion kernels for the two temperatures of panels(b) and (c). As a consequence, the global T -dependenceof χ ( q , z ) is essentially set by the static energy suscepti-bility. As the insets show, the latter exhibits a maximumin the vicinity of T /J z ∼ . For T /J z (cid:29) , χ q ap-proaches its classical limit, decaying ∝ T − , which canalso be read of from Eqs. (11,9). Such behavior is typicalalso for other static susceptibilities of spin systems.In closing, we provide some measure of the finitesize effects involved in our calculations. To thisend we consider both, the absolute and the relativedifference between the diffusion kernels, ∆( q , ω ) = | Re[ D × ( q , ω + ) − D × ( q , ω + )] | and Λ( q , ω ) =2∆( q , ω ) / | Re[ D × ( q , ω + ) + D × ( q , ω + )] | , respec-tively, on N = 30 × and × site systems, for iden-tical wave vectors. Regarding the latter, this implies afactor of difference between their integer representa-tion in terms of G , . The differences are shown in Fig.7. They display statistical noise from the finite num-ber of random gauge realizations and remain acceptablysmall for all ω . Only for ≈ ω (cid:28) . , where Λ( q , ω ) is of O (10%) , the error is not of finite size, or statisticalorigin, but rather stems from the systematic numerical ω/ J z . . . ∆ ( q , ω + ) / ( a J z ) J x , y = J z T = 0.5 J z
200 realz i + ≡ i ω min = 0.005 J z L = L = (a) Λ ( q , ω + ) q =[0,1], L =30 ( q i / L =60)) q =[0,2] q =[0,3] q =[1,2] q =[1,3] (b) Figure 7. (a) Absolute and (b) relative difference of real partof the diffusion kernel Re D ( q , ω + ) in the random gauge state,with 200 realizations, versus ω ≥ . J z at T /J z = 0 . ,for various small momenta q = 2 π/L (cid:80) j =1 , q ij G j on L =30 × and q i → q i on L = 60 × , within the small- q region of insets panel (a), of the irreducible wedge of the BZand for imaginary broadening i + . inaccuracies, mentioned in the preceding, of the denomi-nator in Eq. (12) with χ ( q , ω + ) obtained from Eq. (9) as ω → . In turn, D ( q , ω + ) at very low ≈ ω (cid:28) . maybe inacurate by ∼ . As to be expected, the actualfinite size errors are largest for the smallest wave vector. V. SUMMARY
In summary, above an intermediate temperature scale T ∼ . J , which is still well below the classical limit, theenergy density in Kitaev magnets at finite momentumrelaxes remarkably similar to diffusion in random me-dia, with, however, a clearly notable difference. Namely,while the momentum scaling is practically hydrodynamic ∝ q , the diffusion kernel is not completely energy in-dependent, i.e., it displays some retardation within itssupport. The origin of the latter can be traced backto the presence of two distinct relaxation channels forthe energy density, comprising particle-hole and particle-particle excitations of the Dirac fermions. Both propa-gate in a strongly disordered landscape, created by ther-mally induced gauge excitations. Their combined effect,however, does not lead to a constant diffusion rate. Atextremely low energies we observe a dip in the diffusionkernel. This is consistent with similar claims for thedynamical thermal conductivity, to which we find thatour results connect consistently via generalized Einsteinrelations. Future analysis, focusing on the real space,instead of the momentum dependence of energy-densityrelaxation could be interesting in order to predict finitetemperature quench dynamics. Acknowledgments : This work has been supported inpart by the DFG through project A02 of SFB 1143(project-id 247310070), by Nds. QUANOMET (project . C ( q , ω ) ≡ ω − e − ω / T χ q R e D ( q , ω + ) / ( a J z ) ( a ) T = 0.1 J x , y = J z L = 60200 realz i + ≡ i ω min = 0.005 J z ω/ J z . . . ( b ) T = 0.5 q = [0,2] q = [0,4] q = [0,6] q = [2,4] q = [2,6] Figure 8. Momentum evolution of Einstein relation for cur-rent correlation function in the random gauge state, with200 realizations, versus ω ≥ . J z for T /J z = 0 . , and0.5 in panels (a) and (b), for various small momenta q =2 π/L (cid:80) j =1 , q ij G j within the small- q region identical to Fig.5, for L = 60 × sites and for imaginary broadening i + . NP-2), and by the National Science Foundation underGrant No. NSF PHY-1748958. W.B. also acknowledgeskind hospitality of the PSM, Dresden.
Appendix A: Current correlation functions
In the limit of q → one may speculate, that the dy-namical energy-density diffusion kernel is related to thedynamical energy-current correlation function via a gen-eralized Einstein relation. The zero momentum currentcorrelation function has been considered in Ref. . Forcompleteness sake, we now clarify a relation of the lat-ter quantity to the present work. From Mori-Zwanzig’smemory function method we have z ( χ q − χ ( q , z )) = 1 z − M ( q , z ) χ q χ q , (A1)where M ( q , z ) = (cid:104) Lh q | ( z − QL ) − QLh q (cid:105) is the mem-ory function. L is the Liouville operator LA = [ H, A ] and (cid:104) A | B (cid:105) = (cid:82) β (cid:104) A + ( λ ) B (cid:105) dλ − β (cid:104) A + (cid:105)(cid:104) B (cid:105) is Mori’s scalarproduct with A ( λ ) = e λH Ae − λH = e λL A and β = 1 /T is the inverse temperature. χ q is the isothermal energysusceptibility χ q = (cid:104) h q | h q (cid:105) and Q is a projector perpen-dicular to the energy density, which is formulated usingMori’s product as Q = 1 − | h q (cid:105) χ − q (cid:104) h q | . We emphasize,that Eq. (A1) is not a “high-frequency”, or “slow-mode”approximation. It is a rigorous statement. Due to time-reversal invariance, QLh q = Lh q . Moreover, using thecontinuity equation in the hydrodynamic regime, i.e., dis-carding the lattice structure, we have Lh q = − q · j q ,where j q is the energy current. Altogether iχ q D ( q , z ) = (cid:88) µν e q µ e q ν (cid:104) j q µ | z − QL j q ν (cid:105) , (A2)where e q µ are the components of the unit vector into q -direction. While for arbitrary q the right hand siderefers to a so-called current relaxation-function witha dynamics governed by a projected Liouville opera-tor QL , for q → , one finds that lim q → (cid:104) j q µ | ( z − QL ) − j q ν (cid:105) = (cid:104) j µ | ( z − L ) − j ν (cid:105) , which is the genuinecurrent relaxation-function comprising the complete Li-ouvillean dynamics. This turns Eq. (A2) into an Einsteinrelation for q → . Finally, the spectrum of the currentrelaxation-function can be related to that of a standardcurrent correlation-function C µν ( t ) = (cid:104) j µ ( t ) j ν (cid:105) by theKubo relation and the fluctuation dissipation theorem ω − e ω/T lim q → χ q Re D ( q , z ) = C ( ω ) . (A3) Here we have discarded questions of anisotropy. Whilethe present work’s focus is on q (cid:54) = 0 , it is now very tempt-ing to evaluate the left hand side of Eq. (A3) using, e.g.,the two temperatures of Fig. 5 (a,b) and to considerits evolution with momentum. This is shown in Fig. 8,which should be compared with Fig. 5 (b,d) of Ref. .For this, and because of a different energy unit and nor-malization of spectral densities in the latter Ref., T has torescaled by and the y-axis by /π . While the rescaledtemperatures T = 0 . and . of Ref. are not abso-lutely identical to the ones we use, it is very satisfyingto realize, that the limit q → of Eq. (A3), which canbe anticipated from Fig. 8 is completely in line withFig. 5 of Ref. , including the dip at very low ω . This iseven more remarkable in view of the numerical represen-tation and treatment of the Majorana fermions used inthe present work and in Ref. being decisively different. ∗ [email protected] † [email protected] A. V. Sologubenko, E. Felder, K. Giannò, H. R. Ott, A.Vietkine, and A. Revcolevschi, Phys. Rev. B , R6108(2000). T. Kawamata, N. Takahashi, T. Adachi, T. Noji, K. Kudo,N. Kobayashi, Y. Koike, J. Phys. Soc. Jp. , 034607(2008). N. Hlubek, P. Ribeiro, R. Saint-Martin, A. Revcolevschi,G. Roth, G. Behr, B. Büchner, C. Hess, Phys. Rev. B ,020405 (2010). A. V. Sologubenko, K. Giannò, H. R. Ott, U. Ammerahl,and A. Revcolevschi, Phys. Rev. Lett. , 2714 (2000). C. Hess, C. Baumann, U. Ammerahl, B. Büchner, F.Heidrich-Meisner, W. Brenig, and A. Revcolevschi, Phys.Rev. B , 184305 (2001). X. Zotos, F. Naef, and P. Prelovsek, Phys. Rev. B ,11029 (1997). Heidrich-Meisner, F., A. Honecker, D. C. Cabra, and W.Brenig, Phys. Rev. B , 134436 (2003). Heidrich-Meisner, F., A. Honecker, D. C. Cabra, and W.Brenig, Phys. Rev. Lett. , 069703 (2004). C. Hess, Eur. Phys. J. Spec. Top. , 73 (2007). F. Heidrich-Meisner, A. Honecker, and W. Brenig, Eur.Phys. J. Special Topics , 135 (2007). C. Hess, Phys. Rep. , 1 (2019). B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen,R. Steinigeweg, and M. Znidaric, arXiv:2003.03334 G. S. Dixon, Phys. Rev. B , 2851 (1980). C. Hess, B. Büchner, U. Ammerahl, L. Colonescu, F.Heidrich-Meisner, W. Brenig, and A. Revcolevschi, Phys.Rev. Lett. , 197002 (2003). S.P. Bayrakci, B. Keimer, and D.A. Tennant, arXiv:1302.6476 A. L. Chernyshev and Wolfram, Phys. Rev. B , 054409(2015). L. Balents, Nature , 199-208 (2010). L. Savary and L. Balents, Rep. Prog. Phys. , 016502(2017). M. Yamashita, N. Nakata, Y. Senshu, M. Nagata, H. M.Yamamoto, R. Kato, T. Shibauchi, and Y. Matsuda, Sci-ence , 1246 (2010). P. Bourgeois-Hope, F. Laliberté, E. Lefrançois, G. Grisson-nanche, S. R. de Cotret, R. Gordon, S. Kitou, H. Sawa, H.Cui, R. Kato, L. Taillefer, and N. Doiron-Leyraud, Phys.Rev. X , 041051 (2019). J. M. Ni, B. L. Pan, B. Q. Song, Y. Y. Huang, J. Y. Zeng,Y. J. Yu, E. J. Cheng, L. S. Wang, D. Z. Dai, R. Kato,and S. Y. Li, Phys. Rev. Lett. , 247204 (2019). G. Kolland, O. Breunig, M. Valldor, M. Hiertz, J. Frielings-dorf, and T. Lorenz, Phys. Rev. B , 060402(R) (2012). W. H. Toews, S. S. Zhang, K. A. Ross, H. A. Dabkowska,B. D. Gaulin, and R. W. Hill, Phys. Rev. Lett. , 217209(2013). Y. Tokiwa, T. Yamashita, D. Terazawa, K. Kimura, Y.Kasahara, T. Onishi, Y. Kato, M. Halim, P. Gegenwart,T. Shibauchi, S. Nakatsuji, E.-G. Moon, and Y. Matsuda,J. Phys. Soc. Jpn. , 064702 (2018). M. Hirschberger, R. Chisnell, Y. S. Lee, and N. P. Ong,Phys. Rev. Lett. , 106603 (2015). D. Watanabe, K. Sugii, M. Shimozawa, Y. Suzuki, T. Ya-jima, H. Ishikawa, Z. Hiroi, T. Shibauchi, Y. Matsuda, andM. Yamashita, PNAS , 8653 (2016). M. Yamashita, M. Akazawa, M. Shimozawa, T. Shibauchi,Y. Matsuda, H. Ishikawa, T. Yajima, Z. Hiroi, M. Oda,H. Yoshida, H.-Y. Lee, J. H. Han, and N. Kawashima, J.Phys.: Condens. Matter , 074001 (2019). M. Hirschberger, J. W. Krizan, R. J. Cava, and N. P. Ong,Science , 106 (2015). A. Kitaev, Ann. Phys. (N.Y.) , 2 (2006). X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. , 087204 (2007). H.-D. Chen and Z. Nussinov, J. Phys. A: Math. Theor. ,075001 (2008). Z. Nussinov and G. Ortiz, Phys. Rev. B , 214440 (2009). S. Mandal, R. Shankar and G. Baskaran, J. Phys. A: Math.Theor. , 335304 (2012). 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 , 041112 (2014). S. Trebst,
Kitaev Materials , Lecture Notes of the 48th IFFSpring School 2017, S. Blügel, Y. Mokrousov, T. Schäpers,Y. Ando (Eds.), ISBN 978-3-95806-202-3 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. , 733 (2016). A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone,M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moess-ner, and S. E. Nagler, Science , 6342 (2017). A. Banerjee, P. Lampen-Kelley, 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. Lums-den, D. A. Tennant, R. Moessner, D. G. Mandrus, and S.E. Nagler, Nat. Part. J. Quantum Mater. , 8 (2018). J. Knolle, G.W. Chern, D.L. Kovrizhin, R. Moessner, andN.B. Perkins, Phys. Rev. Lett. , 187201 (2014). D. Wulferding, Y. Choi, S.-H. Do, C. H. Lee, P. Lemmens,C. Faugeras, Y. Gallais, and K.-Y. Choi, Nat Commun ,1 (2020). 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. , 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. , 227208 (2017). D. Hirobe, M. Sato, Y. Shiomi, H. Tanaka, and 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. , 187203 (2017). R. Hentrich, A. U. B. Wolter, X. Zotos, W. Brenig, D.Nowak, A. Isaeva, T. Doert, A. Banerjee, P. Lampen-Kelley, D. G. Mandrus, S. E. Nagler, J. Sears, Y.-J. Kim,B. Büchner, C. Hess, Phys. Rev. Lett. , 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. ,067202 (2018). Alexandros Metavitsiadis and Wolfram Brenig, Phys. Rev.B , 035103 (2020). 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 , 227 (2018). T. Yokoi, S. Ma, Y. Kasahara, S. Kasahara, T. Shibauchi,N. Kurita, H. Tanaka, J. Nasu, Y. Motome, C. Hickey,S. Trebst, and Y. Matsuda,
ArXiv:2001.01899 [Cond-Mat] (2020). A. Metavitsiadis and W. Brenig, Rev. B , 041115(R)(2017). J. Nasu, J. Yoshitake, and Y. Motome, Phys. Rev. Lett. , 127204 (2017). A. Metavitsiadis, A. Pidatella, and W. Brenig, Phys. Rev.B , 205121 (2017). A. Pidatella, A. Metavitsiadis, and W. Brenig Phys. Rev.B , 075141 (2019). A. Metavitsiadis, C. Psaroudaki, and W. Brenig, Phys.Rev. B , 205129 (2019). A. Metavitsiadis and W. Brenig, arXiv:2009.04467 [Cond-Mat] (2020). J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. B ,115122 (2015). D. Forster,