High-temperature expansion of the viscosity in interacting quantum gases

TL;DR

This paper calculates the high-temperature shear and bulk viscosity spectral functions of an interacting Fermi gas using a quantum virial expansion and diagrammatic field theory, revealing broad spectra, bound states, and nonanalytic bulk viscosity behavior.

Contribution

It provides a detailed high-temperature expansion of viscosity spectral functions with analytic continuation, including bound-state effects and nonanalytic bulk viscosity dependence.

Findings

Spectral functions show a broad continuous spectrum at all frequencies.

Bound-state contributions appear above the dimer-breaking energy.

Bulk viscosity exhibits nonanalytic dependence on scattering length.

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 $z=e^{\beta \mu}$, which is small at high temperatures. Calculations are carried out using a diagrammatic finite-temperature field-theoretic framework, in which the analytic continuation from Matsubara to real frequencies is carried out in closed analytic form. Besides a possible 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 nonanalytic dependence on the scattering length, at variance with kinetic theory.

Figures

Equations

$$\eta (\omega )$$

$$\zeta (\omega )$$

$$⟨\hat{\mathcal{A}}⟩$$

$${\mathrm{Tr}}_{1,c}[e^{-\beta \hat{H}}\hat{\mathcal{A}}]$$

$${\mathrm{Tr}}_{2,c}[e^{-\beta \hat{H}}\hat{\mathcal{A}}]$$

$$n{\lambda }_T^d$$

$$\underset{a\to -\infty }{\lim }\frac{{\lambda }_T^3}{\hslash z^2}\zeta {|}_{\omega \to 0}$$

$$\hat{H}$$

$${\hat{\Pi }}_{ij}$$

$$\displaystyle\quad-\delta_{ij}\Bigl{[}\frac{\hbar^{2}}{4m}\nabla^{2}(\hat{\psi}_{\sigma}^{\dagger}\hat{\psi}_{\sigma})+g\hat{\psi}_{\uparrow}^{\dagger}\hat{\psi}_{\downarrow}^{\dagger}\hat{\psi}_{\downarrow}\hat{\psi}_{\uparrow}\Bigr{]}.$$

$${\hat{\Pi }}_{ii}$$

$$\zeta (\omega )$$

$$\hat{\mathcal{O}}$$

$$C$$

$$\eta (\omega )=\underset{q\to 0}{\lim }\frac{m{\omega }^2}{q^2}\mathrm{Im}{\chi }_T(\omega ,q),$$

$$\zeta (\omega )+\frac{2(d-1)}{d}\eta (\omega )=\underset{q\to 0}{\lim }\frac{m{\omega }^2}{q^2}\mathrm{Im}{\chi }_L(\omega ,q),$$

$$\delta H$$

$$\frac{dE}{dt}$$

$${\chi }_{\mathcal{A},\mathcal{B}}(i{\omega }_n)$$

$$G(\tau ,\mathbf{k})$$

$$T_{xy}^0(\mathbf{k})=\frac{{\hslash }^2}{m}{k}_x{k}_y.$$

$$G^{(m)}(\tau ,\mathbf{k})$$

$${\int }_0^{\infty }d(t_1,\dots ,t_n)$$

$$\times [\Theta (t_1)f_1(t_1)]\cdots [\Theta (t_n)f_n(t_n)]\Theta (\tau -t_1-\dots -t_n)$$

$$={\int }_{\mathrm{Bw}}\frac{ds}{2\pi i}{e}^{-\tau s}{f}_1(s)\cdots f_n(s),$$

$$T_2^{-1}(s,\mathbf{P})$$

$$L_{\mathbf{P}}(\tau )$$

$$L_{\mathbf{P}}(s)$$

$$T_2^{-1}(s)$$

$${\tilde{T}}_2(s)$$

$${\lambda }_T^{4}C$$

$$\mathrm{Im}{\chi }_{xy,xy}^{\mathrm{R},\mathrm{ni}}(\omega )=\delta (\beta \hslash \omega )(1-e^{-\beta \hslash \omega })$$

$$\displaystyle\qquad\times 2\int\frac{d^{d}k}{(2\pi)^{d}}\biggl{(}-\frac{\partial f(\varepsilon_{{\bf k}})}{\partial\varepsilon_{\bf k}}\biggr{)}T^{0}_{xy}({\bf k})T^{0}_{xy}({\bf k}),$$

$${\chi }_{xy,xy}^{\mathrm{SE},\mathrm{a}}(i{\omega }_n)=2z^2\int \frac{d^{d}p}{(2\pi {)}^d}\int \frac{d^{d}k}{(2\pi {)}^d}{T}_{xy}^0(\mathbf{k})T_{xy}^0(\mathbf{k})$$

