Bulk viscosity and contact correlations in attractive Fermi gases

TL;DR

This paper investigates how pair fluctuations influence bulk viscosity in attractive Fermi gases, revealing deviations from scale invariance and effects near the superfluid transition, with numerical results across the BEC-BCS crossover.

Contribution

It provides the first numerical Luttinger-Ward calculations of contact correlations and bulk viscosity throughout the BEC-BCS crossover, highlighting the role of pair fluctuations.

Findings

Bulk viscosity is linked to contact density correlations.

Pair fluctuations break scale invariance, affecting viscosity.

Bulk viscosity peaks near the superfluid transition.

Abstract

The bulk viscosity determines dissipation during hydrodynamic expansion. It vanishes in scale invariant fluids, while a nonzero value quantifies the deviation from scale invariance. For the dilute Fermi gas the bulk viscosity is given exactly by the correlation function of the contact density of local pairs. As a consequence, scale invariance is broken purely by pair fluctuations. These fluctuations give rise also to logarithmic terms in the bulk viscosity of the high-temperature nondegenerate gas. For the quantum degenerate regime I report numerical Luttinger-Ward results for the contact correlator and the dynamical bulk viscosity throughout the BEC-BCS crossover. The ratio of bulk to shear viscosity $\zeta/\eta$ is found to exceed the kinetic theory prediction in the quantum degenerate regime. Near the superfluid phase transition the bulk viscosity is enhanced by critical fluctuations and has observable effects on dissipative heating, expansion dynamics and sound attenuation.

Equations

$$\zeta (\omega )=-\frac{1}{\omega d^2}\mathrm{Im}{\int }_0^{\infty }dt{e}^{i\omega t}\int d^{d}x⟨[{\hat{\Pi }}_{ii}(\mathbf{x},t),{\hat{\Pi }}_{jj}(0,0)]⟩,$$

$$\displaystyle\hat{\mathcal{H}}=\sum_{\sigma}\psi_{\sigma}^{\dagger}\bigl{(}-\frac{\nabla^{2}}{2m}\bigr{)}\psi_{\sigma}+g_{0}\psi_{\uparrow}^{\dagger}\psi_{\downarrow}^{\dagger}\psi_{\downarrow}\psi_{\uparrow}.$$

