Viscosity spectral functions of resonating fermions in the quantum virial expansion

TL;DR

This paper calculates the spectral functions of bulk and shear viscosities for resonating fermions in two and three dimensions using the quantum virial expansion, revealing new insights into their high-temperature behavior and non-analyticities.

Contribution

It provides the first systematic evaluation of viscosity spectral functions in the quantum virial expansion framework, including a novel expression for the bulk viscosity via contact response functions.

Findings

Bulk viscosity spectral function expressed with contact-contact response.

Static shear viscosity matches kinetic theory; bulk viscosity shows unexpected non-analyticity.

Discrepancy in static bulk viscosity suggests challenges to the crossover hypothesis.

Abstract

We consider two-component fermions with a zero-range interaction both in two and three dimensions and study their spectral functions of bulk and shear viscosities for an arbitrary scattering length. Here the Kubo formulas are systematically evaluated up to the second order in the quantum virial expansion applicable to the high-temperature regime. In particular, our computation of the bulk viscosity spectral function is facilitated by expressing it with the contact-contact response function, which can be measured experimentally under the periodic modulation of the scattering length. The obtained formulas are fully consistent with the known constraints on high-frequency tail and sum rule. Although our static shear viscosity agrees with that derived from the kinetic theory, our static bulk viscosity disagrees. Furthermore, the latter for three dimensions exhibits an unexpected non-analyticity of $\zeta\sim(\ln a^2)/a^2$ in the unitarity limit $a\to\infty$, which thus challenges the "crossover" hypothesis.

Equations

$$\hat{H}=-\sum _{\sigma =\uparrow ,\downarrow }\int d\mathbf{x}{\hat{\psi }}_{\sigma }^{†}(\mathbf{x})\frac{{\nabla }^2}{2m}{\hat{\psi }}_{\sigma }(\mathbf{x})+\sum _{\sigma ,\tau }\int d\mathbf{x}{\hat{\psi }}_{\sigma }^{†}(\mathbf{x}){\hat{\psi }}_{\tau }^{†}(\mathbf{x})\frac{g}{2}{\hat{\psi }}_{\tau }(\mathbf{x}){\hat{\psi }}_{\sigma }(\mathbf{x}),$$

$$Z=\mathrm{Tr}[e^{-\beta (\hat{H}-\mu \hat{N})}]=\sum _{N=0}^{\infty }{z}^N\underset{Z_N}{\underbrace{{\mathrm{tr}}_N[e^{-\beta \hat{H}}]}},$$

$$|{\sigma }_1{\mathbf{p}}_1,\dots ,{\sigma }_N{\mathbf{p}}_N⟩=\frac{1}{\sqrt{N!}}\prod _{n=1}^N{\hat{\psi }}_{{\sigma }_n{\mathbf{p}}_n}^{†}|0⟩,$$

$${\mathrm{tr}}_N[\hat{\mathcal{O}}]=\sum _{\{\sigma \}}\sum _{\{\mathbf{p}\}}⟨{\sigma }_1{\mathbf{p}}_1,\dots ,{\sigma }_N{\mathbf{p}}_N|\hat{\mathcal{O}}|{\sigma }_1{\mathbf{p}}_1,\dots ,{\sigma }_N{\mathbf{p}}_N⟩.$$

$${\text{openone}}_N=\sum _{\{\sigma \}}\sum _{\{\mathbf{p}\}}|{\sigma }_1{\mathbf{p}}_1,\dots ,{\sigma }_N{\mathbf{p}}_N⟩⟨{\sigma }_1{\mathbf{p}}_1,\dots ,{\sigma }_N{\mathbf{p}}_N|,$$

$$\hat{H}=\sum _{N=0}^{\infty }\underset{{\hat{H}}_N}{\underbrace{{\text{openone}}_N\hat{H}{\text{openone}}_N}}.$$

$${\hat{H}}_1=\sum _{\sigma }\sum _{\mathbf{p}}|\sigma \mathbf{p}⟩\frac{{\mathbf{p}}^2}{2m}⟨\sigma \mathbf{p}|$$

$$Z_1={\mathrm{tr}}_1[e^{-\beta {\hat{H}}_1}]=2\sum _{\mathbf{p}}{e}^{-\beta \frac{{\mathbf{p}}^2}{2m}}.$$

$${\hat{H}}_2=\underset{{\hat{T}}_2}{\underbrace{\sum _{\sigma ,\tau }\sum _{\mathbf{p},\mathbf{q}}|\sigma \mathbf{p},\tau \mathbf{q}⟩\frac{{\mathbf{p}}^2+{\mathbf{q}}^2}{2m}⟨\sigma \mathbf{p},\tau \mathbf{q}|}}+\underset{{\hat{V}}_2}{\underbrace{\sum _{\sigma ,\tau }\sum _{\mathbf{p},\mathbf{q}}\sum _{{\mathbf{p}}^{\prime },{\mathbf{q}}^{\prime }}|\sigma \mathbf{p},\tau \mathbf{q}⟩{\delta }_{\mathbf{p}+\mathbf{q},{\mathbf{p}}^{\prime }+{\mathbf{q}}^{\prime }}\frac{g}{L^d}⟨\sigma {\mathbf{p}}^{\prime },\tau {\mathbf{q}}^{\prime }|}}.$$