Full text

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 $\eta$ and bulk viscosity $\zeta$ Schäfer and Chafin (2012); Schäfer (2014), the experimental determination of which has been an active topic over the past decade Cao et al. (2011a, b); Vogt et al. (2012); Elliott et al. (2014); Joseph et al. (2015). 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 for the shear viscosity to entropy ratio $\eta/s\sim\hbar/k_{B}$ 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 scale-invariant systems Son (2007).

Formally, the shear and bulk viscosities are given as the zero-frequency limit of corresponding viscosity spectral functions $\eta(\omega)$ and $\zeta(\omega)$, which are defined in terms of Kubo relations for the retarded stress tensor correlation function

[TABLE]

where $d$ is the dimension of the gas. The retarded response function $\chi_{ij,kl}^{R}(\omega)$ is the Fourier transform of the real-time response $\chi_{ij,kl}^{R}(t)=-i\hbar^{-1}\Theta(t)\langle[\hat{\Pi}_{ij}(t),\hat{\Pi}_{kl}(0)]\rangle$, with $\langle\cdots\rangle$ the thermal average and $\hat{\Pi}_{ij}$ the stress tensor. In Eq. (2) there is a summation convention for the doubly occurring space indices $i$ and $j$. 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). The frequency-dependent viscosity correlation functions are not only important in theoretical studies to obtain the static limit, however. Linear combinations are also related to the long-wavelength limit of the dynamic structure factor Taylor and Randeria (2010), which can be measured using Bragg spectroscopy Veeravalli et al. (2008); Carcy et al. (2019). Moreover, as discussed further in this paper, the bulk viscosity spectral function describes the energy absorption in response to an oscillating scattering length Son and Thompson (2010).

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 $\eta(\omega),\zeta(\omega)\stackrel{{\scriptstyle\omega\to\infty}}{{\longrightarrow}}C\omega^{-\frac{4-d}{2}}$ Taylor and Randeria (2010); Hofmann (2011); Goldberger and Khandker (2012), with a magnitude set by the contact parameter $C$ 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 spectral functions.

An 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 some limiting cases, such as 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). For the spin-balanced gas, kinetic theory applies at both very low Landau and Khalatnikov (1949) and very high temperatures Yan et al. (2019). In this paper we study the high-temperature or nondegenerate limit, in which the thermal wavelength $\lambda_{T}=\sqrt{2\pi\hbar^{2}\beta/m}$ is much smaller than the interparticle distance, $n\lambda_{T}^{d}\ll 1$, where $n$ is the density, $\beta$ inverse temperature, and $m$ the mass of an atom. A drawback of kinetic theory is that the high-frequency tails of the viscosity spectral functions are not correctly captured and the extrapolation to smaller temperatures is not controlled. While a formal correspondence between a kinetic description and the classical limit of the microscopic theory is 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 $\hat{\cal A}$ as a sum over expectation values restricted to the $N$-particle sector

[TABLE]

where the expansion parameter $z=e^{\beta\mu}$ is called the fugacity. Here ${\rm Tr}_{N,c}$ denotes the connected trace restricted to the $N$-particle sector, i.e.,

[TABLE]

and so on, where ${\rm Tr}_{N}$ sums all $N$-particle states. Since the $N$th-order term of Eq. (3) involves only matrix elements of $N$-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. In particular, for the particle density, we obtain

[TABLE]

with $b_{i}$ the virial coefficients Liu (2013). Truncating the expansion (3) or (6) after the first few terms holds for $z\ll 1$, which corresponds to the nondegenerate or high-temperature regime $n\lambda_{T}^{d}\ll 1$. Note that the quantum virial expansion is valid for any interaction strength, and in this sense it is a nonperturbative 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 is the exact viscosity spectral functions $\eta(\omega)$ and $\zeta(\omega)$ 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 nonanalytic 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 a power law $\lim_{\omega\to 0}\eta(\omega)\sim 1/\omega^{2}$, 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 low-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 function 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 in 3D near the unitary limit

[TABLE]

where the zero-frequency limit of the spectral function was taken before the unitary limit, with $a$ the scattering length and $\gamma_{E}$ the Euler-Mascheroni constant. This result differs from kinetic theory, which predicts a quadratic dependence $\lambda_{T}^{3}\zeta/\hbar z^{2}\sim(\lambda_{T}/a)^{2}$ Dusling and Schäfer (2013), but it is consistent with a predicted critical scaling near unitarity Maki et al. (2018); Maki and Zhou (2019).

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 order in the fugacity $z$. 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 Appendix 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. Furthermore, we demonstrate that 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

[TABLE]

where $\sigma=\uparrow,\downarrow$ and $g$ is the bare interaction strength. The stress tensor $\hat{\Pi}_{ij}$ in the zero-range model (8) takes the form (the summation over $\sigma$ is implied)

[TABLE]