$$d\cdot \hat{P}={\hat{\Pi }}_{ii}=[\hat{\mathcal{H}},i\hat{D}]=2\hat{\mathcal{H}}+\{\begin{array}{ll}\frac{\hat{\mathcal{C}}}{2\pi m}& \text{(2D)},\\ \frac{\hat{\mathcal{C}}}{4\pi ma}& \text{(3D),}\end{array}$$

$${\hat{\Pi }}_{\text{an}}\equiv {\hat{\Pi }}_{ii}-2\hat{\mathcal{H}}=c_d\hat{\mathcal{C}},$$

$$\zeta (\omega )=-\frac{c_d^2}{\omega d^2}\mathrm{Im}{\int }_0^{\infty }dt{e}^{i\omega t}\int d^{d}x⟨[\hat{\mathcal{C}}(\mathbf{x},t),\hat{\mathcal{C}}(0,0)]⟩.$$

$$\displaystyle\chi(\bm{x},t)\equiv\langle[\hat{\mathcal{C}}(\bm{x},t),\hat{\mathcal{C}}(0,0)]\rangle=\begin{cases}2\pi m\bigl{(}\frac{\partial\langle\hat{\mathcal{C}}(\bm{x},t)\rangle}{\partial\ln a(0,0)}\bigr{)}_{s}&\text{(2D)},\\ -4\pi m\bigl{(}\frac{\partial\langle\hat{\mathcal{C}}(\bm{x},t)\rangle}{\partial a^{-1}(0,0)}\bigr{)}_{s}&\text{(3D)}\end{cases}$$

$$\displaystyle S\equiv\frac{2}{\pi}\int_{0}^{\infty}d\omega\,\zeta(\omega)=\begin{cases}-\frac{1}{8\pi m}\bigl{(}\frac{\partial\mathcal{C}}{\partial\ln a}\bigr{)}_{s}&\text{(2D)},\\ \frac{1}{36\pi ma^{2}}\bigl{(}\frac{\partial\mathcal{C}}{\partial a^{-1}}\bigr{)}_{s}&\text{(3D)}.\end{cases}$$

$${\zeta }_{\text{2D,vac}}(\omega )=\frac{{\mathcal{C}}_0}{4m\omega }\frac{\Theta (\omega -{\epsilon }_B)}{{\ln }^2(\omega /{\epsilon }_B-1)+{\pi }^2}$$

$${\zeta }_{\text{3D,vac}}(\omega )=\frac{{\mathcal{C}}_0\Theta (a)}{36\pi ma}\frac{\sqrt{{\epsilon }_B(\omega -{\epsilon }_B)}\Theta (\omega -{\epsilon }_B)}{{\omega }^2}.$$

$$\zeta (\omega \to \infty )=\{\begin{array}{ll}\frac{\mathcal{C}}{4m\omega {\ln }^2(\omega /{\epsilon }_B)}& \text{(2D)},\\ \frac{\mathcal{C}}{36\pi a^2(m\omega {)}^{3/2}}& \text{(3D)}.\end{array}$$

$$\zeta_{\text{3D,vir}}(\omega>0)=\frac{2\sqrt{2}}{9}z^{2}\lambda^{-3}v^{2}\frac{1-e^{-\beta\omega}}{\beta\omega}\\ \times\Bigl{[}\Theta(v)2ve^{v^{2}}\frac{\sqrt{\beta\omega-v^{2}}\,\Theta(\beta\omega-v^{2})}{\beta\omega}\\ +\frac{1}{\pi}\int_{0}^{\infty}dy\,\frac{e^{-y}\sqrt{y(y+\beta\omega)}}{(y+v^{2})(y+\beta\omega+v^{2})}\Bigr{]}.$$

$$S_{\text{3D,vir}}=\frac{2\sqrt{2}}{9}z^{2}T\lambda^{-3}v^{2}\Bigl{[}(1+2v^{2})e^{v^{2}}(1+\operatorname{erf}(v))\\ +\frac{2v}{\sqrt{\pi}}-\Theta(v)4v^{2}e^{v^{2}}\Bigr{]}.$$

$${\zeta }_{\text{3D,vir}}^{\text{unitary}}(\omega )=\frac{2\sqrt{2}}{9\pi }{z}^2{\lambda }^{-3}{v}^2\frac{\sinh (\beta \omega /2)}{\beta \omega /2}{K}_0(\beta \omega /2).$$

$$\displaystyle\zeta_{\text{3D,vir}}(v)=\frac{2\sqrt{2}}{9\pi}z^{2}\lambda^{-3}v^{2}\bigl{[}-1-(1+v^{2})e^{v^{2}}\operatorname{Ei}(-v^{2})\bigr{]}.$$

$$\zeta_{\text{2D,vir}}(\omega)=2\pi z^{2}\lambda^{-2}\frac{1-e^{-\beta\omega}}{\beta\omega}\Bigl{[}\frac{\beta\varepsilon_{B}e^{\beta\varepsilon_{B}}\;\Theta(\omega-\varepsilon_{B})}{\ln^{2}(\omega/\varepsilon_{B}-1)+\pi^{2}}\\ +\int_{0}^{\infty}dy\,\frac{e^{-y}}{[\ln^{2}(yT/\varepsilon_{B})+\pi^{2}][\ln^{2}((yT+\omega)/\varepsilon_{B})+\pi^{2}]}\Bigr{]}.$$

$${\zeta }_{\text{2D,vir}}({\epsilon }_B/T)\simeq \frac{2\pi z^2{\lambda }^{-2}}{[{\ln }^2(T/2{\epsilon }_B)+{\pi }^2{]}^2}.$$

$$\zeta /\eta \propto ({\Pi }_{\text{an}}/{\Pi }_{ii}{)}^2=[(P-2\mathcal{E}/3)/P{]}^2,$$

$$\chi (\mathbf{x},\tau )=⟨{\Delta }^{†}\Delta (\mathbf{x},\tau ){\Delta }^{†}\Delta (0,0)⟩$$

$$\chi (\mathbf{x},\tau )$$

$$=m^4\Gamma (\mathbf{x},\tau )\Gamma (-\mathbf{x},-\tau )+\text{vertex corrections}.$$

$$\begin{gathered} X(i{\omega }_m)=\int d^{d}x\chi (\mathbf{x},i{\omega }_m) \\ =\frac{m^4}{\beta V}\sum _{\mathbf{q}{\epsilon }_n}\Gamma (\mathbf{q},i{\epsilon }_n)\Gamma (\mathbf{q},i{\epsilon }_n+i{\omega }_m)+\text{vtx. corr.}, \end{gathered}$$

$$\begin{gathered} X^R(\omega )=\frac{m^4}{V}\sum _{\mathbf{q}}\int d\epsilon A(\mathbf{q},\epsilon )b(\epsilon ) \\ \times [{\Gamma }^R(\mathbf{q},\epsilon +\omega )+{\Gamma }^A(\mathbf{q},\epsilon -\omega )]+\text{vtx. corr.} \end{gathered}$$

$$\begin{gathered} \frac{d^2}{c_d^2}\zeta (\omega )=-\frac{\mathrm{Im}{X}^R(\omega )}{\omega } \\ =\frac{\pi m^4}{\omega V}\sum _q{\int }_{-\infty }^{\infty }d\epsilon A(\mathbf{q},\epsilon )A(\mathbf{q},\epsilon +\omega )[b(\epsilon )-b(\epsilon +\omega )] \\ +\text{vtx. corr.} \end{gathered}$$

$${\Gamma }^R(\mathbf{q},\epsilon {)}^{-1}$$

$$A_{\text{2D}}(\mathbf{q},\epsilon )$$

$$\displaystyle\qquad\qquad+\frac{\Theta(\varepsilon+2\mu-\omega_{q})}{\ln^{2}[(\varepsilon+2\mu-\omega_{q})/\varepsilon_{B}]+\pi^{2}}\Bigr{]},$$

$$A_{\text{3D}}(\mathbf{q},\epsilon )$$

$$\displaystyle\qquad\qquad+\frac{1}{\pi}\,\frac{\sqrt{\varepsilon+2\mu-\omega_{q}}\,\Theta(\varepsilon+2\mu-\omega_{q})}{\varepsilon+2\mu+(ma^{2})^{-1}-\omega_{q}}\Bigr{]}.$$

$$\chi (\mathbf{x},\tau )$$

$$\simeq ⟨\Delta (\mathbf{x},\tau ){\Delta }^{†}(0,0)⟩⟨{\Delta }^{†}\Delta (0,0)⟩$$

$$=m^2\Gamma (\mathbf{x},\tau )\mathcal{C}(0,0).$$

$$\begin{gathered} \frac{d^2}{c_d^2}\zeta (\omega )=\frac{\pi m^4}{\omega }\int \frac{d^{3}q}{(2\pi {)}^3}\int d\epsilon A_{\text{3D}}(\mathbf{q},\epsilon )A_{\text{3D}}(\mathbf{q},\epsilon +\omega ) \\ \times z^2{e}^{-\beta (\epsilon +2\mu )}[1-e^{-\beta \omega }]. \end{gathered}$$

$$\frac{d^{2}}{c_{d}^{2}}\zeta(\omega)=16\pi^{3}mz^{2}\,\frac{1-e^{-\beta\omega}}{\omega}\int\frac{d^{3}q}{(2\pi)^{3}}\,e^{-\beta\omega_{q}}\int dy\,e^{-y}\\ \times\Bigl{[}2\sqrt{\beta\varepsilon_{B}}\delta(y+\beta\varepsilon_{B})\Theta(v)+\frac{1}{\pi}\frac{\sqrt{y}\,\Theta(y)}{y+v^{2}}\Bigr{]}\Bigl{[}y\mapsto y+\beta\omega\Bigr{]},$$

Full text

