High-temperature expansion of the viscosity in interacting quantum gases

# High-temperature expansion of the viscosity in interacting quantum gases

Johannes Hofmann Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Cambridge CB3 0WA, United Kingdom TCM Group, Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom
July 1, 2019
###### Abstract

We compute the frequency-dependent shear and bulk viscosity spectral functions of an interacting Fermi gas in a quantum virial expansion up to second quadratic order in the fugacity parameter , which is small at high temperatures. Calculations are carried out using a diagrammatic finite-temperature field-theory framework, in which the analytic continuation from Matsubara to real frequencies is carried out in closed analytic form. Besides the zero-frequency Drude peak, our results for the spectral functions show a broad continuous spectrum at all frequencies with an additional bound state contribution for frequencies larger than the dimer-breaking energy. Our results are consistent with various sum rules and universal high-frequency tails. In the low-frequency limit, the shear viscosity spectral function is recast as a collision integral, which reproduces known results for the static shear viscosity from kinetic theory. Our findings for the static bulk viscosity of a Fermi gas near unitarity, however, show a non-analytic dependence on the scattering length, at variance with kinetic theory.

## I Introduction

Strongly interacting quantum gases form a paradigmatic quantum many-body system in which the interaction strength and trap parameters are tuned at will Bloch et al. (2008); Giorgini et al. (2008). Properties of cold gases are probed accurately by studying their collective oscillations and their expansion dynamics. The damping or dissipation in these processes is set by the shear viscosity and bulk viscosity  Schäfer and Chafin (2012); Schäfer (2014), the experimental determination of which has been an active topic over the past years Cao et al. (2011a, b); Vogt et al. (2012); Elliott et al. (2014); Joseph et al. (2015a, b). Indeed, strongly interacting Fermi gases are said to be “perfect fluids” in which both bulk and shear viscosity are anomalously small Chafin and Schäfer (2013); Schäfer (2014). The shear viscosity comes close to a conjectured lower limit on the shear viscosity to entropy ratio  Kovtun et al. (2005), which should apply to strongly interacting quantum systems close to a scale-invariant point. Likewise, the bulk viscosity will vanish in a scale-invariant system Son (2007).

Formally, the shear and bulk viscosities are the zero-frequency limit of corresponding viscosity spectral functions and , which are defined in terms of Kubo relations for the retarded stress tensor correlation function:

 η(ω) =ImχRxy,xy(ω)ω (1) ζ(ω) =ImχRii,jj(ω)d2ω, (2)

where is the dimension of the gas. The retarded response function is the Fourier transform of the real-time response with the thermal average and the stress tensor. Quite generally, the extrapolation of numerical results for the response functions to zero frequency is an ill-posed problem that requires numerical extrapolation methods. A first-principles calculation of the shear and bulk viscosity from Kubo relations is thus intricate Enss et al. (2011); Wlazłowski et al. (2012, 2013, 2015).

Exact constraints on the spectral functions (1) and (2) exist in the form of sum rules and universal high-frequency tails Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012). At large frequencies, the spectral functions have universal power-law high-frequency tails that decay as  Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012), with a magnitude set by the contact parameter  Tan (2008a, b, c). Moreover, the total integrated spectral weight depends on the derivative of the contact with respect to the scattering length Taylor and Randeria (2010, 2012). These exact results should be obeyed by any calculation of the spectral functions.

A often-used framework to compute transport coefficients in quantum gases is kinetic theory Bruun and Smith (2005); Massignan et al. (2005); Bruun and Smith (2007); Schäfer (2012); Enss et al. (2012); Bruun (2012); Dusling and Schäfer (2013). However, strictly speaking, kinetic theory only applies in two limiting cases: the first is the low-temperature limit of a strongly spin-imbalanced Fermi gas, which is described by Fermi liquid theory Pines and Nozières (1966); Abrikosov et al. (1975); Recati and Stringari (2012). Here, a kinetic description exists for long-lived quasiparticle excitations derived from an underlying density functional Recati and Stringari (2010); Hofmann et al. (2018). The second one is the high-temperature or non-degenerate limit Yan et al. (2019), in which the thermal wavelength is much smaller than the interparticle distance, , where is the density, inverse temperature, and the mass of an atom. This is the limit we study in this paper. A drawback of kinetic theory is that the high-frequency tails of the viscosity spectral functions are not correctly captured and the extrapolation to low temperatures is not controlled. While a formal correspondence between a kinetic description and the classical limit of the microscopic theory can be established within the Keldysh formalism Kamenev (2013), a direct link between the results for transport coefficients obtained within kinetic theory and microscopic calculations is not apparent.