This expression is derived, for example, from the Euler equation in operator form, i.e., the Heisenberg equation of motion for the current operator, $-\frac{im}{\hbar}[\hat{j}_{i},\hat{H}]=-\partial_{j}\hat{\Pi}_{ij}$, where $\hat{j}_{i}=\frac{-i\hbar}{2m}[\psi^{\dagger}_{\sigma}(\partial_{i}\psi)-(\partial_{i}\psi^{\dagger}_{\sigma})\psi]$ and the Hamiltonian is given in Eq. (8). 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:

[TABLE]

We use a summation convention where the sum over the space index $i$ is implied. It was established in Refs. Martinez and Schäfer (2017); Czajka and Jeon (2017); Fujii and Nishida (2018) that the bulk viscosity is also defined in terms of the response function

[TABLE]

with ${\cal O}=\frac{1}{d}(\hat{\Pi}_{ii}-2\hat{\cal E})$, where $\hat{\cal E}$ is the energy density with $\hat{H}=\int d^{d}x\,\hat{\cal E}$. The operator $\hat{\cal O}$ is linked to the Tan contact operator and is a measure of the breaking of scale invariance in the quantum gas Hofmann (2012)

[TABLE]

where $a$ is the effective scattering length in $d$ dimensions and ${\hat{C}}$ is the contact operator Tan (2008a, b, c), ${\hat{C}}=\frac{m^{2}g^{2}}{\hbar^{4}}\hat{\psi}_{\uparrow}^{\dagger}\hat{\psi}_{\downarrow}^{\dagger}\hat{\psi}_{\downarrow}\hat{\psi}_{\uparrow}$ in the zero-range model Braaten and Platter (2008); Braaten et al. (2008). The contact parametrizes 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 correlation 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)

[TABLE]

Indeed, the contact also governs the nonanalytic short-time correlation of the stress correlator $\chi_{ij,kl}^{R}(t)$, 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 ${\hat{O}}$ sets the deviation of the pressure ${\cal P}$ from the scale-invariant value: ${\cal P}-\frac{2}{d}{\cal E}=\langle\hat{O}\rangle$. In two dimensions, the operator $\hat{\cal O}$ indicates a quantum scale anomaly, where formally it describes an anomalous contribution that opens up the nonrelativistic 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 $\chi_{{\cal O},{\cal O}}^{R}$ 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)

[TABLE]

where $\chi_{L}$ and $\chi_{T}$ are the longitudinal and transverse parts of the current response function, $\chi^{JJ}_{ij}=\chi_{L}\frac{q_{i}q_{j}}{q^{2}}+\chi_{T}(\delta_{ij}-\frac{q_{i}q_{j}}{q^{2}})$. However, this definition requires the computation of correlation functions not only at finite external frequency but also nonzero momentum, which complicates the calculation in the present case.

Both the shear viscosity spectral function and the bulk viscosity spectral function can be measured in experiments. Since the longitudinal current response function is related to the density response function (and thus the dynamic structure factor) by the continuity equation, the linear combination $\zeta(\omega)+\frac{2(d-1)}{d}\eta(\omega)$ is measured in Bragg spectroscopy experiments Veeravalli et al. (2008); Carcy et al. (2019). In particular, in scale-invariant systems (such as the unitary Fermi gas), where the bulk viscosity vanishes, the dynamic structure factor is proportional to the shear viscosity. Moreover, the contact correlation function used in Eq. (11) to define the bulk viscosity determines the response to an oscillating modulation of the inverse scattering length, i.e., a perturbation of the form $\delta\alpha(t)=\delta a^{-1}(t)=\alpha_{0}\cos\omega t$ Son and Thompson (2010). Such a perturbation changes the Hamiltonian as

[TABLE]

The rate at which energy is absorbed by the system in response to this perturbation, for example, is proportional to the contact correlation function

[TABLE]

II.1 Diagrammatic virial expansion

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

[TABLE]

where $i\omega_{n}$ is a bosonic Matsubara frequency and ${\cal T}_{\tau}$ denotes time ordering with operators at later times placed to the left. Retarded response functions as needed for the spectral functions (1) and (11) are obtained by analytically continuing from Matsubara to real frequencies $i\omega_{n}\to\omega+i0$ Mahan (2000). When working in Matsubara space, the fugacity parameter $z=e^{\beta\mu}$ is contained in Bose and Fermi factors used to convert Matsubara frequency summations to contour integrals Mahan (2000). An expansion in $z$ 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 Green’s function

[TABLE]

where ${\bf k}$ is the wave vector, $\varepsilon_{\bf k}=\frac{\hbar^{2}k^{2}}{2m}$ the single-particle energy, and $n_{F}(\varepsilon_{\bf k}-\mu)$ the Fermi factor, which can be expanded directly as $n_{F}(\varepsilon_{\bf k}-\mu)=-\sum_{n=1}^{\infty}(-1)^{n}z^{n}e^{-n\beta\varepsilon_{\bf k}}$. 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 this section we collect the Feynman rules as used in the remainder of this paper.

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