The bulk viscosity is a fundamental transport property which determines friction and dissipation in fluids during hydrodynamic expansion Landau and Lifshitz (1987); Forster (1975). In particular, scale invariant fluids can expand isotropically without dissipation and therefore have zero bulk viscosity Son (2007). In a generic interacting fluid, instead, a nonzero value of the bulk viscosity quantifies the breaking of scale invariance in physical systems ranging from QCD Karsch et al. (2008); Moore and Saremi (2008); Romatschke and Son (2009); Akamatsu et al. (2018) to condensed matter Taylor and Randeria (2010); Enss et al. (2011); Taylor and Randeria (2012); Dusling and Schäfer (2013); Chafin and Schäfer (2013); Elliott et al. (2014); Fujii and Nishida (2018). An intriguing example is the two-dimensional dilute Fermi gas, where the classical model is scale invariant but a quantum scale anomaly breaks this symmetry Pitaevskii and Rosch (1997); Olshanii et al. (2010); Hofmann (2012); Vogt et al. (2012); this has recently been observed via breathing dynamics in cold-atom experiments Holten et al. (2018); Peppler et al. (2018); Murthy et al. (2019).

The bulk viscosity is necessary to understand and predict the real-time evolution and hydrodynamic modes of dissipative quantum fluids and to quantitatively interpret current experiments. However, measurements of the bulk viscosity remain challenging even for classical fluids Dukhin and Goetz (2009). Now a novel experimental probe via the dissipative heating rate due to a change in scattering length has been proposed for atomic gases Fujii and Nishida (2018). It is therefore important to compute the bulk viscosity theoretically for quantum gases, which moreover includes predictions for the classical gas in the high-temperature limit.

The bulk viscosity is defined as the correlation function of local pressure (the trace of the stress tensor). Since it vanishes in a scale invariant system, only the scale breaking part of the pressure contributes, the so-called trace anomaly Martinez and Schäfer (2017); Czajka and Jeon (2017); Fujii and Nishida (2018). This provides a formal link between the breaking of scale invariance and bulk viscosity. The bulk viscosity of the nonrelativistic, strongly interacting Fermi gas has been calculated from kinetic theory in the nondegenerate high-temperature limit Dusling and Schäfer (2013); Chafin and Schäfer (2013) and in the low-temperature superfluid state Escobedo et al. (2009); Hou et al. (2013). Its value is largest in the strongly coupled region of the BEC-BCS crossover Zwerger (2012) near unitarity, but not precisely at unitarity where is must vanish by scale invariance Werner and Castin (2006); Son (2007); Enss et al. (2011). Furthermore, hydrodynamic fluctuations give rise to nonanalytic corrections to the bulk viscosity at small frequencies Onuki (2002); Martinez and Schäfer (2017). However, key open questions include the bulk viscosity in degenerate Fermi gases at strong interaction, the relative importance of bulk and shear viscosity, and critical scaling near the superfluid phase transition.

In this work, I rewrite the bulk viscosity of the dilute Fermi gas as a correlation function of the contact density of local fermion pairs. This exact mapping explicitly links the bulk viscosity to pairing fluctuations as the relevant degrees of freedom and provides a genuine strong-coupling formulation which is valid in the whole BEC-BCS crossover including the quantum critical regime Nikolić and Sachdev (2007); Sachdev (2011); Enss (2012). New results include (a) dominant logarithmic corrections to the bulk viscosity at high temperature, (b) numerical Luttinger-Ward results for the quantum degenerate gas throughout the BEC-BCS crossover predict a large bulk viscosity well observable with current experimental technology, (c) the transport ratio of bulk to shear viscosity deviates from the kinetic theory prediction in the quantum degenerate regime, and (d) critical scaling near the superfluid transition is less singular than predicted Kadanoff and Swift (1968); Onuki (2002), but pairing fluctuations dynamically enhance the scale anomaly.

Bulk viscosity.—The bulk viscosity $\zeta$ is defined as the stress correlation function Kadanoff and Martin (1963); Landau and Lifshitz (1981); Taylor and Randeria (2010)

[TABLE]

where the trace of the stress tensor $\hat{\Pi}_{ii}(\bm{x},t)=d\cdot\hat{P}$ determines the pressure operator $\hat{P}$ in dimension $d$. The two-component dilute Fermi gas is described by the Hamiltonian density Zwerger (2012)

[TABLE]

The first term denotes the kinetic energy with fermion operators $\psi_{\sigma}(\bm{x},t)$. The attractive contact interaction in the second term is characterized by the $s$-wave scattering length $a$. For a given value of $a$, the bare coupling strength $g_{0}$ is determined according to $g_{0}^{-1}=-(m/2\pi)\ln(a\Lambda)$ in two dimensions (2D) and $g_{0}^{-1}=(m/4\pi)(1/a-2\Lambda/\pi)$ in 3D, with ultraviolet momentum cutoff $\Lambda$. The trace of the stress tensor is given by the scale variation of the Hamiltonian Hofmann (2012),

[TABLE]

where the dilatation operator $\hat{D}=\int d^{d}x\,\bm{x}\cdot m\bm{j}(\bm{x})$ generates scale transformations. The first term on the right-hand side is the scale invariant result $[\hat{\mathcal{H}},i\hat{D}]=2\hat{\mathcal{H}}$. If only this is present, the pressure is proportional to the Hamiltonian and commutes with itself in (1), hence the bulk viscosity $\zeta(\omega)\equiv 0$ vanishes identically in the scale invariant case Werner and Castin (2006); Son (2007); Enss et al. (2011); Zwerger (2016).

The second term, in turn, is proportional to the local pair contact density $\hat{\mathcal{C}}=-m^{2}(\partial\hat{\mathcal{H}}/\partial g_{0}^{-1})=m^{2}g_{0}^{2}\psi_{\uparrow}^{\dagger}\psi_{\downarrow}^{\dagger}\psi_{\downarrow}\psi_{\uparrow}$ Tan (2008). Scale invariance is recovered for the ideal quantum gas where $\hat{\mathcal{C}}=0$, and also for the 3D unitary Fermi gas where $1/a=0$ at the scattering resonance. A nonzero bulk viscosity therefore quantifies the breaking of scale invariance, which is generally expected in the interacting Fermi gas, except at unitarity.

Contact correlation.—By conservation of energy, the Hamiltonian in (3) does not contribute to the pressure commutator (1), and the bulk viscosity is given by the correlator of the scale breaking term. The scaling violation in the trace of the stress tensor is the so-called trace anomaly Martinez and Schäfer (2017); Fujii and Nishida (2018)

[TABLE]

where $c_{d}=-(\partial g_{0}^{-1}/\partial\ln\left\lvert a\right\rvert)/m^{2}$ denotes the scale variation of the bare coupling (beta function). For the dilute gas, the equilibrium bulk viscosity is thus exactly given by the contact correlator,

[TABLE]

The contact operator is the term in the Hamiltonian which couples to the scattering length. In linear response, the bulk viscosity thus captures how the local pair contact density at time $t$ changes in response to a variation of the scattering length at earlier time $t=0$ Fujii and Nishida (2018),

[TABLE]

at constant entropy per particle $s=S/N$. The time dependent contact response captures how quickly the contact adjusts to a change in scattering length; this directly determines the dynamical bulk viscosity according to Eq. (5). This makes the contact correlation, and hence the dynamical bulk viscosity, directly accessible in cold atom experiments where the scattering length can be controlled in time by the magnetic field near a Feshbach resonance and the time evolution of the contact has already been measured using RF spectroscopy Bardon et al. (2014); Luciuk et al. (2017).

Viscosity sum rule.—Since the pressure operator is hermitean, the dynamical bulk viscosity is an even and positive function of frequency, $\zeta(\omega)\geq 0$ Taylor and Randeria (2010). The integral over all frequencies in Eqs. (5), (6) immediately yields the bulk viscosity sum rule Taylor and Randeria (2010, 2012) with $\mathcal{C}=\langle\hat{\mathcal{C}}\rangle$,

[TABLE]

Using the Tan adiabatic relation to express the contact $\Pi_{\text{an}}=c_{d}\,\mathcal{C}=(\partial\mathcal{E}/\partial\ln\left\lvert a\right\rvert)_{s}$ as the scale variation of the energy density $\mathcal{E}$ Tan (2008); Hofmann (2012); Werner and Castin (2012), the sum rule is given by the scale “susceptibility” $S=-(1/d^{2})(\partial^{2}\mathcal{E}/\partial(\ln\left\lvert a\right\rvert)^{2})_{s}\geq 0$ in $d$ dimensions. The sum rule is taken at constant entropy per particle to ensure that the bulk viscosity of the ideal gas is zero Taylor and Randeria (2010).

*Pair fluctuations.—*The local contact density can equivalently be interpreted as the density operator $\hat{\mathcal{C}}=\hat{\Delta}^{\dagger}\hat{\Delta}(\bm{x})$ of the local fermion pair field $\hat{\Delta}(\bm{x})=mg_{0}\psi_{\downarrow}(\bm{x})\psi_{\uparrow}(\bm{x})$. The bulk viscosity thus depends directly, and only, on pairing fluctuations within the attractive Fermi gas; it is given exactly by the four-point pair correlation function $\chi(\bm{x},t)=\langle[\hat{\Delta}^{\dagger}\hat{\Delta}(\bm{x},t),\hat{\Delta}^{\dagger}\hat{\Delta}(0,0)]\rangle$. One can anticipate that the bulk viscosity has a strong signature at the superfluid phase transition which is driven by pair fluctuations (see below). While pair fluctuations are strong also at unitarity, the prefactor $c_{d}^{2}\sim 1/a^{2}$ ensures that $\zeta$ vanishes in this case.

To summarize, the bulk viscosity is the response function of the trace anomaly and is therefore sensitive to scaling violation. For the dilute quantum gas, the trace anomaly is proportional to the contact density of local pairs and depends only on the pairing properties. This establishes the link between pairing Murthy et al. (2018) and the quantum scale anomaly Murthy et al. (2019) suggested by recent experiments in 2D Fermi gases.

Analytical results.—The contact correlations and bulk viscosity can be computed exactly in several limiting cases: at (i) zero density (two-body), (ii) high frequency, and (iii) high temperature (virial expansion).

The zero-density case (i) is determined solely by two-body physics. In this limit, the only source of dissipation is the dissociation of a bound molecule at the two-body binding energy $\varepsilon_{B}=\hbar^{2}/ma^{2}$; this yields a high-frequency tail above the threshold $\omega>\varepsilon_{B}$ to break a pair Taylor and Randeria (2012); SM ,

[TABLE]

in 2D, where $\mathcal{C}_{0}$ denotes the two-body contact. In 3D, a two-body bound state exists only on the BEC side for $a>0$, and

[TABLE]

This two-body result serves to disentangle the dissipation due to two-body pair breaking from the genuine many-body bulk viscosity below Taylor and Randeria (2012).

In the limit (ii) of high frequency $\omega\gg\varepsilon_{F},T$, the contact correlator is evaluated at small times where it factorizes as $\chi(\bm{x},t\to 0)\simeq m^{2}\Gamma(\bm{x},t)\mathcal{C}(0,0)$; at large frequency, the pair propagator $\Gamma(\bm{x},t)$ approaches the zero-density form SM . It follows immediately that for $\omega\to\infty$ the bulk viscosity is proportional to the contact density and decays with a characteristic frequency dependence,

[TABLE]

This derivation reproduces earlier results Hofmann (2011); Taylor and Randeria (2012); Goldberger and Khandker (2012) in a dramatically simpler calculation. The zero-density results (8) and (9) approach the high-frequency limit with two-body contact density $\mathcal{C}_{0}$. However, the exact high-frequency limit is more general and holds at arbitrary density, temperature and interaction in terms of the total contact density $\mathcal{C}(n,T,a)$. This asymptotic behavior is important because it guarantees convergence of the sum rule (7).

Finally, the dynamical bulk viscosity can be computed exactly in the high-temperature limit (iii) by virial expansion Dusling and Schäfer (2013); Chafin and Schäfer (2013). To second order in fugacity $z=e^{\beta\mu}$, the pair distribution $b(\varepsilon)=z^{2}e^{-\beta(\varepsilon+2\mu)}$ is combined with the zero-density spectral function to yield SM

[TABLE]

Here, $v=(\lambda/a)/\sqrt{2\pi}$ denotes the dimensionless interaction parameter as the inverse scattering length in units of the thermal length $\lambda=\sqrt{2\pi/mT}$. The dynamical viscosity has two terms as illustrated in Fig. 1: the first, bound-continuum contribution occurs only on the BEC side $v>0$ and arises from breaking up bound states at high frequency $\left\lvert\omega\right\rvert>\varepsilon_{B}$, which leads to strong damping as seen before in the two-body limit (9). The second term is the continuum-continuum contribution of dissociated pairs, which extends over all frequencies but has most of its spectral weight at small frequencies $\omega\lesssim\varepsilon_{B}$. At this order there is no bound-bound contribution because an ideal Bose gas of bound pairs is scale invariant; corrections arise from atom-dimer scattering at order $\mathcal{O}(z^{3})$. Both contributions in (11) are necessary to exhaust the sum rule (cf. Fig. 2)

[TABLE]

This agrees with the adiabatic derivative (7) of the contact Yu et al. (2009); Dusling and Schäfer (2013) \mathcal{C}_{\text{3D,vir}}=16\pi z^{2}\lambda^{-4}\bigl{[}1+\sqrt{\pi}\,ve^{v^{2}}(1+\operatorname{erf}(v))\bigr{]}.

At unitarity $v\to 0$, the analytical dynamical viscosity

[TABLE]

At this order, the unitary contact correlation has a logarithmic singularity $a^{2}\zeta\sim\ln(T/\omega)z^{2}$ for small frequencies from the modified Bessel function $K_{0}(\beta\omega/2)$, as shown in Fig. 1. The logarithmic singularity for small frequencies corresponds via Fourier transform to the logarithmic singularity of the bulk viscosity at long times, $\zeta(t)\sim\ln(t)/(a^{2}t)$ Maki and Zhou (2019). Precisely at unitarity, the bulk viscosity vanishes for all frequencies due to the $v^{2}\propto a^{-2}$ factor. Throughout the BEC-BCS crossover, the dc bulk viscosity is then given by (see Fig. 2)

[TABLE]

The exponential integral $\operatorname{Ei}(x)$ yields a logarithmic singularity in scattering length $a^{2}\zeta\sim\ln(a^{2}/\lambda^{2})z^{2}$ shown in the inset of Fig. 2. The singular coefficient of the virial expansion is regularized by higher-order terms $\mathcal{O}(z^{3})$ from the fermionic self-energy Enss et al. (2011); Dusling and Schäfer (2013); these are resummed in the Luttinger-Ward computation and yield a finite dc limit in Fig. 3(b) below.

In 2D, there is always a bound state with binding energy $\varepsilon_{B}>0$ even for arbitrarily weak attractive interaction. The dynamical bulk viscosity is obtained as SM

[TABLE]

The dc bulk viscosity is then approximately given by

[TABLE]

This result for the bulk viscosity based on contact correlations is similar in structure to the fermionic Boltzmann calculation Chafin and Schäfer (2013) but larger by a factor $4\pi^{2}$, which is necessary to satisfy the sum rule SM and the high-frequency asymptotics with the contact density Ngampruetikorn et al. (2013); Barth and Hofmann (2014) \mathcal{C}_{\text{2D,vir}}=16\pi^{2}z^{2}\lambda^{-4}\bigl{[}\beta\varepsilon_{B}e^{\beta\varepsilon_{B}}+\int_{0}^{\infty}dy\,\frac{e^{-y}}{\ln^{2}(yT/\varepsilon_{B})+\pi^{2}}\bigr{]}.

Luttinger-Ward results.—The Luttinger-Ward (LW) technique is a diagrammatic strong-coupling approach to fermions in the BEC-BCS crossover Haussmann et al. (2007); Bauer et al. (2014) which treats fermions $\psi_{\sigma}$ and the pair field $\Delta$ on equal footing. Its predictions for the unitary shear viscosity Enss et al. (2011) agree well with recent data Bluhm et al. (2017), and similarly for spin diffusion Enss and Haussmann (2012); Enss and Thywissen (2019). In this work, I extend the previous LW approach to compute the bulk viscosity (5) via the contact correlation function (6). It uses the self-consistent pair propagator $\Gamma$ and includes vertex corrections which represent the scattering between pairs, resummed to arbitrary order SM . While contact vertex corrections are subleading in the high-temperature limit and could be neglected, they are crucial in the quantum degenerate regime and need to be included for an accurate numerical solution.

The dynamical bulk viscosity $\zeta(\omega)$ determines the dissipation when the scattering length in Eq. (6) is modulated at frequency $\omega$; the hydrodynamic limit is obtained for $\omega\to 0$. While $\zeta(\omega)$ vanishes at unitarity as $1/a^{2}$, the contact correlations $a^{2}\zeta(\omega)$ are nonzero at unitarity as shown in Fig. 3(a). At low temperature there is a pronounced peak at low frequencies $\omega\lesssim T$ that crosses over into the universal high-frequency tail $\zeta(\omega)\sim\mathcal{C}\omega^{-3/2}$ (dashed) for $\omega\gtrsim\varepsilon_{F}$. At higher $T\gtrsim T_{F}$ the thermal peak for $\omega\lesssim T$ leads directly into the tail. The peak width $\propto T$ is consistent with quantum critical scaling.

The dc bulk viscosity $a^{2}\zeta(T)$ shown in Fig. 3(b) is one of the central results: it is largest near the superfluid transition and decreases toward high temperature where pair fluctuations become weaker, as discussed below.

*Bulk/shear ratio.—*At high temperature kinetic theory predicts the ratio of bulk viscosity $\zeta$ to shear viscosity $\eta$,

[TABLE]

to be proportional to the squared pressure deviation from scale invariance Bluhm et al. (2012); Dusling and Schäfer (2013). Using LW bulk/shear and thermodynamic data Enss et al. (2011), this is tested by comparing $\zeta$ to the kinetic theory prediction $\zeta_{\text{kin}}\equiv\eta[(P-2\mathcal{E}/3)/P]^{2}$, which is shown in Fig. 3(b) as the dashed line. There is very good agreement with unit proportionality factor at high temperature $T\geq T_{F}$, where a quasiparticle picture is expected to hold. Consequently, the shear viscosity at high temperature is fully determined by scale breaking pair fluctuations as reflected in $\zeta$ and in the contact.

In the quantum degenerate regime, the bulk viscosity grows monotonously as the temperature is lowered toward the superfluid phase transition and can reach large values $\zeta\gtrsim\hbar n$ near $T_{c}$. At low temperature, $\zeta>\zeta_{\text{kin}}$ and also $\zeta/\eta>1$ can exceed unity since pair fluctuations near the superfluid phase transition affect the bulk viscosity more strongly than the shear viscosity.