In this paper, we consider the high-temperature limit of the viscosity spectral functions in a microscopic calculation using a quantum virial expansion. The virial expansion applies in the grand canonical ensemble and expands the thermal expectation value of an operator as a sum over expectation values restricted to the -particle sector:

 ⟨^A⟩ =Tr[e−β(^H−μN)^A]=∑NzNTrN[e−β^H^A], (3)

where the expansion parameter is called the fugacity. In particular, for the particle density, we obtain

 nλdT =z+2b2z2+3b3z3+…, (4)

with the virial coefficients Liu (2013). Truncating the expansion (3) or (4) after the first few terms holds for , which corresponds to the non-degenerate or high-temperature regime . Since the -th order term of Eq. (3) involves only matrix elements of -particle states, the virial expansion provides a link between the high-temperature properties of the quantum gas and few-particle solutions, which are often known. Note that the quantum virial expansion is valid for any interaction strength, and in this sense it is a non-perturbative method. In order to make contact with many-body theory, in this paper we use a diagrammatic method to compute the viscosity spectral functions.

The main result of this paper are the exact viscosity spectral functions and up to second order in the fugacity. An advantage of our calculation is that the analytical continuation to real frequencies is performed exactly and does not require extrapolation schemes. For illustration, Fig. 1 shows the general form of the spectral functions obtained in this paper. There is a broad continuous spectrum at all frequencies, which arises from interactions between scattering atoms. If two-particle bound states are present, there is additional weight from bound-free transitions, which contribute at frequencies larger than the bound-state energy with a non-analytic frequency-dependence right above threshold. In particular, for the bulk viscosity, the bound state part has a very steep onset and dominates the spectral function. At small nonzero frequencies, only the continuum part remains. For the shear viscosity, this continuum part diverges as as a power-law , whereas for the bulk viscosity, it saturates to a finite value. There is an additional delta-peak contribution at zero frequency for both shear and bulk viscosity. For the bulk viscosity, this peak arises from bound-bound transitions and has nonzero weight only if a bound state is present. Both the continuum and bound state parts contribute at high frequencies.

It turns out that the interacting contribution to the small-frequency divergence of the shear viscosity spectral function takes the form of a collision integral, which allows us to obtain a static shear viscosity using a memory functions resummation. Our results agree with calculations of the shear viscosity within kinetic theory. Interestingly, in the absence of bound states in the two-particle sector, the static limit of the bulk viscosity spectral function can be taken directly. We quote a result of our calculation near the unitary limit

 lima→−∞λ3Tℏz2ζ (5)

with the scattering length and the Euler-Mascheroni constant. This result differs from kinetic theory, which predicts a quadratic dependence  Dusling and Schäfer (2013).

This paper is structured as follows: At the beginning of Sec. II, we introduce the zero-range model and discuss the Kubo formula that define the spectral functions. The remainder of Sec. II then presents the main calculation of this paper: In Sec. II.1, we set up a diagrammatic framework to compute the high-temperature expansion of correlation functions exactly to second in the fugacity . The diagrammatic calculation for the stress correlator and the contact correlator used to compute the viscosity spectral functions is presented in Secs. II.2 and II.3, with details of the analytical continuation from Matsubara space to real frequencies relegated to App. A. The results for the shear and bulk viscosity spectral functions are presented in Sec. III. We discuss the general form of the shear viscosity spectral function in Sec. III.1, with particular attention to the universal high-frequency behavior. We show how a memory function approach reproduces the kinetic results in the static limit. The bulk viscosity spectral function is discussed in Sec. III.2. Our results are consistent with universal high-frequency tails and sum rule constraints. In addition, we provide analytic results for the static bulk viscosity and comment on the discrepancy with kinetic theory. The paper is concluded by a summary and outlook in Sec. IV.

## Ii High-temperature expansion of the viscosity spectral function

We consider a Fermi quantum gas with zero-range interactions described by the Hamiltonian

 ^H =∫ddx{∑σ^ψ†σ−ℏ2∇22m^ψσ−g^ψ†↑^ψ†↓^ψ↓^ψ↑}, (6)

where and is the bare interaction strength. The stress tensor in the zero-range model (6) takes the form (the summation over is implied)

 ^Πij =ℏ22m[(∇i^ψ†σ)(∇j^ψσ)+(∇j^ψ†σ)(∇i^ψσ)] (7)