[TABLE]

Likewise, the contact operator insertion carries a factor of $\frac{m^{2}g^{2}}{\hbar^{4}}$. As shown in Fig. 2, we indicate the ${\cal O}(z^{m})$ term of the single-particle Green’s function (19) by a line that is slashed $m$ times. It contributes factors of

[TABLE]

Contributions ${\cal O}(z^{N})$ to some correlation function are represented by Feynman diagrams that contain $N$ slashes. A key observation is that the leading order in $z$ of the propagator is purely retarded, i.e., it runs from the left to right in any Feynman diagram. 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

[TABLE]

where we define the Laplace transform $f(s)=\int_{0}^{\infty}dt\,e^{st}f(t)$ with inverse $f(t)=\int_{\rm Bw}\frac{ds}{2\pi i}\,e^{-st}f(t)$. Here 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 subblocks. It turns out that to leading order in the fugacity ${\cal O}(z^{0})$, 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 $T_{2}$, which we represent by a box as shown in Fig. 2. Using the convolution theorem and performing the Laplace transform, $T_{2}$ is given by a geometric series Barth (2015)

[TABLE]

where ${\bf P}$ is the center-of-mass momentum and we define the bare two-particle propagator

[TABLE]

The momentum integration diverges in $d=2,3$, which is canceled by the renormalization of the coupling constant $g$. For reference, we note the result for the $T$ matrix

[TABLE]

where, by Galilean invariance, $T_{2}(s,{\bf P})=T_{2}(s-\frac{1}{2}\varepsilon_{\bf P})$. Recall the structure of the imaginary part of the $T$ matrix, ${\rm Im}\,T_{2}(s+i0)$, which consists of a bound-state peak (where present) at the bound-state energy $E_{b}=-\frac{\hbar^{2}}{ma^{2}}$, and a continuum of scattering states for positive energies. It will be convenient at a later point to use dimensionless scattering matrices $\tilde{T}_{2}(s)$ defined via

[TABLE]

where $\Omega_{d}=2\pi^{d/2}/\Gamma(d/2)$ is the $d$-dimensional surface element, i.e., we separate a factor $\frac{2\hbar^{2}}{m\lambda_{T}}$, $\frac{2\pi\hbar^{2}}{m}$, and $\frac{4\pi\hbar^{2}\lambda_{T}}{m}$ in $d=1,2$, and $3$, respectively.

For further reference, we note a useful representation of the contact to leading order ${\cal O}(z^{2})$ in the fugacity:

[TABLE]

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.

Note that in this paper, we focus on the viscosity of Fermi gases, but all calculations apply with minor modifications to the Bose gas as well. The zero-range model for Bose gases replaces the Fermi by Bose fields in Eq. (8) and conventionally sets $g\to g/4$. The scattering $T$ matrix in Eq. (23) then has an additional symmetry factor of $1/2$ in front of the loop integral, and the overall Bose $T$ matrix is larger by a factor of $2$ compared to the Fermi case Barth and Hofmann (2015). The inclusion of a three-body term in the zero-range model will not affect any results to leading order ${\cal O}(z^{2})$ in the fugacity, which is fully determined by two-body correlations Barth (2015).

II.2 Stress correlator

In this section we discuss the Feynman diagram expansion of the stress correlator $\chi_{xy,xy}^{R}(\omega)$ used in the computation of the shear viscosity (1) up to quadratic order in the fugacity ${\cal O}(z^{2})$. A Feynman diagram representation of the imaginary-time correlator $\chi_{xy,xy}(\tau)$ is obtained by drawing all connected diagrams with a stress tensor insertion $\hat{\Pi}_{xy}$ at $\tau=0$ and at $0<\tau<\beta$. All such diagrams up to quadratic order are shown in Fig. 3. We then take the discrete Fourier transform to obtain the response as a function of Matsubara frequencies. The leading-order ${\cal O}(z)$ contribution [Fig. 3(a)] is the same as for noninteracting gas. Indeed, the noninteracting contribution can be summed to all orders in $z$. Using the Feynman rules introduced in the preceding section, we obtain the standard kinetic result

[TABLE]

where the prefactor $2$ accounts for the two spin species. This term will contribute a $\delta$ function to the shear viscosity spectral function $\eta(\omega)$.

