Signatures of fractionalization in spin liquids from interlayer thermal transport

Quantum spin liquids (QSLs) are intriguing phases of matter possessing fractionalized excitations. Several quasi-two dimensional materials have been proposed as candidate QSLs, but direct evidence for fractionalization in these systems is still lacking. In this paper, we show that the inter-plane thermal conductivity in layered QSLs carries a unique signature of fractionalization. We examine several types of gapless QSL phases - a $Z_2$ QSL with either a Dirac spectrum or a spinon Fermi surface, and a $U(1)$ QSL with a Fermi surface. In all cases, the in-plane and $c-$axis thermal conductivities have a different power law dependence on temperature, due to the different mechanisms of transport in the two directions: in the planes, the thermal current is carried by fractionalized excitations, whereas the inter-plane current is carried by integer (non-fractional) excitations. In layered $Z_2$ and $U(1)$ QSLs with a Fermi surface, the $c-$axis thermal conductivity is parametrically smaller than the in-plane one, but parametrically larger than the phonon contribution at low temperatures.


I. INTRODUCTION
Quantum spin liquids (QSLs) are phases of matter with intrinsic topological order, which cannot be characterized by local order parameters as typically used in symmetry-breaking phases. Instead, their primary characteristic is the emergence of excitations with fractional quantum numbers [1][2][3][4][5][6]. The presence of these excitations is related to the existence of long-range entanglement in ground states of such systems [7,8]. In addition, the excitations are accompanied by an emergent gauge field leading to a low-energy description in terms of gauge theories. The relevant gauge group can be discrete (e.g., Z 2 ) or continuous (e.g., U (1)). The matter excitation spectrum may be gapped (as in a gapped Z 2 phase [9-15]) or gapless (as in a gapless Z 2 [16,17] or U (1) [18][19][20][21][22][23][24][25][26] QSL).
The excitations of QSLs can carry fractional quantum numbers corresponding to global symmetries possessed by the system [13,47,48] and also possess fractional (anyonic) statistics [10, [48][49][50]. There have been numerous proposals to detect these fractional quantum numbers and statistics in QSL materials [48,[51][52][53][54][55][56][57]. The presence of fractionalization itself has primarily been de-duced through a diffuse scattered intensity seen in inelastic neutron scattering experiments on various candidate spin liquids at temperatures much smaller than the relevant exchange coupling [43,58]. The absence of sharp features in the neutron scattering intensity is attributed to the presence of a multi-particle continuum [59]. However, such broadening can also arise due to other factors such as disorder and it would be useful to have additional signatures of fractionalization.
In this work, we propose the inter-plane thermal conductivity κ c as a probe for fractionalization in a system of weakly coupled layers of two dimensional gapless QSLs [60]. The in-layer thermal conductivity κ ab in these materials is dominated by the low-energy fractionalized excitations pertinent to the type of QSL in question; in contrast, the thermal current between the planes must be carried by a gauge invariant excitation with integer quantum numbers. This is because the emergent gauge charge carried by fractionalized excitations is conserved separately in each layer, and therefore a single spinon cannot move from one layer to the next. Moreover, a non-gauge invariant fractionalized excitation, such as a spinon, is highly non-local in space (it is composed of a long "string" of local spin operators). This implies that the matrix element of a local operator to transfer pairs of spinons from one layer to another decays exponentially with the spatial separation between the two spinons. Therefore, only pairs of nearby spinons can hop between adjacent layers.
The situation is depicted schematically in Fig. 1, where a single spinon is deconfined and may propagate freely in each plane, while only pairs of spinons may hop between planes. Therefore, κ c in a gapless QSL is expected to be qualitatively different from the in-layer thermal conductivity, and obey a different power law at low temper- In-plane and c-axis thermal conductivity for several types of QSL. Z2 Dirac refers to a Z2 QSL with a Dirac spectrum of fermionic fractional excitations. Z2 FS is a Z2 QSL with a Fermi surface of fractional excitations. U (1) refers to a spinon Fermi surface coupled to a U (1) fluctuating gauge field. α = 6∆A/(π + ∆A), with ∆A the (dimensionless) time-reversal preserving disorder strength [see Eq. (11)]. The result for the clean Z2 FS case is correct up to logarithmic factors.
In-plane c-axis Clean Disordered Clean Disordered Z 2 Dirac T [62] T T 1/3 [63] T [63] T 5/3 T 2 atures [61]. An experimental detection of such a parametrically large anisotropy in ratio κ ab /κ c at low temperatures will be a strong indication of the existence of fractionalized excitations and hence a QSL state. Our findings are summarized in Table I. We have considered three cases: a gapless Z 2 QSL with either a Dirac spectrum or a spinon Fermi surface, and U (1) QSL with a spinon Fermi surface. In all cases, the in-plane and c-axis thermal conductivity follow qualitatively different behavior as a function of temperature, for both clean and mildly disordered systems. In all QSLs we consider, the inter-plane thermal conductivity follows a power law behavior in temperature, with an exponent which is larger than for the corresponding intra-plane behavior. Interestingly, in some cases, the exponent of the inter-plane thermal conductivity is smaller than 3, and therefore it is parametrically larger than the phonon contribution (proportional to T 3 ) at sufficiently low temperatures.

II. CLEAN Z2 QUANTUM SPIN LIQUID
We begin by considering a layered system where each layer forms a Z 2 QSL with gapless fermionic excitations. The fermions may either have a Dirac spectrum, or form a Fermi surface. As a concrete example of the gapless Z 2 QSL, one may consider the gapless phase of the Kitaev honeycomb model [16], which consists of spin-1/2s interacting in an anisotropic manner on a two-dimensional hexagon lattice. We will use this model to facilitate our discussion; our conclusions are generic to any gapless Z 2 QSL.
The low energy theory of the Kitaev QSL phase may be described either as two linearly dispersing Majorana fermions, or equivalently as a single complex Dirac theory. Here, we consider a three-dimensional layered generalization of the Kitaev model. The low-energy effective Hamiltonian of each layer is given by FIG. 1: A schematic representation of the difference between in-plane and interplane transport. In each QSL plane, the spinons are deconfined and may travel freely. However, transport between the planes is only possible via gauge-invariant excitations, such as spinon pairs.
where l is the layer index, and ψ † l (k) = ψ A † l (k), ψ B † l (k) is a spinor of complex fermionic spinon creation operators in layer l, with A, B denoting the sublattice. σ is a vector of Pauli matrices, k = (k x , k y ) is measured relative to the corner of the honeycomb lattice (the K point), and v the Fermi velocity. m and ∆ describe a mass gap and an effective chemical potential, respectively. Throughout the paper we have set = 1. In Appendix A we show an explicit microscopic spin Hamiltonian that leads to low-energy effective Hamiltonian in Eq. (1). The m and ∆ terms arise from time reversal-breaking three spin interactions. On the honeycomb lattice, in the presence of time reversal (TR) symmetry, m = ∆ = 0, and the Fermi energy is at the Dirac point. Breaking time reversal symmetry [17], or considering generalizations of the Kitaev model to other lattices [64][65][66][67][68][69], allows for a stable Fermi surface. In all Z 2 QSLs we consider, the fermionic excitations ("spinons") are gapless, which corresponds to |∆| ≥ |m|. The fluxes of the Z 2 gauge field ("visons") are gapped.
Although the in-plane theory is described by fractional excitations, inter-plane transport must be mediated by gauge-invariant excitations. The most relevant interlayer coupling terms which are allowed by symmetry are given by where l, l are neighboring layers, σ α are the Pauli matrices (with σ 0 the identity matrix), J ⊥ is the strength of the inter-plane coupling, and F 0,...,4 are dimensionless coupling constants. In Sec. A 4 we argue that generically, J ⊥ is proportional to the microscopic spin-spin inter-layer interactions. For simplicity, we will mostly focus on the case where only F 0 = F is nonzero. A derivation of such a coupling term from a microscopic spin-spin interaction is given in the Appendix A. We believe that the particular form of the interlayer coupling is not important; the contribution to the thermal conductivity from other terms give the same parametric dependence on temperature. The crucial point is that the inter-plane coupling term must contain an even number of fermion operators from each layer, as a single fractional excitation may not hop from one layer to another.
In a generic Z 2 QSL, there are also short-range intraplane interactions between the fermionic spinons. However, for most of the following discussion we may ignore such interactions, as they are irrelevant in the Dirac case, and lead to a Landau Fermi liquid state with well defined quasiparticles in the Fermi surface case.
For a clean Z 2 QSL, whose low energy theory is described by weakly interacting fermions, the interlayer thermal conductivity may be calculated to lowest order in J ⊥ using Fermi's golden rule. We work in the basis of the eigenvalues of the in-plane Hamiltonian; we therefore revert from the sublattice (α = A, B) to the band (λ = ±1) basis, and consider the transformed function F λ1...λ4 k,k ,q in this basis. In the case of a Dirac spectrum, the eigenstates of the in-plane Hamiltonian are given by a λ=± . Energy is transported between layers by the excitation of spinon pairs; thus, if a temperature difference δT is applied between two adjacent layers l and l , the rate with which energy transfer occurs, for the specific momenta k, k + q, k , k − q, is where |i l ,|f l are the initial and final many-body states of layer l (which are eigenstates of the J ⊥ = 0 Hamiltonian), with energies E i,l and E f,l , respectively, and sim-ilarly for layer l . Z is the partition function. The thermal conductivity is then given by (here J Q is the thermal current) where n F ( ) is the Fermi function.
For the case of a Z 2 QSL with a Dirac spectrum, the dependence of the integral on temperature can be evaluated easily by rescaling ⇒ /T and {k, k , q} ⇒ {k/vT, k /vT, q/vT }. This gives the result The case of a Z 2 QSL with a Fermi surface corresponds to ∆ = 0 in Eq. (1). To simplify the calculation, we set the mass term in (1) such that ∆ > m but |∆ − m| m. In this limit, the eigenstates of the band which crosses the Fermi energy simplify to a(k) = ψ A (k), with a nonrelativistic dispersion k = k 2 /2m * − µ, with µ = ∆ − m, and m * = m/v 2 . The result should not depend on this choice.
The evaluation of the integrals in Eq. (4) for the case of a Fermi surface is described in Appendix B. After integrating over k, k , and 1,2,3 , κ c has the form: where ν = m * /2π is the density of states on the Fermi energy, and k F = √ 2m * µ. This integral is logarithmically divergent; this is similar to the divergence of the electronic self energy in a Fermi liquid in two dimensions [70]. As in a Fermi liquid, intralayer short-range interaction between the spinons lead to a finite spinon lifetime τ ∝ 1/T 2 . The associated broadening of the spinon spectral function provides an infra-red cutoff for the logarithm [71], giving with Λ a high-energy cut-off, of the order of the Fermi energy (which is proportional to the exchange coupling between the original spins). The in-plane thermal conductivity of the Z 2 QSL with a Fermi surface is given, using the Einstein relation, by κ ab ∼ c V v 2 F τ /2, where c V = π 2 νT /3 is the specific heat of the system at low temperatures, v F = k F /m * is the Fermi velocity, and τ is the spinon lifetime. In a perfectly clean crystal, the lifetime comes from weak shortrange interaction between the spinons mediated by the gapped gauge field (assuming that Umklapp processes are available to relax the total momentum of the scattering spinons). The lifetime is given by τ −1 ∼ T 2 log(Λ/T ) as discussed earlier, and therefore we have: III. DISORDERED Z2 QUANTUM SPIN LIQUID As we shall now show, quenched disorder changes the low-temperature inter-plane transport in a qualitative way. The effects of disorder depend crucially on the type of disorder, which is subject to the symmetry of the problem. Consider, for example, the case of the honeycomb Kitaev model with time reversal symmetry. Then, disorder can take the form of a random bond strength, that translates to a random vector potential [72] in the lowenergy Dirac Hamiltonian, Eq. (1). Breaking time reversal symmetry can induce random scalar potential and mass terms, as well (see Appendix C 1 for a demonstration of how such terms arise in a disordered version of the Kitaev model).
Here, we will focus on random vector and scalar potentials; a random mass term is important at the transition between different gapped spin liquid states, a case we will not consider in the present work. The disordered part of the low-energy effective Hamiltonian in layer l is given by where V k−k and A k−k are random scalar and vector potentials, respectively. We assume that the disordered potentials in different layers are statistically independent. First, we study the case of a Dirac QSL with time reversal symmetry, in which only a random vector potential term is allowed, V k−k = 0. The effects of a vector potential disorder on a system with a Dirac dispersion were studied extensively in Ref. [73], where it was shown that such a term leads to a line of fixed points, characterized by scaling exponents which depend continuously on the disorder strength. Using the methods introduced in Ref. [73], we can find the scaling form of correlation functions at this fixed point, as described in detail in Appendix C 2 a. This allows us to show that vector potential disorder results in a modification of the exponent of the thermal conductivity, which is given by κ z ∼ T 5−α (disordered Z 2 with Dirac spectrum), (10) with α = 6∆ A /(π+∆ A ), ∆ A being the disorder strength: and the average is over disorder configurations. We consider smooth disorder, such that ψ † ψ † terms (corresponding to intervalley scattering in the Majorana model) are negligible. Next, we consider the effect of disorder on a Z 2 QSL with a Fermi surface, corresponding to ∆ = 0 in Eq. (1). In this case, since time reversal symmetry is broken, both scalar and vector disorder potentials are allowed. To simplify the computation, we will neglect the vector potential in this case, and assume that the scalar potential is short range correlated in space: V q V q dis = δ(q + q )/(2πντ ), where ν is the density of states at the Fermi level, and τ is the mean free time of quasi-particles at the Fermi surface. Moreover, we will again set the mass term in Eq.
(1) such that ∆ > m but |∆ − m| m. We expect none of the qualitative aspects of the solution to depend on these choices.
In the presence of disorder, the calculation of the thermal conductivity is most conveniently done using the Luttinger prescription [74][75][76]. The thermal conductivity is written as where Π(ω) is the retarded thermal current-thermal current correlation function, The c-axis thermal current operator can be derived using the energy continuity equation: , where h(q c ) is the energy density operator at wavevector q c . An explicit calculation to leading order in J ⊥ using Eq. (2) gives (see Appendices D 1, D 2) Here, we have suppressed the eigenstate indices λ 1,...,4 , since in the non-relativistic limit |∆ − m| m, the wavefunctions of states at the Fermi surface are confined to a single sublattice. Similarly, we have suppressed the eigenstate indices in F 0 , which is now momentumindependent. Note that, similarly to the inter-plane coupling, the thermal current operator in our model is quartic in the fermionic operators, corresponding to the fact that energy is carried between the plane by the hopping of fermion pairs.
The diagrams describing the leading-order contribution to κ c are shown in Fig. 2. The computation is lengthy but straightforward, and we will only describe the main steps here, deferring the details to Appendix D 5. We assume that the disorder is weak, such that k F 1, where k F is the Fermi momentum and = k F τ /m is the mean free path. Under these conditions, we may use the self-consistent Born approximation [77], equivalent to summing only non-crossed diagrams [78].
A key object is the disorder averaged four-point correlator within a single layer, Υ, depicted in Fig. 2 The thermal current correlation function, Eq. (C12), is then given as a convolution of two four-point correlation functions of two adjacent layers: . The clean, free fermion limit of this expression, with the correlator Υ(k 1 , k 1 , q; iν n , iν m ) = δ(k 1 − k 1 )G(k 1 , iν n )G(k 1 + q, iν m ), reproduces the Fermi golden rule calculation, Eq. (4). In the presence of disorder, the computation of Υ for small q (such that q 1) involves a summation over a ladder series (see Appendix D 5); this results in where D the diffusion constant D = v 2 τ /2, with τ the disorder induced single particle lifetime, and Note the appearance of the diffusion kernel in Eq. (17); this is related to the diffusive behavior of the dynamical charge correlation function in a disordered system. The computation of the sums in Eq. (16) is described in Appendix D 5. The dominant contribution comes from low frequencies and momenta, where the four-point correlator takes the form (17). At low temperatures, T < 1/τ , the result is At higher temperatures, T > ∼ 1/τ , κ c crosses over to the clean form, Eq. (7). Eq. (19) can also be derived from scaling arguments, assuming that the intra-plane densitydensity correlation function has a diffusion form; see Appendix C 2 b.
In Appendix D 6, we show that the F 4 pair hopping The disorder averaged four point correlator within each layer, Υ(k, k , q). The black lines denote fully dressed fermion propagators, dashed lines represent the effects of disorder, the squiggly lines are the bare thermal current vertex, while the green area stands for the fully renormalized two-particle vertex. We work in the self-consistent Born approximation, applicable for kF 1, where only ladder diagrams are taken into account. We suppress the frequency dependence for clarity.
inter-layer term results in the same power law, κ c ∼ T 2 .

IV. U (1) QUANTUM SPIN LIQUID-
We further study the case of a layered U (1) QSL with a spinon Fermi surface. In addition to the fermionic spinons, there exist gapless gauge field photons, which also contribute to transport. The low energy sector of each layer is described by the Lagrangian density [20][21][22][23][24][25][26] where ψ † l,σ creates a spinon at layer l with spin σ, (a 0 , a) is the U (1) gauge field, µ is a chemical potential that sets the size of the spinon Fermi surface and m is the spinon effective mass. A "Maxwell" term for a ν , 1 2g νλ f νλ f νλ where g is a coupling constant and f νλ = ∂ ν a λ − ∂ λ a ν , is also allowed by symmetry; however, it gives rise to subleading contributions at low momenta and frequencies, and hence we will drop it in the following.
Under the random phase approximation (RPA), the clean system is described by a strong-coupling fixed point, with the retarded gauge boson and spinon propagators (D R (q, ω) and G R (k, ω), respectively) given by Fermi momentum. The use of the RPA has been formally justified in a large-N expansion, where N is the number of fermion flavors [26], but this has been shown to be problematic [79]. Additional expansion parameters have been proposed, that essentially reproduce the RPA results [80, 81]. We shall use the RPA approximation, assuming it is pertinent to at least some area in parameter space.
In a layered U (1) QSL, heat may be transferred between the layers both by spinon and photon excitations. The most relevant inter-layer interaction term of each sector is given by where a T is the transverse part of the gauge field. The coupling functions F sp and F ph depend on the spatial structure of the inter-layer coupling; their explicit form is unimportant. In real space, the gauge invariant term ∇ × a is related the chirality of the underlying spin degrees of freedom [25], and therefore the J ph ⊥ term corresponds to an interaction between the chiralities of the spin textures in the two layers. Micropically, this term may be small compared to J sp ⊥ , since it is of higher order in the inter-plane Heisenberg exchange coupling. However, as we shall see below, in a clean case, it gives a dominant contribution to κ c at asymptotically low temperatures.
The calculation of the spinon-mediated inter-plane thermal conductivity proceeds in a similar fashion as in the Z 2 QSL case. κ c is given by a similar expression to Eq. (D20) (with the replacement F → F sp ). It is given by where W is an appropriate UV cut-off (see Appendix E). However, in the clean case, the dominant source of low-T thermal transport turns out to be the exchange of gauge fluctuations; this contribution may also be calculated by the Kubo formula, and is given by (see Appendix E for details) χ 2 k 6 +γ 2 2 is the photon spectral function.
Thus, at sufficiently low temperature, κ c,ph κ c,sp . Note that the thermal conductivity can be written as κ c,ph ∼ T 2−1/z , where z = 3 is the dynamical critical exponent of the fixed point described by RPA.
The introduction of disorder to the U (1) theory is likely to destabilize the z = 3 fixed point, leading instead to diffusive behavior, similar to that of a disordered Fermi liquid. In the RPA approximation, the propagators of the disordered theory are given by [82] with D a diffusion constant and τ the disorder-induced finite lifetime. The calculation of the c-axis thermal conductivity is then similar to the disordered Z 2 QSL case. Inserting Eq. (25) in the Kubo formula for the c-axis conductivity leads to for both spinon and photon contribution.

V. EXPERIMENTAL CONSIDERATIONS
In this section, we discuss possible experimental candidate systems where thermal conductivity provides a gateway to observing QSL physics. In order to observe the magnetic contribution to the inter-layer thermal conductivity, one has to be able to separate it from the phonon contribution. Since the phonon contribution scales as T 3 , the magnetic contribution in certain QSLs dominates at sufficiently low temperatures. This happens in QSLs with a disordered spinon Fermi surface and in strongly disordered Dirac QSLs (see Table I). Below, provide a rough order-of-magnitude estimate for the temperature T * at which the magnetic contribution exceeds the phonon one, as a function of system parameters (such as the strength of the inter-plane coupling, the Debye temperature, and the disorder strength). As we elaborate below, this estimate indicates that at least in some material candidates, the crossover to magnetically dominated thermal transport may occur at accessible temperatures.
We base our estimate of T * on the case of a QSL with a spinon FS, whose magnetic c-axis thermal conductivity is given by Eq. (19). We set the unit of length to be the lattice spacing a, and estimate ν ∼ 1/J, D = 1 2 v F sp ∼ 1 2 J sp , where J is the in-plane exchange coupling, and sp is the spinon mean-free path in the plane. This gives Next, we estimate the contribution of the phonons. The acoustic phonon specific heat is where c s is the sound velocity. Therefore, by the Einstein relation, The temperature below which the spinon contribution to the thremal conductivity becomes larger than the phonon contribution is given by equating (27) to (28). The result is Eq. (29) highlights the parameters that control T * : T * is higher the stronger the disorder, the higher is Θ D , and the smaller is J. [Note that Eq. (27) is only valid for T J; therefore, T * in Eq. (29) cannot exceed J].
As an illustrative example, we roughly estimate the crossover temperature T * for kapellasite, a kagome gapless QSL candidate [83]. This is a polymorph of Herbertsmithite; however, the in-plane exchange coupling is about an order of magnitude smaller. The exchange couplings of kapellasite have been estimated from from firstprinciple calculations [84]: J ≈ 10K, J ⊥ ≈ 0.5K. We assume that kapellasite has a spinon Fermi surface, and that the Debye temperature is Θ D ∼ 300K. The mean free paths of the spinons and the phonons are not known. However, disorder in the planes is believed to be substantial. To get a rough estimate of the order of magnitude of T * , let us assume a strongly disordered sample, such that sp = 20a and ph = 200a. This gives T * ≈ 6 0.5 2 · 300 2 10 3 · 20 · 200 ≈ 35 mK.
kapellasite does not order magnetically at least down to 20 mK [83]. Thus, for sufficiently strong disorder, we get that the crossover temperature is within experimental reach.
Let us discuss other QSL candidate materials where the spinon contribution to κ c may be measurable. A promising candidate material is the recently discovered 2d spin-orbit coupled iridate H 3 LiIr 2 O 6 , which has been observed to be paramagnetic to very low temperatures, and hosts gapless excitations [85,86]. Compared to other similar compounds like Na 2 LiO 3 and Li 2 IrO 3 (which or-der at low temperatures), in H 3 LiIr 2 O 6 the interlayer distance is smaller due to replacement of Li by smaller H atoms in between layers, which increases J ⊥ . Further, the in-plane bond length is also larger which reduces the scale of in-plane exchange interactions J. As per Eq. (29), both these factors are conducive to a larger crossover temperature T * where the magnetic contribution becomes large.
Other candidate materials are magnetic insulators with strong spin-orbit coupling, where the Kitaev interaction is the dominant term. Some of these materials, like α−RuCl 3 , are believed to be proximate to a QSL phase [87]. Further, the magnetic order can be suppressed by doping, making such materials an interesting playground for observing spin-liquid physics [88], although the nature of the field-induced QSL phase is still unclear.
In the layered organic insulators [6], the inter-plane exchange coupling is estimated to be three orders of magnitude below the intra-plane coupling [89], and therefore it is likely that phonons dominate the c-axis thermal transport at accessible temperatures. Herbertsmithite [6] is believed to have a gapped QSL ground state [46], although the spin gap seems to be quite small (∆ gap < ∼ 10K [46,90]). An applied magnetic field can induce a finite spinon density of states at zero energy, opening the way to measure the spinon contribution to κ c . However, the in-plane exchange coupling J is about an order of magnitude larger larger than in kapellasite, while the ratio J ⊥ /J is comparable in the two systems [84]. Therefore, we expect T * in Herbertsmithite to be smaller than in kapellasite.
Finally, we discuss a few techniques can be used to isolate the magnetic contribution to the thermal conductivity from that of phonons.
(i) In gapless spin liquid candidates where the magnetic contribution is a power law of the form T θ with θ < 3, one can isolate the magnetic contribution from the phononic one (which scales as T 3 ), since the magnetic contribution is dominant at low sufficiently low temperature. Plotting a curve of κ/T θ vs. T 3−θ , the slope of the curve gives us the phonon contribution, while the intercept gives us the magnetic contribution to the thermal conductivity. This is possible as long as the sample temperature is not much higher than T * .
(ii) In addition, in some materials an applied magnetic field may be used to establish long range order, suppressing the spinon contribution to the thermal conductivity, while weakly affecting the phonon contribution. Contrasting the measurements of the c-axis thermal conductivity in the presence and absence of such a field may enable us to isolate the spinon contribution.

VI. CONCLUSIONS
We have studied the thermal conductivity in layered, gapless QSLs. The key observation is that the mechanisms of in-plane and out-of-plane thermal transport are qualitatively different: the former is carried by fractionalized excitations, while the latter is carried by gaugeneutral, non-fractionalized excitations. Thus, in all the cases we have studied, κ ab and κ c follow different power law dependences at low temperature; in particular, the anisotropy κ ab /κ c diverges in the limit T → 0. This property is a clear hallmark of a fractionalized, layered system. A large number of layered QSL candidates have been proposed in the last few years, and inter-plane thermal conductivity can serve as an unambiguous probe for fractionalization in these experimental candidates.   [86] Kevin Slagle, Wonjune Choi, Li Ern Chern, and Yong Baek Kim, "Theory of a quantum spin liquid in hydrogen-intercalated honeycomb iridate, h3liir2o6," arXiv preprint arXiv:1710.01307 (2017). We model the layered Z 2 QSL system as layers of the Kitaev honeycomb model, coupled by a weak inter-layer interaction. The Kitaev honeycomb model [16] is an exactly solvable model of interacting spin-1/2s. It is composed of a honeycomb lattice of spins interacting via direction-dependent exchange interactions, where j, k are nearest neighbors on the hexago lattice, and S α jk are the x, y, or z component of the spin operator, depending on the type of link between j and k. The links are denoted x, y, or z, based on their orientation, as shown in Fig. 3

. Each of the spins is represented in terms of Majorana fermions
however, the representation in terms of these fermions spans a larger Fock space, and must be restricted to the physical Hilbert space of the spins by the gauge On each α-direction link, u α ij = ib α i b α j is conserved, and a theorem by Lieb [91] guarantees that the ground state is in the sector where it is possible to set u α ij = 1 (the flux-free sector). Thus, the ground state manifold is described by the free Majorana Hamiltonian with A i and B i+δ the c i Majorana on the A and B sublattice, respectively. The δs are the three vectors connecting the even and odd sublattices, as shown in Fig. 3.
In order to probe the Z 2 QSL with a Fermi surface, we consider further two specific time-reversal breaking terms, which result in a Fermi surface without creating vison excitations that would take us out of the ground state manifold: Here, the sites labeled 1...6 on each plaquette are shown in Fig. 3. In the ground state manifold, these terms are given by We comment that on different lattices, one can also obtain a QSL with a spinon Fermi surface even in presence of TRS [64][65][66][67][68][69].
Lastly, we consider also the effect of disorder in the system via a term H dis , the exact form of which will be given later. Thus, the spin liquid we consider is described by with H T RB = 0 in the TR invariant case.

Inter-layer coupling
The most relevant interlayer coupling terms are those that leave each layer in its ground state; that is, the interlayer coupling must commute with all the u α ij s. A general form of such a tunneling term in the language of the original spins will consist of spin operators from one layer coupled to spin operators from another. In order to maintain exact solvability, we consider the interlayer coupling term H ⊥ = J ⊥ l,l =l±1 plaquettes S z 1,l S y 2,l S x 3,l − S x 2,l S z 3,l S y 4,l + S y 3,l S x 4,l S z 5,l − S z 4,l S y 5,l S x 6,l + S x 5,l S z 6,l S y 1,l − S y 6,l S x 1,l S z 2,l × S z 1,l S y 2,l S x 3,l − S x 2,l S z 3,l S y 4,l + S y 3,l S x 4,l S z 5,l − S z 4,l S y 5,l S x 6,l + S x 5,l S z 6,l S y 1,l − S y 6,l S x 1,l S z 2,l .
In the Majorana representation this term is given by (again settingû j,δ = 1) We have chosen this term as it allows for simple calculations, and in particular it reduces to a "density-density" interlayer interactions [the F 0 term in Eq.
(2)] in the continuum limit. We expect the exact form of the coupling term to be unimportant; the important point is that it must have at least two spinon operators from each layer. This is because only a spinon pair excitation can be transferred between different layers, and this requires at least two spinon operators from each. We will comment on the case of a more generic form of the inter-layer Hamiltonian, having on spin operator in each layer, in Sec. A 4 below.

Continuum Hamiltonian
Using f k = δ e ik·δ , g k = i sin(k · n i ), ψ A l (k) = e iπ/4 j A l j e ik·Rj , ψ B l (k) = e −iπ/4 j B l j e ik·Rj , ∆ k = 4J 1 T RB g k , and m k = 4J 2 T RB g k (where the vectors n i=1,2,3 are defined in Fig. 3), the disorder free Hamiltonian of each layer is given by The low energy theory is centered near the Dirac points K, K = 2π 3 ± 1 √ 3 , 1 . Near these points, Jf k ≈ v(kx±ik y ), with v = 3J/2. The low energy in-plane Hamiltonian can thus be written as two Majorana theories with a Dirac dispersion. It is convenient to consider an equivalent system, of complex fermions which reside only on half the Brillouin zone; this makes use of the equivalence ψ A (k) = ψ A † (−k). The low energy theory of this system is given by a single Dirac cone centered at K, and its in-plane Hamiltonian is given by Here ψ †l k = (ψ A † l (k), ψ B † l (k)) is a spinor of complex fermions, σ = (σ x , σ y ) a vector of Pauli matrices, ∆ = ∆ k=K , and m = m k=K . This is Eq. (1) in the main text.
In terms of the continuum theory, the inter-plane term in Eq. (A7) is given by (neglecting the variation of ∆ k around the Dirac points) We expand the g k form factors for small deviations away from the Dirac point K; we set F 0 = g 2K , while the lowest order term in the expansion of g k−k away from the K point vanishes, introducing additional factors of momentum. This will introduce additional factors of temperature T in the contribution to the thermal conductivity, and we therefore neglect the pair hopping term.
We work in the basis of the eigenstates of H l 0 + H l T RB . In the TR symmetric case, where ∆ = m = 0, the eigenstates are given by a λ (k) = [ψ A (k) + e iφ k ψ B (k)]/ √ 2, where φ k is the angle between k x and k y , and their energies are λ k = λvk. For simplicity, in the analysis of the Z 2 QSL with a Fermi surface, we consider the regime ∆ > m > 0, m ∆ − m. In this limit, the eigenstates with energy close to the Fermi surface are located almost entirely on the A sublattice, and we may ignore the sublattice degree of freedom; these eigenstaes are denoted by a(k) ∼ ψ A (k). In this basis, the single layer Green's functions are G(k, iν n ) = [iν n − k 2 /2m * + µ] −1 , with µ = ∆ − m, and m * = v 2 /m.

Generic inter-layer Hamiltonian
The inter-layer coupling term (A6) is designed to maintain the exact solvability of the model, and for computational convenience. However, generically, we expect the largest components of the inter-layer coupling Hamiltonian to be quadratic in the spin operators. Within the Kitaev model, a quadratic inter-layer coupling term (such as a Heisenberg term, J ⊥ α S α l S α l+1 with α = x, y, z) does not merely create a pair of spinon excitations in each layer. Rather, it creates a pair of spinons and a pair of gapped vison (flux) excitations. In order to annihilate the pair of flux excitations and return to the low-energy subspace, we have to apply the J ⊥ term again. Thus, it appears that the effective inter-plane interaction in the low-energy effective Hamiltonian (2) is proportional to J ⊥ ∼ (J ⊥ ) 2 /∆ v , where J ⊥ is the "microscopic" strength of the inter-layer coupling, and ∆ v is the vison gap. One may then worry that the effective inter-plane interaction J ⊥ is too small to contribute significantly to the inter-layer thermal conductivity.
However, we argue that for a generic intra-layer Hamiltonian that contains also non-Kitaev terms (such that even the intra-layer Hamiltonian is not exactly solvable), this is not the case; in fact, J ⊥ ∝ J ⊥ . This is because in the generic case, the vison excitations are not static, even within the single-layer Hamiltonian. A pair of vison excitations can annihilate each other without the need for another application of the inter-layer J ⊥ term.
To illustrate this, consider the case where there is an additional intra-plane interaction: where i, j denotes two nearest-neighbor sites i, j on the honeycomb lattice, and Γ αβ is a 3 × 3 symmetric matrix. Such terms are present in real "Kitaev materials" [92]. Consider a quadratic inter-plane coupling term of the form J ⊥ S z 1,l S z 1,l+1 , where the position of the site 1 is indicated in Fig. 3. We can now derive the effective inter-plane coupling J ⊥ (that creates a pair of spinons in each layer and does not create any visons) perturbatively in both J ⊥ and Γ αβ . One can check explicitly that acting with the following sequence of operators: amounts, in a certain gauge, to acting c 1,l c 2,l c 1,l+1 c 2,l+1 and not changing the number of visons in either layer. The Thus, the term in the effective Hamiltonian that creates a pair of fermionic spinons in each of the adjacent layers l, l + 1 is proportional to J ⊥ . Generically, there is no reason to expect Γ αβ to be much smaller in magnitude than ∆ v , since both energy scales characterize the intra-plane Hamiltonian, and do not involve inter-layer coupling. We conclude that in a generic situation, J ⊥ and J ⊥ are of the same order of magnitude. k = vδk x sin θ + vδk y cos θ, k+q = −vδk x sin θ + vδk y cos θ, k = −vδk x sin θ − vδk y cos θ, k −q = vδk x sin θ − vδk y cos θ. (B2) Here, sin θ = q 2k F (see Fig. 4), and v = k F m * is the Fermi velocity. The integrals over δk, δk can now be performed easily, giving (B3) The integral over 1,2,3 can now be performed; by scaling, this integral is proportional to T 5 . Substituting θ with q, we get that . (B4) Using ν = k F 2πv F , we arrive at Eq. (4) in the text. The logarithmic divergence comes from small q scattering, as well as from scattering with momentum transfer close to q = 2k F . The inelastic lifetime of the quasi-particles in each layer, τ ∝ 1/T 2 , is needed in order to cut off this divergence. There is another contribution from q > 2k F ; this contribution is parametrically the same as Eq. (B4). Appendix C: Disordered Z2 QSL

Possible forms of disorder
In a TR invraiant system, we consider disorder in the spin-spin couplings J: where again j, k are nearest neighbors and α jk = x, y, z, according to the type of link. In the Majorana fermion representation, this becomes with A k,k x = 1 J j,δ [e −ik δ δJ j,δ e i(k−k )Rj ], A k,k y = 1 J j,δ [e −ik δ δJ j,δ e i(k−k )Rj ]; thus, the low energy, longwavelength disorder is of the form of a random vector potential.
Disorder which affects the next nearest neighbor hopping results from three-spin interaction terms in the original spin model; these terms break time reversal symmetry. Depending on the relative sign of the disorder between the A and B sublattices, just as in Eq. (A4), these terms will give a rise to a scalar potential term or a random mass term

Inter-layer thermal conductivity
To compute the c-axis conductivity in a disordered, layered Z 2 QSL, we compute the rate at which pairs of spinon excitations tunnel between planes. This can done using the Fermi golden rule, analogously to Eq. (4), replacing the momentum eigenstates with eigenstates of the disordered intra-plane Hamiltonian.
The general expression for the thermal conductivity is Here, H ⊥ (l, l ) is the part of the inter-plane Hamiltonian that couples layers l and l . |i, f are the initial and final many-body eigenstates of the system with J ⊥ = 0 (decoupled planes), with corresponding energies E i,l and E f,l (at layer l). Since we neglect intra-plane interactions, we can expand the fermionic spinon operators in the basis of the single-particle eigenstates of each layer (which include the effects of disorder): where η = A, B labels the sublattice, and f λ,l annihilates an eigenstate with energy ξ λ,l , which has the wavefunction ϕ l,λ,η (x). As in the main text, we will mostly work with an inter-layer Hamiltonian of the "density-density" form, H ⊥ (l, l ) = η,η d 2 x ψ η † l (x)ψ η l (x)ψ η † l (x)ψ η l (x), commenting along the way about other forms of H ⊥ . Using Eq. (C6), we can write the disorder-averaged c-axis thermal conductivity as where . . . dis represents disorder averaging. We denote such that the thermal conductivity is given by The remaining task is to compute the function g l (q, ε, ε ) for a Z 2 QSL with either a Dirac spectrum or a Fermi surface.

a. Disordered Dirac
The properties of two-dimensional Dirac fermions coupled to a random vector potential, that corresponds to a disordered Z 2 Dirac QSL with time reversal symmetry, has been studied extensively in Ref. [73]. Here, we will briefly review some of the results of Ref. [73], and use them to determine the scaling of the c-axis thermal conductivity with temperature.
Since the problem is non-interacting, the actions for different frequency modes decouple (before disorder averaging). The ω = 0 system is described by a fixed line of interacting theories in d = 1 + 1 dimensions [73]. The frequency ω corresponds to a relevant operator with scaling dimension 2 − z, where z = 1 + ∆ A /π is the dynamical critical exponent, and ∆ A is the disorder strength, defined in Eq. (11). I.e., under scaling, q → q = q/b, ω → ω = ω/b z . The fixed line is characterized by ω/T scaling.
The function g l (q, ε, ε ) satisfies the scaling relation where y is a critical exponent related to the scaling dimension of the fermion density operator, which we compute below, and b is a rescaling factor. Choosing b = |ε| −1/z , we get that g l can be written as where Φ is a universal scaling function.
To determine y, we notice that g l is related to the density-density correlator: where n F (ε) is the Fermi function. Using Eq. (C11), this can be written as Here, we have used Eq. (C11) and performed a change of variables, ε = T ξ, ε = T ξ . On the other hand, χ(q, ω n = 0) can be expressed as Here, n l,ν (x) is the Matsubara frequency ν component of the density operator in layer l [93]. In the second to last line we have applied scaling to the correlation function n l,ν (x)n l,ν (0) , using the fact that the scaling dimension of n l,ν (x) is 2 − z. [73] Choosing b = T −1/z in Eq. (C14), we get that χ(q, ω n = 0, T ) = T 2−z z Ψ(q/T 1/z ), where Ψ is a scaling function. Comparing this to Eq. (C12), we can extract the exponent y: Now, we are in a position to find the scaling of the c-axis thermal conductivity with temperature. Inserting Eq. (C11) into Eq. (C9) results in Rescaling the integral,ω = ω/T ,ε 1,2 = ε 1,2 /T , andq = q/T 1/z , we get that This result coincides with that of the clean case in the limit ∆ A → 0 (i.e. z = 1 + ∆ A /π → 1).
The analysis above has been done for an inter-plane interaction of the density-density form. A similar analysis can be done for any quartic inter-plane interaction. The only difference is the scaling dimension of the fermion bilinear operator that appears in the interaction term, that can be determined using the methods of Ref. [73]. It turns out, however, that for any ∆ A > 0, the density operator is the fermion bilinear with the smallest scaling dimension. Hence a density-density interaction gives the dominant contribution to κ c at low temperatures.

b. Disordered FS
In the disordered FS case, we know that the density-density correlation function takes a diffusive form at small q, ω n : where D is the diffusion constant, and ν is the density of states at the Fermi level. Comparing this to Eq. (C12), we deduce that g l (q, ε, ε ) should satisfy the following scaling relation: Hence, g l (q, ε, ε ) can be written as where Ω(ξ, ξ ) is a dimensionless scaling function.
Starting from Eq. (D20), we perform the Sommerfeld expansion with respect to 1 . The first non-vanishing contribution, in powers of T , occurs for the term: