Quantum critical thermal transport in the unitary Fermi gas

Strongly correlated systems are often associated with an underlying quantum critical point which governs their behavior in the finite temperature phase diagram. Their thermodynamical and transport properties arise from critical fluctuations and follow universal scaling laws. Here, we develop a microscopic theory of thermal transport in the quantum critical regime expressed in terms of a thermal sum rule and an effective scattering time. We explicitly compute the characteristic scaling functions in a quantum critical model system, the unitary Fermi gas. Moreover, we derive an exact thermal sum rule for heat and energy currents and evaluate it numerically using the nonperturbative Luttinger-Ward approach. For the thermal scattering times we find a simple quantum critical scaling form. Together, this determines the heat conductivity, thermal diffusivity, Prandtl number and sound diffusivity from high temperatures down into the quantum critical regime. The results provide a quantitative description of recent sound attenuation measurements in ultracold Fermi gases.

Strongly correlated systems are often associated with an underlying quantum critical point which governs their behavior in the finite temperature phase diagram. Their thermodynamical and transport properties arise from critical fluctuations and follow universal scaling laws. Here, we develop a microscopic theory of thermal transport in the quantum critical regime expressed in terms of a thermal sum rule and an effective scattering time. We explicitly compute the characteristic scaling functions in a quantum critical model system, the unitary Fermi gas. Moreover, we derive an exact thermal sum rule for heat and energy currents and evaluate it numerically using the nonperturbative Luttinger-Ward approach. For the thermal scattering times we find a simple quantum critical scaling form. Together, this determines the heat conductivity, thermal diffusivity, Prandtl number and sound diffusivity from high temperatures down into the quantum critical regime. The results provide a quantitative description of recent sound attenuation measurements in ultracold Fermi gases.

I. INTRODUCTION
Thermal transport caused by temperature gradients is ubiquitous in nature and typically occurs in a diffusive manner. A calculation of the corresponding thermal conductivity κ and the associated thermal diffusion constant D T = κ/c p is often based on a kinetic theory description like the Boltzmann equation. This works well, e.g., in metals at low temperature and allows one to understand the origin of universal laws like the Lorenz ratio L = κ/σT → L 0 = π 2 k 2 B /3e 2 between the thermal and the electrical conductivity σ as predicted by Wiedemann and Franz. In strongly correlated systems, sometimes called bad metals [1], the underlying Fermi liquid description does not apply, however, and L deviates substantially from its ideal value L 0 [2,3]. Developing a microscopic theory for thermal transport in non-Fermi liquids has been a major challenge for many years. It has been approached using a number of different techniques like the memory function formalism [4]. In a number of cases, a possible and phenomenologically often successful strategy to describe transport in the absence of welldefined quasiparticles is to assume the existence of an underlying quantum critical point (QCP). Transport in the quantum critical regime above the QCP may then be analyzed in terms of critical fluctuations where decay and scattering rates typically scale linearly with temperature according to a Planckian law τ −1 ∝ k B T / [5][6][7], a behavior which has been observed recently in the thermal diffusivity of near optimally doped cuprates above the superconducting transition [8]. The aim of our present work is to develop a microscopic theory for thermal transport in a much simpler system with a quantum critical point, namely the unitary Fermi gas (UFG) [9,10]. This system has a QCP at zero density which is both scale and conformally invariant [11][12][13][14]. In the quantum critical regime above this point, the thermal wavelength gases have realized homogeneous Fermi gases [24,25] and opened the possibility to access local thermal transport via the diffusive spreading of density and thermal wave packets propagating in a sufficiently large box [23,[26][27][28]. These experiments are considerably more sensitive than previous global transport measurements from trap collective modes. For instance, the measurements of the hydrodynamic sound dispersion ω q = c s q − iD sound q 2 /2 + · · · in a homogeneous unitary Fermi gas at MIT [23] provide both the speed of sound c s and the sound diffusivity D sound = (4/3)D η + (c p /c V − 1)D T [29]. Knowledge of the kinematic viscosity D η = η/(mn) [16,17] and the Landau-Placzek ratio LP = c p /c V −1 (Fig. 7 below) gives then access to the thermal diffusivity D T in the quantum degenerate gas, see Fig. 1. Theoretical results for thermal transport are so far available only at high temperature from the virial expansion [30]. It is the goal of this work to compute thermal transport at low temperature and in particular in the quantum critical regime.
In the following, we compute thermal transport in the quantum critical region of the unitary Fermi gas based on a decomposition of the thermal conductivity as a product of a nontrivial, thermodynamic sum rule χ T qq for the heat current response and a thermal scattering time τ κ which can formally be derived within a memory function approach, cf. Sec. II. We show that both factors of this decomposition can be described by universal scaling forms which smoothly connect the quantum critical to the high-temperature regime, where a virial expansion for the thermodynamic properties and a Boltzmann equation for the associated scattering time is applicable. In Sec. III we derive an exact expression for the thermal sum rule χ T qq in terms of Green's functions with the help of Ward identities for energy and particle number conservation. Based on nonperturbative results for the Green's functions from a fully self-consistent Luttinger-Ward computation [31][32][33] we evaluate χ T qq numerically. We find a strong enhancement of spectral weight in the quantum critical regime compared to the noninteracting gas which reaches two orders of magnitude in the quantum critical regime just above the superfluid transition. In Sec. IV we compute the thermal scattering time τ κ of order T / using a large-N expansion. Quite unexpectedly, the time τ κ extrapolates in a simple manner from the Boltzmann gas limit down into the quantum critical regime. In Sec. V, we combine the results for the sum rule with the scattering times in Eq. (1) to predict the thermal transport coefficient κ, the diffusivity D T shown in Fig. 1, and the Prandtl number Pr. In particular, we find good agreement with the experimentally observed values in the quantum critical regime. We conclude with a discussion in Sec. VI.

Nondegenerate gas
FIG. 2. Phase diagram of the spin-balanced, unitary Fermi gas at finite temperatures [34]. The QCP at T = 0, µ = 0 is the starting point for the phase boundary of the homogeneous superfluid at Tc = 0.4µ (solid line). The dashed lines mark the crossover to the quantum critical regime above the QCP.

II. QUANTUM CRITICAL THERMAL TRANSPORT
In this section we first define the quantum critical regime of the unitary Fermi gas in part II A, and discuss the crossover to classical critical behavior close to the finite-temperature superfluid transition. In part II B, we discuss the formal structure of how to compute thermal transport in linear response from the Kubo formula and its evaluation within the memory function formalism.

A. Quantum critical regime
Dilute ultracold Fermi gases interact via a short-range attractive interaction between different spin components. At low temperature, atoms scatter predominantly in the s-wave channel with scattering amplitude f (k) = −1/(a −1 + ik), which is fully characterized by the s-wave scattering length a. Here we focus on the unitary limit 1/a = 0 that gives rise to a strongly interacting system as the standard perturbative expansion in a small gas parameter n|a| 3 1 breaks down. The phase diagram shown in Fig. 2 exhibits a quantum critical point at vanishing chemical potential and temperature µ = T = 0, which separates the vacuum state at µ < 0 from a homogeneous superfluid (SF) state at µ > 0 [10,11,34]. Here, all energies are expressed in terms ofĒ, which is of the order of the van der Waals energy that sets the cutoff scale beyond which details of the interaction potential start to matter. The universal description based on the model Hamiltonian (10) below is thus applicable only for µ, T Ē . In the absence of a finite effective Zeeman field h = (µ − µ )/2 which may lead to nontrivial phases with a finite spin population imbalance [32,33] the phase diagram is characterized by a single dimensionless parameter βµ. The superfluid state remains stable for temperatures below the critical curve T c 0.4µ or equivalently (βµ) c 2.5 [35]. Instead, for high temperature or small fugacity z = e βµ 1 the system forms a dilute, non-degenerate gas which can be described in terms of the virial expansion. Increasing the fugacity z to values of order unity one enters the QCR, as shown in Fig. 2. In this regime, both thermodynamic and transport properties follow universal scaling laws associated with the zero density fixed point at T = µ = 0, with βµ as the single relevant scaling variable [11].
The quantum critical scaling is replaced by the one characteristic for a classical d = 3 XY model close to the superfluid phase transition at µ c (T ) 2.5T . This crossover occurs when the Gaussian correlation length ξ G of the quantum model-defined by the quadratic term 1/ξ 2 G = 2m(µ − µ c )/ 2 in the Ginzburg-Landau free energy-becomes of the same order as the characteristic length ξ 1 . The length ξ 1 1/u 0 is associated with the coefficient u 0 of the quartic term (u 0 /4!) φ 2 of the classical φ 4 theory for a complex scalar field φ(x) = φ 1 (x) + iφ 2 (x) that depends only on space. This term may in principle be derived from the usual complex order parameter ψ(x, τ ) for the superfluid transition by integrating out all nonzero Matsubara frequencies Ω n = 0. In explicit form, this has been worked out for a generalization of the proper N = 2 component model of a weakly interacting Bose gas to a large number N , which yields u BEC 0 (a) = 96π 2 a/λ 2 T in the N = ∞ limit [36]. In the case of the unitary Fermi gas at 1/a = 0, simple dimensional analysis requires that u 0 1/λ T , however the value of the numerical prefactor is unknown. Qualitatively, the crossover condition ξ 1 ξ G thus gives the simple relation µ − µ c (T ) k B T , which entails a Ginzburg parameter of order unity and a very large Ginzburg region that extends up to about 2T c as discussed by Debelhoir and Dupuis [37].
In the vicinity of the superfluid transition, the quantum critical scaling of dynamical quantities is replaced by the classical dynamical scaling. In particular, the thermal conductivity of the UFG is described within Model F [38] with dynamical critical exponent z = 3/2 for the superfluid transition in the universality class of the d = 3 XY model. As shown by Ferrell et al. [39], this implies a divergent thermal conductivity as T → T + c , which diverges with an exponent close to 1/3 since ν ≈ 2/3.

B. Linear response and memory function formalism for the thermal conductivity
A formally exact expression which in principle allows to calculate transport coefficients for an arbitrary form of the underlying Hamiltonian is based on linear response theory and the resulting Kubo formula. In the special case of the thermal conductivity at external momentum q = 0, it is convenient to consider the heat current density j q [40], which is defined as the energy current j E at constant particle number, i.e., with the enthalpy per particle w/n times the number current density j subtracted: Here, we have used the Gibbs-Duhem relation w = ε + p = µn + T s and defined the entropy per particlẽ s = s/n. In standard hydrodynamic terms this corresponds to the definition of the thermal conductivity via Fourier's law j E = −κ∇T in the absence of a particle current. Microscopically, the effect of a weak temperature gradient is encoded in the equilibrium retarded heat current response function from linear response theory, where we suppress the argument q = 0 from now on. The retarded commutator in Eq. (4) defines a positive and even spectral representation of the dynamic thermal conductivity Since a fully microscopic evaluation of the frequencydependent response function in a strongly interacting many-body system is impossible, it is necessary to reduce the problem by restricting attention to the dc-response and a simplified procedure to evaluate the characteristic time scale τ κ defined in Eq. (1). Such a procedure is provided by the memory function formalism. It has been used to determine the dynamical charge conductivity of simple metals some time ago by Götze and Wölfle [41] and it provides a systematic and unified description for the derivation of hydrodynamic equations of motion in fluids, see, e.g., the textbook by Forster [29]. More recently, the method has been applied successfully to calculate transport coefficients in systems without well-defined quasiparticles [4,6,42]. It is based on a formal expression for the Laplace transform of the relaxation function in terms of a matrix χ T AB of static thermodynamic susceptibilities of slow variables A, B [43] and an associated frequency-dependent memory matrix M AB (z) (we assume that the operators A and B have the same sign under time reversal, otherwise an additional contribution appears in the denominator). Provided that this matrix has a finite limit M (0) at vanishing frequency, this leads to an expansion of the dynamical response function at low frequencies, which defines a matrix of relaxation times Identifying κ AB T = φ AB (z = 0) as the dc-transport coefficient, this leads to κ AB T = χ T AC τ CB , which is precisely of the form given in Eq. (1). In principle, therefore, the memory function formalism determines transport coefficients in quantum many-body systems in terms of the matrix χ T AB of associated static thermodynamic susceptibilities and the zero-frequency limit M AB (z = 0) of the memory matrix. The formal expression for M AB (z) shows that it is again a relaxation function but now for operatorsQȦ in which the dynamics of the slow variables A, B is projected out byQ = 1 − AB (χ T ) −1 AB |A) (B|. In practice, the memory matrix can hardly be determined exactly. In systems without long-lived quasiparticles, however, even approximate results for the scattering times are often sufficient because the only exactly or approximately conserved quantities are then particle number, momentum and energy while all other variables relax on microscopic time scales. In fact, much of the nontrivial structure of transport coefficients near quantum critical points is determined by the associated thermodynamic susceptibilities, which is behind the success of the memory function formalism in this context. This turns out to be the case also for the unitary Fermi gas studied here. Indeed, as will be shown in Sec. III, the relevant susceptibility χ T qq exhibits a rather strong dependence on the scaling variable βµ (see Fig. 4 below), while the scattering time in Fig. 5 evolves rather smoothly from the high temperature limit down into the quantum critical regime, essentially extrapolating the result obtained from a Boltzmann equation calculation. A similar situation also applies to other transport coefficients, such as the shear viscosity η = pτ η , where the sum rule is given by the pressure p [16], or particle transport, where an analog of the product form (1) for the thermal conductivity has also been found to hold.
In the following, we will determine the thermodynamic susceptibility χ T qq by a direct Green function approach, using an extension of exact Ward identities first derived by Polyakov [44] in the context of transport in the vicinity of a thermal critical point. Since the heat current is an ergodic variable, the result must coincide with the associated dc-susceptibility, which is given by the standard Kramers-Kronig relation, as the integral of the frequency dependent heat conductivity κ(ω, q = 0) times the temperature, including a proper regularization of the divergences which arise as a result of the assumption of a zero-range interaction in Eq. (10) below.

III. THERMAL SUM RULE
In this section we first introduce the model for the interacting Fermi gas in part III A, and we then express the linear response theory for thermal transport in the field theoretical formulation based on Green's functions (part III B). In particular, we derive a new Ward identity for the interaction part of the heat current, which gives rise to a novel exact expression for the thermal sum rule (25)- (27) in terms of one-and two-particle Green's functions. Next we discuss in part III C the necessary regularization of the high-momentum asymptotics. Finally, we numerically evaluate the sum rule in the quantum critical regime using the nonperturbative Luttinger-Ward approach in part III D.

A. Model
The many-body physics of ultracold Fermi gases can be described by the Hamiltonian for two-component fermions with contact interaction [10], with a scale dependent coupling constantḡ(Λ). It is related to the physical s-wave scattering length a viā where the high-momentum cutoff Λ takes the effects of a finite effective range into account. In the experimentally relevant case of open-channel dominated Feshbach resonances the zero-range limit Λ → ∞ is realized, which leads to characteristic power laws in the high-momentum asymptotics of correlation functions. For instance, the momentum distribution satisfies n σ (k → ∞) = Ck −4 where the Tan contact density C accounts for the physics at short distances [45]. The unitary limit 1/a = 0 can be reached by tuning the interaction directly to the Feshbach resonance, which is controlled by an external magnetic field. As a result, there is no small interaction parameter available and a nonperturbative treatment is mandatory to obtain quantitative results. The Luttinger-Ward approach results in single-particle Green's functions G σ at finite temperature with selfconsistently resummed interaction effects and is in good agreement with thermodynamic measurements in the strong-coupling regime around the unitary limit [35,46]. In addition to the fermionic Green's function, the Luttinger-Ward theory also allows to determine the pair propagator where T τ denotes time ordering in imaginary time τ . At the superfluid transition temperature T c , the pair propagator Γ(Q = 0, Ω n = 0) diverges according to the Thouless criterion. Furthermore, the Tan contact is obtained from the short-distance limit [47] 2 C where the anomalous contribution from the superfluid order parameter ∆ vanishes in the normal phase considered here. In the following, both G σ and Γ form the basis for the evaluation of the thermal sum rule.

B. Linear response
In order to determine the thermal conductivity of the UFG we first evaluate the thermal sum rule χ T qq . In contrast to other transport coefficients such as the viscosity [16] or the spin diffusivity [20], χ T qq cannot be directly attributed to standard thermodynamic quantities but requires an additional thermal operator [48].
In general, the heat current response χ qq (ω) is obtained within linear response by adding the perturbation δĤ(t) =´x q (x, t) · h(x, t) to the Hamiltonian. Rather than working in real time, the problem is more conveniently treated in imaginary time τ ∈ [0, β). Furthermore, we consider only homogeneous source terms h(x, τ ) = h(τ ) since we are interested in the q = 0 response. We express the grand canonical partition function in the presence of the external field h as a coherent state path integral with fermionic action S F , Then log Z is a generating functional for connected heat current correlations From the latter function the retarded response in real frequency is obtained by Fourier transforming τ to the bosonic Matsubara frequency ω n and subsequent analytic continuation iω n → ω + i0 + .
For the Hamiltonian (10) the particle and energy current operators read [29,49] The bare energy current operator E has a kinetic and an interaction contribution. Considering the corresponding operators in momentum space, shows that the prefactor of the interaction part is only sensitive to the center-of-mass momentum Q of the pair of fermions participating in the interaction. Therefore, this term is most easily discussed in two-channel variables with a bosonic pair field ∆(x) =ḡψ ↓ (x)ψ ↑ (x). The latter can be easily introduced by decoupling the action (14b) in the pairing channel by a Hubbard-Stratonovich transformation. We notice that the presence of h(τ ) leads to the shift ε p → ε p + (ε p − µ σ − Ts)p/m · h(τ ) of the bare fermionic dispersion relation and to the rescalinḡ g(Λ) →ḡ(Λ)(1 + Q/m · h(τ )) in S F . With these substitutions we obtain the path integral Z[h] within the two-channel formulation in momentum space, From Eq. (15a), we thus find the expectation value of the heat current where we have defined the bare fermionic and bosonic heat current vertices and furthermore recovered the single-particle Green's function G σ as well as the pair propagator Γ defined above. A vanishing perturbation h(τ ) = 0 implies  q (q = 0, τ ) = 0 by rotation invariance. For nonzero perturbation, we obtain the susceptibility by taking a functional derivative of the current (19) according to Eq. (15b), Here, the dressed current vertices for spin component σ and the pairs are defined as the time-ordered expectation values while the last contribution to χ qq arises from the second derivative of the∆∆ prefactor in the action (18b). The thermal sum rule follows by Fourier transformation to the external bosonic Matsubara frequencies ω n = 0, which yields This form of the sum rule in terms of current vertex functions is analogous to the sum rules for momentum [16,50,51] and spin currents [52]. IntroducingT q σ,pair as the amputated counterparts of T q σ,pair allows one to represent the Kubo formula for χ qq and thus also the sum rule in a diagrammatic manner, as depicted in Fig. 3, except for the last line.
Quite crucially, the exact heat current vertexT q satisfies a Ward identity [44]. Extending the latter from the fermionic to the bosonic sector, it reads in momentum space (cf. Appendix A) at vanishing external arguments ω = 0, q → 0 relevant for the sum rule. The first line contains the fermionic part expressed via the single-particle Green's function G σ , while the second line denotes the bosonic contribution in terms of the pair propagator Γ. The bare vertices T q(0) σ andT q(0) pair in Eq. (20) are obtained simply by using the noninteracting Green function G −1 0,σ (p, ε) = ε − ε p + µ and the inverse bare coupling Γ −1 0 (Q, Ω) =ḡ(Λ) −1 inside the Ward identity (see Appendix A).
In order to obtain the sum rule as the static limit of the current response function given by the Kubo formula (Fig. 3) we insert the Ward identities into Eq. (23). As a result, we obtain the exact thermal sum rule expressed in terms of the Green's and vertex functions, with two contributions: a fermionic part (26) and a new interaction part arising from the bosonic pairs of the form Both terms can be evaluated by inserting the previously computed Luttinger-Ward results for G σ (p, iε n ) and Γ(Q, iΩ n ) [31][32][33] as functions of momentum p (Q) and Matsubara frequency iε n (iΩ n ). Note that the full fermionic and bosonic energy current vertices as defined by the Ward identity (24) provide an exact solution of the Luttinger-Ward transport equations formulated in terms of fermionic and bosonic transport vertices [16,51]. This proves that the Luttinger-Ward approach implements exact energy conservation, even when fermionic and bosonic Green functions are obtained within the self-consistent T-matrix approximation. This was indeed the goal of constructing a conserving approximation, which in our case furthermore satisfies the exact Tan relations [34]. However, as indicated by the bar, these terms still depend explicitly on the momentum cutoff Λ, which is manifest for the second term. Moreover, a finite value of Λ is necessary to render the momentum integrals in the fermionic part finite. Therefore, we first have to discuss how to extract the universal results for the sum rule before presenting the numerical results.

C. Short-distance asymptotics
Due to the contact interaction, several terms in the sum rule (26,27) diverge in the zero-range limit Λ → ∞. This is apparent from the pair contribution, which in the unitary limit 1/a = 0 is directly proportional to Λ. Indeed, one quite generally expects a cutoff dependence of the static sum rules for these quantum critical systems. This can be attributed to high-frequency tails of the dynamic transport coefficients such as κ(ω) defined in Eq. (5) above [6]. For instance, in case of the shear viscosity the full sum rule reads [16,50] Π xyΠxy ω=0 = p + 4 2 CΛ 15π 2 m , where the second term arises from the high-frequency tail η(ω → ∞) = 3/2 C/15π √ mω. The static transport coefficient p is given instead by the regularized form of the sum rule with the Λ-dependent terms subtracted. As a result for the shear viscosity of the unitary gas, one has η = pτ η in analogy to Eq. (1) in the thermal case.
These divergences arise from the asymptotic largemomentum behavior of the fermionic momentum distribution, The appearance of the two leading contributions in asymptotic power-law decay of the momentum distribution arises from the nonanalytic contributions proportional to |x| and |x| 3 , respectively, in the short-distance operator product expansion Rψ † ↑ψ † ↓ψ ↓ψ↑ (R) + · · · (30) of the one-particle density matrix [53]. The coefficients C and D 1 in Eq. (29) are defined through the expectation values of the contact operatorĈ(R) = −4 m 2ḡ2 (Λ)ψ † ↑ψ † ↓ψ ↓ψ↑ (R) and its second derivative ∇ 2 RĈ (R); note that in d = 3 the Fourier transform of |x| is −8π/p 4 while |x| 3 gives 96π/p 6 . The presence of a subleading term D 1 /p 6 in the momentum distribution of two-component Fermi gases has been discussed in detail by Werner and Castin [54]. In general, the coefficient D 1 also contains a contribution which involves the derivative of the energy with respect to the effective range of the interaction. In our model, no such contribution appears and the full expression for D 1 is given in terms of the firstorder time and second-order spatial derivative of the pair propagator, see Eqs. (B5) and (B7) in Appendix B.
Within the self-consistent T-matrix approximation to the Luttinger-Ward functional the powers of momentum are correctly reproduced, whereas the contact coefficients that characterize the short-distance correlations in the many-body medium are obtained approximately. The asymptotic behavior of the numerical data of the fermionic momentum distribution is consistent both with the OPE and Ref. [54] up to k −6 , but to our knowledge the k −7 contribution has not been discussed before. The latter arises from an anomalous contribution to the pair In the fermionic part (26), the leading divergence O(Λ 3 ), which could arise from the C/p 4 tail of the momentum distribution, cancels between the first and last term in the square brackets, hence there is no Λ 3 divergence. According to Eq. (29), this leaves terms of order O(Λ) and O(log(Λ/k)), wherek denotes the momentum scale beyond which the algebraic power laws of the terms in χ T qq dominate; in practice, one hask 10/λ T . The coefficients of these subleading divergences depend on C, D 1 , and D 2 (for the log term). Similarly, for the pair momentum distribution we find the asymptotic expansion n pair (Q → ∞) = 64π 2 nC/3Q 6 + · · · , see Eq. (B12) in App. B. This implies that the momentum sum in the interaction term (27) is finite, while the inverse bare coupling in the prefactor diverges as O(Λ). In the numerical evaluation, we subtract all divergent terms to obtain the regularized sum rule (9), as has been done for the shear viscosity [16]. At unitary, in particular, the interaction term does not contribute to the regularized sum rule as its contribution scales like 1/a. Away from unitarity 1/a = 0, in turn, the bosonic part gives rise to a new contact correlation contribution to the thermal conductivity similar to what has been found in the bulk viscosity [51,55,56]. Regarding the dynamic thermal conductivity, the O(Λ) contribution implies a tail κ(ω → ∞) ∼ ω −1/2 , while O(ln Λ) causes a subleading contribution to the high-frequency behavior proportional to ω −1 , in analogy to the discussion below Eq. (28).

D. Numerical Results
After subtracting from Eq. (25) all terms that diverge in the zero-range limit we find the exact result for the thermal conductivity sum rule at unitarity where the pair contribution (27) vanishes, This defines the dimensionless quantum critical scaling function f χ T qq (βµ). In the high-temperature regime the sum rule is given analytically by the result − 10(βµ +s (0) ) Li 5/2 (−e βµ ) (32) of a free Fermi gas. The entropy per particle reads s (0) = 5 Li 5/2 (−e βµ )/[2 Li 3/2 (−e βµ )] − βµ is expressed in terms of the polylogarithm Li s (z). In terms of density, the result approaches mχ T qq → (5/2)nT 2 for βµ → −∞. Using the Luttinger-Ward thermodynamic data allows us to extend the sum rule from the high-temperature regime into the quantum degenerate regime down to the critical temperature of the superfluid transition at (βµ) 2.5 [35]. The result is shown in Fig. 4: while it agrees with χ T (0) qq in the virial limit, it shows large deviations in the quantum degenerate regime from strong pairing fluctuations, which lead to an enhancement of up to two orders of magnitude close to the superfluid transition.

IV. QUANTUM CRITICAL SCATTERING TIMES
Within kinetic theory, the thermal scattering time τ κ is obtained as the collision time in the Boltzmann equa- tion in response to temperature gradients [57,58]. The collision integral I[f ] is evaluated for a generic distribution function f p , which deviates from the thermal equilibrium distribution f 0 p as f p = f 0 p + δf p . For small variations δf p , the collision integral can be linearized as where the linearized collision operator H[f 0 p ] acts on δf p but itself only depends on the equilibrium distribution f 0 p . The solution δf p = f 0 p (1 − f 0 p )U p of the Boltzmann equation minimizes the scattering rate, hence the particles choose a distribution U p to best avoid scattering. Within a family of trial functions U p , an upper bound to the true scattering rate is found in variational kinetic theory as [57,59] The scalar products (A, B) =´dΓ p f 0 p (1 − f 0 p )A p B p are defined with respect to the equilibrium distribution function f 0 p . The system is driven out of equilibrium by the perturbation X p : it determines which transport channel is considered, e.g., X p = pz m (ε p − (Ts + µ)) for thermal and X p = pxpy m for shear transport. The variational functions U p are arbitrary functions of momentum that have the same angular dependence as the perturbation X p .
The linear collision operator H[f 0 ] for 2 → 2 scattering between fermionic (quasi)particles (p 1 , p 2 → p 1 , p 2 ) reads where momentum conservation p 1 + p 2 = p 1 + p 2 and energy conservation ε p1 + ε p2 = ε p 1 + ε p 2 are satisfied in elastic scattering, and Ω denotes the angle between the incoming and outgoing scattering planes. The scattering cross section is given as dσ/dΩ = |f | 2 in terms of the swave scattering amplitudef ; in the strongly interacting Fermi gas, the medium scattering amplitude reads [33,34] to leading order in the systematic large-N expansion. While the first two terms reproduce the s-wave scattering amplitude at the two-particle level, the integral takes corrections caused by the presence of a finite density medium into account. At unitarity 1/a = 0 the constant offset vanishes, and the dimensionless scattering amplitudẽ f /λ T depends on βµ alone. The properties at high temperature are obtained to leading order in the virial expansion in small fugacity z = e βµ 1, with f 0 p = e −β(εp−µ) the Boltzmann distribution. The resulting scattering times are and [59] already from the first variational basis function U p ∝ X p , and corrections from higher basis functions are less than 1.5% for the shear viscosity [60]. Note that the high-temperature results at unitarity already satisfy the quantum critical scaling form τ x T / = f x (βµ) with x = κ, η.
In the quantum degenerate regime, one instead has to use the Fermi-Dirac distribution f 0 p = [e β(εp−µ) + 1] −1 . Two competing effects thus modify the scattering times τ x : Pauli blocking in Eq. (34) reduces the phase space for scattering and strongly increases the scattering time, while medium scattering in Eq. (35) has the opposite effect and reduces the scattering time. In case of the shear viscosity the scattering cross section dσ/dΩ even diverges at the superfluid transition due to gapless pairing fluctuations if only a single variational basis function U p ∝ X p is considered. This would lead to the unphysical result η → 0 at T c , which arises from the divergence of the T-matrix Γ ∼ Q −2 for small energies. However, an improved variational solution in a larger basis set yields a finite result, as is expected for the viscosity near the superfluid transition [38].
The full results are shown in Fig. 5: the surprising and remarkable observation is that the scattering time f x (βµ) is nearly the same for the Boltzmann distribution ("Boltzmann") and for the Fermi-Dirac distribution ("large-N medium"), not only for viscous [61] but also for thermal transport. Changing the distribution f 0 p from Boltzmann to Fermi-Dirac modifies the calculation in three places: (i) in the scalar product in the variational expression (33), (ii) in the occupation numbers of the collision integral (34), and finally (iii) in the medium scattering amplitude (35). In the quantum critical regime, the subtle interplay between these effects leads to an almost perfect cancellation between the Pauli blocking and medium scattering corrections in the large-N medium result. We find a similar coincidence also for spin diffusion (see App. C). Hence, there appears to be a general mechanism at work that does not depend on the angular, spin or energy weight of the driving term X p .
What has not been appreciated before is that, even more remarkably, also the strong-coupling Luttinger-Ward computations [16] (red dots) confirm this result for the scattering time as a function of βµ for the whole quantum critical regime βµ 1 (T 2T c ) within a 15% error bound, where the scattering time has been extracted from the relation η = pτ η in analogy to Eq. (1). We thus conjecture that the large-N expansion is similarly accurate for the thermal scattering time τ κ in the quantum critical regime, and we use the large-N result (36) henceforth. Closer to the phase transition, however, the quantum critical scaling crosses over into the classical critical scaling of the 3D XY universality class near the superfluid phase transition (see Sec. II A above).
At unitarity, the scattering times thus satisfy the quantum critical scaling form [5,34] τ x = f x (βµ)( /T ), where the dimensionless scaling function f x (βµ) depends only on the value of βµ, not only in the quantum critical regime but also in the high-temperature nondegenerate gas; in the quantum degenerate region βµ ≥ 0, the scaling function attains values of order unity (Fig. 5). For spin transport (App. C), this is consistent with the experimental observation of quantum critical spin drag [19] and Planckian dissipation for spin [21,22,62]. , red above Tc [33]) and from experiment [35]. The dashed line indicates the high-temperature limit LP = 2/3.