The off-diagonal terms of the stress tensor as used in the shear viscosity spectral function (1) only involve bilinear operators. For direct calculations of the bulk viscosity, however, the expression (2) is inconvenient as the stress tensor trace involves both bilinear and quartic operators:

 ^Πii =ℏ2m[(∇i^ψ†σ)(∇i^ψσ)−d4∇2(^ψ†σ^ψσ)]−dg^ψ†↑^ψ†↓^ψ↓^ψ↑. (8)

It was established in Refs. Martinez and Schäfer (2017); Czajka and Jeon (2017) that the bulk viscosity can also be defined in terms of the response function

 ζ(ω) =ImχRO,O(ω)ω (9)

with , where is the energy density with . The operator is linked to the Tan contact operator and is a measure of the breaking of scale invariance in the quantum gas Hofmann (2012):

 ^O =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩aℏ2^C2m(1D)ℏ2^C4πm(2D)ℏ2^C12πma(3D), (10)

where is the effective scattering length in dimensions and is the contact operator Tan (2008a, b, c), in the zero-range model Braaten and Platter (2008); Braaten et al. (2008). The contact parameterizes universal short-range correlations and thermodynamic properties and is thus a central experimental observable Sagi et al. (2012); Hoinka et al. (2013); Carcy et al. (2019); Mukherjee et al. (2019). The expectation value of the contact sets the magnitude of the leading short-distance divergence of the pair correlations function, hence intuitively it describes the number of pairs with opposite spin in close proximity Braaten (2012). The contact is linked to the partition function of the gas by the adiabatic relation Tan (2008c); Barth and Zwerger (2011); Valiente et al. (2011); Braaten (2012)

 C =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩2mℏ2∂(Ω/V)∂a(1D)2πmℏ2a∂(Ω/V)∂a(2D)4πmℏ2a2∂(Ω/V)∂a(3D). (11)

Indeed, the contact also governs the non-analytic short-time correlation of the stress correlator , and thus sets the magnitude of a power-law high-frequency tail of the viscosity spectral functions Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012), which will be discussed later in this paper.

The expectation value of sets the deviation of the pressure from the scale-invariant value: . In two dimensions, the operator indicates a quantum scale anomaly, where formally, it describes an anomalous contribution that opens up the non-relativistic conformal operator algebra of a scale-invariant quantum gas Hofmann (2012). Indeed, the bulk viscosity spectral function vanishes identically in a scale-invariant system and is nonzero only when scale-invariance is broken Son (2007). On a technical level, the contact response does not contain bilinear operators, which simplifies calculations of the bulk viscosity considerably. Note that the viscosity spectral functions can also be defined in terms of the transverse and longitudinal current response functions Forster (1990); Taylor and Randeria (2010); Hofmann (2011). However, this definition requires the computation of correlation functions at finite external frequency and momentum, which complicates the calculation in the present case.

### ii.1 Diagrammatic virial expansion

Within finite-temperature field theory, we compute the discrete Fourier transform of the time-ordered correlation function

 χA,B(iωn) =∫β0dτeiℏωnτ⟨TτA(0)B(τ)⟩, (12)

where is a bosonic Matsubara frequency and denotes time-ordering with operators at later times placed to the left. Retarded response functions as needed for the spectral functions Eqs. (1) and (9) are obtained by analytically continuing from Matsubara to real frequencies  Mahan (2000). When working in Matsubara space, the fugacity parameter is contained in Bose and Fermi factors used to convert Matsubara frequency summations to contour integrals. An expansion in is not apparent in this framework as it changes the analytic structure of the integrand and hence does not commute with the frequency integration. Little is gained by computing the full result first and then expanding in the fugacity.

It turns out that a fugacity expansion can be developed systematically when working in imaginary time as opposed to Matsubara space. There, the fugacity only enters through the imaginary-time single-particle Greens function

 G(τ,k) =−e−τ(εk−μ)[Θ(τ)−nF(εk−μ)], (13)