*Critical pair fluctuations.—*The fact that the bulk viscosity is the dynamical correlator of order-parameter fluctuations $\Delta(\bm{x},t)$ suggests that $\zeta$ might diverge at $T_{c}$ Onuki (2002); instead, vertex corrections in the LW calculation substantially reduce the contact vertex at low momenta and render the bulk viscosity large but finite SM . The absence of divergent critical scaling might depend on how the critical point is approached, as found in QCD Martinez et al. (2019).

Finally, Fig. 3(c) shows the viscosity sum rule $S$. It is large in the quantum degenerate regime and decreases toward high temperature as $T^{-3/2}$ (12) (dot-dashed), i.e., faster than the contact $\mathcal{C}\sim T^{-1}$ itself (dashed) Yu et al. (2009); Enss et al. (2011); Mukherjee et al. (2019).

To conclude, the bulk viscosity identifies the breaking of scale invariance with the strength of pair fluctuations, which become very large near $T_{c}$ and on the BEC side. This provides a strong signature in cold atom experiments, either directly in the response of the contact Bardon et al. (2014); Luciuk et al. (2017); Fujii and Nishida (2018) to a change in scattering length, or by modulating the scattering length periodically and measuring the dissipative heating rate $\dot{\mathcal{E}}=d^{2}a^{2}\zeta\cdot(\partial_{t}a^{-1})^{2}$ Fujii and Nishida (2018) proportional to $a^{2}\zeta(\omega,T)$ shown in Fig. 3(a,b), which is nonzero also at unitarity. Further signatures of enhanced dissipation $\zeta$ can be found in the hydrodynamic description of scaling or breathing dynamics Vogt et al. (2012); Taylor and Randeria (2012); Elliott et al. (2014); Murthy et al. (2019) and sound attenuation $D_{s}=[\frac{4}{3}\eta+\zeta+\kappa(c_{v}^{-1}-c_{p}^{-1})]/mn$ Forster (1975); Patel et al. (2019).

Appendix A Supplemental material

A.1 Contact correlations

The contact correlation can be written in terms of the pair field $\Delta(\bm{x},\tau)$ in imaginary time $\tau$ as

[TABLE]

with imaginary time ordering understood. The pair propagator $\langle\Delta(\bm{x},\tau)\Delta^{\dagger}(0,0)\rangle=m^{2}\Gamma(\bm{x},\tau)$ can be expressed in terms of the T matrix $\Gamma(\bm{x},\tau)$, and one can write the contact correlation as the scale variation of the T matrix,

[TABLE]

In imaginary Matsubara frequency $i\omega_{m}$, the spatially integrated contact correlator is given by

[TABLE]

while the contact density itself is given by $\mathcal{C}=(m^{2}/\beta V)\sum_{\bm{q}\varepsilon_{n}}\Gamma(\bm{q},i\varepsilon_{n})$. The vertex corrections are important in the quantum degenerate case and are computed below using the Luttinger-Ward approach. At high temperature or low density, the vertex corrections are subleading and the first term can be computed analytically. After analytical continuation to real frequency, the retarded correlator reads

[TABLE]

in terms of the retarded/advanced pair propagators and the pair spectral function $A(\bm{q},\varepsilon)=-(1/\pi)\operatorname{Im}\Gamma^{R}(\bm{q},\varepsilon)$. This determines the dynamical bulk viscosity as the contact correlator spectral function

[TABLE]

A.2 Zero-density limit

At zero density, the vertex corrections in (19) vanish and the contact correlator is completely determined by the particle-hole excitations of pairs (21). The pair propagator (T matrix), in turn, is given diagrammatically by repeated particle-particle scattering,

[TABLE]

Here, $f(\varepsilon)$ denotes the Fermi distribution and $\xi_{p}=\varepsilon_{p}-\mu$ measures the dispersion $\varepsilon_{p}=p^{2}/2m$ from the chemical potential $\mu$. At zero density there is only a single up and a single down fermion, such that the Fermi functions $f(\xi_{\bm{p}})$ vanish: the momentum integral is then performed analytically and yields the pair spectral functions $A(\bm{q},\varepsilon)=-(1/\pi)\operatorname{Im}\Gamma^{R}(\bm{q},\varepsilon)$ given as

[TABLE]

The first term is the bound-state peak, which appears always for attractive interaction in 2D and for $a>0$ (BEC side) in 3D, followed by the scattering continuum; $\omega_{q}=q^{2}/2M$ denotes the dispersion of fermion pairs (molecular bound states) of mass $M=2m$.

In order to compute the bulk viscosity in the zero-density limit, one has to set $T=0$ and chemical potential $\mu=-\varepsilon_{B}/2$ at the threshold of the two-body binding energy $\varepsilon_{B}=\hbar^{2}/ma^{2}$ ($\mu=0$ when there is no two-body bound state on the BCS side in 3D). At zero density, dissociating a bound state at $\varepsilon=\omega_{q}=0$ in 2D yields Eq. (8) in the main text. This satisfies the sum rule (7), $S_{\text{2D,vac}}=\mathcal{C}_{0}/4\pi m=\varepsilon_{B}/V=-(1/4)(\partial^{2}\mathcal{E}_{0}/\partial(\ln\left\lvert a\right\rvert)^{2})_{s}$ where $\mathcal{E}_{0}=-\mathcal{C}_{0}/4\pi m$. In 3D, a vacuum bound state exists only for $a>0$, and one finds Eq. (9), which again satisfies the sum rule (7) with $S_{\text{3D,vac}}=\mathcal{C}_{0}/36\pi ma=-(1/9)(\partial^{2}\mathcal{E}_{0}/\partial(\ln\left\lvert a\right\rvert)^{2})_{s}$ and $\mathcal{E}_{0}=-\mathcal{C}_{0}/8\pi ma$.

A.3 High-frequency limit

The response at high frequencies is only sensitive to the behavior of the contact correlation at short times, which factorizes in the limit $\tau\to 0$ as

[TABLE]

In this limit, the T matrix $\Gamma(\bm{x},\tau)$ is unaffected by finite density and approaches the zero-density form (23), (24). Therefore, the high-frequency limit $\omega\to\infty$ of the bulk viscosity is proportional to the contact density and decays with the asymptotic frequency dependence (10) quoted in the main text for $\omega\gg\varepsilon_{B}$.

A.4 Virial expansion