V. RESULTS AND QUANTUM CRITICAL TRANSPORT RATIOS
Based on the hydrodynamic arguments from above, we arrive at the first prediction for the thermal conductivity (1), κT = χ T qq τ κ , in the quantum critical regime, as shown in Fig. 6. Here, χ T qq is evaluated within Luttinger-Ward theory (Fig. 4) and combined with the thermal scattering time in the Boltzmann limit (36). In the limit βµ → −∞ one finds the Boltzmann value for the thermal conductivity κ B = 225/(64 √ 2) T /( λ T ) [30] by using the noninteracting sum rule (32).
Weighting the thermal diffusion D T with the thermodynamic Landau-Placzek ratio LP = c p /c V − 1 (Fig. 7) yields the thermal contribution LP×D T to the sound dif- fusion D sound shown in Fig. 1. At low temperatures above T c the decrease of LP seems to suggest that heat diffusion becomes less important for sound attenuation near T c , but this is more than compensated by the increase of D T which makes heat diffusion rather more important. A fluid is characterized by the relative importance of different transport channels, which is quantified by transport ratios. Here, we consider the Prandtl number, which is defined as the ratio of shear and thermal diffusivities (D T is reported in Fig. 1), As the last term shows, the transport ratio is a product of a thermodynamic term that incorporates nontrivial temperature scaling from the full equation of state, and a ratio of transport times which we have found to remain nearly constant at τ η /τ κ = 2/3 throughout the quantum critical regime. Therefore, in the unitary gas the transport ratios derive their temperature dependence predominantly from the equation of state, and we use the best available Luttinger-Ward equation of state [31,33] to obtain the theory prediction for the Prandtl number in Fig. 8. Note that Pr starts from a value of 2/3 in the high-temperature limit and then grows to about 0.7 near T ≈ T F , before it falls to much smaller values below 0.2 near the superfluid transition. This nonmonotonic behavior results from the Landau-Placzek ratio [26,33] shown in Fig. 7 for the unitary Fermi gas, and is consistent with the virial expansion [30]. At the classical superfluid phase transition (Model F) [38] one expects that η remains finite while κ diverges according to Eq. (2), suggesting a vanishing Pr → 0. This nonmonotonic behavior is very well confirmed by a recent measurement of sound attenuation in the unitary gas [23]. The value of the Prandtl number also has an important interpretation in terms of possible nonrelativistic gravity duals, which, however, predict Pr = 1 [63] and can therefore be excluded as a model for the unitary Fermi gas. Another important transport ratio is the bulk-to-shear viscosity ratio ζ/η computed in [51], which shows that viscous transport occurs via quasiparticles only for T T F but deviates in the quantum degenerate regime. The Schmidt number Sc = D η /D s comparing shear with spin transport is shown in Fig. 10, see Appendix C.