Two-particle interactions enter at second order ${\cal O}(z^{2})$. The Feynman diagrams that represent the interacting ${\cal O}(z^{2})$ contribution to the shear viscosity spectral function (1) are shown in Figs. 3(b) and 3(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, focusing on the main steps of the calculation. To illustrate the calculations in detail, Appendix B contains an in-depth discussion of the self-energy contribution.

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, the first one evaluates to

[TABLE]

The overall prefactor $2$ accounts for the spin summation. The Bromwich integral and the analytic continuation from Matsubara to real frequency $i\omega_{n}\to\omega+i0$ are performed as outlined in Appendix A. This gives the self-energy contribution to the retarded response:

[TABLE]

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

[TABLE]

The imaginary part of the self-energy contribution is antisymmetric in the frequency, as expected for a retarded response function. Appendix B contains a detailed discussion of the above derivation.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The diagram now takes a simpler form

[TABLE]

and we split off a second term

[TABLE]

II.3 Contact correlator

We now turn to the contact correlator $\chi_{C,C}(\omega)$ 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 ${\cal O}(z^{2})$. Corresponding diagrams are shown in Fig. 4. They involve the $T$ matrix coefficients and the bare two-particle propagator $L_{\bf P}(\tau)$ [Eq. (24)] as subblocks. The calculation proceeds in a similar way as in the preceding section. After analytic continuation, the sum of all four diagrams is

[TABLE]

where $L_{\bf P}(s)$ is the Laplace transform of the bare two-particle propagator defined in Eq. (25). Equation (40) is simplified using the definition of the $T$ matrix (23), which gives

[TABLE]

This expression is antisymmetric in $\omega$ as can be seen by shifting the integration variable $s\to s-\omega$. The integral (41) is convenient for a direct numerical evaluation.

III Results

In this section we combine the results of the preceding 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 noninteracting ${\cal O}(z)$ contribution [Eq. (29)]. At finite frequency, the spectral function is ${\cal O}(z^{2})$ and given by the sum of self-energy, Maki-Thompson, and Aslamazov-Larkin contribution

[TABLE]

where we have symmetrized the stress tensor matrix element $T_{xy}^{0}$. For a numerical evaluation, it is convenient to transform to relative and center-of-mass momenta ${\bf P}=({\bf k}+{\bf p})/2$ and ${\bf q}={\bf k}-{\bf p}$, which restricts the angular integration to the matrix elements $T^{0}_{xy}({\bf q})T^{0}_{xy}({\bf q})$. The center-of-mass momentum integration is then performed directly.

The imaginary part of the scattering $T$ matrix $T_{2}(s+i0)$ [Eq. (26)] has a continuous spectrum at positive energies and (if present) a state pole at the bound-state energy $E_{b}=-\frac{\hbar^{2}}{ma^{2}}$. From the structure of the integrand in Eq. (42), we deduce the general form of the shear viscosity spectral function shown in Fig. 1(a). First, the free-particle part of the $T$ 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 (42) gives rise to a bound-state contribution with a threshold behavior $\sim(\hbar\omega+E_{b})^{(d+2)/2}$, which is nonanalytic in 3D. At small frequencies, only the continuum part contributes to the spectrum. As is apparent from the integral (42), the spectral function diverges as $1/\omega^{2}$ 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 $\lambda_{T}/a$. 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. (42) with $\omega\to-\omega$. Transforming to a relative-energy integration variable $s=2\varepsilon_{\bf q}$, we shift $s\to s+\omega$ and extend the $s$ integration to run from $-\infty$ to $+\infty$. Separating the leading power in $\hbar\omega$ from the integrand, we recognize the remaining integral as the leading term in the virial expansion of the contact [Eq. (28)]. This gives

[TABLE]

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

In the three-dimensional unitary limit, we can compare with a Keldysh calculation Enss et al. (2011), which was performed using a finite lifetime $\gamma/2$ in the single-particle self-energy. Here, $C=16\pi z^{2}/\lambda_{T}^{4}=4m^{2}z^{2}/\pi\beta^{2}\hbar^{4}$ and we state the separate contributions of the self-energy, Maki-Thompson, and Aslamazov-Larkin terms to the high-frequency tail

[TABLE]

which agrees with Ref. Enss et al. (2011). This indicates 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)

[TABLE]

which has a finite limit for $\omega\to 0$. The static shear viscosity $\eta=\lim_{\omega\to 0}\frac{{\rm Im}\,\chi(\omega)}{\omega}=W\tau$ is then expressed as the ratio of a Drude weight $W$ and a viscous scattering rate $\tau^{-1}$. The scattering rate $\tau^{-1}$, which vanishes for noninteracting systems, admits a separate fugacity expansion. A direct virial expansion of the correlation function thus applies in the limit $\omega\tau\gg 1$ and does not connect with the hydrodynamic expression (45). This is the origin of the $1/\omega^{2}$ divergence of the shear viscosity spectral function (42) for $\beta\hbar\omega\ll 1$.

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 $M(\omega)$ defined through

[TABLE]

with $\eta_{0}$ a real constant. Note that the relation between the correlation function and the memory function is similar to the relation between the single-particle Green’s function and self-energy. Since $\chi(-\omega)=\chi^{*}(\omega)$, we have $M(-\omega)=-M^{*}(\omega)$. Using the low-frequency expansion $M(\omega)=iM_{0}+\omega M_{1}+{\cal O}(\omega^{2})$ with $M_{0}$ and $M_{1}$ real, we find