The virial expansion to second order $\mathcal{O}(z^{2})$ in the fermionic fugacity $z=e^{\beta\mu}$ correctly describes the properties of the interacting Fermi gas as long as the pair fugacity $z_{\text{pair}}=z^{2}e^{\beta\varepsilon_{B}}\ll 1$ remains small. Because vertex corrections in Eq. (21) appear only at higher order in $z$, the second-order virial result is fully determined by the first term in that equation. Since the occupation factor $b(\varepsilon)$ is already order $z^{2}$, it suffices to use the zero-density form of the pair spectral function (24) and the Boltzmann distribution to obtain

[TABLE]

With the explicit form of $A_{\text{3D}}$ (24) and in the new frequency variable $y=\beta(\varepsilon+2\mu-\omega_{q})$ one finds

[TABLE]

where $v=(\lambda/a\sqrt{2\pi})$ denotes the dimensionless 3D interaction parameter. Since $\zeta(\omega)$ is even in frequency it suffices to consider $\omega>0$; then there is the bound-continuum contribution from the bound state $y=-\beta\varepsilon_{B}$ in the first bracket and the continuum $y+\beta\omega>0$ in the second. On the other hand, for $y>0$ one obtains the continuum-continuum contribution, and both terms are combined to yield

[TABLE]

Note that there appears no bound-bound contribution $\delta(\omega)$ because an ideal Bose gas of pairs is scale invariant. Multiplication with the coefficient $(c_{d}/d)^{2}=T^{2}\lambda^{2}v^{2}/(288\pi^{3})$ in 3D yields the dynamical bulk viscosity (11).

Frequency integration determines the spectral weight of the terms in brackets as $S_{\text{bc}}=[2(1+2v^{2})e^{v^{2}}\operatorname{erf}(v)-4v^{2}e^{v^{2}}+4v/\sqrt{\pi}]\Theta(v)$ for the bound-continuum contribution, which becomes large in the BEC limit. Furthermore, the continuum-continuum contribution has weight $S_{\text{cc}}=[(1+2v^{2})e^{v^{2}}(1-\operatorname{erf}(\left\lvert v\right\rvert))-2\left\lvert v\right\rvert/\sqrt{\pi}]$, and by combining both one exhausts the sum rule (12).

The contact density is given to the same order in the virial expansion by the occupied spectral function of pairs,

[TABLE]

This is equivalent to the derivative of the second virial coefficient, $\mathcal{C}=8\pi\sqrt{2\pi}z^{2}\lambda^{-4}\partial b_{2}(v)/\partial v$ with $b_{2}(v)=-1/4\sqrt{2}+(1/\sqrt{2})(1+\operatorname{erf}v)e^{v^{2}}$, $b_{2}^{\prime}=\sqrt{2/\pi}+2vb_{2}$, and $b_{2}^{\prime\prime}=2(b_{2}+vb_{2}^{\prime})$. The contact grows monotonously with the interaction parameter $v$, as shown in Fig. 2; this generalizes earlier results Yu et al. (2009); Dusling and Schäfer (2013).

In order to compute the adiabatic derivative of the contact, one has to keep the entropy per particle $s=\beta(h-\mu)$ fixed. This is given in terms of the enthalpy per particle $h=(p+\mathcal{E})/n$ and the chemical potential $\mu$, and at second order virial expansion one finds

[TABLE]

in extension of the Sackur-Tetrode entropy formula $s=5/2-\ln(z)$. The adiabatic derivative with respect to $v$ is related to the grand canonical derivative keeping $\mu$ and $T$ fixed, plus an additional term adjusting $z^{\prime}/z$ to compensate for the change in entropy: $(\partial z^{2}b_{2}^{\prime}/\partial v)_{s}=(2\frac{z^{\prime}}{z}b_{2}^{\prime}+b_{2}^{\prime\prime})z^{2}$. In the BEC limit, $z^{\prime}/z=-v$ and we find the adiabatic derivative $(\partial z^{2}b_{2}^{\prime}/\partial v)_{s}=2z^{2}b_{2}$. The corresponding two-body contact is $\mathcal{C}_{0}=4\pi n/a=16\pi z^{2}\lambda^{-4}\sqrt{\pi}2ve^{v^{2}}$, and the adiabatic derivative $(\partial\mathcal{C}_{0}/\partial v)_{s}=\mathcal{C}_{0}/v=16\pi z^{2}\lambda^{-4}\sqrt{\pi}2e^{v^{2}}$ yields the bound-continuum contribution to the sum rule (12) without a bound-bound contribution. Similarly, the adiabatic derivative of the continuum part of the contact at fixed fermionic entropy $s=5/2-\ln(z)$ yields the continuum-continuum contribution to the sum rule (12), thus confirming Eq. (7) by explicit computation in the high-temperature limit.

The dynamical bulk viscosity in two dimensions is computed in a completely analogous fashion to obtain the bulk viscosity quoted in the main text. Its total spectral weight is given by the sum rule (7) with

[TABLE]

The sum rule agrees with the adiabatic derivative of the contact density Ngampruetikorn et al. (2013); Barth and Hofmann (2014)

[TABLE]

Appendix B Luttinger-Ward approach

The Luttinger-Ward approach Haussmann et al. (2007) to the attractive Fermi gas is based on two-component fermions $\psi_{\sigma}(\bm{x},\tau)$ which interact by forming (virtual) pairs $\Delta(\bm{x},\tau)$. It is constructed to be a conserving approximation which exactly conserves not only particle number and momentum current but also the dilatation current (scale invariance) Enss et al. (2011) and the Tan relations Enss (2012). The pair propagator is given by the T matrix (cf. Eq. (22))

[TABLE]

with bare coupling $g_{0}$, and the pair self-energy is given by the fermion particle-particle bubble $M(\bm{q},i\epsilon_{n})$. This is a convolution of two dressed fermion propagators $G(\bm{q},i\epsilon_{n})$ in momentum $\bm{q}$ and Matsubara frequency $i\epsilon_{n}$. However, numerically it is more conveniently computed by a Fourier transform to real space $\bm{r}$ and imaginary time $\tau$, where the particle-particle bubble is a local product

[TABLE]

the resulting $M(\bm{r},\tau)$ is then Fourier transformed back to momentum and frequency $M(\bm{q},i\epsilon_{n})$ to be inserted in the pair Dyson equation (25). In turn, the fermion propagator for each spin component of the balanced gas is given by the Dyson equation