where is the wave vector, the single-particle energy, and the Fermi factor, which can be expanded directly as . This expansion commutes with any internal integration over imaginary time or momenta. Diagrammatic calculations of the equation of state date back to the early days of quantum many-body theory with works by Montroll and Ward Montroll and Ward (1958) and Vedenov and Larkin Vedenov and Larkin (1958) for the screened plasma. For quantum gases with short-range interactions, a diagrammatic framework was put forward in Refs. Kaplan and Sun (2011); Leyronas (2011), which compute the third virial coefficient of a strongly interacting Fermi gas. These techniques are not restricted to thermodynamic quantities but can also be applied to correlation functions as well. Indeed, besides virial coefficients Kaplan and Sun (2011); Leyronas (2011); Ngampruetikorn et al. (2015); Barth and Hofmann (2015), diagrammatic calculations have been performed for the spectral function Barth and Hofmann (2014); Ngampruetikorn et al. (2013); Sun and Leyronas (2015); Sun et al. (2017), the momentum distribution of Bose gases including three-body Efimov effects Barth and Hofmann (2015), and also the electron gas Hofmann et al. (2013). Most of these works employ a framework similar to the one introduced in Ref. Leyronas (2011). In the remainder of this section, we collect the Feynman rules as used in this paper.

The Feynman rules for the model (6) in imaginary time are as follows: Imaginary time runs from the left to the right of the Feynman diagram in an interval . Every vertex is assigned an imaginary-time index and contributes a factor . In addition, we denote the bare off-diagonal stress tensor insertion (7) on a propagator at a specific imaginary time by

 T0xy(k)=ℏ2mkxky. (14)

Likewise, the contact operator insertion carries a factor of . As shown in Fig. 2, we indicate the term of the single-particle Green’s function (13) by a line that is slashed times. It contributes factors of

 G(m)(τ,k) (15)

Contributions of order to some correlation function are represented by Feynman diagrams that contain slashes. A key observation is that the leading order in of the propagator is purely retarded. Leading-order diagrams thus contain a minimum number of backward-propagating (advanced) propagators. We impose momentum conservation and integrate over internal momenta and imaginary-time labels. Imaginary-time integrals take the form of convolution integrals which are simplified using the convolution theorem for the Laplace transform

 ∫∞0d(t1,…,tn) ×[Θ(t1)f1(t1)]…[Θ(tn)fn(tn)]Θ(τ−t1−…−tn) =∫BWds2πie−τsf1(s)…fn(s), (16)

where we define the Laplace transform with inverse . BW is the Bromwich contour which runs from negative to positive infinity with all singularities of the integrand to the right of the contour.

Typical Feynman diagrams constructed as described above involve repeated few-body scattering processes as sub-blocks. It turns out that to leading order in the fugacity , the sum of these processes gives the free-particle scattering matrices. These blocks link the high-temperature properties of an interacting many-body system to known few-body solutions. We will use the two-particle scattering matrix , which we represent by a box as shown in Fig. 2. Using the convolution theorem and performing the Laplace transform, is given by a geometric series Barth (2015)

 T−12(s,P) =−1g−LP(s), (17)

where is the center-of-mass momentum and we define the bare two-particle propagator

 LP(τ) =Θ(τ)∫ddl(2π)de−τ(εl−εP−l) (18) LP(s) =∫ddl(2π)d−1s−εl−εP−l. (19)

The momentum integration diverges in , which is canceled by the renormalization of the coupling constant . For reference, we note the result for the -matrix

 T−12(s) =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩m2ℏ2(−a+√−ℏ2ms)(1D)m2πℏ2lna√−mℏ2s(2D)m4πℏ2(−a−1+√−mℏ2s)(3D), (20)

where, by Galilean invariance, . It will be convenient at a later point to use dimensionless scattering matrices , where we separate a factor , , and in , and , respectively.

For further reference, we note a useful representation of the contact to leading order in the fugacity:

 λ4TC =4π2z2β22d/2λdT∫∞−∞dse−βsπImT2(s+i0). (21)

The integral can be evaluated in closed analytical form in 1D Hoffman et al. (2015) and 3D. This expression will be used in Sec. III when extracting the high-frequency tails of the viscosity spectral function, the magnitude of which is set by the contact.

### ii.2 Stress correlator

In this section, we discuss the Feynman diagram expansion of the stress correlator used in the computation of the shear viscosity (1) up to quadratic order in the fugacity . A Feynman diagram representation of the imaginary-time correlator is obtained by drawing all connected diagrams with a stress tensor insertion at and at . All such diagrams up to quadratic order are shown in Fig. 2. We then take the discrete Fourier transform to obtain the response as a function of Matsubara frequencies. The leading order contribution, Fig. 2(a), is the same as for non-interacting gas. Indeed, the non-interacting contribution can be summed to all orders in . Using the Feynman rules introduced in the previous section, we obtain the standard kinetic result

 ImχR,NIxy,xy(ω)=δ(βℏω)(1−e−βℏω) ×2∫ddk(2π)d(−∂f(εk)∂εk)T0xy(k)T0xy(k), (22)

where the prefactor accounts for the two spin species. This term will contribute a delta function to the shear viscosity spectral function .

Two-particle interactions enter at second order . The Feynman diagrams that represent the interacting contribution to the shear viscosity spectral function (1) are shown in Fig. 2(b)-(c). These diagrams are commonly divided into three separate classes Larkin and Varlamov (2005): (a) Self-energy diagrams, which affect the spin-symmetric part of the response function and contain a self-energy insertion on one propagator; (b) Maki-Thompson diagrams, which contribute to the spin-antisymmetric response and contain a vertex correction; and (c) Aslamazov-Larkin diagrams. In the following, we compute each of these contributions separately.

#### ii.2.1 Self-energy

There are two self-energy diagrams as shown in Fig. 3(b). Using the imaginary-time Feynman rules and applying the convolution theorem for the Laplace transform gives, the first one evaluates to

 χSE,axy,xy(iωn)=2z2∫ddp(2π)d∫ddk(2π)dT0xy(k)T0xy(k) ×∫BWds2πie−βsT2(s+iℏωn,k+p)(s−εk−εp)(s+iℏωn−εk−εp)2. (23)

The overall prefactor accounts for the spin summation. The Bromwich integral and the analytic continuation from Matsubara to real frequency are performed as outlined in App. A. This gives the self-energy contribution to the retarded response:

 ImχR,SE,axy,xy(ω)1−e−βℏω=−2z2∫ddp(2π)d∫ddk(2π)de−β(εk+εp) ×Im[T2(ℏω+i0+12εk−p)(ℏω+i0)2]T0xy(k)T0xy(k). (24)

The second self-energy diagram in Fig. 3(b) contributes an identical term up to the sign with :

 ImχR,SE,bxy,xy(ω) =−ImχR,SE,axy,xy(−ω). (25)

The imaginary part of the self-energy contribution is antisymmetric in the frequency, as expected for a retarded response function.

#### ii.2.2 Maki-Thompson

The Maki-Thompson diagrams are shown in Fig. 3(c). The evaluation proceeds in a similar way as for the self-energy part. For the first diagram, we obtain:

 ImχR,MT,axy,xy(ω)1−e−βℏω=−2z2∫ddp(2π)d∫ddk(2π)de−β(εk+εp) ×Im[T2(ℏω+i0+12εk−p)(ℏω+i0)2]×T0ij(k)T0kl(p). (26)

Again, the second diagram in Fig. 3(c) contributes a term of equal magnitude with the replacement :

 ImχR,MT,bxy,xy(ω) =−ImχR,MT,axy,xy(−ω). (27)

#### ii.2.3 Aslamazov-Larkin

The computation of the Aslamazov-Larkin diagram is more involved compared to self-energy and Maki-Thompson contributions. Corresponding Feynman diagrams are shown in Fig. 3(d). They evaluate to

 χALxy,xy(iωn)=4z2∫ddp(2π)d∫ddk(2π)dT0xy(k) ×∫BWds2πie−βsT2(s,k+p)T2(s+iℏωn,k+p)(s−εk−εp)(s+iℏωn−εk−εp) ×∫ddl(2π)dT0xy(l)(s−εl−εk+p−l)(s+iℏωn−εl−εk+p−l). (28)

There are two distinct diagrams in Fig. 3(d) and two spin species, hence the overall prefactor . The integral in the last line of Eq. (28) corresponds to the internal three-propagator loop at the center of the diagram (sometimes called a triangle diagram). We evaluate the integral using Eq. (17):

 ∫ddl(2π)dT0xy(l)(s−εl−εk+p−l)(s+iℏωn−εl−εk+p−l) =T0xy(k+p)4(iℏωn)[T−12(s,k+p)−T−12(s+iℏωn,k+p)]. (29)

Note that Eq. (29) is a Ward identity which holds for the full -matrix. Furthermore, we simplify the product of the two -matrices using

 T2(s,k+p)T2(s+iℏωn,k+p) =T2(s+iℏωn,k+p)−T2(s,k+p)T−12(s,k+p)−T−12(s+iℏωn,k+p). (30)

The diagram now takes a simpler form:

 ImχR,AL,axy,xy(ω)1−e−βℏω=−z2∫ddp(2π)d∫ddk(2π)de−β(εk+εp) ×Im[T2(ℏω+i0+12εk−p)(ℏω+i0)2]T0xy(k)T0xy(k+p), (31)