[TABLE]

which implies, for the static limit of the viscosity spectral function,

[TABLE]

We can extract the coefficients $\eta_{0}$ and $M_{0}$ from the virial expansion results by matching the low-frequency terms of the virial expansion with the expansion of Eq. (46),

[TABLE]

In this way, we obtain a static limit of the shear viscosity using (48). Intuitively, the memory kernel result for the static viscosity (48) is obtained from a geometric resummation of the leading two low-frequency terms in Eq. (49). Note that a similar approach has 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 three-dimensional case in the following. First from Eq. (29) we obtain $\eta_{0}$,

[TABLE]

valid to leading order in the fugacity. Second, the coefficient $\eta_{0}M_{0}$ follows from the low-frequency divergence of Eq. (42), which is obtained by expanding the prefactor and setting $\omega=0$ in the integrand. While at first sight this expression might seem unwieldy, we will show that it takes the form of a (linearized) collision integral such that $M_{0}$ can be identified with a viscous scattering time $\tau_{\eta}^{-1}$. In this way, using Eq. (48), the virial expansion combined with the memory kernel approach reproduces the calculation of the viscosity within kinetic theory.

Within kinetic theory, the static shear viscosity takes the form Massignan et al. (2005); Bruun and Smith (2005, 2007)

[TABLE]

where we use the definition ion of $\eta_{0}$ in Eq. (50) and following the literature we introduce a viscous relaxation time $\tau_{\eta}$ Massignan et al. (2005); Bruun and Smith (2005, 2007),

[TABLE]

with a normalization factor

[TABLE]

which has dimension of inverse volume. In Eq. (52), $\frac{d\sigma}{d\Omega}$ is the differential scattering cross section for the scattering of two particles with opposite spin and wave vectors ${\bf k}$ and ${\bf p}$ to states ${\bf k}^{\prime}$ and ${\bf p}^{\prime}$, respectively, and ${\bf v}_{\rm rel}=\frac{2\hbar{\bf q}}{m}$ is the relative velocity.

We now show that $\tau_{\eta}^{-1}$, given in Eq. (52), and $M_{0}$, given by Eq. (42), are identical at high temperatures, i.e., $\lim_{T\to\infty}\tau_{\eta}^{-1}=M_{0}$. At high temperatures, Pauli blocking is not effective and the Fermi factors in Eq. (52) reduce to a Maxwell-Boltzmann distribution

[TABLE]

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

[TABLE]

The equality between $\tau_{\eta}^{-1}$ and $M_{0}$ then follows directly by comparing with Eq. (42). From Eq. (48) we then obtain the static shear viscosity. 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). 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

[TABLE]

Our result agrees with Ref. Enss et al. (2011) if we identify $\hbar\omega\to 2i\gamma$.

III.2 Bulk viscosity

This section presents results for the bulk viscosity spectral function $\zeta(\omega)$. From the definitions (11) and (12) and our result for the contact correlation function (41), we obtain the integral representation of the spectral function:

[TABLE]

where the dimensionless scattering matrix $\tilde{T}_{2}$ is defined in Eq. (27). This expression is a symmetric function of the frequency $\omega$. 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. (57), we directly infer the general form of the spectrum shown in Fig. 1(b). The 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. (57) is finite at $\omega=0$; hence the spectral function approaches a constant value for small frequencies. The second contribution appears for frequencies larger than the bound-state energy, for which we pick up the pole of the first scattering matrix. Notably, it has a nonanalytic frequency dependence right above the threshold energy which follows from the nonanalytic form of the $T$ matrix (26). Finally, a third contribution at $\omega=0$ arises from the pinching of the bound-state poles in the integrand (57). This part gives rise to a $\delta$-function contribution $\sim\delta(\omega)$. 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 $\lambda_{T}/a$ are shown in Fig. 6 for all dimensions and various values of the scattering length. The entire spectral function vanishes identically for $a\to\infty$. This point marks the unitary or noninteracting scale-invariant limit, where the bulk viscosity is expected to vanish on general grounds Son (2007). Away from this limit, we 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 nonanalytic 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 $a\sqrt{\frac{m\omega}{\hbar}}$ is kept fixed, the dependence on the integration variable of the second term in the integrand of Eq. (57) can be neglected. The remaining integral is identified with the contact expectation value at this order in the virial expansion [Eq. (28)]. We obtain the high-frequency tails

[TABLE]

These results agree with the OPE result in 3D Hofmann (2011); Goldberger and Khandker (2012) and the calculations in 2D Hofmann (2011). 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)

[TABLE]

These sum rules are obeyed by our numerical results.

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

[TABLE]

where ${\rm Ei}$ is the exponential integral. In the perturbative limit of small negative scattering length, this becomes