[TABLE]

with the fermionic self-energy $\Sigma(\bm{q},i\epsilon_{n})$ determined by scattering fermions off virtual pairs. The particle-hole bubble of a pair (particle) and a fermion (hole) is again concisely written as a local multiplication in Fourier space,

[TABLE]

These four equations (25)–(28) form a closed set of equations, which is solved self-consistently by iteration until convergence is reached. The propagators $\Gamma(\bm{q},i\epsilon_{n})$ and $G(\bm{q},i\epsilon_{n})$ are first initialized with the bare propagators, then Fourier transformed to real space. The self-energies (26) and (28) are computed in real space, then Fourier transformed back to momentum space. These are then inserted into the Dyson equations (25) and (27), which provide the starting point for the next iteration.

The Luttinger-Ward approach has previously been used for fermionic shear and spin transport Enss et al. (2011); Enss and Haussmann (2012); by considering the response to variations of an external field (shear strain or spin gradient), a new set of self-consistent transport equations is obtained for the renormalized current vertices which include vertex corrections to satisfy the Ward identities exactly. For the bosonic contact correlations, instead, one has to consider the response to a time-dependent scale variation. Specifically, the scale breaking variation $\delta=-\partial/\partial(4\pi ma)^{-1}$ at external drive frequency $i\omega_{m}$ of each of the above equations (25)–(28) yields a new set of self-consistent transport equations for the renormalized trace anomaly response vertices:

[TABLE]

In the first line (29), $\delta\Gamma^{-1}$ denotes the interaction scale variation of the bosonic pair propagator: it consists of the bare contact vertex $\delta g_{0}^{-1}=-m^{2}$, which is obtained from the scale variation $c_{d}$ of the bare coupling, and the contact vertex correction $\delta M$, which arises from the scale variation of the particle-particle bubble (30). While the bare fermionic propagator has no scale dependence, $\delta G_{0}^{-1}=0$, the dressed fermion propagator (27) acquires a scale dependence from the self-energy term (28). This gives rise to a fermionic trace anomaly vertex $\delta G^{-1}$ in (B) with two distinct types of vertex corrections. The first, so-called Maki-Thompson vertex correction $\delta\Sigma_{\text{MT}}$ in (32) is built from a fermionic anomaly vertex $\delta G$, while the second, Aslamazov-Larkin vertex corrections $\delta\Sigma_{\text{AL}}$ in (33) arises from the renormalized bosonic contact vertex $\delta\Gamma$. The scale variation of the propagators is computed from the inverse propagators as

[TABLE]

The self-consistent transport equations are solved separately for each value of the external frequency $i\omega_{m}$; as before, a Fourier transform converts the anomaly vertices between Fourier space $(\bm{q},i\epsilon_{n})$ and real space $(\bm{r},\tau)$ at fixed parameter $i\omega_{m}$.

In the first iteration of the transport equations, only the bare pair propagator depends on scale and yields the bare contact vertex $\delta g_{0}^{-1}=-m^{2}$, while the bare fermions are independent of the interaction scale, $\delta G_{0}^{-1}(\bm{q},i\epsilon_{n};i\omega_{m})=0$. In subsequent iterations, the fermionic anomaly vertex $\delta G^{-1}$ picks up an interaction scale dependence from the MT and AL vertex corrections (32) and (33). The scale dependence of the fermionic vertex, in turn, generates vertex corrections $\delta M$ which renormalize the contact vertex $\delta\Gamma^{-1}$. Once convergence is reached, the Maki-Thompson and Aslamazov-Larkin vertex corrections include contributions of arbitrarily high perturbative order in the bare coupling.

The spatially integrated contact correlation function (19) is computed by the Kubo formula

[TABLE]

for each value of the external Matsubara frequency $i\omega_{m}$. Finally, analytical continuation to real frequencies $\omega+i0$ is performed using Padé approximants for the first few tens of Matsubara frequencies to obtain the retarded contact correlator $X^{R}(\omega)$ in (20), which yields the bulk viscosity via Eq. (21).

At high temperature, the vertex corrections are subleading and it suffices to use the bare contact vertex $\delta\Gamma^{-1}=\delta g_{0}^{-1}$, such that (34) simplifies to

[TABLE]

in agreement with the first term in Eq. (19). Hence, the virial limit is contained within the LW approach.

In the BEC limit $a\to 0^{+}$, the pair propagator has a large gap of order $\varepsilon_{B}\gg\varepsilon_{F},T$ between the bound state branch and the continuum of dissociated pairs. However, the bound state branch is no longer a $\delta$ peak as in Eq. (24) but is broadened due to atom-dimer and dimer-dimer scattering contained within the dressed propagators and the vertex corrections in the LW equations. This broadening gives rise to a finite bound-bound contribution to the bulk viscosity, which should approach the dc bulk viscosity of a weakly repulsive Bose gas as $a\to 0^{+}$.

Appendix C Importance of contact vertex corrections

While in the high-temperature limit the fully dressed contact vertex $\delta\Gamma^{-1}(\bm{q},i\epsilon_{n};i\omega_{m})$ is close to the bare contact vertex $\delta g_{0}^{-1}=-m^{2}$, it is found to be substantially renormalized in the quantum degenerate regime approaching $T_{c}$. In particular, at the superfluid phase transition the pair propagator $\Gamma\sim 1/(\xi^{-2}+q^{2})$ becomes gapless and gives rise to divergent critical fluctuations as the correlation length $\xi\to\infty$ diverges. One might expect that the contact correlations (34) would also diverge with a positive power of $\xi$. However, the contact vertex corrections $\delta M$ which are included within the Luttinger-Ward approach strongly suppress the static contact vertex at small momenta and result in a scaling form $\delta\Gamma^{-1}(\bm{q},0;0)\sim q^{2}$, which renders the contact correlation less singular when approaching $T_{c}$. This is illustrated in Fig. 4, which shows the fully renormalized contact vertex in units of the bare contact vertex $\delta g_{0}^{-1}$: at large momenta $q\gg k_{F}$ it remains unrenormalized, but is strongly suppressed as $\sim q^{2}$ for small momenta $\xi^{-1}\ll q\ll k_{F}$ before it eventually saturates for $q\ll\xi^{-1}$, where $\xi$ depends on the distance from $T_{c}$.