$$e^{-\beta {\hat{H}}_2}$$

$$={\int }_{-\infty }^{\infty }\frac{dE}{\pi }{e}^{-\beta E}\mathrm{Im}\left[\frac{1}{E-{\hat{T}}_2-i{0}^{+}}+\sum _{n=0}^{\infty }\frac{1}{E-{\hat{T}}_2-i{0}^{+}}{\hat{V}}_2{\left(\frac{1}{E-{\hat{T}}_2-i{0}^{+}}{\hat{V}}_2\right)}^n\frac{1}{E-{\hat{T}}_2-i{0}^{+}}\right].$$

$$\sum _{n=0}^{\infty }{\hat{V}}_2{\left(\frac{1}{E-{\hat{T}}_2-i{0}^{+}}{\hat{V}}_2\right)}^n$$

$$=\frac{1}{L^d}\sum _{\sigma ,\tau }\sum _{\mathbf{p},\mathbf{q}}\sum _{{\mathbf{p}}^{\prime },{\mathbf{q}}^{\prime }}|\sigma \mathbf{p},\tau \mathbf{q}⟩{\delta }_{\mathbf{p}+\mathbf{q},{\mathbf{p}}^{\prime }+{\mathbf{q}}^{\prime }}\frac{1}{\frac{1}{g}-\frac{1}{L^d}{\sum }_{\mathbf{k}}\frac{1}{E-\frac{(\mathbf{p}+\mathbf{q}{)}^2}{4m}-\frac{{\mathbf{k}}^2}{m}-i{0}^{+}}}⟨\sigma {\mathbf{p}}^{\prime },\tau {\mathbf{q}}^{\prime }|,$$

$$\displaystyle e^{-\beta\hat{H}_{2}}=e^{-\beta\hat{H}_{2}}\big{|}_{\mathrm{free}}+e^{-\beta\hat{H}_{2}}\big{|}_{\mathrm{int}},$$

$$\displaystyle e^{-\beta\hat{H}_{2}}\big{|}_{\mathrm{free}}=\sum_{\sigma,\tau}\sum_{{\bm{p}},{\bm{q}}}|\sigma{\bm{p}},\tau{\bm{q}}\rangle e^{-\beta\frac{{\bm{p}}^{2}+{\bm{q}}^{2}}{2m}}\langle\sigma{\bm{p}},\tau{\bm{q}}|$$

$$\displaystyle e^{-\beta\hat{H}_{2}}\big{|}_{\mathrm{int}}$$

$$\times {\int }_{-\infty }^{\infty }\frac{d\epsilon }{\pi }{e}^{-\beta \epsilon }\mathrm{Im}\left[\frac{{\mathcal{T}}_2(\epsilon -i{0}^{+})}{(\epsilon -\frac{(\mathbf{p}-\mathbf{q}{)}^2}{4m}-i{0}^{+})(\epsilon -\frac{({\mathbf{p}}^{\prime }-{\mathbf{q}}^{\prime }{)}^2}{4m}-i{0}^{+})}\right]⟨\sigma {\mathbf{p}}^{\prime },\tau {\mathbf{q}}^{\prime }|,$$

$${\mathcal{T}}_2(\epsilon -i{0}^{+})\equiv \frac{1}{\frac{1}{g}-\frac{1}{L^d}{\sum }_{\mathbf{k}}\frac{1}{\epsilon -\frac{{\mathbf{k}}^2}{m}-i{0}^{+}}}$$

$$Z_2={\mathrm{tr}}_2[e^{-\beta {\hat{H}}_2}]=\underset{Z_1^2/2}{\underbrace{2\sum _{\mathbf{p},\mathbf{q}}{e}^{-\beta \frac{{\mathbf{p}}^2+{\mathbf{q}}^2}{2m}}}}-\sum _{\mathbf{p},\mathbf{q}}{e}^{-\beta \frac{{\mathbf{p}}^2+{\mathbf{q}}^2}{2m}}{\delta }_{\mathbf{p}\mathbf{q}}+\frac{1}{L^d}\sum _{\mathbf{K},\mathbf{k}}{e}^{-\beta \frac{{\mathbf{K}}^2}{4m}}{\int }_{-\infty }^{\infty }\frac{d\epsilon }{\pi }{e}^{-\beta \epsilon }\mathrm{Im}\left[\frac{{\mathcal{T}}_2(\epsilon -i{0}^{+})}{(\epsilon -\frac{{\mathbf{k}}^2}{m}-i{0}^{+}{)}^2}\right],$$

$$\frac{p}{T}=\frac{\ln Z}{L^d}$$