VI. DISCUSSION
In conclusion, we have found that transport scattering times τ κ and τ η in the quantum critical regime follow a remarkably simple scaling law, which extends to the vicinity of the superfluid transition where pairing fluctuations become dominant. We have chosen specifically the unitary gas where the quantum critical regime extends to high temperature to demonstrate this point. This information is combined with a new exact sum rule for thermal transport, which depends on the equation of state and thermal operators beyond, to predict the thermal conductivity κ in the quantum degenerate regime. For κ and the Prandtl number Pr we find good agreement with recent experiments [23].
The remarkable coincidence of the quantum critical scattering times from the high-temperature Boltzmann calculation and the strong-coupling large-N and Luttinger-Ward results is a unique feature of the quantum critical point at infinite scattering length 1/a = 0: the scattering times must follow the quantum critical scaling form, which in the particular case of the unitary Fermi gas must extend up to high temperature by dimensional analysis, in contrast to lattice models. At high temperature, the scattering times are reliably obtained from kinetic theory as τ T / ∝ z −1 proportional to the inverse fugacity. Now quantum critical scaling predicts that this form continues from the dilute gas throughout the QCR until near T c . This remarkable observation is supported by the fact that it leads to good agreement with recent experimental data in the regime where quantum critical scaling can be applied. It will be interesting to see if our approach can be extended to other types of QCPs.
While at unitarity the bosonic part of the exact sum rule (25) provides only a regularization, away from unitarity (1/a = 0) it gives a new regular contribution that arises from local pair fluctuations, the so-called contact correlations [51]. This new contribution to thermal transport is not captured by fermionic kinetic theory and is particularly important at low temperatures near the superfluid phase transition.
We find that not only shear and spin diffusion, but also the thermal diffusion D T in units of /m exhibit quantum limited diffusion near T c . For thermal transport, the diffusion minimum D T 4.2 /m occurs well in the quantum critical region at T 0.7 T F (see Fig. 1). Hence, the quantum degenerate unitary Fermi gas is a nearly perfect fluid not only regarding momentum transport but also for thermal transport.
With current sound propagation measurements in box traps reaching into the superfluid regime [23,64], it will be particularly interesting to study critical scaling of the transport coefficients and observe the increase of D T shown in Fig. 1. This, as well as the related monotonic decrease of the Prandtl number indicated in Fig. 8, is due to the growing thermal conductivity associated with the crossover to classical critical fluctuations as expressed asymptotically in Eq. (2). For the sound diffusion D sound (T ) (Fig. 1) both our quantum critical prediction and the experimental data indicate a monotonic decrease, while at even lower temperatures T /T F 0.2 an increase of D sound is again theoretically expected from critical fluctuations. In the strongly interacting 2D Fermi gas the recently observed quantum scale anomaly [65] will have a large effect on the transport coefficients.  The spin scattering time τ s shown in Fig. 9 also exhibits the quantum critical scaling that we observed already for shear and thermal transport: the medium scattering time is, within our numerical resolution, identical to the Boltzmann scattering time τ s T / = 3π 8 √ 2 e −βµ . The quantum critical scattering time is now combined with the Luttinger-Ward equation of state for density n and spin susceptibility χ s to obtain the spin diffusivity D s . Now, the Schmidt number [57] Sc = D η D s = (p/mn)τ η (n/mχ s )τ s = pχ s n 2 × τ η τ s (C4) is defined as the dimensionless transport ratio of shear and spin diffusion and characterizes the relative importance of momentum and spin relaxation. As shown in Fig. 10, the Schmidt number starts from a value of Sc = τ η /τ s = 5/2 in the high-temperature limit and drops to around 0.3 near T c , indicating that momentum diffusion is suppressed by a factor of almost 10 relative to spin diffusion. This is physically expected because viscosity is carried both by single fermions and pairs and therefore strongly affected by pair fluctuations near T c , whereas pairs carry no spin current.