[TABLE]

In the unitary limit $a\to\infty$, we obtain the result given in the Introduction, Eq. (7). The prefactor of $(\lambda_{T}/a)^{2}$ in Eq. (7) is consistent with an expected scaling with $({\cal P}-\frac{2}{3}{\cal E})^{2}/{\cal P}^{2}\sim\langle\hat{\cal O}\rangle^{2}/{\cal P}^{2}$ that vanishes in the unitary limit and measures the deviation from scale invariance Dusling and Schäfer (2013); Martinez and Schäfer (2017). Note that for the bulk viscosity in high-temperature QCD, a similar logarithmic-quadratic scaling with interaction strength is found Arnold et al. (2006). Our findings are at variance with a previous kinetic theory calculation near the unitary limit, which does not include logarithmic renormalization effects and gives a different result Dusling and Schäfer (2013)

[TABLE]

While the exact high-frequency tail and the sum rule constraints provide very compelling checks of our virial expansion and indicate the consistency of our results, a disagreement with kinetic theory is unexpected, and further work warranted.

Note that the finite result (7) for the static bulk viscosity at unitarity is obtained only when the unitary limit of infinite scattering length ($a\to\infty$) is taken after the zero-frequency limit ($\omega\to 0$). However, these limits do not commute. By contrast, the low-frequency limit of the bulk viscosity spectral function near unitarity, i.e., taking the low-frequency limit after the unitary limit, diverges logarithmically with frequency:

[TABLE]

This result appears consistent with recent work on the nonequilibrium dynamics of a nearly unitary quantum gas Maki et al. (2018); Maki and Zhou (2019); Maki and Ordóñez (2019), which shows that the unitary conformal fixed point governs the long-time quench dynamics in these systems at times $1\ll\ln t\ll a$. The scaling for the bulk viscosity in the time domain, $\zeta(t)\sim(\ln t)/a^{2}t$ [Eq. (77) of Ref. Maki and Zhou (2019)], is consistent with the logarithmic low-frequency singularity (63). The finite result for the static bulk viscosity quoted in the Introduction, Eq. (7), then corresponds to evaluating the correlation function outside of the conformal critical window.

IV Conclusion

In this paper, we have presented exact nonperturbative results for the shear and bulk viscosity spectral functions calculated to quadratic order in the fugacity. The field-theoretical methods employed here allow us to make contact with many-body techniques such as the Keldysh formalism and can be compared with numerical calculations. A significant advantage of our method is that the analytic continuation from discrete Matsubara space to real frequencies is performed in closed analytical form and does not require numerical extrapolation schemes. The form of the viscosity spectral functions is linked to the two-particle scattering structure: Spectra consist of a broad continuous band superimposed with an additional bound-state contribution above the dimer-breaking threshold frequency. The results of the virial expansion are internally consistent in that they obey universal Tan relations for the high-frequency tail and, in the case of the bulk viscosity, saturate exact sum rules.

It is interesting to note that the correspondence between the exact virial expansion in the static limit and kinetic theory is quite nontrivial. For the static shear viscosity, a memory kernel resummation of the leading low-frequency divergence in the spectral function is required to match with kinetic theory. A similar resummation scheme for the bulk viscosity is not apparent as the spectral function saturates to a finite value, albeit with an additional dimer peak where present. A direct extrapolation to the static limit for three-dimensional fermions with negative scattering length, i.e., without bound states in the two-body spectrum, gives a nonanalytic logarithmic-quadratic result near unitarity [Eq. (7)], which differs from kinetic theory. It would be interesting to resolve this discrepancy in further work.

There are several additional applications of the methods used in this paper. For example, at the next-to-leading (third) order in the fugacity, the spectral functions include three-body effects, which provide an additional way to break scale invariance (besides an explicit finite two-body scattering length or a quantum anomaly). For the bulk viscosity of the unitary Bose gas, which is scale invariant in the two-body sector, this third-order term will be the leading nonzero order. Likewise, calculations of thermal transport coefficients will be performed in a similar way as outlined in this paper.

Note added. Recently, Refs. Enss (2019); Nishida (2019) appeared, which discuss the virial expansion of the viscosity spectral functions using other methods and have overlap with this work. The overlapping results are in agreement.

Acknowledgements.

I thank Sergej Moroz and Wilhelm Zwerger for discussions. I also thank Tilman Enss, Jeff Maki, Yusuke Nishida, and Fei Zhou for pointing our their work and discussions. This work was supported by Peterhouse, Cambridge.

Appendix A Analytic continuation

The quantum cluster expansion gives expressions for response functions in imaginary time or, after discrete Fourier transformation, in Matsubara space for frequencies $i\hbar\omega_{n}=2\pi n/\beta$. Frequency-dependent retarded response functions are then obtained by analytic continuation $i\omega_{n}\to\omega+i0$. A considerable advantage of the cluster expansion compared to numerical methods is that this analytic continuation is performed in closed analytic form and does not require an interpolation or other numerical techniques. However, some care must be taken Mahan (2000). This appendix outlines the analytic continuation for the integrals discussed in the main text.