$$=\frac{2z}{L^d}\sum _{\mathbf{p}}{e}^{-\beta \frac{{\mathbf{p}}^2}{2m}}-\frac{z^2}{L^d}\sum _{\mathbf{p}}{e}^{-\beta \frac{{\mathbf{p}}^2}{m}}+\frac{z^2}{L^{2d}}\sum _{\mathbf{K},\mathbf{k}}{e}^{-\beta \frac{{\mathbf{K}}^2}{4m}}{\int }_{-\infty }^{\infty }\frac{d\epsilon }{\pi }{e}^{-\beta \epsilon }\mathrm{Im}\left[\frac{{\mathcal{T}}_2(\epsilon -i{0}^{+})}{(\epsilon -\frac{{\mathbf{k}}^2}{m}-i{0}^{+}{)}^2}\right]+O(z^3).$$

$$\frac{p}{T}=\frac{2z}{{\lambda }_T^d}-\frac{z^2}{2^{d/2}{\lambda }_T^d}+\frac{2^{d/2}{z}^2}{2{\Omega }_{d-1}{\lambda }_T^d}{\int }_{-\infty }^{\infty }\frac{d\epsilon }{\pi }{e}^{-\epsilon /T}\mathrm{Im}\left[\frac{m^2{\mathcal{T}}_2(\epsilon -i{0}^{+})}{(-m\epsilon +i{0}^{+}{)}^{2-d/2}}\right]+O(z^3),$$

$${\mathcal{T}}_2(\epsilon -i{0}^{+})=\frac{{\Omega }_{d-1}}{m}\frac{d-2}{a^{2-d}-(-m\epsilon +i{0}^{+}{)}^{d/2-1}},$$

$$\frac{1}{g}=\frac{m{a}^{2-d}}{(d-2){\Omega }_{d-1}}$$

$$\mathcal{C}=-{\Omega }_{d-1}m{a}^{d-1}\frac{\partial p}{\partial a}=\frac{2^{d/2}{z}^2}{{\lambda }_T^d}{\int }_{-\infty }^{\infty }\frac{d\epsilon }{\pi }{e}^{-\epsilon /T}\mathrm{Im}[m^2{\mathcal{T}}_2(\epsilon -i{0}^{+})]+O(z^3),$$

$$a^{d-1}\frac{\partial {\mathcal{T}}_2(\epsilon -i{0}^{+})}{\partial a}=-\frac{2(-m\epsilon +i{0}^{+}{)}^{2-d/2}}{m}\frac{\partial {\mathcal{T}}_2(\epsilon -i{0}^{+})}{\partial \epsilon }$$

$$\mathrm{Im}[{\mathcal{T}}_2(\epsilon -i{0}^{+})]=\theta (a)\frac{2{\Omega }_{d-1}}{m^2{a}^{4-d}}\pi \delta (\epsilon +\frac{1}{m{a}^2})+\theta (\epsilon )\times \{\begin{array}{ll}\frac{2\pi }{m}\frac{2\pi }{{\ln }^2(m{a}^2\epsilon )+{\pi }^2}& (d=2),\text{vskip}6.0ptplus2.0ptminus2.0pt\\ \frac{4\pi }{m}\frac{\sqrt{m\epsilon }}{m\epsilon +\frac{1}{a^2}}& (d=3),\end{array}$$

$$\mathcal{C}=\frac{16\pi z^2}{{\lambda }_T^4}{\left[1+\sqrt{\pi }\frac{1+\mathrm{erf}({\tilde{a}}^{-1})}{\tilde{a}}\exp ({\tilde{a}}^{-2})\right]}_{\tilde{a}=\sqrt{2\pi }a/{\lambda }_T}+O(z^3).$$

$${\chi }_{ij,kl}(\omega )=\frac{i}{Z{L}^d}{\int }_0^{\infty }dt{e}^{i(\omega +i{0}^{+})t}\mathrm{Tr}[[{\hat{\Pi }}_{ij}(t),{\hat{\Pi }}_{kl}(0)]e^{-\beta (\hat{H}-\mu \hat{N})}],$$

$$\frac{\mathrm{Im}[{\chi }_{ij,kl}(\omega )]}{\omega }=\left[\zeta (\omega )-\frac{2}{d}\eta (\omega )\right]{\delta }_{ij}{\delta }_{kl}+\eta (\omega )\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right).$$

$$\zeta (\omega )=\frac{1}{d^2}\sum _{i,j=1}^d\frac{\mathrm{Im}[{\chi }_{ii,jj}(\omega )]}{\omega }.$$

$${\hat{\Pi }}_{ij}=-\sum _{\sigma =\uparrow ,\downarrow }\int d\mathbf{x}{\hat{\psi }}_{\sigma }^{†}(\mathbf{x})\frac{{\partial }_i{\partial }_j}{m}{\hat{\psi }}_{\sigma }(\mathbf{x})+\sum _{\sigma ,\tau }\int d\mathbf{x}{\hat{\psi }}_{\sigma }^{†}(\mathbf{x}){\hat{\psi }}_{\tau }^{†}(\mathbf{x})\frac{g}{2}{\delta }_{ij}{\hat{\psi }}_{\tau }(\mathbf{x}){\hat{\psi }}_{\sigma }(\mathbf{x}).$$

$$\sum _{i=1}^d{\hat{\Pi }}_{ii}=2\hat{H}+\frac{\hat{C}}{{\Omega }_{d-1}m{a}^{d-2}},$$

Full text

1. I Introduction, summary, and conclusion

2. II Preliminaries

3. II.1 Formulation

4. II.2 Thermodynamics

5. III Bulk viscosity

6. III.1 Formulation

7. III.2 Evaluation

8. IV Shear viscosity

9. IV.1 Formulation

10. IV.2 Evaluation

11. A Bulk viscosity formula at $O(z^{3})$

12. B Derivation of Eq. ()

I Introduction, summary, and conclusion

Transport coefficients provide valuable insights into strongly correlated systems complementary to thermodynamic quantities. In particular, the shear viscosity has attracted significant interest across diverse subfields in physics Schafer:2009 ; Adams:2012 , motivated by the universal lower bound on the shear viscosity to entropy density ratio Kovtun:2005 . Although the proposed conjecture was eventually violated Brigante:2008a ; Brigante:2008b ; Kats:2009 ; Buchel:2009 , the quest for perfect fluids in strongly correlated quantum systems keeps on going.

Ultracold atoms offer ideal grounds to investigate strongly correlated quantum fluids because atomic interactions are tunable via Feshbach resonances Chin:2010 . Here the shear viscosity was measured for a strongly correlated Fermi gas of two components whose interaction saturates the unitarity bound Cao:2011a ; Cao:2011b . This measurement was later extended to the BCS-BEC crossover and the shear viscosity was found to be lower on the BEC side Elliott:2014b . On the other hand, the bulk viscosity in the unitarity limit vanishes identically because of the conformal invariance Son:2007 , which was also confirmed experimentally Elliott:2014a . Furthermore, the bulk and shear viscosities were measured in the BCS-BEC crossover for a two-dimensional Fermi gas Vogt:2012 .

The purpose of this paper is to present comprehensive analyses of viscosity spectral functions with full frequency dependence throughout the BCS-BEC crossover both in two and three dimensions. We consider two-component fermions with a zero-range interaction parametrized only by the scattering length and evaluate the Kubo formulas systematically up to the second order in the quantum virial expansion, which is applicable to the high-temperature regime where the fugacity diminishes as $z\sim n/T^{d/2}$ at a fixed density. To this end, the quantum virial expansion is reviewed in Sec. II for thermodynamic quantities, which is then applied to the spectral functions of bulk and shear viscosities in Secs. III and IV, respectively. In particular, our computation of the bulk viscosity spectral function is facilitated by expressing it with the contact–contact response function, which allows us to derive its formula even at $O(z^{3})$ as shown in Appendix A. We set $\hbar=k_{B}=1$ throughout this paper.

Because our analyses below appear to be involved, we first summarize our results and conclusions to guide readers who may not be interested in the intermediate technical details. The spectral functions of bulk and shear viscosities up to the second order in fugacity are presented in Eqs. (40) and (IV.1), respectively, which are valid for an arbitrary scattering length both in two and three dimensions. We confirm in Eqs. (41), (44) for bulk and in Eqs. (71), (74) for shear that the obtained formulas are fully consistent with the known constraints on high-frequency tail and sum rule Taylor:2010 ; Enss:2011 ; Hofmann:2011 ; Goldberger:2012 ; Taylor:2012 , which involves the contact density presented in Eq. (22).

The static limit of the shear viscosity spectral function requires the resummation, which leads to the static shear viscosity presented in Eq. (77) and agrees with that derived from the kinetic theory Massignan:2005 ; Bruun:2005 ; Bruun:2012 ; Schafer:2012 . On the other hand, the static limit of the bulk viscosity spectral function is presented in Eq. (47) and turns out to disagree with that derived from the kinetic theory Dusling:2013 ; Chafin:2013 . Furthermore, our static bulk viscosity for three dimensions exhibits an unexpected non-analyticity of $\zeta\sim(\ln a^{2})/a^{2}$ in the unitarity limit $a\to\infty$, which thus challenges the “crossover” hypothesis Zwerger:2012 . We also find that the static shear viscosity for three dimensions exhibits a weaker non-analyticity starting with $\eta\sim(\ln a^{2})/a^{6}$ in the same limit.

The origin of these discrepancies between the Kubo formalism and the kinetic theory for the bulk viscosity is currently unknown and needs to be elucidated in a future study. To this end, it may be useful to investigate higher order corrections to the bulk viscosity spectral function, which at the third order in fugacity is presented in Eq. (A). Also, as shown in Eqs. (32), (52), and (53), the bulk viscosity spectral function can be extracted experimentally by measuring the contact, energy, or entropy density under the periodic modulation of the scattering length. Therefore, ultracold atom experiments may be qualified to discriminate the two contradicting predictions from the Kubo formalism and the kinetic theory.

II Preliminaries

II.1 Formulation

In this paper, we consider two-component fermions with a zero-range interaction in $d$ spatial dimensions described by

[TABLE]

where the annihilation and creation operators obey the usual anti-commutation relation of $\{\hat{\psi}_{\sigma}({\bm{x}}),\hat{\psi}_{\tau}^{\dagger}({\bm{y}})\}=\delta_{\sigma\tau}\delta({\bm{x}}-{\bm{y}})$. Before we work on the quantum virial expansion for viscosity spectral functions, let us warm up by reviewing that for the partition function.

The quantum virial expansion is a systematic expansion in terms of fugacity, which utilizes a simple fact that the grand canonical trace can be written as a sum of canonical traces Liu:2013 . In particular, the partition function reads

[TABLE]

where $z\equiv e^{\beta\mu}$ is the fugacity and obviously $Z_{0}=1$. It is convenient to choose the complete set of $N$-body states to be momentum eigenstates,

[TABLE]

where $\hat{\psi}_{\sigma{\bm{p}}}^{\dagger}\equiv L^{-d/2}\int\!d{\bm{x}}\,e^{i{\bm{p}}\cdot{\bm{x}}}\,\hat{\psi}_{\sigma}^{\dagger}({\bm{x}})$ is the Fourier transform of the creation operator in a periodic box of volume $L^{d}$. The canonical trace in the $N$-body sector is then provided by

[TABLE]

Also, with the use of the identity operator,

[TABLE]

the Hamiltonian (1) can be decomposed into distinct particle number sectors according to

[TABLE]

The one-body sector is trivial because the Hamiltonian is

[TABLE]

and does not involve the interaction. Consequently, the partition function at $O(z)$ is simply provided by

[TABLE]

On the other hand, the Hamiltonian in the two-body sector is

[TABLE]

In order to evaluate the corresponding Boltzmann operator, we write the Hamiltonian as a sum of kinetic and potential energy operators and utilize the identity of

[TABLE]

Because the summation over $n$ can be performed according to

[TABLE]

the Boltzmann operator is found to be

[TABLE]

respectively. Here the integration variable is changed to $\varepsilon=E-({\bm{p}}+{\bm{q}})^{2}/4m$ and

[TABLE]

is the two-body scattering $T$-matrix in the center-of-mass frame. We note that Eq. (14) is valid for any complex $\beta$ and the interacting part (14) vanishes when $\sigma=\tau$ for fermions. Consequently, the partition function at $O(z^{2})$ is provided by

[TABLE]

where ${\bm{K}}={\bm{p}}+{\bm{q}}$ and ${\bm{k}}=({\bm{p}}-{\bm{q}})/2$ are the center-of-mass and relative momenta, respectively.

II.2 Thermodynamics

From the above quantum virial expansion for the partition function, we obtain that for the pressure according to

[TABLE]

Finally, by taking the thermodynamic limit, $L^{-d}\sum_{\bm{p}}\to\int\!d{\bm{p}}/(2\pi)^{d}$, and performing the resulting momentum integrations, we find

[TABLE]

where $\lambda_{T}\equiv\sqrt{2\pi/mT}$ is the thermal de Broglie wavelength and $\Omega_{d-1}\equiv(4\pi)^{d/2}/2\Gamma(2-d/2)=2,2\pi,4\pi$ coincides with the surface area of the unit $(d-1)$-sphere for $d=1,2,3$. Similarly, the two-body scattering $T$-matrix in the same limit is found to be

[TABLE]

where $a$ is the scattering length introduced via

[TABLE]

in the dimensional regularization Fujii:2018 .

The partial derivative of the pressure with respect to the scattering length then leads to the so-called contact density Tan:2008a ; Tan:2008b ; Tan:2008c ,

[TABLE]

where the integration by parts following

[TABLE]

is utilized. Because the imaginary part of the two-body scattering $T$-matrix is provided by

[TABLE]

the contact density for $d=3$ has the analytic expression of Ho:2004 ; Yu:2009

[TABLE]

In particular, it is smooth even in the unitarity limit at $\lambda_{T}/a=0$ so as to be consistent with the “crossover” hypothesis Zwerger:2012 .

III Bulk viscosity

III.1 Formulation

We now turn to the quantum virial expansion for viscosity spectral functions. In terms of the stress–stress response function,

[TABLE]

the spectral functions of bulk and shear viscosities are defined according to Taylor:2010

[TABLE]

Here $\hat{\Pi}_{ij}(t)=e^{i\hat{H}t}\hat{\Pi}_{ij}e^{-i\hat{H}t}$ is the stress tensor operator in the Heisenberg picture and its trace is directly related to the bulk viscosity spectral function via

[TABLE]

For two-component fermions described by the Hamiltonian (1), the stress tensor operator assuming the dimensional regularization is provided by Goldberger:2012 ; Fujii:2018

[TABLE]

Its trace is then simplified into

[TABLE]

where the last term quantifies the conformal symmetry breaking with