and we split off a second term

 ImχR,AL,bxy,xy(ω)=−ImχR,AL,axy,xy(−ω). (32)

### ii.3 Contact correlator

We now turn to the contact correlator used in the computation of the bulk viscosity spectral function. The leading-order contribution to the correlation function is of second order in the fugacity . Corresponding diagrams are shown in Fig. 4. They involve the -matrix coefficients and the bare two-particle propagator (Eq. (18)) as subblocks. The calculation proceeds in a similar way as in the previous section. After analytic continuation, the sum of all four diagrams is

 χC,C(ω)=m4z2g4ℏ8∫ddP(2π)d∫∞−∞dsπe−βs ×Im[LP(s+ω+i0) ×(1+T2(s+ω+i0,P)LP(s+ω+i0))] ×Im[LP(s+i0)(1+T2(s+i0,P)LP(s+i0))], (33)

where is the Laplace transform of the bare two-particle propagator defined in Eq. (19). Equation (33) is simplified using the definition of the -matrix (17), which gives

 χC,C(ω)1−e−βℏω=m4z2ℏ82d/2λdT ×∫∞−∞dsπe−βsIm[T2(s+i0)]Im[T2(s+ω+i0)]. (34)

This expression is antisymmetric in as can be seen by shifting the integration variable . The integral (34) is convenient for a direct numerical evaluation.

## Iii Results

In this section, we combine the results of the previous section to compute, in turn, the shear viscosity and bulk viscosity spectral functions in Secs. III.1 and III.2, respectively. We present integral representations of the spectra, from which the general three-part structure of the spectral function as shown in Fig. 1 — a Drude or dimer zero-frequency peak, a broad free-particle continuum, and a bound-state contribution above threshold — becomes apparent. Numerical results for the spectral functions are presented for various dimensions and scattering lengths. Universal high-frequency tails and sum rules are obeyed in our calculations, which demonstrates the internal consistency of the virial expansion up to this order. We discuss the static limit in some detail and make connections with kinetic theory.

### iii.1 Shear viscosity

We begin by computing the shear viscosity spectral function. The zero-frequency Drude peak has a leading non-interacting contribution, Eq. (22). At finite frequency, the spectral function is of order and given by the sum of self-energy, Maki-Thompson, and Aslamazov-Larkin contribution:

 η(ω)=Im[χR,SExy,xy(ω)+χR,MTxy,xy(ω)+χR,ALxy,xy(ω)]ω =z24(1−e−βℏω)ℏ2ω3∫ddp(2π)d∫ddk(2π)de−β(εk+εp) ×Im[T2(ℏω+i0+12εk−p)]T0xy(k−p)T0xy(k−p) +(ω→−ω), (35)

where we have symmetrized the stress tensor matrix element . For a numerical evaluation, we transform to relative and center-of-mass momenta and , which restricts the angular integration to the matrix elements . The centre-of-mass momentum integration is then performed directly.

The imaginary part of the scattering -matrix , Eq. (20), has a continuous spectrum at positive energies and (if present) a state pole at the bound-state energy . From the structure of the integrand in Eq. (35), we deduce the general form of the shear viscosity spectral function shown in Fig. 1(a): First, the free-particle part of the -matrix contributes at all frequencies and forms the continuum indicated in the figure. Second, if the frequency exceeds the dimer binding energy, the second term in (35) gives rise to a bound state contribution with a threshold behavior , which is non-analytic in 3D. At small frequencies, only the continuum part contributes to the spectrum. As is apparent from the integral (35), the spectral function diverges as in this limit. This divergence will be discussed in more detail at the end of the section. Figure 5 shows numerical results for the shear viscosity spectral function in two and three dimensions for various values of the dimensionless inverse scattering length . The general structure of Fig. 1(a) is clearly reflected in these results.

#### iii.1.1 High-frequency tails

The shear viscosity spectral function has a power-law high-frequency tail with a magnitude set by the contact Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012). The leading-order term at large frequencies stems from the second term in Eq. (35) with . Transforming to a relative-energy integration variable , we shift and extend the -integration to run from to . Separating the leading power in from the integrand, we recognize the remaining integral as the leading term in the virial expansion of the contact, Eq. (21). This gives

 limω→∞η(ω) =⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ℏ2C8mω(2D)ℏ3/2C15π√mω(3D). (36)

These high-frequency tails agree with general results obtained using the operator product expansion (OPE) Hofmann (2011); Goldberger and Khandker (2012).