After performing the Laplace transformation for the Feynman diagram integrals, we obtain expressions of the form

[TABLE]

where $f$ and $g$ are products of the single-particle propagator or the two-body $T$ matrix, which have branch cuts or poles on the real axis (we omit momentum integrals for brevity). Here Bw is the Bromwich contour, which by definition avoids every singularity of the integrand. This is shown in Fig. 7, where the black solid lines indicate the branch cuts or poles of the integrand for ${\rm Im}\,s=0$ and ${\rm Im}(s+i\omega_{n})=0$. Deforming the contour of integration to run along the cuts ($Bw^{\prime}$ in Fig. 7), we obtain

[TABLE]

where the integration range is extended to run along the complete real axis. We can now safely perform the analytic continuation $i\omega_{n}\to\omega+i0$ to obtain the retarded response function. The imaginary part reads

[TABLE]

The prefactor $(1-e^{-\beta\hbar\omega})$ is expected on general grounds from the fluctuation-dissipation theorem Forster (1990). For the stress correlator self-energy and Maki-Thompson diagrams, one factor $f$ or $g$ involves only simple poles from single-particle propagators, so the integration is performed using the residue theorem. The same is true for the Aslamazov-Larkin diagrams when using the Ward identity (36).

Appendix B Example calculation

To illustrate the diagrammatic framework discussed in the main text, this appendix presents the calculation of the self-energy contribution to the shear viscosity spectral function in detail. The corresponding Feynman diagrams are shown in Fig. 3(b) and also in Fig. 8 including the imaginary-time axis. The calculation of other diagrams proceeds in a similar way.

The Feynman diagrams discussed in this paper lie on an imaginary-time axis, with imaginary time running from the left-hand side to the right-hand side. We partition the interval $[0,\tau]$ in three parts with interval lengths $t_{1}$,$t_{2},$ and $t_{3}$ as shown in Fig. 8, where $\tau=t_{1}+t_{2}+t_{3}$. Both diagrams in Fig. 8 contain two retarded propagators, both with momentum ${\bf k}$. Using the Feynman rules in Eq. (21), they contribute factors of $[-\Theta(t_{1})\,e^{-t_{1}\varepsilon_{\bf k}}]$ and $[-\Theta(t_{3})\,e^{-t_{3}\varepsilon_{\bf k}}]$ for the first diagram and $[-\Theta(t_{1})\,e^{-t_{1}\varepsilon_{\bf k}}]$ and $[-\Theta(t_{2})\,e^{-t_{2}\varepsilon_{\bf k}}]$ for the second diagram (we neglect an additional factor $e^{\tau\mu}$, the argument of which will add to zero as all fermion lines are closed). The diagrams contain two-body scattering matrices, which sum all repeated scattering events shown in Fig. 2 and which contribute terms of $[\Theta(t_{2})T_{2}(t_{2},{\bf k}+{\bf p})]$ and $[\Theta(t_{3})T_{2}(t_{3},{\bf k}+{\bf p})]$, respectively. The overall leading order in the fugacity of these diagrams is determined by the backward propagating lines, which are ${\cal O}(z)$, such that the total diagram contributes at order ${\cal O}(z^{2})$ to the response function. The advanced Green’s functions contribute factors $[ze^{-(\beta-t_{1}-t_{2}-t_{3})\varepsilon_{\bf k}}]$ and $[ze^{-(\beta-t_{2})\varepsilon_{\bf p}}]$ or $[e^{-(\beta-t_{2})\varepsilon_{\bf p}}]$, respectively. Having evaluated these factors, we integrate over the undetermined loop momenta ${\bf k}$ and ${\bf p}$ and take the discrete Fourier transform to obtain the response in Matsubara space.

Collecting all terms, the first Feynman diagram evaluates to

[TABLE]

where the prefactor of $2$ accounts for a spin summation. Rearranging the factors gives

[TABLE]

We now apply the convolution theorem (22), which gives the expression (30). The second line of Eq. (30) takes the form (64) with

[TABLE]

Performing the analytic continuation as outlined in Appendix A, we end up with an expression of the form (66). As mentioned above, this expression is further simplified using ${\rm Im}[f(s+i0)]=-\pi\delta(s-\varepsilon_{{\bf k}}-\varepsilon_{{\bf p}})$, which removes the $s$ integration. This then yields the expression (31).

Repeating the same calculation for the second diagram in Fig. 8, the Bromwich integral takes a similar form as before with the functions $f$ and $g$ interchanged:

[TABLE]

Evaluating the residue of $g(s)$ then gives the expression (32).