[TABLE]

being the contact operator Braaten:2008 ; Barth:2011 ; Hofmann:2012 . Because the commutator of the Hamiltonian with any operator in the grand canonical average vanishes, the bulk viscosity spectral function turns into the favorite form of

[TABLE]

where

[TABLE]

is the contact–contact response function Martinez:2017 ; Fujii:2018 .

The quantum virial expansion for the bulk viscosity spectral function is thus to evaluate

[TABLE]

Here the contact operator (31) can be decomposed into distinct particle number sectors according to

[TABLE]

where nontrivial matrix elements first appear from the two-body sector in the form of

[TABLE]

On the other hand, the time evolution operator such as $e^{-i\hat{H}_{2}t}$ is obtained from Eq. (14) with the simple replacement of $\beta\to it$. Because the bare coupling constant vanishes as $mg\to-2\pi/\ln\Lambda$ for $d=2$ and as $mg\to-2\pi^{2}/\Lambda$ for $d=3$ in the cutoff regularization with $|{\bm{k}}|<\Lambda$, the free part of the time evolution operator $e^{-i\hat{H}_{2}t}|_{\mathrm{free}}$ has vanishing contributions in the limit of $\Lambda\to\infty$ when sandwiched by the contact operators. Therefore, only the interacting part of the time evolution operator $e^{-i\hat{H}_{2}t}|_{\mathrm{int}}$ has non-vanishing contributions to the contact–contact response function. By utilizing the identity of

[TABLE]

we obtain

[TABLE]

Consequently, the contact–contact response function at $O(z^{2})$ is found to be

[TABLE]

from which the bulk viscosity spectral function in the thermodynamic limit reads

[TABLE]

The higher order correction at $O(z^{3})$ is presented in Eq. (A).

III.2 Evaluation

Our concise formula for the bulk viscosity spectral function is even under $\omega\to-\omega$ and has the high-frequency tail of

[TABLE]

where $\mathcal{C}$ is the contact density from Eq. (22). This asymptotic form agrees with Refs. Hofmann:2011 ; Goldberger:2012 for $d=3$ and with Ref. Hofmann:2011 for $d=2$.111Although Eq. (62) of Ref. Hofmann:2011 presents a different result for $d=2$, the substitution of Eqs. (45), (47) into Eqs. (57), (58) in Ref. Hofmann:2011 indeed leads to the same result as our Eq. (41). Also, the integration over $\omega$ leads to

[TABLE]

With the use of the Kramers–Kronig relation followed by the identity of

[TABLE]

we obtain

[TABLE]

which is consistent with the sum rule derived in Ref. Taylor:2010 for $d=3$ and in Ref. Taylor:2012 for $d=2$. We note that the partial derivative with respect to the scattering length at fixed entropy and particle number is equivalent to that at fixed $T$ and $\mu$ to the lowest order in fugacity.

After confirming our formula (40) in light of the known constraints, we are ready to evaluate it more closely. The substitution of Eq. (24) into Eq. (40) leads to

[TABLE]

where the imaginary part of the two-body scattering $T$-matrix is normalized according to

[TABLE]

The three terms in Eq. (III.2) correspond to bound–bound, bound–continuum, and continuum–continuum transitions, which contribute to the bulk viscosity spectral function at $\omega=0$, $|\omega|>1/ma^{2}$, and $\forall\omega$, respectively. Therefore, the first two terms contribute only when the two-body bound state exists, $a>0$, which is always the case for $d=2$. The resulting dynamic bulk viscosity is plotted in Fig. 1 as well as its static limit, $\zeta\equiv\lim_{\omega\to 0^{+}}\zeta(\omega)$, which excluding the singular term of $\delta(\omega)\theta(a)$ is provided by

[TABLE]

Here the integration variable is changed to $\tilde{\varepsilon}=\varepsilon/T$ and the dimensionless scattering length is introduced via $\tilde{a}\equiv\sqrt{mT}a=\sqrt{2\pi}a/\lambda_{T}$.

The static bulk viscosity for $d=2$ is a smooth function of $\lambda_{T}/a>0$ and has the asymptotic forms of

[TABLE]

On the other hand, the static bulk viscosity for $d=3$ has the asymptotic forms of

[TABLE]

and

[TABLE]

which turns out to exhibit the non-analyticity in the unitarity limit as shown in the right panel of Fig. 1. Our results disagree with the static bulk viscosity derived from the kinetic theory, where reported are $\lim_{a\to\infty}\lambda_{T}^{3}\zeta=z^{2}/(12\sqrt{2}\tilde{a}^{2})$ for $d=3$ Dusling:2013 and $\lim_{a\to\infty}\lambda_{T}^{2}\zeta=z^{2}/(2\pi\ln^{4}\tilde{a}^{2})$ for $d=2$ Chafin:2013 . The origin of these discrepancies is currently unknown. We also note that the dynamic bulk viscosity for $d=3$ toward the unitarity limit is provided by

[TABLE]