In the 3D unitary limit, we can compare with a Keldysh calculation Enss et al. (2011), which was performed using a finite lifetime in the single-particle self-energy. Here, and we state the separate contributions of the self-energy, Maki-Thompson, and Aslamazov-Larkin terms to the high-frequency tail

 limω→∞ηunitary(ω) =4m2z2πβ2ℏ41+1+030ℏ3/2π√mω, (37)

which agrees with Ref. Enss et al. (2011). This indicate that the set of diagrams considered in Ref. Enss et al. (2011) contains all the diagrams in Fig. 3 that yield the correct high-temperature limit.

#### iii.1.2 Static shear viscosity and memory function approach

Our results for the shear viscosity spectral function diverge at small frequencies, such that a direct extrapolation to the static limit is not possible. In the following, we discuss how this divergence is resummed using the memory function formalism Götze and Wölfle (1972); Forster (1990); Enss et al. (2011), which yields a finite static shear viscosity that reproduces the kinetic theory result.

Hydrodynamics dictates a low-frequency form of the retarded response Forster (1990)

 χ(ω)ω =Wω−i/τ=Wωτ21+(ωτ)2+Wiτ1+(ωτ)2, (38)

which has a finite limit for . The static shear viscosity is then expressed as the ratio of a Drude weight and a viscous scattering rate . The scattering rate, which vanishes for non-interacting systems, admits a separate fugacity expansion. A direct virial expansion of the correlation function thus applies in the limit and does not connect with the hydrodynamic expression (38). This is the origin of the divergence of the shear viscosity spectral function (35) for .

Progress can be made using the memory function formalism Götze and Wölfle (1972); Forster (1990); Enss et al. (2011), which instead of the full correlation function proposes to expand the memory function defined through

 χ(ω)ω =η0ω−M(ω), (39)

with a real constant. The memory function thus has the role of a self-energy of the viscous spectral function. Since , we have . Using the low-frequency expansion with and real, we find

 W =η01−M1,1τ=M01−M1, (40)

such that

 η=η0M0. (41)

By comparing the virial expansion results with the expansion

 χ(ω)ω =η0ω+η0ω2M(ω)+…, (42)

we extract the coefficients and and obtain a static limit of the shear viscosity. A similar approach as been used to determine the line shift in the dynamic structure factor of interacting quantum gases Hofmann and Zwerger (2017). To be definite, we consider the 3D case in the following. From Eq. (22), we obtain for :

 η0 =2∫d3k(2π)3(−∂f(εk)∂εk)[T0xy(k)]2 =(2m)3/2z2π3/2ℏ3β5/2+O(z2), (43)

valid to leading order in the fugacity. We make contact with kinetic theory by recasting the coefficient of the quadratic divergence of Eq. (35) at small frequencies as a (linearized) collision integral. Equation (35) is expressed as

 limω→0+η(ω) =η0M0ω2=η0τ−1ηω2, (44)