so that the limits of $\omega\to 0$ and $a\to\infty$ do not commute mutually.

The bulk viscosity spectral function can be measured experimentally via the contact–contact response function (33). When the scattering length is periodically modulated as $a(t)=a+\delta a\sin(\omega t)$, the linear response theory predicts that the contact density responds according to Fujii:2018

[TABLE]

where $\mathcal{C}_{\mathrm{eq}}[a]\equiv\operatorname{Tr}[\hat{C}e^{-\beta(\hat{H}-\mu\hat{N})}]/(ZL^{d})$ is the contact density in thermodynamic equilibrium for the scattering length $a$ and $\chi_{CC}(-\omega)=\chi_{CC}^{*}(\omega)$ is utilized. Then, the thermodynamic identity together with the adiabatic and dynamic sweep theorems Tan:2008a ; Tan:2008b ; Tan:2008c leads to the energy and entropy densities at the production rates of Fujii:2018

[TABLE]

Therefore, by measuring the contact, energy, or entropy density under the periodic modulation of the scattering length, it is possible in principle to extract the contact–contact response function and thus the bulk viscosity spectral function via Eq. (32).

IV Shear viscosity

IV.1 Formulation

The shear viscosity spectral function in Eq. (27) is related to the off-diagonal element of the stress tensor operator via

[TABLE]

where the stress–stress response function has the quantum virial expansion of

[TABLE]

The off-diagonal element of the stress tensor operator (29) then reads

[TABLE]

and its decomposition into one-body and two-body sectors is

[TABLE]

and

[TABLE]

respectively.

The stress–stress response function at $O(z)$ is simply provided by

[TABLE]

and its contribution to the shear viscosity spectral function is thus

[TABLE]

It is also straightforward to evaluate the stress–stress response function at $O(z^{2})$,

[TABLE]

with the use of the time evolution operator such as $e^{-i\hat{H}_{2}t}$ obtained from Eq. (14). The free part with only $e^{-i\hat{H}_{2}t}|_{\mathrm{free}}$ involved is simply provided by

[TABLE]

Its contribution to the shear viscosity spectral function is thus

[TABLE]

On the other hand, the interacting part with at least one $e^{-i\hat{H}_{2}t}|_{\mathrm{int}}$ involved is found to be

[TABLE]

Its contribution to the shear viscosity spectral function is thus

[TABLE]

where ${\bm{K}}={\bm{p}}+{\bm{q}}$ and ${\bm{k}}^{(\prime)}=({\bm{p}}^{(\prime)}-{\bm{q}}^{(\prime)})/2$ are the center-of-mass and relative momenta, respectively, and the first two terms correspond to the self-energy and Maki–Thompson diagrams and the last term corresponds to the Aslamazov–Larkin diagram. The last term can be further simplified by utilizing the identity of

[TABLE]

so that we obtain

[TABLE]

Finally, by combining Eqs. (60), (63), and (IV.1), taking the thermodynamic limit, and performing the resulting momentum integrations, we find that the shear viscosity spectral function is provided by

[TABLE]

where $p$ in the first term is the pressure from Eq. (19) and the integration variable is changed to $\varepsilon={\bm{k}}^{2}/m$.

IV.2 Evaluation

Our concise formula for the shear viscosity spectral function is even under $\omega\to-\omega$ and has the high-frequency tail of

[TABLE]

where $\mathcal{C}$ is the contact density from Eq. (22). This asymptotic form agrees with Refs. Enss:2011 ; Hofmann:2011 ; Goldberger:2012 for $d=3$ and with Ref. Hofmann:2011 for $d=2$. Also, the integration over $\omega$ leads to

[TABLE]

By performing the integration over $\varepsilon$ with the divergent piece transposed to the left hand side and utilizing the identity of

[TABLE]

we obtain

[TABLE]

which for $d=3$ is consistent with the sum rule derived in Ref. Taylor:2010 up to the factors in front of $\mathcal{C}$ and in Refs. Enss:2011 ; Hofmann:2011 .

After confirming our formula (IV.1) in light of the known constraints, we are ready to evaluate its static limit, $\omega\to 0$, where Eq. (IV.1) is formally ill-defined and thus requires the resummation. To this end, we write $\pi\delta(\omega)=\operatorname{Im}(\omega-i0^{+})^{-1}$ in the first term and regard the remainder as corrections in a geometric series Gotze:1972 ; Enss:2011 . Consequently, the shear viscosity spectral function is brought into the Drude form of

[TABLE]

where the memory function is found to be

[TABLE]

Because of $\operatorname{Re}[\mathcal{M}(0)]=0$, the static limit is now provided by $\eta\equiv\lim_{\omega\to 0}\eta(\omega)=p/\operatorname{Im}[\mathcal{M}(0)]$ and the substitution of Eqs. (19) and (24) therein leads to

[TABLE]