where following the literature we introduce a viscous relaxation time  Massignan et al. (2005); Bruun and Smith (2005, 2007)

 1τη =1N∫d3p(2π)3∫d3k(2π)3∫dΩdσdΩ|vrel|β2ℏ2m2q2xq2y ×f(εp)f(εk)(1−f(εp′))(1−f(εk′), (45)

with a normalization

 N =2β∫d3k(2π)3(−∂f(εk)∂εk)[T0xy(k)]2=βη0, (46)

which has dimension of inverse volume. Here, is the differential scattering cross section for the scattering of two particles with opposite spin and wave vector and to a state and , and is the relative velocity. The identity (44) is understood to hold at high temperatures, where Pauli blocking is not effective and the Fermi factors in Eq. (45) reduce to a Maxwell-Boltzmann distribution

 f(εp)f(εk)(1−f(εp′))(1−f(εk′) =z2e−β(εp+εk). (47)

Furthermore, the optical theorem relates the -wave scattering cross section to the imaginary part of the -matrix

 dσdΩ|vrel| =12πℏImT2(2εq). (48)

The expression (44) then follows directly by comparing with Eq. (35). From Eq. (41), we obtain the static shear viscosity

 η =2τη∫d3p(2π)3[T0xy(k)]2(−∂fσ∂εk). (49)

This is the result from kinetic theory Massignan et al. (2005); Bruun and Smith (2005, 2007). A discussion of this result is provided in the literature and will not be repeated here. Our calculations provide a rigorous correspondence between the high-temperature limit and kinetic theory for all scattering lengths.

A similar result was obtained for the unitary Fermi gas Enss et al. (2011), and we note an additional check of our results for the shear viscosity spectral function in the unitarity limit: writing the divergent contribution of the self-energy, Maki-Thompson, and Aslamazov-Larkin terms separately, we find

 limω→0ηunitary(ω) =z2(2m)3/215π5/2β7/2ℏ4ω2×43+3−302√2. (50)

Our result agrees with Ref. Enss et al. (2011) if we identify .

### iii.2 Bulk viscosity

This section presents results for the bulk viscosity spectral function . From the definitions (9), (10), and our result for the contact correlation function (34), we obtain the following integral representation of the spectral function:

 λdTℏζ(ω)=1−e−βℏωd2ℏω(λ2Ta2)d−2 ×∫∞−∞dsπe−βsIm[~T2(s+i0)]Im[~T2(s+ℏω+i0)], (51)

where the dimensionless scattering matrix is defined after Eq. (20). In the following, we discuss the structure of the spectral function and provide numerical results. We also discuss the high-frequency tail and sum rules, and obtain analytical results for the static limit.

From the structure of the integrand in Eq. (51), we directly infer the general form of the spectrum shown in Fig. 1(b): A first contribution arises from the integration over the free spectrum at positive energies. This part contributes at all frequencies and forms a continuous band. The integral in Eq. (51) is finite at , hence the spectral function approaches a constant value for small frequencies. A second contribution appears for frequencies larger than the bound state energy, for which we can pick up the pole of the first scattering matrix. Notably, it has a non-analytic frequency-dependence right above the threshold energy which follows from the non-analytic form of the -matrix (20). Finally, a third contribution at arises from the pinching of the bound-state poles in the integrand (51), which, on the two-particle level, describes the spectral function of a single dimer. This part gives rise to a delta-function contribution . The bound-state and delta contributions can be evaluated in closed analytical form.

Numerical results for the bulk viscosity spectral function as a function of the dimensionless inverse scattering length are shown in Fig. 6 for all dimensions and various values of the scattering length. The entire spectral function vanishes identically for . This point marks the unitary or non-interacting scale-invariant limit, where the bulk viscosity is expected to vanish on general grounds Son (2007). Away from this limit, we can clearly distinguish the continuum and the bound state part of the spectral function. As is apparent from the figure, the bound state contribution, when present, vastly dominates the continuum part. This observation is consistent with results for the zero-temperature spectral function Taylor and Randeria (2012).

#### iii.2.1 High-frequency tails and sum rules

The high-frequency tail of the bulk viscosity spectral function is determined by the non-analytic short-time behavior of the quantum gas and is expected to decay as a power-law Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012). Imposing a simultaneous scaling form, where the scattering length is assumed to scale with frequency (i.e, the dimensionless scaling variable is kept fixed), the dependence on the integration variable of the second term in the integrand of Eq. (51) can be neglected. The remaining integral is identified with the contact expectation value at this order in the virial expansion, Eq. (21). We obtain the high-frequency tails

 limω→∞ζ(ω) =⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ℏ5/2C2(mω)3/21(a−1√ℏ/mω)2+1(1D)ℏ2C4mω1π2+ln2(a−1√ℏ/mω)2(2D)ℏ3/2C36π√mω(a−1√ℏ/mω)2(a−1√ℏ/mω)2+1(3D). (52)

These results agree with the OPE result in 3D Hofmann (2011); Goldberger and Khandker (2012) and correct a result in 2D Hofmann (2011). To the best of our knowledge, the result in 1D is new. The analytic high-frequency tails match our numerical calculations of the full bulk viscosity spectral function, and provide a check of our results.

Unlike the shear viscosity spectral function, the bulk viscosity is finite for all nonzero frequencies, and hence so is the total spectral weight. An additional strong check of our result is provided by computing sum rules for the total spectral weight, which are related to the derivative of the contact Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012); Taylor and Randeria (2012)

 1π∫∞0dωζ(ω) =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ℏ216πma∂C∂a−1(2D)ℏ272πma2∂C∂a−1(3D). (53)

These sum rules are obeyed by our numerical results. Note that the dimer delta-peak (which contributes with half its weight at the integral boundary) is essential to saturate the sum rules.

#### iii.2.2 Static bulk viscosity

In three dimensions for negative scattering length, where no bound state exists, the zero-frequency limit of the bulk viscosity is finite. We state the exact result

 limω→0λ3Tℏz2ζ(ω) (54)

where is the exponential integral. In the perturbative limit of small negative scattering length, this becomes

 lima→−0λ3Tℏz2ζ =4√29(aλT)2.