Here the integration variable is changed to $\tilde{\varepsilon}=\varepsilon/T$ and the dimensionless scattering length is introduced via $\tilde{a}\equiv\sqrt{mT}a=\sqrt{2\pi}a/\lambda_{T}$. Our results agree with the static shear viscosity derived from the kinetic theory both for $d=3$ Massignan:2005 ; Bruun:2005 and for $d=2$ Bruun:2012 ; Schafer:2012 . Although it is a smooth function of $\lambda_{T}/a>0$ for $d=2$, we find that the static shear viscosity for $d=3$ exhibits the weak non-analyticity starting with $\lambda_{T}^{3}\eta\sim(\ln\tilde{a}^{2})/\tilde{a}^{6}$ in the unitarity limit.

Acknowledgements.

The author thanks Thomas Schäfer for the valuable discussions. He also thanks Tilman Enss and Johannes Hofmann for the correspondences regarding their recent preprints Enss:2019 ; Hofmann:2019 , which overlap with this work. This work was supported by JSPS KAKENHI Grant Nos. JP15K17727 and JP15H05855.

Appendix A Bulk viscosity formula at $O(z^{3})$

Let us derive the formula for the contact–contact response function at $O(z^{3})$, which is to evaluate

[TABLE]

The three-body sector is spanned by $|\sigma_{1}{\bm{p}}_{1},\sigma_{2}{\bm{p}}_{2},\sigma_{3}{\bm{p}}_{3}\rangle$, which is further decomposed into four subsectors where $\{\sigma_{1},\sigma_{2},\sigma_{3}\}=\{\uparrow,\uparrow,\uparrow\},\,\{\uparrow,\uparrow,\downarrow\},\,\{\uparrow,\downarrow,\downarrow\},\,\{\downarrow,\downarrow,\downarrow\}$. The first and last subsectors do not involve the interaction and thus do not contribute to the contact–contact response function. Because the second and third subsectors contribute identically, it is sufficient to consider only the second subsector. By introducing the abbreviated notation,

[TABLE]

and utilizing the identity operator,

[TABLE]

the projection of the Hamiltonian (1) and the contact operator (31) onto this subsector leads to

[TABLE]

and

[TABLE]

In order to evaluate the corresponding Boltzmann operator, we write the Hamiltonian as a sum of kinetic and potential energy operators and utilize the identity of

[TABLE]

Because the summation over $n$ can be performed according to

[TABLE]

the Boltzmann operator is found to be

[TABLE]

Here an infinite series,

[TABLE]

is the three-body scattering $T$-matrix between atom and dimer from incoming $({\bm{p}}_{1},{\bm{K}}-{\bm{p}}_{1})$ to outgoing $({\bm{p}}_{1}^{\prime},{\bm{K}}-{\bm{p}}_{1}^{\prime})$ momenta and solves the integral equation of Braaten:2006 ; Barth:2015

[TABLE]

With the use of the above Boltzmann operator valid for any complex $\beta$ as well as the identity of Eq. (37), we obtain

[TABLE]

where ${\bm{K}}={\bm{p}}_{1}+{\bm{p}}_{2}+{\bm{q}}$ and ${\bm{k}}^{(\prime)}={\bm{p}}_{1}^{(\prime)}-{\bm{K}}/3$ are the center-of-mass and relative momenta, respectively, and the integration variables are changed to $\varepsilon^{(\prime)}=E^{(\prime)}-{\bm{K}}^{2}/6m$. We also note that $\mathcal{T}_{3}(\varepsilon-i0^{+};{\bm{k}},{\bm{k}}^{\prime})\equiv\tilde{\mathcal{T}}_{3}(\varepsilon+\frac{{\bm{K}}^{2}}{6m}-i0^{+},{\bm{K}};\frac{{\bm{K}}}{3}+{\bm{k}},\frac{{\bm{K}}}{3}+{\bm{k}}^{\prime})$ is the three-body scattering $T$-matrix in the center-of-mass frame and is actually independent of ${\bm{K}}$. Consequently, the contact–contact response function at $O(z^{3})$ is found to be

[TABLE]

where the overall factor of 2 accounts for the existence of two interacting subsectors and the first term corresponds to the disconnected diagram, the second and third terms correspond to the self-energy diagrams, and the last term corresponds to the Aslamazov–Larkin diagram. Finally, by combining Eq. (39), taking the thermodynamic limit, and performing the resulting momentum integrations, we find that the bulk viscosity spectral function is provided by

[TABLE]

where the three-body scattering $T$-matrix in the center-of-mass frame solves the integral equation of

[TABLE]

The higher order correction at $O(z^{3})$ may be useful to better understand the static bulk viscosity in the Kubo formalism, which is to be investigated in a future study.

Appendix B Derivation of Eq. (IV.1)

We first decompose the left-hand side of Eq. (IV.1) as

[TABLE]

Then, the definition of the two-body scattering $T$-matrix,

[TABLE]

leads to

[TABLE]

Finally, with the use of

[TABLE]

and

[TABLE]

we obtain

[TABLE]

which is the right-hand side of Eq. (IV.1) under $\omega\delta^{\prime}(\omega)=-\delta(\omega)$.