# Spectral functions of the strongly interacting 3D Fermi gas

Christian H. Johansen  Max-Planck-Institut für Physik komplexer Systeme, 01187
Dresden, Germany  Institute of Physics, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University in Toruń, Grudzia̧dzka 5, 87-100
Toruń, Poland  Bernhard Frank  Institut für Theoretische Physik and Würzburg-
Dresden Cluster of Excellence ct.qmat, Technische Universität Dresden, 01062
Dresden, Germany  Johannes Lang<EMAIL_ADDRESS>Institut für Theoretische
Physik, Universität zu Köln, Zülpicher Straße 77, 50937 Cologne, Germany  Max-
Planck-Institut für Physik komplexer Systeme, 01187 Dresden, Germany

###### Abstract

Computing dynamical properties of strongly interacting quantum many-body
systems poses a major challenge to theoretical approaches. Usually, one has to
resort to numerical analytic continuation of results on imaginary frequencies,
which is a mathematically ill-defined procedure. Here, we present an efficient
method to compute the spectral functions of the two-component Fermi gas near
the strongly interacting unitary limit directly in real frequencies. To this
end, we combine the Keldysh path integral that is defined in real time with
the self-consistent T-matrix approximation. The latter is known to predict
thermodynamic and transport properties in good agreement with experimental
observations in ultracold atoms. We validate our method by comparison with
thermodynamic quantities obtained from imaginary time calculations and by
transforming our real-time propagators to imaginary time. By comparison with
state-of-the-art numerical analytic continuation of the imaginary time
results, we show that our real-time results give qualitative improvements for
dynamical quantities. Moreover, we show that no significant pseudogap regime
exists in the self-consistent T-matrix approximation above the critical
temperature $T_{c}$, an issue that has been under significant debate. We close
by pointing out the versatile nature of our method as it can be extended to
other systems, like the spin- or mass-imbalanced Fermi gas, other Bose-Fermi
models, 2D systems as well as systems out of equilibrium.

††preprint: APS/123-QED

##  I Introduction

The two-component Fermi gas, realized in experiments with ultra-cold atoms
both in 2D and 3D, nowadays plays a central role in the study of strongly
correlated systems [1, 2, 3]. This is a consequence of the simple form of the
interactions that only depend on the scattering length. However, at the same
time the system exhibits a highly non-perturbative quantum many-body regime
close to the quantum mechanical unitary scattering limit reached in the
vicinity of a Feshbach resonance. The presence of only a single relevant
interaction parameter gives rise to universal physics (with corrections in 2D
[4, 5]). Despite this simplicity the low-temperature phase diagram contains a
variety of neutral superfluid phases depending on the spin polarization [6,
7]. Equally important for the scientific interest in the Fermi gas is the fact
that it can be investigated experimentally with a high degree of control. One
of the first examples was the experimental determination of the phase diagram
as a function of the spin-polarization [8]. Later, the thermodynamic
properties of the spin-balanced gas at the Feshbach resonance were measured
with high precision in both the superfluid and normal phase [9]. More
recently, further technological advances have allowed the direct measurement
of the Tan contact [10, 11, 12] for the balanced Fermi gas at unitarity with a
homogeneous box trap [13]. The latter also paved the way for the determination
of local transport coefficients like sound diffusion [14].

The remarkably high precision reached in experiments in combination with the
universal and non-perturbative nature of the system makes it ideal for
benchmarking theoretical approaches. The self-consistent T-matrix
approximation is a theoretical method that has proven surprisingly accurate
throughout the BCS-BEC crossover. In particular, its predictions for
thermodynamic quantities including critical temperature, Bertsch parameter,
and the Tan contact agree extremely well with experimental observations [15,
16, 17]. Furthermore, combining the self-consistent T-matrix approximation
with suitable scattering time approximations even describes the transport
properties on a quantitative level [18].

Experiments not only focus on thermodynamic and transport quantities but also
aim for the investigation of dynamical properties. The latter is achieved
using radio-frequency (rf) spectroscopy which provides access to the spectral
functions of the system [19]. With this technique, experiments have been able
to measure the superfluid gap of an imbalanced gas [20]. The spectral
measurements have also revealed the existence of an intermediary regime in a
spin-imbalanced Fermi gas where the quasiparticle picture breaks down,
manifesting through the increasing width of the spectral features [21]. In a
2D gas, another interesting spectral feature has been measured, namely the
existence of a pseudogap above the critical temperature [22]. While the
existence of a pseudogap was expected in 2D, its fate in 3D is less obvious.
Although a typical back-bending of the dispersion has been observed well above
the critical temperature [23] the expected suppression of the spin
susceptibility is absent [24].

Having direct experimental access to the spectral information motivates
theoretical computation of the spectrum with equally high reliability as
thermodynamic quantities. It is expected, that self-consistent approximations
produce reasonable spectra across all scales without obvious artifacts.
Combined with the good agreement between experimental thermodynamic
observables, we thus expect similar quantitative agreement for the spectral
properties obtained through the self-consistent T-matrix approximation.
However, reliably extracting spectral properties with high accuracy has proven
difficult, as the procedure to date still involves the numerical analytic
continuation (NAC) of data in imaginary frequencies to the real frequency
axis, a mathematically ill-defined process, that results in uncontrolled
errors [25, 26]. As a result, the comparison of dynamical information gained
from rf-spectra lacks the triumphs of determining the thermodynamic and
universal parameters [27, 28, 29, 30].

In this work, we develop a method based on the Keldysh path integral that
allows us to calculate the self-consistent T-matrix approximation to very high
precision directly in real frequencies. We exploit the convolution structure
of the T-matrix approximation to express it through Fourier transforms. The
key to the success of our approach lies in a piecewise interpolation scheme
that allows one to analytically treat the fast oscillations that inevitably
appear in such a real-time formulation. The method is validated against
previous results for thermodynamic quantities obtained from Matsubara
propagators in imaginary time and differs by less than $1\%$. This level of
precision is remarkable because the previous results are known to be very
accurate in imaginary time, demonstrating that our method is capable of
reliably computing thermodynamic properties. For systems in thermal
equilibrium, the main advantage of our technique is the direct access to
dynamical quantities as opposed to the imaginary-time approach, which yields
uncontrolled errors when analytically continued to real frequencies. As a
result, the formulation of the self-consistent T-matrix within Keldysh field
theory allows us to quantitatively describe both spectral features and
thermodynamics.

The paper is structured as follows, in Sec. II, we introduce the two-component
Fermi gas as realized in experiments with ultra-cold atoms. Readers familiar
with this subject can skip ahead to Eq. (7). The Keldysh formalism is
recapitulated and applied to the Fermi gas in Sec. III. Our central
technological development is presented in Sec. IV and carefully validated in
Secs. V and VI. We then apply the method to spectral functions near the
superfluid transition in the BCS-BEC crossover, addressing the issue of the
pseudogap at unitarity (Sec. VII) before we provide a short outlook in Sec.
VIII.

##  II The interacting Fermi gas

A gas with density $n$, by definition, describes a dilute system, meaning that
the average distance between particles $n^{-1/3}$ is much larger than the size
of the particles. Under these conditions, neutral particles quite generally
interact via induced dipole-dipole interactions resulting in an attractive
two-body potential generated by the van der Waals force [31]. To describe the
two-body scattering potential one can use a simple model, which cuts off the
attraction at some small atomic distance $R_{0}$ [3]

$V(r)=\begin{cases}-C_{6}/r^{6},&r>R_{0},\\\ \infty,&r\leq R_{0},\end{cases}$
(1)

where $C_{6}$ is a positive constant that is determined by the specific atoms
used in the experiments.

At low temperatures, the gas becomes degenerate, meaning that the particles’
wave functions overlap. This behavior manifests when $n\lambda_{T}^{3}\gg 1$
with $\lambda_{T}=\hbar\sqrt{2\pi/mT}$ being the thermal de Broglie
wavelength, implying that scattering needs to be treated quantum mechanically
[3].

To simplify the theoretical model one can exploit a separation of length
scales in the ultra-cold dilute quantum gas

$R_{0}\ll l_{vdW}\ll n^{-1/3}\ll\lambda_{T},$ (2)

where $l_{vdW}=(mC_{6}/\hbar^{2})^{1/4}/2$ is the length scale associated with
the van der Waals interaction in Eq. 1. Due to the small kinetic energy of the
gas, it is sufficient to consider the s-wave scattering channel with the
isotropic scattering amplitude [32]

$f(k)=\frac{-1}{a^{-1}-r_{e}k^{2}/2+...+ik},$ (3)

with $a$ being the scattering length and $r_{e}$ the effective range. For the
potential in Eq. 1 the scattering length and effective range can be
analytically computed [3, 33, 34]

$\displaystyle a$
$\displaystyle=\tilde{a}\bigg{(}1-\tan\left(\Phi-3\pi/8\right)\bigg{)},$ (4a)
$\displaystyle r_{e}$
$\displaystyle=2.92\tilde{a}\left(1-2\frac{\tilde{a}}{a}+2\left(\frac{\tilde{a}}{a}\right)^{2}\right).$
(4b)

The results can be parametrized through the mean scattering length
$\tilde{a}=0.956l_{vdW}$ and the WKB-phase $\Phi=2l_{vdW}^{2}/R_{0}^{2}$.
Looking at Eq. 4a it is seen that the scattering length diverges for specific
values of $\Phi$. Such a divergence appears when the collision of the incoming
particles is resonant with a bound state. Even though these resonances are
determined by the atomic scales $l_{vdW}$ and $R_{0}$, Feshbach resonances
[35, 36] allow one to tune the scattering properties via an external magnetic
field.

Here we only briefly describe the main concept in the experimentally relevant
case of alkali atoms with one valence electron and refer the interested reader
to the detailed review [31]. A Feshbach resonance emerges because the external
magnetic field polarizes the valence electrons of the atoms. For low energies,
the physics of a collision is therefore mainly described by the triplet state.
However, the wave function for the singlet state of two fermions is symmetric,
which means that the singlet potential at small distances is much more
attractive, compared to the triplet potential, and therefore supports many
bound states. Because the singlet and triplet states have different magnetic
moments, it is possible to create a resonance between the incoming triplet
states (open channel) and a bound singlet state (closed channel), by tuning
the magnetic field. As the magnetic field is finite there will be interactions
between the close-to-resonance open channel and the closed channel. To
describe two fermionic atoms of the same type scattering into the open
channel, one can use the two-channel model [37, 38]

$\displaystyle\begin{split}H_{2c}=&\int\differential^{3}r\Bigg{[}\sum_{\sigma}\psi_{\sigma}^{\dagger}(\textbf{r})\left(-\frac{\hbar^{2}\nabla^{2}}{2m}\right)\psi_{\sigma}(\textbf{r})\\\
&+\Delta^{\dagger}(\textbf{r})\left(-\frac{\hbar^{2}\nabla^{2}}{4m}+\nu(B)\right)\Delta(\textbf{r})\\\
&+\\!\tilde{g}\\!\int\\!\\!\differential^{3}r^{\prime}\chi(\absolutevalue{\textbf{r}-\textbf{r}^{\prime}})\\\
&\qquad\quad\times\\!\left(\\!\Delta^{\dagger}\\!\left(\frac{\textbf{r}+\textbf{r}^{\prime}}{2}\right)\psi_{\uparrow}(\textbf{r})\psi_{\downarrow}(\textbf{r}^{\prime})+\text{h.c.}\\!\right)\\!\Bigg{]},\end{split}$
(5)

where $\psi_{\sigma}$ and $\psi_{\sigma}^{\dagger}$ are the fermionic
annihilation and creation operators with spin $\sigma$ and mass $m$. The
closed channel is described by the annihilation operator of the bosonic
pairing field $\Delta$, as well as the corresponding creation operator
$\Delta^{\dagger}$. The magnetic field $B$ leads to a detuning of the closed
channel by $\nu(B)=\delta\mu\left(B-B_{c}\right)$, with the difference in the
magnetic moment between the triplet and singlet states being $\delta\mu$ and
$B_{c}$ denoting the field strength at which they become resonant. The
scattering of two fermions in the open channel into the pairing field happens
through an isotropic form factor that is normalized
$\int\differential^{3}r\,\chi(\absolutevalue{\textbf{r}})=1$. This means that
the strength of the interaction is given by $\tilde{g}$. By eliminating the
pairing field one can compute the corresponding scattering length and
effective range in Eq. 3 originating from the Feshbach resonance [39]

$\displaystyle a^{-1}$
$\displaystyle=-\frac{mr^{\star}}{\hbar^{2}}\nu(B)+\frac{1}{2\tilde{a}},$ (6a)
$\displaystyle r_{e}$
$\displaystyle=-2r^{\star}+3\tilde{a}\left(1-\frac{4\tilde{a}}{3a}\right),$
(6b)

where $r^{\star}=2\pi\hbar^{4}/(m\tilde{g})^{2}$ is the length scale
associated to the Feshbach coupling. Importantly there exists a critical field
strength $B_{0}=\hbar^{2}(2\tilde{a}\,\delta\mu\,m\,r^{\star})^{-1}+B_{c}$ at
which the scattering length $a$ diverges.

The behavior of the system qualitatively differs depending on the relation
between the van der Waals length scale and the strength of the Feshbach
resonance. This motivates introducing the resonance strength
$s_{res}=\tilde{a}/r^{\star}$ [31], which compares the two scales. In the
following, we will consider a so-called open-channel dominated Feshbach
resonance where $s_{res}\gg 1$, such that $\tilde{g}$ is large. Furthermore,
such a resonance directly leads to $r_{e}$ being on the order of $l_{vdW}$.
Since the Fermi momentum $k_{F}\approx n^{1/3}$, the separation of length
scales in an ultra-cold gas in Eq. 2 means that $r_{e}k_{F}\ll 1$. The system
is therefore well described by considering a contact interaction where
$r_{e}\rightarrow 0$. By rescaling the field operator
$\Delta\rightarrow\Delta/\tilde{g}$ the Hamiltonian Eq. 5 can be rewritten as

$\displaystyle H=\int\differential^{3}r\Bigg{[}$
$\displaystyle\sum_{\sigma}\psi_{\sigma}^{\dagger}(\textbf{r})\left(-\frac{\nabla^{2}}{2m_{\sigma}}-\mu_{\sigma}\right)\psi_{\sigma}(\textbf{r})$
(7)
$\displaystyle-\Delta^{\dagger}(\textbf{r})\frac{1}{g_{0}}\Delta(\textbf{r})$
$\displaystyle+\Delta^{\dagger}\left(\textbf{r}\right)\psi_{\uparrow}(\textbf{r})\psi_{\downarrow}(\textbf{r})+\text{h.c.}\Bigg{]},$

where we have now switched to units of $\hbar=1$ and $k_{B}=1$, which will be
used throughout the rest of the paper. Furthermore, the system has also been
generalized by including different chemical potentials of the two different
fermion species $\mu_{\sigma}$ and mass imbalance. We have also defined
$g_{0}=-\tilde{g}^{2}/\nu(B)$ similarly to [40] and neglected the kinetic
energy of the bare $\Delta$-field due to $\tilde{g}$ being large. This
Hamiltonian quantitatively describes experiments with an open-channel
dominated Feshbach resonance, which for example is the case for several of the
experiments that have investigated the BCS-BEC crossover [41, 9]. The zero-
range interactions lead to a divergence and the interactions therefore have to
be cut off at the length scale $l_{vdW}$. This renormalization connects
$g_{0}$ to the experimental scattering length, which is discussed in App. A.1.

##  III Field theoretical method

With the effective Hamiltonian for the system defined by Eq. 7, we now seek to
compute observables. The system is theoretically challenging because it is
strongly interacting, three dimensional and at finite density. For computing
thermodynamic properties imaginary-time quantum field theory has been very
successful. From this theory dynamical quantities can be computed using
numerical analytical continuation (NAC). To avoid this inherently uncontrolled
method, while still taking advantage of the successes of quantum field theory,
we use the real-time Keldysh path integral approach.

In Sec. III.1 we will give a brief introduction to the generic structure of
the Keldysh path integral. For the reader already familiar with the theory, it
will be sufficient to consider the conventions fixed in Eqs. 15, 17 and 22, as
these will be used in the remaining parts of the paper. For a more in-depth
discussion of the construction of real-time quantum field theories and its
different formulations, the interested reader is referred to existing
textbooks and reviews [42, 43, 44, 45, 46, 47].

Having described the generic structure of the real-time path integral we then,
in Sec. III.2, apply it to the strongly interacting Fermi gas and derive the
real-time formulation of the self-consistent T-matrix approximation.

###  III.1 Keldysh path integral

To construct a quantum field theory in real time, one defines the partition
function as a trace over the density matrix

$\mathcal{Z}=\Tr{\rho(t)}=\Tr{U(t,t_{0})\rho(t_{0})U(t,t_{0})^{\dagger}},$ (8)

where $U(t,t_{0})$ is the time-evolution operator. With the partition function
defined through the density matrix, this formalism allows to consider systems
in and out of equilibrium as well as time evolution. For the current paper we
focus on an equilibrium system and use the real-time formulation to avoid the
NAC. Because $U(t,t_{0})$ and $U(t,t_{0})^{\dagger}$ contain opposite time-
ordering, it is necessary to differentiate between forward and backward
propagation. This can conveniently be done using the time contour
$\mathfrak{C}$ shown in Fig. 1.

Figure 1: Time contour $\mathfrak{C}$ used to construct a real-time action.

Using this contour, the system is evolved from some initial time $t_{0}$ to
some time $t$ on the $+$ branch and then back to $t_{0}$ along the $-$ branch
[42]. To evaluate the partition function using the contour $\mathfrak{C}$, one
uses coherent Grassmann fields to describe fermionic degrees of freedom and
coherent complex fields to describe bosonic degrees of freedom, similar to the
imaginary time procedure [46]. The differentiation between forward and
backward propagation is achieved through the use of two different fields
$\phi_{+/-}$, which respectively have support exclusively on the $+$ or $-$
branch in Fig. 1. The partition function can then be written as a functional
path integral over the fields $\phi_{\pm}$

$\mathcal{Z}=\int\mathcal{D}\left\\{\phi,\bar{\phi}\right\\}\text{e}^{i\mathcal{S}}.$
(9)

For complex numbers, the bar denotes complex conjugation while for Grassmann
numbers it indicates a different, independent Grassmann number. The measure in
Eq. 9 should be understood as a functional measure over the independent
fields, $\phi_{\alpha}(\textbf{r},t)$ and $\bar{\phi}_{\alpha}(\textbf{r},t)$
for Grassmann fields, and $\text{Re}\,\phi_{\alpha}(\textbf{r},t)$ and
$\text{Im}\,\phi_{\alpha}(\textbf{r},t)$ for complex fields. In both cases
$\alpha\in\\{+,-\\}$ specifies the branch of the Keldysh contour. Following
this procedure, the continuous action in time $\mathcal{S}$ appearing in Eq. 9
takes the form

$\mathcal{S}=\int\differential^{4}x\sum_{\alpha=+,-}\alpha\left(\bar{\phi}_{\alpha}i\partial_{t}\phi_{\alpha}-H\left[\bar{\phi}_{\alpha},\phi_{\alpha}\right]\right),$
(10)

where $H\left[\bar{\phi}_{\alpha},\phi_{\alpha}\right]$ is the normal-ordered
Hamiltonian density with all field creation (annihilation) operators replaced
by $\bar{\phi}_{\alpha}\,(\phi_{\alpha})$. As the system is time-independent
and only has local interactions, all fields are evaluated at the same space-
time point $x=(\textbf{r},t)$. Generalizing the construction to time-dependent
Hamiltonians and non-local interactions is straightforward. For this action
the contour has been extended to the entire real axis and interactions are
switched on adiabatically in the infinite past [46]. For interacting systems
it is also convenient to write the action in the form

$\mathcal{S}=\mathcal{S}_{0}+\mathcal{S}_{I},$ (11)

where $\mathcal{S}_{0}$ is quadratic in the fields and $\mathcal{S}_{I}$
contains the interactions.

The continuous representation Eq. 10 does not explicitly contain the two
boundaries at $\pm\infty$ of the contour in Fig. 1. These can be conveniently
included after a unitary linear transformation on the contour degrees of
freedom [42]. The linear transformation then maps to the symmetric and anti-
symmetric superpositions
$\phi_{c/q}=\left(\phi_{+}\pm\phi_{-}\right)/\sqrt{2}$. When $\phi$ is a
complex number, therefore representing a bosonic degree of freedom, one can
derive mean-field equations where $\phi_{c}$ will represent the possibly
finite mean-field. Conversely, in a mean-field description of a bosonic
system, $\phi_{q}$ can be connected to quantum fluctuations. To this extent
$\phi_{c}$ is often denoted the classical field and $\phi_{q}$ is denoted as
the quantum field. For fermionic degrees of freedom, where $\phi$ is a
Grassmann number, we only use these conventions to make the equivalence of the
microscopic formulation for fermions and bosons apparent.

Considering a spatially homogeneous equilibrium theory, means that the non-
interacting action is diagonalized by a transformation to reciprocal space,
represented by the four-vector $p=(\textbf{k},\omega)$. It is related to real
space by the Fourier transformation

$f(x)=\int\frac{d^{4}p}{(2\pi)^{4}}\text{e}^{ip\cdot x}f(p),$ (12)

where the four-vector dot product is defined as $p\cdot
x=\textbf{k}\cdot\textbf{r}-\omega t$ and
$d^{4}p=\differential\omega\,\differential\textbf{k}^{3}$. In reciprocal
space, the bare part of the action is given by

$\mathcal{S}_{0}=\int\frac{d^{4}p}{(2\pi)^{4}}\matrixquantity(\bar{\phi}_{c}(p)\\\
\bar{\phi}_{q}(p))^{T}\mathcal{G}_{0}^{-1}(p)\matrixquantity(\phi_{c}(p)\\\
\phi_{q}(p)),$ (13)

where the inverse contour propagator takes the form

$\mathcal{G}_{0}^{-1}(p)=\matrixquantity(0&(G_{0}^{A})^{-1}(p)\\\
(G_{0}^{R})^{-1}(p)&P^{K}(p)).$ (14)

The lower off-diagonal element, $(G_{0}^{R})^{-1}(p)$, is the inverse bare
retarded propagator which is causal in the sense that
$G_{0}^{R}(x)\sim\theta(t)$. Consequently, it is connected to the inverse,
bare advanced propagator $(G_{0}^{A})^{-1}(p)$ by complex conjugation. The
lower-diagonal component, $P^{K}$, is a purely imaginary and infinitesimally
small term that ensures the normalization of the partition function.

The effect of interactions, on correlation functions, are included through the
self-energy matrix $\underline{\Sigma}$, which has a similar causality
structure as the inverse propagator $\mathcal{G}_{0}^{-1}$, such that the
interacting inverse contour propagator is

$\mathcal{G}^{-1}=\matrixquantity(0&(G_{0}^{A})^{-1}-\Sigma^{A}\\\
(G_{0}^{R})^{-1}-\Sigma^{R}&-\Sigma^{K}).$ (15)

Note, that the bare infinitesimal $P^{K}$ has been dropped since interactions
generally makes $\Sigma^{K}$ finite. Inverting Eq. 15 one arrives at the
contour propagator

$\mathcal{G}=\matrixquantity(G^{K}&G^{R}\\\ G^{A}&0).$ (16)

This inversion leads to the Dyson equation for the retarded propagator

$G^{R}(p)=-i\left\langle\phi_{c}(p)\bar{\phi}_{q}(p)\right\rangle=\frac{1}{\left(G_{0}^{R}\right)^{-1}(p)-\Sigma^{R}(p)},$
(17)

and a similar equation for $G^{A}$. Transforming back to real-space and the
$\pm$-basis, one recovers the operator definition of the retarded propagator

$G^{R}(x,x^{\prime})=-i\theta(t-t^{\prime})\left\langle\left[\hat{\phi}(x),\hat{\phi}^{\dagger}(x^{\prime})\right]_{\mp}\right\rangle,$
(18)

where $\hat{\phi}$ denotes the field operator and the commutator (minus-sign)
is used for bosonic fields while the anti-commutator (plus sign) is used for
fermionic fields. In this basis the relation between the advanced and retarded
propagators are $G^{A}(x,x^{\prime})=\bar{G}^{R}(x^{\prime},x)$. The retarded
propagator contains information about the spectrum of the system, as its
imaginary part is directly linked to the spectral function

$\mathcal{A}(p)=-2\imaginary G^{R}(p).$ (19)

Probability conservation $\mathcal{Z}\equiv 1$ guarantees that the lower
diagonal entry of Eq. 16 vanishes. Explicitly, one has
$\left\langle\phi_{q}(p)\bar{\phi}_{q}(p)\right\rangle=0$. The upper-diagonal
element in $\mathcal{G}$ is referred to as the Keldysh propagator and given by

$G^{K}(p)=-i\left\langle\phi_{c}(p)\bar{\phi}_{c}(p)\right\rangle=\absolutevalue{G^{R}(p)}^{2}\Sigma^{K}(p),$
(20)

which is anti-Hermitian and related to how the spectrum is occupied [42]. In
analogy to Eq. 18 for the retarded propagator, the Keldysh propagator can be
connected to the operator definition by transforming to real-space and
$\pm$-basis

$G^{K}(x,x^{\prime})=-i\left\langle\left[\hat{\phi}(x),\hat{\phi}^{\dagger}(x^{\prime})\right]_{\pm}\right\rangle,$
(21)

where the upper (lower) sign represents bosonic (fermionic) fields. It is
often convenient to separate out the vacuum term of $G^{K}$ and define the
occupied part of the spectrum as

$\delta
G^{K}(p)=G^{K}(p)+i\mathcal{A}(p)=\absolutevalue{G^{R}(p)}^{2}\delta\Sigma^{K}(p),$
(22)

where $\delta\Sigma^{K}(p)=\Sigma^{K}(p)-2i\imaginary\Sigma^{R}(p)$. With this
decomposition one can compute the momentum distribution by integrating $\delta
G^{K}$ over frequency

$n(\textbf{k})=\pm i\int\frac{\differential\omega}{4\pi}\delta
G^{K}_{\sigma}(p),$ (23)

where the sign is determined by the statistics of the degree of freedom as in
Eq. 21. For degrees of freedom with well-defined particle number the density
$n$ follows directly by integrating over momentum.

Out of equilibrium, $\delta G^{K}$ is independent of $G^{R}$, but in thermal
equilibrium, they are linked by the fluctuation-dissipation relation (FDR)

$i\delta G^{K}(p)=\pm 2n_{B/F}(\omega)\mathcal{A}(p),$ (24)

with the same sign conventions as in Eq. 23 and the thermal distribution
functions $n_{B/F}(\omega)=\left(\text{e}^{\beta\omega}\mp 1\right)^{-1}$.
Note that in equilibrium $\delta\Sigma^{K}$ and $\imaginary\Sigma^{R}$ also
satisfy the FDR (24). This means that in thermal equilibrium, similar to the
formulation in imaginary time, all interaction effects are encoded in a single
function: the retarded self-energy. However, as discussed in the appendix,
computing the Keldysh self-energy directly can be advantageous numerically as
it is generally less sensitive to numerical cutoffs in frequency and momentum
space than $G^{R}$. The FDR can therefore be used to remove numerical noise
from the spectral function.

###  III.2 Self-consistent T-matrix

Having outlined the generic construction of the Keldysh field theory we now
apply it to the strongly interacting Fermi gas. The Hamiltonian in Eq. 7
contains two different fermionic degrees of freedom, which we represent with
the Grassmann fields $\psi_{c/q,\sigma}$ with
$\sigma\in\\{\uparrow,\downarrow\\}$. The bosonic pairing field on the other
hand is represented with the complex fields $\Delta_{c/q}$. Following the
procedure in the previous section the action can be split into three parts

$\mathcal{S}=\mathcal{S}_{\psi}+\mathcal{S}_{\Delta}+\mathcal{S}_{I}.$ (25)

The action of free fermions $\mathcal{S}_{\psi}$ is described by the inverse
retarded/advanced propagators

$(G_{0,\sigma}^{R/A})^{-1}(p)=\omega-\bm{\mathbf{k}}^{2}/2m_{\sigma}+\mu_{\sigma}\pm
i0^{+},$ (26)

while the inverse free propagators for the pairing field are

$\left(\Gamma_{0}^{R/A}\right)^{-1}=1/g_{0}\pm i0^{+}.$ (27)

In both cases the infinitesimal imaginary part ensures the correct causal
form. When deriving the free action for the pairing field
$\mathcal{S}_{\Delta}$, the dynamic part has been neglected. This is
consistent with the contact interaction approximation which requires a
renormalization of the coupling $g_{0}$. In that procedure $g_{0}\rightarrow
0$ simultaneously with a UV cutoff that goes to $\infty$ in such a way, that
the experimental binding energy is recovered. This procedure means $\omega\ll
1/g_{0}$ for all physically relevant values of $\omega$, and the
$\omega$-dependence can therefore be neglected.

The final term of the action in Eq. 25 describes interactions between the
fermions and the pairing field

$\displaystyle S_{I}=$
$\displaystyle-\frac{1}{\sqrt{2}}\int\\!\prod_{i=1}^{3}\frac{d^{4}p_{i}}{(2\pi)^{4}}\Big{[}M^{\gamma}_{\alpha,\beta}\bar{\psi}_{\alpha,\uparrow}(p_{1})\bar{\psi}_{\beta,\downarrow}(p_{2})\Delta_{\gamma}(p_{3})$
(28)
$\displaystyle\qquad\qquad\qquad\qquad\times\delta(p_{3}-p_{1}-p_{2})+c.c\Big{]},$

where the sum over Keldysh indices $\\{\alpha,\beta,\gamma\\}\in\\{c,q\\}$ is
implicit but the four-momentum structure has been kept explicit. The
components of the tensor $M$ are given by the first Pauli matrix
$M^{c}_{\alpha,\beta}=\sigma^{1}_{\alpha,\beta}$ and
$M^{q}_{\alpha,\beta}=\delta_{\alpha,\beta}$.

The inclusion of interactions necessitates finding a suitable approximation
scheme for the self-energies. To compute thermodynamic quantities that satisfy
all thermodynamic relations, it is necessary to employ a self-consistent
treatment of the pairing channel [16, 48, 30]. To derive such an approximation
one can start from the the 2-particle-irreducible effective action [49, 50] or
a similar formulation using the Luttinger-Ward functional $\Phi$ [51], which
is the interaction part of the grand canonical potential. The formulation in
terms of the former is more common in the context of Keldysh field theory, as
it is more general, applying also when the system is out of equilibrium. There
the self-consistent approximation preserves all conservation laws of the exact
action, wherefore it is called a conserving approximation in this context.
Explicitly for the Hamiltonian in Eq. 7, this entails the conservation of
particle number, energy, and momentum. We point out that in the present
equilibrium setting both formalisms are identical. To relate to the existing
literature we express our approximation in terms of the $\Phi$-functional
augmented by the Keldysh structure of the propagators.

The $\Phi$-functional is elegantly represented diagrammatically as a sum of
closed diagrams constructed from the vertices of the action, such that in each
diagram more than two lines have to be cut for it to become disconnected. For
the interaction in Eq. 28 the lowest order $\Phi$-functional that can be
constructed is the two-loop diagram shown in Fig. 2. This forms the basis of
the self-consistent T-matrix approximation. Due to the contact nature of the
interaction, the Bethe-Salpeter equation in Fig. 2 has the same form as the
Dyson equation, which is the reason why the pairing field description is
convenient. The self-energies can be derived by taking functional derivatives
of $\Phi$ with respect to the propagators

$\Sigma_{\sigma}=\frac{\delta\Phi}{\delta G_{\sigma}}\text{ and
}\Sigma_{\Gamma}=\frac{\delta\Phi}{\delta\Gamma}.$ (29)

The resulting self-energies are also drawn in Fig. 2, where all propagators
are understood as dressed by self-energies. Hence, Eq. 29 leads to three
coupled self-consistent integral equations for the retarded self-energies,
which together with the three corresponding Dyson equations of the form of Eq.
17 give rise to a closed set of equations.

Figure 2: The $\Phi$-functional, the resulting self-energies, and the Dyson
equations for both the pairing field and spin-$\uparrow$ fermion. The full
fermion propagators are black bold lines, with an arrow indicating the species
and the pair propagator is the grey box labeled with $\Gamma$. The arrows
indicate the direction of propagation and dashed line is the bare contact
interaction for the electrons, while the thin fermion lines represent the bare
fermion propagators. The self-energy and Dyson equation for $\downarrow$ is
obtained by flipping the spins of the spin-$\uparrow$ counterparts.

At this point, the general structure of the theory has been fully defined,
recovering the self-consistent T-matrix approximation [15] but formulated
directly in real frequencies.

By comparison with high-precision experiments [9, 13] and bold diagrammatic
Monte-Carlo results [52, 53] this approximation has proven surprisingly
accurate for thermodynamic quantities in the experimentally relevant regime of
strong interactions $|k_{F}a|\gtrsim 1$. Outside this regime larger deviations
from exactly known results occur: for weak interactions, BCS theory is
recovered, lacking the quantitatively important corrections found by Gor‘kov
and Melik-Barkhudarov [54]. For very deeply bound pairs, the gas is well
described in terms of weakly interacting bosons. However, the scattering
length between these pairs [55] $a_{\text{pp}}=0.6a$ is not correctly
recovered by the T-matrix approximation [15]. Focusing on parameters near
unitarity, several other approximation schemes have been developed. These
include the historic non-selfconsistent approach by Nozières and Schmidt-Rink
[56], an expansion in the inverse number of fermionic components [57], the
Chevy ansatz [58], Gaussian fluctuations [59, 60] as well as partially self-
consistent T-matrix approximations [30]. Comparison across these methods
systematically shows that thermodynamic quantities are most accurately
approximated by the self-consistent T-matrix.

In contrast to an imaginary-time formulation, self-energies are not only
distinguished by particle species but also by their Keldysh index $R/A/K$.
Hence, one has to account for the contour degree of freedom ($q,c$), which we
have suppressed in our discussion of Eq. 29. Its inclusion is straightforward:
One constructs all diagrams that involve an odd number of quantum fields at
each vertex, thereby accounting for the fact that $S_{\text{int}}$ is
antisymmetric under exchange of the $+$ and $-$ branches of the Keldysh
contour $\mathfrak{C}$. The contraction of two legs is then identified as one
of the four possible propagators
$\left\langle\phi_{\alpha}(p)\bar{\phi}_{\beta}(p)\right\rangle$, with
$\alpha,\beta\in\\{c,q\\}$ and $\phi\in\\{\psi_{\sigma},\Delta\\}$ thus
removing all diagrams that contain a $q$-$q$ propagator.

After extracting the vacuum contribution as discussed previously, one finds a
retarded self-energy for the pair propagator given by

$\displaystyle\Sigma_{\Gamma}^{R}(p)=\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}\Big{(}$
$\displaystyle
2G^{R}_{\uparrow}\left(p-p^{\prime}\right)G^{R}_{\downarrow}\left(p^{\prime}\right)$
(30) $\displaystyle+\delta
G^{K}_{\uparrow}\left(p-p^{\prime}\right)G^{R}_{\downarrow}\left(p^{\prime}\right)$
$\displaystyle+G^{R}_{\uparrow}\left(p-p^{\prime}\right)\delta
G^{K}_{\downarrow}\left(p^{\prime}\right)\Big{)}\,.$

The first term determines the vacuum limit as it is independent of $\delta
G^{K}$. The two remaining terms describe corrections due to occupation in
either of the two fermion species. In appendix A.1 we discuss how the vacuum
term can be used to renormalize the bare divergent propagator and connect the
theory to the experimental scattering length. As expected, the Keldysh self-
energy describes the occupation of the fluctuations in the pair propagator and
thus requires finite occupation of both fermion species

$\delta\Sigma_{\Gamma}^{K}(p)=\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}\delta
G^{K}_{\uparrow}\left(p-p^{\prime}\right)\delta
G^{K}_{\downarrow}\left(p^{\prime}\right).$ (31)

The full pair propagator $\Gamma$ is then used in the retarded self-energy for
the fermions given by

$\displaystyle\Sigma^{R}_{\sigma}\left(p\right)=-\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}\Big{(}$
$\displaystyle\Gamma^{R}\left(p+p^{\prime}\right)\delta
G^{K}_{\bar{\sigma}}\left(p^{\prime}\right)$ (32)
$\displaystyle+\delta\Gamma^{K}\left(p+p^{\prime}\right)G^{A}_{\bar{\sigma}}\left(p^{\prime}\right)\Big{)},$

where $\bar{\sigma}$ is the opposite species of $\sigma$. Since pairs are
absent in the vacuum limit, $\Sigma^{R}_{\sigma}$ is non-zero only in the
presence of a finite density of fermions with opposite spin. This is also the
case for the occupation of the fermion spectrum which is determined through
Eq. 22 with

$\displaystyle\delta\Sigma^{K}_{\sigma}\left(p\right)=-\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}\delta\Gamma^{K}\left(p+p^{\prime}\right)\Big{(}$
$\displaystyle\delta G_{\bar{\sigma}}^{K}\left(p^{\prime}\right)$ (33)
$\displaystyle-2i\mathcal{A}_{\bar{\sigma}}\left(p^{\prime}\right)\Big{)}.$

The full theory is then defined by the self-energies Eqs. 30, 31, 32 and 33
dressing the bare propagators of fermions and pairing field and the
corresponding Dyson equations in Eqs. 17 and 22.

##  IV Numerical implementation

Having formulated the theoretical framework, we now focus on the numerical
solution of the coupled self-consistent integral equations for the
propagators.

As the retarded and Keldysh propagators are simultaneously diagonal in
$p$-space, the inversion of the Dyson equations and the computation of the
Keldysh propagator are straightforward. The main challenge lies in computing
the self-energies self-consistently. It is possible to do this efficiently and
accurately thanks to the convolution nature of Eqs. 30, 31, 32 and 33. This
means that the self-energies can be obtained by transforming the propagators
into the position and time domain then computing their products and
transforming these back to momentum and frequency space.

However, performing these transformations accurately is difficult even in
imaginary time, because the correlation functions decay slowly at high momenta
and frequencies. For instance, the momentum distribution, according to the Tan
energy theorem [11], acquires the form
$n(\mathbf{k})=\mathcal{C}/\mathbf{k}^{4}$. This implies that the momentum-
integrated RF spectrum $I(\omega)=i\int_{\mathbf{k}}\delta
G^{K}(\mathbf{k},\mathbf{k}^{2}/2m-\mu-\omega)/2$ scales like $\omega^{-3/2}$
for large frequencies in case of negligible final-state interactions [61] in
agreement with the experimental observations [19]. Furthermore, the fermionic
spectral functions themselves show power-law asymptotics, which are given
below in Eqs. 46 and 47. Even in imaginary time, it turns out that using a
fast Fourier transform and hence an equidistant grid, to discretize
propagators, is unfeasible. Instead it has been shown that using adaptive non-
equidistant grids one can solve the imaginary time problem with high accuracy
[16, 48].

In a real-time formulation the need for long grids is much more challenging
than in imaginary time. This is because the high-momentum behaviour of the
propagators is quickly oscillating with a frequency that grows linearly with
the momentum. This so-called chirped oscillation creates the predicament that
a direct implementation requires an increasingly dense grid for larger
momentum. In real-time one therefore has to use long grids which are
simultaneously dense, making the direct implementation intractable.

Here we solve this problem by constructing a method which, through a series of
different interpolations, analytically separates out the chirped oscillations
such that numerical sampling of the latter is completely avoided. As a result,
similar to the calculation on Matsubara frequencies, one can use grids that
sparsely sample high momenta such that the power-law tails of the correlations
are accurately described. The method, therefore, allows one to compute
convolutions both efficiently and accurately, such that the self-consistent
T-matrix can be computed directly in real time without NAC.

To understand the procedure it is insightful to consider the simplest case
where a chirped oscillation appears, namely the transformation of the bare
retarded propagator from $p$ to $x$

$G^{R}_{0}(x)=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{\exp(ip\cdot
x)}{\omega+\mu-\frac{\textbf{k}^{2}}{2m}+i0^{+}}.$ (34)

The frequency integral can be computed by contour integration and leads to the
propagator in $(\mathbf{k},t)$ space

$G^{R}_{0}(\bm{\mathbf{k}},t)=-i\theta(t)\exp[\left(-i\frac{\textbf{k}^{2}}{2m}+i\mu-0^{+}\right)t],$
(35)

which, due to isotropy only depends on $\absolutevalue{\textbf{k}}$. The first
term in the argument of the exponential in Eq. 35 leads to oscillations with a
linear chirp in momentum, which is exceedingly hard to sample on a finite
numerical grid for larger values of $\absolutevalue{\textbf{k}}$. However, for
the bare propagator, the momentum integral can be performed analytically

$G^{R}_{0}(\textbf{r},t)=\theta(t)\exp\left(i\frac{\textbf{r}^{2}m}{2t}+(i\mu-0^{+})t\right)\\!\left(\frac{im}{2\pi
t}\right)^{\\!3/2},$ (36)

where the oscillations due to $\textbf{k}^{2}t$ translates into oscillations
with an argument $\textbf{r}^{2}/t$. This is equally hard to sample and both
the oscillations and the entire propagator diverge at $t=0$. The challenge to
be solved thus lies in the efficient and accurate computation of Fourier
transforms without relying on sampling these fast oscillations. To explain our
method we will consider the transformation $p\rightarrow x$ and discuss its
relation to the inverse transformation at the end of the section.

Our method to perform the transformations relies on the observation that, in
the ultraviolet, propagators are only weakly modified by the effects of a
finite density. Hence, one can expect well-defined quasiparticles with a
dispersion mostly determined by the bare mass, which in the case of the pair
propagator is $m_{\Gamma}=m_{\uparrow}+m_{\downarrow}$. More formally, we
require that, by changing frequency coordinates to the bare dispersion, the
ultraviolet parts of functions are almost constant. To this extent, the
frequency coordinate is transformed according to
$\omega_{\alpha}=\omega+\textbf{k}^{2}v_{\alpha}+\mu_{\alpha}$, where for the
transformation to real-space we have $v_{\alpha}=1/2m_{\alpha}$ with
$\alpha\in\\{\Gamma,\uparrow,\downarrow\\}$. On the transformed coordinates
our functions $f_{\alpha}(\textbf{k},\omega_{\alpha})$ are now slowly varying
as a function of momentum as illustrated in Fig. 3.

Figure 3: Typical spectral function of an isotropic strongly interacting
system with both broadening and a non-trivial structure at small momenta. The
bare quadratic dispersion leads to a large change of the energy at large
momentum in the original coordinate system shown in a). Performing the
coordinate transformation makes the energy constant for large momentum as
illustrated in b).

Additionally, for long-lived quasiparticles or near the critical point the
functions can become very strongly peaked in $\omega_{\alpha}$ which, due to
the frequency transformation, can be solved by locally and adaptively
increasing the density of the $\omega_{\alpha}$-grid. Having denser grids only
locally is possible, because $f_{\alpha}$ varies slowly with $k$, such that
the same adaptations to the frequency grid are useful across a large range of
momenta. Performing a numerical Fourier transform on an irregular grid, with
$N$ points, scales as $N^{2}$ compared to $N\log\left(N\right)$ for an
equidistant or logarithmic grid [62, 63, 64]. To mitigate having to use
smaller grids, the temporal Fourier transform is done by using a spline
interpolation [16, 48]. For a system with rotational invariance, the angular
directions of the momentum integral can be solved exactly which leaves a one-
dimensional integral for each time (space) grid point $t_{j}$ ($r_{i}$), which
for a three dimensional system takes the form

$\displaystyle f_{\alpha}(r_{i},t_{j})=$
$\displaystyle\text{e}^{i\mu_{\alpha}t_{j}}\int_{0}^{\infty}\frac{\differential
k}{(2\pi)^{2}}\frac{2k\sin\left(kr_{i}\right)}{r_{i}}$ (37)
$\displaystyle\times f_{\alpha}(k,t_{j})\exp(-iv_{\alpha}k^{2}t_{j}),$

where from now on a notation convenient for an isotropic system
$k=\absolutevalue{\textbf{k}}$ and $r=\absolutevalue{\textbf{r}}$ will be
used. As the fast oscillation $v_{\alpha}k^{2}$ has been extracted by the
transformation to $\omega_{\alpha}$, $f_{\alpha}$ can be interpolated with a
spline of order $M$. The $\sin\left(kr_{i}\right)$ factor is not included in
the interpolation for two reasons: The first is that it oscillates quickly for
large positions and therefore requires a dense $k$-grid for it to be
efficiently interpolated with a spline. Secondly, it would lead to the spline
coefficient being a higher-rank tensor and would greatly increase the memory
requirements. To circumvent these issues we only do a low-order spline (cubic
or quintic) on $kf_{\alpha}$ given by

$kf_{\alpha}(k,t_{j})\approx\sum_{l=0}^{M}a^{l}_{\alpha;n,j}(k-k_{n})^{l},\text{
where }k_{n}\leq k<k_{n+1},$ (38)

leading to the real-space representation

$f_{\alpha}(r_{i},t_{j})\approx\sum_{n=0}^{N_{k}-1}\sum_{l=0}^{M}a^{l}_{\alpha;n,j}W_{\alpha;i,j,n}^{l},$
(39)

where $N_{k}$ is the number of grid points in the momentum grid and the tensor
$W$ is given by

$W_{\alpha;i,j,n}^{l}=\int_{k_{n}}^{k_{n+1}}\frac{\differential
k}{\left(2\pi\right)^{2}}\frac{2\sin(kr_{i})}{r_{i}}\text{e}^{-iv_{\alpha}t_{j}k^{2}}\left(k-k_{n}\right)^{l}.$
(40)

To make the calculations tractable this tensor has to be factorized. To keep
the evaluation of the $k^{2}$ oscillation exact we factorize
$\sin(kr_{i})/r_{i}$ by using a 2-point Hermite interpolation to order $P$.
Such an interpolation guarantees that both the function values and the first
$P-1$ derivatives are correct at every point. It leads to the approximation
[65]

$2\sin(kr_{i})/r_{i}\approx\sum_{p=0}^{2P-1}b_{i,n}^{p}\left(k-k_{n}\right)^{p-c(p)}\left(k-k_{n+1}\right)^{c(p)},$
(41)

for $k_{n}\leq k<k_{n+1}$. Here $c(p)=\lfloor
p/(2P)\rfloor\left(p\%(2P-1)\right)$, where $\lfloor\cdot\rfloor$ and $\%$ are
the floor and modulo operators respectively. To ensure that the oscillations
for large values of $r$ are correctly captured, while still keeping the
density of the $k$-grid manageable, one can choose $P$ based on an upper bound
of the error discussed in appendix B. The coefficients $b_{i,n}^{p}$ can be
efficiently computed through the generalized divided differences [66] and are
independent of the frequency grid. With this factorization, the tensor $W$
takes the form

$W^{l}_{\alpha;i,j,n}=\sum_{p=0}^{2P-1}b^{p}_{i,n}I^{l+p}_{\alpha,n,j},$ (42)

where the elements in $I^{l+p}_{\alpha,n,j}$ only depend on the index $\alpha$
describing the quickly oscillating factor of the function to be transformed.
Through these interpolations, we are no longer relying on sampling linear
chirped oscillations because the elements

$\displaystyle
I^{l+p}_{\alpha,n,j}=\sum_{q=0}^{c(p)}z_{n,p,q}\text{e}^{-iv_{\alpha}t_{j}k_{n}^{2}}\int_{0}^{\Delta
k_{n}}\frac{\differential k}{\left(2\pi\right)^{2}}$
$\displaystyle\text{e}^{-iv_{\alpha}t_{j}\left(k^{2}-2kk_{n}\right)}$ (43)
$\displaystyle\times$ $\displaystyle k^{l+p-c(p)+q},$

can be computed exactly. The integrals depend on the time/momentum grids and
the bare mass, while the coefficient

$z_{n,p,q}=(-1)^{c(p)-q}\matrixquantity(c(p)\\\ q)\Delta k_{n}^{c(p)-q},$ (44)

only depend on the grid-spacing in momentum $\Delta k_{n}=k_{n+1}-k_{n}$. With
these elements the real-space representation is

$f_{\alpha}(r_{i},t_{j})\approx\sum_{p=0}^{2P-1}\sum_{l=0}^{M}\bigg{(}\mathbf{B}^{p}\cdot\Big{(}\mathbf{I}^{p+l}_{\alpha}\circ\left(\mathbf{A}^{l}\cdot\mathbf{f}_{\alpha}\right)\Big{)}\bigg{)}_{i,j},$
(45)

where $\cdot$ is a matrix contraction and $\circ$ is the element-wise
(Hadamard) product. The elements of the matrix $\textbf{B}^{p}$ are
$b^{p}_{i,n}$, while the matrix $\textbf{I}^{p+l}_{\alpha}$ is build up of
$I_{\alpha,n,j}^{p+l}$. The matrix $\textbf{A}^{l}$ is independent of $\alpha$
and leads to the spline coefficients $a^{l}_{\alpha;n,j}$ in Eq. 39 once
contracted with the matrix built by
$\left(\textbf{f}\right)_{\alpha;m,j}=f_{\alpha}(k_{m},t_{j})$. The first
contraction scales as $M\times N_{k}^{2}\times N_{t}$, whereas the Hadamard
step scales as $(2P+M)\times N_{k}\times N_{t}$. The last contraction scales
with $2P\times N_{r}\times N_{k}\times N_{t}$. As $N_{i}\gg M,P$ this
transformation scales the same way as a naive DFT apart from an overhead of
$2P$, or $M$ depending on the relations between the grid lengths. We point out
that, since $f(k,\omega_{\alpha})$ is slowly varying, $M=3$ is sufficient for
the present application. Furthermore, typical values for $P$, as determined by
the condition in App. B, fall between $9$ and $18$ for momentum grids with a
maximum length of $k\sqrt{\beta}=50$.

An important aspect of the transformation procedure is that $\textbf{B}^{n}$,
$\textbf{I}^{n}_{\alpha}$ and $\textbf{A}^{n}$ only depend on the $r,k$ and
$t$ grids and $v_{\alpha}$, which means that all tensors can be precomputed at
initialization and are unaffected by the adaptive frequency grid. To reliably
transform all functions to $r$ we need three different coordinate
transformations with the choices
$v_{\alpha}\in\\{1/m_{\Gamma},\,1/m_{\uparrow},\,1/m_{\downarrow}\\}/2$. For
transforming from $r$ to $k$, the transformation is identical apart from the
substitutions $r\leftrightarrow k$, $v_{\alpha}t\rightarrow-1/(v_{\alpha}t)$
and the factor of $\left(2\pi\right)^{2}$ in $\textbf{I}_{\alpha}^{n}$.

###  IV.1 Comparison to conventional FFT

The alternative implementation using a naive fast Fourier transform (FFT)
would require one to be able to sample the fast oscillations from the bare
dispersions. For the transformation to space, one has to sample the
oscillating factor $\text{e}^{-ik^{2}t}$, where for simplicity we have chosen
units with $m=1/2$. By expanding the argument around the largest $k$ and $t$
value it becomes clear that one needs to be able to sample an oscillation with
a frequency of $2k_{max}t_{max}$. To satisfy the Nyquist criterion the
momentum grid has to have a spacing smaller than
$\Delta_{k}=\pi/(2k_{max}t_{max})$. Using a conventional FFT such a grid
spacing requires that $r_{max}=2\pi/\Delta_{k}$. However, the same applies to
the inverse transformation from position to momentum space. There one has to
be able to sample oscillations of the form $\text{e}^{ir^{2}/t}$. Again,
equidistant sampling requires $N_{\text{fft}}=\lceil
2r_{max}^{2}/(t_{min}\pi)\rceil$ grid points, where $\lceil\cdot\rceil$ is the
ceiling operator. Imposing the constraint from the momentum grid one finds
that the necessary grid length is $N_{\text{fft}}=\lceil
32k_{max}^{2}t_{max}^{2}/(t_{min}\pi)\rceil$. For each discrete time value,
this implies $\mathcal{O}(N_{\text{fft}}\log(N_{\text{fft}}))$ operations
compared to $\mathcal{O}(2PN_{k}^{2})$ operations for the method presented in
the previous section. Because we do not have to sample the fast oscillations
and because the $k$ and $r$ grid are only connected by the range they cover,
since $r_{\text{min/max}}\sim\left(k_{\text{max/min}}\right)^{-1}$, we can
work with around 500 grid points. For a typical calculation we use values
$k_{max}\sqrt{\beta}=50,\,r_{max}/\sqrt{\beta}=40,\,t_{min}/\beta=2\times
10^{-5}$ and $t_{max}/\beta=50$ which is accurately sampled with $P=18$ as
discussed in App. B. By comparison, a conventional FFT-based method would
require around $10^{9}$ times more grid points and consequently $10^{5}$ times
the number of operations. The overhead due to the interpolations and the
$N^{2}$ scaling of the transforms is thus vastly overcompensated by the
decrease in the grid lengths. The result is that our method allows us to
compute self-energies efficiently by relying on their convolution structure,
without using grids with a prohibitively large number of points.

To further increase the efficiency of the implementation, it is useful to
subtract the analytically known UV behavior, while also preventing the
appearance of singularities in the IR. A detailed account of the analytical
subtractions and how to apply the transformations is presented in appendix A.

##  V Validation

To test the accuracy of our method, we are now going to compare it with known
analytical results, imaginary-time calculations, and state-of-the-art
numerical analytic continuation. For this comparison, we focus on the balanced
case where $n_{\sigma}=n$ and $m_{\sigma}=m$ and will therefore drop the spin-
dependence of all fermion propagators.

###  V.1 Asymptotics

One of the main challenges of the real-frequency formalism lies in the
accurate description of the asymptotic behavior of propagators at high
momenta, as well as $|\omega|\to\infty$. We focus on quantities that are fully
determined by the effect of interactions, i.e. that vanish for the non-
interacting gas, as these will provide key insights into the precision of our
method.

At high energies, the fermionic spectral function exhibits a scattering
continuum with the analytically known frequency dependence [27]

$\displaystyle\mathcal{A}(k,\omega\to\infty)=4\pi nm^{-3/2}\omega^{-5/2}\,,$
(46)

where $n$ is the density of either species of fermions. Similarly, at low
frequencies, it can be shown to decay as

$\displaystyle\mathcal{A}(k,\omega\to-\infty)\sim(-\omega)^{-9/2}\,,$ (47)

with a prefactor that has no closed analytical form. In the limit of weak
coupling the latter is determined by including the behavior at high
frequencies stated in Eq. (46) in the evaluation of the fermonic self-energy.
There one finds $4\pi nm^{-7/2}\mathcal{C}$ for the prefactor. At stronger
coupling this result is modified self-consistently by a resummation to
infinite order in perturbation theory. As we show in the inset of Fig. 4, both
behaviors are very well captured by our method.

The Tan contact density $\mathcal{C}$ can be calculated from the pair
propagator $\mathcal{C}=im^{2}\delta\Gamma^{K}(r=0,t=0)/2$ or the asymptotic
decay of the momentum distribution [11]
$\lim_{k\to\infty}n(k)=\mathcal{C}/k^{4}$. For a $\Phi$-derivable theory, both
calculations must coincide and therefore provide a critical test for the
convergence and accuracy of our method. In Fig. 4, the momentum distribution
of the balanced, unitary Fermi gas at $T/T_{F}=0.2$ is compared with the
asymptotic behavior $\mathcal{C}/k^{-4}$ with the contact obtained from the
pair propagator. The excellent agreement indicates two important observations.
On one hand, we see that the self-consistent evaluation is indeed well-
converged. On the other hand, the clean asymptotics of $n(k)$, i.e. the
absence of any numerical noise in Fig. 4, demonstrates the high fidelity, with
which self-energies are evaluated, even at high momenta $k\gg k_{F}$ and
frequencies $|\omega|\gg\epsilon_{F}$.

Figure 4: The fermionic spectral function $\mathcal{A}(0,\omega)$ is known to
decay as $\omega^{-5/2}$ and $(-\omega)^{-9/2}$ at very high and low
frequencies respectively. Both limits are well reproduced by our method as is
shown in the inset. Furthermore, the momentum distribution at high momenta
behaves as $n(k\gg k_{F})\approx\mathcal{C}/k^{4}$ with $\mathcal{C}$ the Tan
contact density obtained from the pair propagator.

###  V.2 Comparison with Matsubara formalism

Additional insights can be gained from the comparison with results obtained in
Matsubara formalism. Indeed, there, the quickly oscillating factor
$e^{ik^{2}t}$ of the real-time propagators is replaced by a Gaussian envelope
$\sim e^{-k^{2}\tau}$ in the imaginary time interval $\tau\in[0,\beta)$ and
therefore much easier to handle. Consequently, the self-consistent T-matrix in
Matsubara formalism is a well-established method [15, 16] with well-controlled
errors [67]. We will leverage this to compare our new approach against this
method. For the dimensionless Fermi energy $\epsilon_{F}/T$ and the contact,
we find typical differences no larger than $3\%$. However, these integrated
quantities provide little new insight over the previous tests. Instead, a more
telling comparison can be performed between propagators directly in imaginary
time.

Although the analytic continuation from Matsubara frequencies to real
frequencies or the equivalent transformation from imaginary to real-time is
numerically an ill-posed problem, the reverse operations are unproblematic.
Imaginary time propagators are obtained from spectral functions by the
generalized Laplace transform

$\displaystyle G(k,\tau)=\int\frac{d\omega}{2\pi}\frac{e^{-\omega\tau}}{1\mp
e^{-\beta\omega}}\mathcal{A}(k,\omega)\,,$ (48)

where the upper/lower sign refers to bosons/fermions. It is well known that
the spectrum of the latter transformation possesses only a few singular values
larger than any realistic numerical resolution, which immediately explains the
problems one encounters when inverting the transformation as is the case for
numerical analytic continuations. However, it also highlights that propagators
in real frequencies contain genuinely more information than their Matsubara
counterparts.

Nevertheless, the thermodynamic information contained in the imaginary time
propagators needs to be fully recovered by the presented formalism. To this
extent, both the propagators of the fermions and the pair propagator are
transformed to imaginary times and compared in Figs. 5 and 6 with results
previously obtained in the Matsubara formalism [48] at zero momentum. Both
agree for all values of $\tau/\beta$ to within $1\%$. Furthermore, we show the
integrated relative difference

$\displaystyle\epsilon_{\text{rel}}(k)=\sqrt{\frac{\int
d\tau(G_{\text{imag}}(k,\tau)-G_{\text{Keld}}(k,\tau))^{2}}{\int d\tau
G_{\text{imag}}^{2}(k,\tau)}}$ (49)

as a function of momentum in Fig. 7 to highlight that this level of agreement
is global and no larger differences occur at higher momenta.

We emphasize, that although the Matsubara formalism contains less information
than the method developed here, the comparison is even more critical than the
previous inspection of the asymptotics. This is because the comparison is
sensitive to local quantities, like the linewidth or energy of the interacting
fermions.

Figure 5: Absolute value of fermionic imaginary time propagator at vanishing
momentum $k=0$ (black) compared with the absolute value of the difference
between its evaluation directly in imaginary times and via the generalized
Laplace transform of the spectral function (blue). Figure 6: As for the
fermionic propagator in Fig. 5, the difference between the pairing propagator
in imaginary times obtained from real frequencies and in Matsubara formalism
(blue) is small compared to the absolute value of the paring propagator
(black). Figure 7: The integrated, relative difference between single-
particle, imaginary-time propagators obtained in Matsubara and Keldysh
formalism as defined in Eq. (49) shows very good agreement between the two
independent numerical methods.

##  VI Comparison with numerical analytic continuation

With the numerical implementation of the self-consistent T-matrix
approximation validated, it is now time to turn the tables and review the
accuracy of the numerical analytic continuation as the current state-of-the-
art for the evaluation of spectral functions in strongly interacting systems.

Returning to the exemplary case of the unitary Fermi gas at $T/T_{F}=0.2$, we
calculate the Matsubara propagator in imaginary time as reported in Ref. [48].
We then use Bryan’s algorithm [68] for the maximum entropy method [25, 69] to
perform the analytic continuation (see App. C). The maximum entropy method
uses Bayesian inference to complement the limited information that can be
gained from imaginary times by an entropy term $S[\mathcal{A}_{0}]$ relative
to a default model $\mathcal{A}_{0}$. The latter is chosen as a slowly varying
function with the analytically known behavior for large values of $|\omega|$
stated in Eqs. 46 and 47. The relative weight between the constraint in
imaginary times and the entropy is a free parameter, the choice of which can
be optimized in different ways, leading to somewhat different spectral
functions. The details of this choice and the method behind the numerically
continued spectral function $\mathcal{A}_{\text{NAC}}$ in the comparison in
Fig. 8 are given in App. C. As is shown in said figure for $k=0$, spectral
functions obtained from numerical analytic continuation are in general too
broad and typically shifted to higher frequencies, with the scattering
continuum undersized when compared to
$\mathcal{A}_{\text{Keldysh}}(k,\omega)$. These properties are expected since
the default model has a single broad maximum. A pronounced gap between the two
peaks of the spectrum frustrates the entropy term and is therefore
underestimated with the two maxima broadened and drawn too close together by
the maximum entropy method. We emphasize that although the qualitative
features of the spectrum can be gleaned from $\mathcal{A}_{\text{NAC}}$, it is
quantitatively a bad approximation to $\mathcal{A}_{\text{Keldysh}}$.

Figure 8: Comparison between single particle spectral functions at
$\bm{\mathbf{k}}=0$. $\mathcal{A}_{\text{Keldysh}}$ denotes the spectrum
obtained by the method developed here and $\mathcal{A}_{\text{NAC}}$ the
result from numerical analytic continuation explained in Sec. VI. The spectrum
$\mathcal{A}_{\text{mixed}}$, which agrees almost perfectly with
$\mathcal{A}_{\text{Keldysh}}$ is obtained by analytic continuation with
$\mathcal{A}_{\text{Keldysh}}$ as default model, indicating that the latter is
fully consistent with the imaginary time propagator used for
$\mathcal{A}_{\text{NAC}}$.

Given the large discrepancy between the two spectral functions, we also used
$\mathcal{A}_{\text{Keldysh}}$ as the default model (i.e.
$S[\mathcal{A}_{\text{Keldysh}}]$) to confirm that no additional error sources
exist in our implementation of the maximum entropy method. If the default
model is consistent with the Matsubara propagator in imaginary time, the
Lagrange parameter $\alpha$ has no effect and the maximum entropy spectrum is
identical to the default model. Indeed, the spectrum obtained by the same
method as $\mathcal{A}_{\text{NAC}}$ but with $\mathcal{A}_{\text{Keldysh}}$
used for the default model is shown in Fig. 8 labeled
$\mathcal{A}_{\text{mixed}}$ and agrees exceptionally well with the spectrum
obtained by the newly developed method.

We conclude that although the imaginary time results obtained in Matsubara
formalism are consistent with the new method, the common maximum entropy
method for numerical analytic continuation is inaccurate for the strongly
interacting Fermi gas at low temperatures.

##  VII Spectral function near criticality

The large deviation between the validated method presented here, and the
current standard in the form of the numerical analytic continuation,
necessitates a review of previous numerical results for the spectral function
of the strongly interacting Fermi gas.

Figure 9: Spectral functions at criticality for different values of the
scattering length $(k_{F}a)^{-1}=\\{-1,0,1\\}$ where the effect of the
interactions increases from left to right. The critical temperatures as
determined by our methods from left to right are
$T_{c}/T_{F}=\\{0.068,0.156,0.205\\}$ and the chemical potential
$\mu/\epsilon_{F}=\\{0.743,0.407,-0.803\\}$ is indicated by the white dashed
line. Both agree with the values obtained in imaginary time.

We show the results for the spectral functions at the superfluid transition at
different scattering lengths in Fig. 9. We mention that the corresponding pair
susceptibility $\Gamma^{R}(k,\omega)$ is calculated with similar accuracy as
shown in App. D. For $k_{F}a=-1$ the effect of interactions is limited to a
slight reduction of the spectral weight in the immediate vicinity of the Fermi
surface and a broadening of the dispersion deep inside the Fermi sea. The
latter describes the short lifetime of holes deep in the Fermi sea as no
effects of artificial broadening, discussed in App. A.2, are visible in Fig.
9. Even at the superfluid transition of the unitary gas the dispersion in Fig.
9(b) shows only very weak signs of the nearby instability. There is however a
weak scattering continuum visible at small momenta with its peak in reasonable
agreement with the energy expected for Bogoliubov quasiparticles

$\displaystyle
E^{\pm}(k)=\mu\pm\sqrt{\left(\frac{k^{2}}{2m^{*}}-\mu+U\right)^{2}+\Delta^{2}}\,,$
(50)

where the pseudogap $\Delta$, chemical potential $\mu$ and Hartree shift $U$
are fit parameters [70].

As the scattering length is increased further, the gap in the single-particle
spectrum finally opens as the low-energy effective theory transitions from
interacting fermions to weakly coupled bosons formed by deeply bound pairs of
fermions. In agreement with this picture, the hole spectrum, i.e. the occupied
part of the spectrum below the chemical potential now has its maximum at
$k=0$. Interestingly, the part of the spectrum above the chemical potential,
known as the particle spectrum, at small momenta consists of a clear
quasiparticle peak merged with the broad background formed by the scattering
continuum. We emphasize that the narrow peak in the particle spectrum in Fig.
9(c) is necessitated by the large gap between the particle and hole spectrum.
With the latter several times larger than the temperature, there are very few
thermal excitations. The minimum of the upper branch in Fig. 9(c) is therefore
a good approximation to the lowest state available to a single particle added
to the system. Because there are no empty states with lower energy to decay to
and few thermal excitations to scatter with, the state must be long-lived.
This property is not well recovered in the analytic continuation [27], which
otherwise captures the qualitative properties of the spectra shown in Fig. 9.

Although not visible for the small momenta shown in Fig. 9, in all cases the
linewidth at high momenta decreases as
$\Im{\Sigma^{R}(k,\omega=k^{2}/(2m))}\sim k^{-1}$, in agreement with the
analytic result [17]. This reflects the vanishing role interactions play for
particles with high kinetic energy exploited in our numerical implementation
in Sec. IV.

The argument above regarding the lifetime of the particle spectrum illustrates
the importance of an ordering principle in the energy scales. Especially in
the highly correlated, but disordered state at temperatures above, but close
to $T_{c}$, such simplifications are highly desirable due to the strongly
interacting nature of the system. One therefore argues that already above the
critical temperature, fermions near the Fermi surface pair up such that these
pairs eventually condense at $T=T_{c}$. As breaking these pairs costs energy,
one may expect a suppression of the single-particle density of states

$\displaystyle\rho(\omega)=\int\frac{d^{3}k}{(2\pi)^{3}}\mathcal{A}(k,\omega)$
(51)

already above the critical temperature [71]. The existence of this so-called
pseudogap has been observed for example in underdoped high-$T_{c}$
superconductors, where its origin is, however, hard to interpret [72] and in
2D systems [22]. This has sparked both theoretical [73] and experimental [23]
interest in a potential pseudogap phase in the considerably simpler 3D unitary
Fermi gas. One of the key findings obtained from the numerical analytic
continuation of spectra in the self-consistent T-matrix approximation near the
critical point was a much weaker pseudogap [27] than previously reported using
auxiliary field quantum Monte Carlo [28] and non-selfconsistent approximations
[74]. However, all of these methods rely on numerical analytic continuation,
which raises doubts about their validity.

Figure 10: The single-particle density of states of the unitary Fermi gas
exhibits a very weak pseudogap only for temperatures $T\lesssim 0.19\,T_{F}$,
which is close to the critical temperature $T_{c}=0.156\,T_{F}$.

We reinvestigate the issue of the pseudogap in the normal phase of the unitary
Fermi gas and report the results in Fig. 10. It is found that the unitary gas
develops a weak pseudogap at temperatures below $T\approx 0.19\,T_{F}$.
However, even at the critical temperature $T_{c}=0.156\,T_{F}$ the signature
is not very pronounced as the abundance of thermal fluctuations broadens the
spectrum, thereby precluding the possibility of a marked pseudogap. Our
results, therefore, qualitatively confirm those previously reported using the
maximum entropy method [27, 3] and Padé approximation [30]. Moreover, the
present controlled analysis settles the quantitative discrepancy between these
previous results.

Although the underlying self-consistent T-matrix approximation is
uncontrolled, its good agreement with experimentally observed thermodynamic
properties near unitarity [16, 48, 9] lends credibility also to the spectral
properties reported here.

##  VIII Conclusions

We have demonstrated an efficient method to calculate the spectral function of
the 3D Fermi gas in the self-consistent T-matrix approximation directly in
real frequencies. The extension to the spin- or mass-imbalanced setups is
straightforward by including a spin-dependent chemical potential
$\mu_{\sigma}$ or mass $m_{\sigma}$ as already shown in Eq. (7). Also the
symmetry-broken state has been accessed with the self-consistent T-matrix
[16], which requires the translation to the Keldysh formalism. Furthermore,
the theory can be directly applied to two dimensions.

We emphasize that the underlying method relies only on the convolution theorem
and is therefore rather versatile. Natural applications are Bose-Fermi models,
which appear frequently in other fields of condensed matter. For instance,
metallic quantum criticality is typically described by a Fermi surface coupled
to a bosonic order parameter, which describes instabilities of the former [75,
76]. In two dimensions, these systems commonly show non-Fermi liquid behavior
accompanied by Landau-damped critical modes, which are both primarily observed
in the fermionic and bosonic spectral functions. Moreover, for the well-
established Ising-nematic and spin-fermion models, as well as in the context
of a Kondo heterostructure, it was shown that self-consistent 2PI quantum
field theory reproduces imaginary-time quantum Monte Carlo data with good
accuracy [77, 78]. Since the corresponding field theory is based again on
convolutions, our method, which by construction has access to both the
fermionic and the bosonic spectra (for examples cf. App. D), seems a promising
route to obtain a more quantitative understanding of the spectral functions
and the identification of different scaling regimes indicated by analytic
calculations [79].

Another possible application is the quantitative analysis of quench dynamics
in cold atom experiments [80, 81, 82, 83], which is well described by 2PI
field theory [43]. Notably, the method presented here is directly applicable
since both the convolution structure and the conserving nature of the
approximation are preserved under Wigner expansion [84]. Finally, the quickly
growing field of pump-probe experiments in solid-state materials [85] can be
treated analogously [86]. There, the order-parameter dynamics is again
described by a 2PI effective field theory with dissipation. Our approach will
be useful to develop a deeper understanding of long-lived states and scaling
dynamics far from equilibrium [87, 88, 89].

We have presented a detailed validation of our method on the example of the
Fermi gas near unitarity and outlined its applicability to a broad range of
structurally similar, strongly interacting systems. This makes our approach a
prime tool for the investigation of dynamical quantities without the need to
resort to uncontrolled numerical methods.

An important future step will be the extension of the method to less
restrictive approximations to the Luttinger-Ward functional, including for
example higher-body correlations [90, 91, 92] that lack an efficient
representation in terms of convolutions.

A different approach to the self-consistent spectral function for the unitary
Fermi gas has recently been based on Monte Carlo integration [93]. Although
significantly less efficient than the method presented here, its flexibility
may prove useful for the extension to other classes of diagrams.

##  IX Acknowledgements

J.L. would like to thank Eugen Dizer for many stimulating discussions that
prompted us to write this paper. We also thank Sebastian Diehl for comments on
the manuscript. J.L. acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) CRC 1238 project
number 277146847. C.H.J acknowledges support by the project No.
2021/43/P/ST2/02911 co-funded by the National Science Centre and the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement no. 945339. B.F. acknowledges support by the
DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and
Topology in Quantum Matter – ct.qmat (EXC 2147, project id 390858490). For the
purpose of Open Access, the authors have applied a CC-BY public copyright
licence to any Author Accepted Manuscript (AAM) version arising from this
submission.

##  Appendix A Details of the self-energy computation

To optimize the accuracy of the loop integrals involved in the evaluation of
the diagrammatic expressions that have been derived in Sec. III, it is
important to ensure that the different features of terms in the self-energies
are catered for. As we detail below, we employ several different techniques in
this regard, which either allow us to reduce the domains that need to be
covered numerically or improve the overall accuracy. Given their importance to
the performance claims we make in the main text, we provide a detailed account
of the essential procedures used.

###  A.1 Pair propagator

As stated in the main text, the contact nature of the interaction allows one
to write the retarded Bethe-Salpeter equation as a Dyson equation

$\Gamma^{R}(p)=\frac{1}{\frac{1}{g_{0}}-\Sigma_{\Gamma}^{R}(p)}.$ (52)

Including the Keldysh structure and taking advantage of the causality
structure the retarded self-energy is shown in Eq. 30. To remove the UV
divergence of $\Gamma$ arising from the contact interaction the self-energy is
split into four parts

$\displaystyle\Sigma_{\Gamma}^{R}(p)=\Sigma_{\Gamma,v}^{R}(p)+\Sigma_{\Gamma,v\delta}^{R}(p)+\Sigma_{\Gamma,\uparrow}^{R}(p)+\Sigma_{\Gamma,\downarrow}^{R}(p).$
(53)

The labeling here should be understood as the $v$ referring to a bare vacuum,
$\delta v$ to a dressed vacuum and $\uparrow$ ($\downarrow$) is a contribution
due to a finite density in the $\uparrow$ ($\downarrow$) species. The bare
vacuum contribution is

$\displaystyle\Sigma_{\Gamma,v}^{R}(p)=i\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}G^{R}_{\uparrow,0}(p-p^{\prime})G^{R}_{\downarrow,0}(p^{\prime}).$
(54)

Together with the contact interaction this is exactly equal to the 2-body
problem and using the Lippmann-Schwinger equation [94] one can relate this to
the 2-body T-matrix

$\frac{1}{g_{0}}-\Sigma_{\Gamma,v}^{R}(p)=T_{0}^{-1}\left(\omega-\frac{k^{2}}{2m_{\Gamma}}+\mu_{\Gamma}\right),$
(55)

where $m_{\Gamma}=m_{\uparrow}+m_{\downarrow}$ and
$\mu_{\Gamma}=\mu_{\uparrow}+\mu_{\downarrow}$. In three dimensions $T_{0}$
can also be computed from scattering theory as the 2-body scattering with a
$\delta$-potential [95]

$T_{0}^{-1}(\omega)=\frac{m_{r}}{\sqrt{2}\pi
a_{s}}-\frac{m_{r}^{3/2}\sqrt{-\omega}}{\sqrt{2}\pi},$ (56)

where $m_{r}=m_{\downarrow}m_{\uparrow}/m_{\Gamma}$ is the reduced mass and
$a_{s}$ is the scattering length. In this way, the bare coupling $g_{0}$ is
explicitly renormalized and connected to the experimental scattering length.
Using this procedure the full renormalized retarded pair propagator is given
by

$\Gamma^{R}(p)=\frac{1}{T_{0}^{-1}\left(\omega-\frac{k^{2}}{2m_{\Gamma}}+\mu_{\Gamma}\right)-\Sigma_{\Gamma,v\delta}^{R}-\Sigma_{\Gamma,\downarrow}^{R}-\Sigma_{\Gamma,\uparrow}^{R}},$
(57)

where the self-energy dependence on $p$ has been suppressed. The remaining
three self-energies can be read off from Eq. 30. One term results from
modifications of the fermion propagators due to interactions and is still a
vacuum term

$\displaystyle\Sigma_{\Gamma,v\delta}^{R}(p)=i\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}$
$\displaystyle\Big{(}G^{R}_{\uparrow}(p-p^{\prime})G^{R}_{\downarrow}(p^{\prime})$
(58)
$\displaystyle-G^{R}_{\uparrow,0}(p-p^{\prime})G^{R}_{\downarrow,0}(p^{\prime})\Big{)}.$

As opposed to Eq. 54, the integral is finite and requires no further
renormalization. To transform this term of the self-energy one has to remove a
fast oscillation with a mass given by $m_{\Gamma}$, such that the
transformation procedure in Sec. IV is used with $\alpha=\Gamma$.

The two remaining contributions to the self-energy are caused by finite
densities in either of the two species of fermions

$\displaystyle\Sigma_{\Gamma,\downarrow}^{R}(p)$
$\displaystyle=\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}G^{R}_{\uparrow}(p-p^{\prime})\delta
G^{K}_{\downarrow}(p^{\prime}),$ (59)
$\displaystyle\Sigma_{\Gamma,\uparrow}^{R}(p)$
$\displaystyle=\frac{i}{2}\int\frac{d^{4}p^{\prime}}{(2\pi)^{4}}\delta
G^{K}_{\uparrow}(p-p^{\prime})G^{R}_{\downarrow}(p^{\prime}).$

As $\delta G^{K}_{\alpha}$ is cut off at large momentum by the distribution
function, the self-energy $\Sigma^{R}_{\Gamma,\uparrow(\downarrow)}$ is
transformed using the appropriate fermion mass $m_{\uparrow(\downarrow)}$.
Next, we can use the FDR to minimize the noise of the imaginary part of
$\Sigma_{\Gamma,\uparrow(\downarrow)}^{R}$. This is possible because
$\delta\Sigma_{\Gamma}^{K}$ for bosons satisfies the relation

$\delta\Sigma_{\Gamma}^{K}(k,\omega)=-4in_{B}(\omega)\imaginary\Sigma_{\Gamma}^{R}(k,\omega).$
(60)

Instead of using the FDR one can also compute $\delta\Sigma_{\Gamma}^{K}$
directly using the loop integral in Eq. 31. This integral is well-behaved due
to the fast decay of the distribution functions. Hence we compute
$\delta\Sigma_{\Gamma}^{K}(k,\omega)$ with high accuracy and use it to remove
errors in the retarded self-energy through Eq. 60.

###  A.2 Fermions

As we have seen explicitly in the previous section, the pair propagator
depends on the self-energies of the fermions, which in turn depend on the pair
propagator. These self-consistent equations are solved iteratively, starting
from the non-interacting theory. The fermion self-energy has several different
terms with different features. Just as for the pair propagator it is necessary
to isolate these terms and optimize the transformation of each independently.
This is the focus of the first half of this section of the appendix.

In the latter half, we focus on artificial broadening. This is only necessary
to consider for the fermions as the pair propagator inherits its broadening
from the fermions, whose spectra in turn will be broadened by interactions
with the pairing field. To obtain accurate results, it is thus important to
minimize the artificial broadening required for the numerical Fourier
transformation of the propagators.

Considering first the optimal transformation of the different terms in the
self-energy, we start each iteration with all the propagators from the
previous iteration in the $(r,t)$ basis. In this space, the self-energies are
simple products and the retarded fermion self-energy consists of two distinct
parts. One arises from occupation of the paring field

$\Sigma^{R}_{\sigma,1}(r,t)=\frac{i}{2}G^{R}_{\bar{\sigma}}(r,-t)\delta\Gamma^{K}(r,t),$
(61)

which can be transformed to momentum by removing the bare oscillations in the
retarded propagator by choosing $m_{\alpha}=m_{\bar{\sigma}}$. The other
contribution to the self-energy is due to occupation in the
$\bar{\sigma}$-fermions

$\Sigma^{R}_{\sigma,2}(r,t)=\frac{i}{2}\delta
G^{K}_{\bar{\sigma}}(r,-t)\Gamma^{R}(r,t).$ (62)

To transform this self-energy we first subtract the bare T-matrix which
requires us to transform the term

$\Sigma^{R}_{\sigma,2,0}(r,t)=\frac{i}{2}\delta
G^{K}_{\bar{\sigma}}(r,-t)T_{0}(r,t).$ (63)

This term has oscillations at large momentum set by the bare pair propagator
such that we remove oscillations with a mass $m_{\alpha}=m_{\Gamma}$. The
remaining contribution

$\Sigma^{R}_{\sigma,2,\delta}(r,t)=\frac{i}{2}\delta
G^{K}_{\bar{\sigma}}(r,-t)\left(\Gamma^{R}(r,t)-T_{0}(r,t)\right),$ (64)

is generally varying significantly at short times and with a behavior not
related to the bare oscillations $m_{\alpha}=m_{\Gamma}$, which dominate at
larger times. For this reason, only the parts at times larger than a cutoff of
order $(10\beta)^{-1}$ are transformed using $\alpha\rightarrow\Gamma$, while
the short-time behavior is transformed without removing fast oscillations.
This cutoff is found to be fairly robust and doesn’t require fine-tuning
throughout the iterations. By adding the three contributions in Eqs. 61, 63
and 64 together one arrives at $\Sigma^{R}_{\sigma}(k,t)$.

Transforming from $(k,t)$ to $(k,\omega)$ is difficult due to a divergence at
$t=0$ in $T_{0}$ which scales as $t^{-1/2}$. To remove the errors generated by
this divergence, we fit $T_{0}(k,t=0^{+})$ to
$\Sigma^{R}_{\sigma,2,0}(k,t=0^{+})$. By numerically transforming
$\beta_{k}T_{0}(k,t)=\tilde{T}_{0}(k,\omega)$, with $\beta_{k}$ being the
fitting factor, the numerical errors due to the divergence can be identified
as $\tilde{T}_{0}(k,\omega)-\beta_{k}T_{0}(k,\omega)$. The procedure relies on
computing $T_{0}(k,t)$ with very high accuracy. This is possible and
computationally cheap because of the argument structure of
$T_{0}(\omega-k^{2}/(2m_{\Gamma})+\mu_{\Gamma})$, as seen in Eq. 55, where
momenta only enter via a shift of the frequency argument. Furthermore, this
calculation only has to be done once at initialization.

To decrease numerical errors in the spectral functions of the fermions one
can, as done in the pair propagator, use the FDR in Eq. 24. This is done by
directly computing

$\delta\Sigma_{\sigma}^{K}(r,t)=\frac{i}{2}\delta\Gamma^{K}(r,t)\left(\delta
G^{K}_{\bar{\sigma}}(r,-t)-i2\imaginary\Sigma_{\sigma,1}(r,t)\right).$ (65)

As discussed for the pair propagator, this self-energy can be computed with
high precision due to the fast decay of the distribution functions. After
subtracting the errors one arrives at the retarded self-energy for the
fermions.

Having discussed how to optimally transform the different terms in fermion
self-energy we now turn our attention to minimizing the effect of artificial
broadening. The artificial broadening is necessary because the bare fermion
propagators exhibit a $\delta$-peak in the spectral function. To sample this
numerical it must be made finite through artificial broadening. This
broadening is equally important in the time domain where it ensures that the
retarded propagators keep their causal structure and decay sufficiently within
the length of our numerical time grid. Such a broadening is always introduced
in the theory but should be taken to zero after the calculation is performed
[96]. As this is impossible with a numerical result, it is essential to
minimize the artificial broadening. First, we discuss how the artificial
broadening is included followed by how it is minimized throughout the
iterations.

Conventionally, this broadening is included as a constant and positive
imaginary contribution to the bare retarded propagator as

$\tilde{G}_{0,\eta}^{R}(\omega)=\frac{1}{\omega+i\eta},$ (66)

where the propagator has been written on the transformed frequency grid
discussed in Sec. IV. The corresponding Lorentzian spectral function decays as
$\omega^{-2}$. If such tails are included in the calculation, then they will
have a large impact on the high frequency tails of the spectral functions.
Through the convolutions, this then leads to significant errors, for example
on the Tan contact density extracted from the pair-propagator as well as total
fermion density. It is therefore essential to avoid these tails as much as
possible. To this extent, we add the finite broadening through an artificial
frequency-dependent retarded self-energy of the form

$\Sigma_{\rho}^{R}(\omega)=\frac{1}{\rho\omega-\frac{1}{\omega+i}},$ (67)

where $\rho$ is a positive computational parameter. The broadened bare
propagator is then given by

$G_{0,\eta}^{R}(\omega)=\frac{1}{\omega-\eta\Sigma^{R}_{\rho}(\omega)}.$ (68)

Because the artificial self-energy is retarded, it preserves the causality of
the propagator [42]. Differently from the standard case in Eq. 66, our
artificial broadening also contains a real part, which is necessary to satisfy
the Kramers-Kronig relations [45]. However, if one chooses the two artificial
parameters such that $\rho\gg 1\gg\eta$ and $\eta\rho\ll 1$ then the resulting
spectrum is very well described as

$A_{0,\eta}(\omega)=\frac{2\eta}{\eta^{2}+\omega^{2}+\rho^{2}\left(\omega^{4}+\omega^{6}\right)}.$
(69)

In practice, we satisfy the limit by choosing $\rho=10$ and $\eta\leq
10^{-2}$. With this choice the spectral function close to the bare dispersion
(meaning $\omega=0$), is almost Lorentzian but as one deviates from the bare
dispersion the bare spectral function now decays as $\omega^{-6}$. One can
find analytic expressions for the poles of Eq. 68, which allows us to compute
the Fourier transforms to $(r,t)$ analytically, ensuring no additional errors
are added in our subtraction schemes. Including the additional broadening in
this way makes it possible to keep the artificial width $\eta$, close to the
dispersion, large enough to allow for sampling of the peak without introducing
significant errors in the frequency tails of the spectral function.

Throughout the self-consistent iterations, the fermions acquire a natural
broadening due to the interactions. The result is that less artificial
broadening is needed to reliably perform the transformations and it is thus
preferable to remove all the superfluous artificial broadening at each
iteration. For this procedure, we first compute the broadening of the
quasiparticle pole for each momentum

$\eta_{0}(k)=-\imaginary\Sigma^{R}_{\sigma}(k,\epsilon_{0}(k)),$ (70)

where $\epsilon_{0}(k)$ satisfies

$\left(G^{R}_{0,\sigma}(k,\epsilon_{0}(k))\right)^{-1}-\real\Sigma^{R}_{\sigma}(k,\epsilon_{0}(k))=0.$
(71)

Whenever $\eta_{0}$ is greater than $\eta$ no artificial broadening has to be
added. In the case where $\eta_{0}<\eta$ one can identify the minimal
artificial broadening needed to ensure we can resolve the sharp features in
the spectrum. This minimal broadening $\eta_{m}$ is defined as

$\eta_{m}(k)=\theta\big{(}\eta-\eta_{0}(k)\big{)}(\eta-\eta_{0}(k)).$ (72)

We then compute the smallest imaginary value, $\eta_{s}$, over all $k$ and
$\omega$, of the self-energy when the minimum quasiparticle linewidth has been
renormalized

$\eta_{s}=2\,\text{Min}\,\absolutevalue{\imaginary\left(\Sigma^{R}_{\sigma}(k,\omega)+\eta_{m}(k)\Sigma_{\rho}^{R}\big{(}k^{2}v_{\alpha}+\mu_{\alpha}+\omega\big{)}\right)},$
(73)

where $\Sigma^{R}_{\rho}$ is defined in Eq. 67. By removing this constant and
adding $\eta_{m}(k)$ the maximum amount of artificial broadening has been
removed while retaining the causal structure of self-energy at all points in
$p$-space. This procedure gives a quantitative improvement but is not
essential for the stability of the method. Typical values of $\eta_{s}$ are on
the order of $10^{-3}/\beta$. The final retarded self-energy is then

$\Sigma^{R}_{\sigma}(k,\omega)\rightarrow\Sigma^{R}_{\sigma}(k,\omega)+i\eta_{s}+\eta_{m}\Sigma_{\rho}^{R}(k,\omega),$
(74)

and the retarded fermion propagator can be computed through the Dyson equation
Eq. 17. With this procedure, we have ensured that at each iteration the
smallest possible artificial broadening is added.

For imbalanced systems, one repeats all the same steps for the second fermion
species $\bar{\sigma}$.

##  Appendix B Interpolation order

The upper error bound between a function $f(x)$ and its $q$-point Hermite
interpolation of order $P$ at the point $x$, $h(x)$, is given by [97]

$\absolutevalue{f(x)-h(x)}\leq\absolutevalue{\frac{f^{(qP)}(a)}{qP!}\prod_{i=1}^{q}\left(x-x_{i}\right)^{P}},$
(75)

where $a\in\\{x_{1},x_{q}\\}$ is the value that maximizes $f^{(qP)}(a)$.
Specifying this for the 2-point case where $q=2$ and $f(x)=\sin(xb)/b$

$\displaystyle\max\absolutevalue{\frac{\sin(xb)}{b}-h(x)}=$
$\displaystyle\max\Bigg{|}\left(\frac{\differential^{2P}}{\differential
y^{2P}}\frac{\sin(yb)}{b}\right)_{y=a}\frac{1}{(2P)!}$ (76)
$\displaystyle\times\left(x-x_{n}\right)^{P}\left(x-x_{n+1}\right)^{P}\Bigg{|}$
$\displaystyle<$
$\displaystyle\max\left[\frac{b^{2P-1}}{(2P)!}\left(\frac{\Delta
x_{N}}{2}\right)^{2P}\right],$

where $\Delta x_{N}=x_{N+1}-x_{N}$ and it has been assumed that the largest
grid spacing is at the end of the grid. As the backward and forward
transformations require two different interpolations one can put an upper
limit on the interpolation error given by

$\epsilon=\max\left[\frac{r_{N}^{2P-1}}{(2P)!}\left(\frac{\Delta
k_{N}}{2}\right)^{2P},\frac{k_{N}^{2P-1}}{(2P)!}\left(\frac{\Delta
r_{N}}{2}\right)^{2P}\right].$ (77)

From this relation, one can choose the appropriate interpolation order. This
error gives an upper bound and therefore generally gives a larger than
necessary $P$.

##  Appendix C Maximum Entropy method

In this appendix, we provide the details of the maximum entropy method used
for the numerical analytic continuation in Sec. VI.

As explained in the main text, numerical propagators with a finite accuracy in
real frequencies contain significantly more information than their
counterparts in imaginary time or Matsubara frequency. Specifically, NAC aims
to minimize the dimensionless distance from the Matsubara propagator in
imaginary times

$\displaystyle|\chi(k)|^{2}=\sum_{i}\left(G(k,\tau_{i})-\int
d\omega\frac{e^{-\tau_{i}\omega}}{1+e^{-\beta\omega}}\mathcal{A}(k,\omega)\right)^{2}\,,$
(78)

where in our case the imaginary times $\tau_{i}$ are logarithmically spaced
near 0 and $\beta$. The local relative numerical accuracy of the Matsubara
propagator is $\sim 10^{-8}$. Consequently, upon discretization, only a few of
the eigenvectors of the integral transformation are indeed fixed by the
minimization of $|\chi(k)|^{2}\lesssim 10^{-8}$. Using an optimized grid with
400 points that take into account the known asymptotic decay of $\mathcal{A}$
at small and high frequencies and using cubic splines for the integration,
only around 80 of the eigenvectors of the integral transformation in (78)
provide meaningful information. Consequently, at the given accuracy an
infinite set of spectral functions $\mathcal{A}$ produce the same imaginary-
time propagator.

To resolve this ambivalence one needs to introduce some prior information.
Using Bayesian statistics [25, 69] one argues that instead of minimizing
$|\chi(k)|^{2}$ the spectral function should maximize

$\displaystyle Q(k)=\alpha S(k)-\frac{1}{2}|\chi(k)|^{2}$ (79)

where

$\displaystyle\begin{split}S[\mathcal{A}_{0}](k)=&\int
d\omega\left(\mathcal{A}(k,\omega)-\mathcal{A}_{0}(k,\omega)\right.\\\
&\left.-\mathcal{A}(k,\omega)\ln{(\mathcal{A}(k,\omega)/\mathcal{A}_{0}(k,\omega))}\right)\end{split}$
(80)

is the entropy relative to a smooth default model $\mathcal{A}_{0}$ containing
analytically known properties. Since $|\chi(k)|^{2}$ provides little
information in the high and low-frequency tails, we use the default model

$\displaystyle\mathcal{A}_{0}(k,\omega-k^{2}/(2m))=\frac{\gamma\left(\sqrt{\frac{\delta}{\gamma}+\omega^{2}}+\omega\right)}{\left(\rho+\omega^{2}\right)^{7/4}}\,,$
(81)

with $\gamma=2\pi nm^{-3/2}$ and $\delta=8\pi nm^{-7/2}\mathcal{C}$, which
satisfies the known asymptotic behavior from Eq. (46) for $\omega\to\infty$
and from Eq. (47) for $\omega\to-\infty$. Furthermore, $\rho$ is determined by
the normalization condition $\int_{\omega}\mathcal{A}_{0}(k,\omega)=1$ and on
the typical energy scales set by $T$ and $\epsilon_{F}$ the default model is a
smooth, slowly varying function.

Figure 11: Spectral function of the pair susceptibility for the unitary
balanced system at the critical temperature $T/T_{F}=0.156$.

Due to the large information gap left by $|\chi(k)|^{2}$ the choice of the
Lagrange parameter $\alpha$ is very important. We follow Bryan’s method [68]
and obtain the optimal spectral function $\mathcal{A}_{\alpha}(k)$ for each
value of $\alpha$. These are then averaged over $\alpha$ to obtain the best
estimate for the spectral function

$\displaystyle\mathcal{A}_{\text{NAC}}(k)=\mathcal{N}\\!\int\\!d\alpha\mathcal{A}_{\alpha}(k)f(\alpha)e^{Q(\mathcal{A}_{\alpha}(k))}\prod_{i}\left(\frac{\alpha}{\alpha+\lambda_{i}}\right)^{\\!1/2}\,.$
(82)

Here $\mathcal{N}$ denotes the normalization and $\lambda_{i}$ are the
eigenvalues of the curvature matrix of $|\chi(k)|^{2}$ in the space of
spectral functions. Finally, $f(\alpha)$ is an undetermined, dimensionless
scaling function, for which we choose $f(\alpha)=\alpha^{\delta}$. This choice
determines the relative importance of entropy and constraints from imaginary
times. For too large values of $\delta$ the entropy is weighted too heavily
such that eventually
$\mathcal{A}_{\text{NAC}}(\omega,k)\underset{\delta\to\infty}{\longrightarrow}\mathcal{A}_{0}(\omega,k)$.
For too small values of $\delta$ on the other hand,
$\mathcal{A}_{\text{NAC}}(\omega,k)$ develops oscillatory instabilities.
Hence, the choice of the scaling function $f(\alpha)$ has a significant impact
on the spectral function $\mathcal{A}_{\text{NAC}}$. We fix $\delta=-1.74$,
which is the smallest value, that induces no oscillations in the asymptotic
behavior at large positive or negative frequencies. Since this condition is
somewhat arbitrary, we emphasize that the qualitative difference with the
spectral function obtained directly in real frequencies
$\mathcal{A}_{\text{Keldysh}}(\omega,k)$, which we report in Fig. 8, remains
for all stable values of $\delta$.

##  Appendix D Pair susceptibility

A recent theoretical proposal has suggested that it might be possible to probe
the pair propagator using Raman spectroscopy techniques [98]. Here, the pair
propagator for the imbalanced system was argued to contain information about
the FFLO phase transition giving it experimental relevance. Another example
where the spectrum of the pair propagator can be of use is in questions
regarding transport [18]. Here computation of dynamical transport
coefficients, like the dynamical shear viscosity, is simplified if one has
access to frequency-resolved $\Gamma^{R}(k,\omega)$. For bosonic propagators,
the NAC is even more challenging [99], making reliable results difficult to
acquire. With our method the full frequency dependence of $\Gamma^{R}$ is
computed at no additional cost and bosonic spectra are in general computed
with the same high accuracy as their fermionic counterparts. To highlight this
we plot, in Fig. 11, the spectral function of the pair-propagator for the
balanced case at unitarity and at the critical temperature. Here the
singularity at $\omega=0$ for small momentum and the slowly decaying frequency
tails observed in imaginary time in Fig. 6, also appear in real frequencies.

## References

  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Zwerger [2012] W. Zwerger, ed.,  _The BCS-BEC Crossover and the Unitary Fermi Gas_ (Lecture Notes in Physics, Vol. 836, Springer, Berlin, Heidelberg, 2012).
  * Zwerger [2016] W. Zwerger, Strongly interacting Fermi gases, in  _Quantum matter at ultralow temperatures, Proceedings of the International School of Physics ’Enrico Fermi’, Course 191, Varenna, 7-15 July 2014_ , edited by M. Inguscio, W. Ketterle, S. Stringari, and G. Roati (IOS Press, Amsterdam, 2016) pp. 63–141.
  * Hofmann [2012] J. Hofmann, Quantum anomaly, universal relations, and breathing mode of a two-dimensional fermi gas, Phys. Rev. Lett. 108, 185303 (2012).
  * Gao and Yu [2012] C. Gao and Z. Yu, Breathing mode of two-dimensional atomic fermi gases in harmonic traps, Phys. Rev. A 86, 043609 (2012).
  * Son and Stephanov [2006] D. T. Son and M. A. Stephanov, Phase diagram of a cold polarized fermi gas, Phys. Rev. A 74, 013614 (2006).
  * Pilati and Giorgini [2008] S. Pilati and S. Giorgini, Phase separation in a polarized fermi gas at zero temperature, Phys. Rev. Lett. 100, 030401 (2008).
  * Shin  _et al._ [2008] Y.-i. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, Phase diagram of a two-component fermi gas with resonant interactions, Nature 451, 689 (2008).
  * Ku  _et al._ [2012] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, Revealing the superfluid lambda transition in the universal thermodynamics of a unitary fermi gas, Science 335, 563 (2012), https://www.science.org/doi/pdf/10.1126/science.1214987 .
  * Tan [2008a] S. Tan, Large momentum part of a strongly correlated fermi gas, Annals of Physics 323, 2971 (2008a).
  * Tan [2008b] S. Tan, Energetics of a strongly correlated fermi gas, Annals of Physics 323, 2952 (2008b).
  * Tan [2008c] S. Tan, Generalized virial theorem and pressure relation for a strongly correlated fermi gas, Annals of Physics 323, 2987 (2008c).
  * Mukherjee  _et al._ [2019] B. Mukherjee, P. B. Patel, Z. Yan, R. J. Fletcher, J. Struck, and M. W. Zwierlein, Spectral response and contact of the unitary fermi gas, Phys. Rev. Lett. 122, 203402 (2019).
  * Patel  _et al._ [2020] P. B. Patel, Z. Yan, B. Mukherjee, R. J. Fletcher, J. Struck, and M. W. Zwierlein, Universal sound diffusion in a strongly interacting fermi gas, Science 370, 1222 (2020), https://www.science.org/doi/pdf/10.1126/science.aaz5756 .
  * Haussmann [1994] R. Haussmann, Properties of a fermi liquid at the superfluid transition in the crossover region between bcs superconductivity and bose-einstein condensation, Phys. Rev. B 49, 12975 (1994).
  * Haussmann  _et al._ [2007] R. Haussmann, W. Rantner, S. Cerrito, and W. Zwerger, Thermodynamics of the bcs-bec crossover, Phys. Rev. A 75, 023610 (2007).
  * Enss  _et al._ [2011] T. Enss, R. Haussmann, and W. Zwerger, Viscosity and scale invariance in the unitary fermi gas, Annals of Physics 326, 770 (2011).
  * Frank  _et al._ [2020] B. Frank, W. Zwerger, and T. Enss, Quantum critical thermal transport in the unitary fermi gas, Phys. Rev. Res. 2, 023301 (2020).
  * Stewart  _et al._ [2008] J. T. Stewart, J. P. Gaebler, and D. S. Jin, Using photoemission spectroscopy to probe a strongly interacting fermi gas, Nature 454, 744 (2008).
  * Schirotzek  _et al._ [2008] A. Schirotzek, Y.-i. Shin, C. H. Schunck, and W. Ketterle, Determination of the superfluid gap in atomic fermi gases by quasiparticle spectroscopy, Phys. Rev. Lett. 101, 140403 (2008).
  * Yan  _et al._ [2019] Z. Yan, P. B. Patel, B. Mukherjee, R. J. Fletcher, J. Struck, and M. W. Zwierlein, Boiling a unitary fermi liquid, Phys. Rev. Lett. 122, 093401 (2019).
  * Feld  _et al._ [2011] M. Feld, B. Fröhlich, E. Vogt, M. Koschorreck, and M. Köhl, Observation of a pairing pseudogap in a two-dimensional fermi gas, Nature 480, 75 (2011).
  * Gaebler  _et al._ [2010] J. P. Gaebler, J. T. Stewart, T. E. Drake, D. S. Jin, A. Perali, P. Pieri, and G. C. Strinati, Observation of pseudogap behaviour in a strongly interacting fermi gas, Nature Physics 6, 569 (2010).
  * Sommer  _et al._ [2011] A. Sommer, M. Ku, G. Roati, and M. W. Zwierlein, Universal spin transport in a strongly interacting fermi gas, Nature 472, 201 (2011).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data, Physics Reports 269, 133 (1996).
  * Goulko  _et al._ [2017] O. Goulko, A. S. Mishchenko, L. Pollet, N. Prokof’ev, and B. Svistunov, Numerical analytic continuation: Answers to well-posed questions, Phys. Rev. B 95, 014102 (2017).
  * Haussmann  _et al._ [2009] R. Haussmann, M. Punk, and W. Zwerger, Spectral functions and rf response of ultracold fermionic atoms, Phys. Rev. A 80, 063612 (2009).
  * Magierski  _et al._ [2009] P. Magierski, G. Wlazłowski, A. Bulgac, and J. E. Drut, Finite-temperature pairing gap of a unitary fermi gas by quantum monte carlo calculations, Phys. Rev. Lett. 103, 210403 (2009).
  * Schmidt and Enss [2011] R. Schmidt and T. Enss, Excitation spectra and rf response near the polaron-to-molecule transition from the functional renormalization group, Phys. Rev. A 83, 063620 (2011).
  * Pini  _et al._ [2019] M. Pini, P. Pieri, and G. C. Strinati, Fermi gas throughout the bcs-bec crossover: Comparative study of $t$-matrix approaches with various degrees of self-consistency, Phys. Rev. B 99, 094502 (2019).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Gottfried and Yan [2003] K. Gottfried and T. Yan,  _Quantum Mechanics: Fundamentals_, Graduate Texts in Contemporary Physics (Springer New York, 2003).
  * Gribakin and Flambaum [1993] G. F. Gribakin and V. V. Flambaum, Calculation of the scattering length in atomic collisions using the semiclassical approximation, Phys. Rev. A 48, 546 (1993).
  * Flambaum  _et al._ [1999] V. V. Flambaum, G. F. Gribakin, and C. Harabati, Analytical calculation of cold-atom scattering, Phys. Rev. A 59, 1998 (1999).
  * Feshbach [1958] H. Feshbach, A Unified Theory of Nuclear Reactions, Annals of Physics 5, 357 (1958).
  * Fano [1961] U. Fano, Effects of Configuration Interaction on Intensities and Phase Shifts, Phys. Rev. 124, 1866 (1961).
  * Timmermans  _et al._ [2001] E. Timmermans, K. Furuya, P. W. Milonni, and A. K. Kerman, Prospect of creating a composite fermi–bose superfluid, Physics Letters A 285, 228 (2001).
  * Holland  _et al._ [2001] M. Holland, S. J. J. M. F. Kokkelmans, M. L. Chiofalo, and R. Walser, Resonance superfluidity in a quantum degenerate fermi gas, Phys. Rev. Lett. 87, 120406 (2001).
  * Schmidt, R.  _et al._ [2012] Schmidt, R., Rath, S.P., and Zwerger, W., Efimov physics beyond universality, Eur. Phys. J. B 85, 386 (2012).
  * Rath and Schmidt [2013] S. P. Rath and R. Schmidt, Field-theoretical study of the bose polaron, Phys. Rev. A 88, 053632 (2013).
  * Navon  _et al._ [2010] N. Navon, S. Nascimbène, F. Chevy, and C. Salomon, The Ground State of a Fermi Gas with Tunable Interactions, Science 328, 729 (2010).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Berges [2004] J. Berges, Introduction to Nonequilibrium Quantum Field Theory, AIP Conference Proceedings 739, 3 (2004), https://pubs.aip.org/aip/acp/article-pdf/739/1/3/11882676/3_1_online.pdf .
  * Rammer and Smith [1986] J. Rammer and H. Smith, Quantum field-theoretical methods in transport theory of metals, Rev. Mod. Phys. 58, 323 (1986).
  * Morawetz [2018] K. Morawetz,  _Interacting Systems far from Equilibrium_ (Oxford University Press, 2018).
  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_, 2nd ed. (Cambridge University Press, 2010).
  * Haug and Jauho [2007] H. Haug and A.-P. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer, 2007).
  * Frank  _et al._ [2018] B. Frank, J. Lang, and W. Zwerger, Universal phase diagram and scaling functions of imbalanced fermi gases, Journal of Experimental and Theoretical Physics 127, 812 (2018).
  * Baym [1962] G. Baym, Self-consistent approximations in many-body systems, Phys. Rev. 127, 1391 (1962).
  * Cornwall  _et al._ [1974] J. M. Cornwall, R. Jackiw, and E. Tomboulis, Effective action for composite operators, Phys. Rev. D 10, 2428 (1974).
  * Luttinger and Ward [1960] J. M. Luttinger and J. C. Ward, Ground-state energy of a many-fermion system. ii, Phys. Rev. 118, 1417 (1960).
  * Van Houcke  _et al._ [2012] K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek, and M. W. Zwierlein, Feynman diagrams versus fermi-gas feynman emulator, Nature Physics 8, 366 (2012).
  * Rossi  _et al._ [2018] R. Rossi, T. Ohgoe, K. Van Houcke, and F. Werner, Resummation of diagrammatic series with zero convergence radius for strongly correlated fermions, Phys. Rev. Lett. 121, 130405 (2018).
  * Gor’kov and Melik-Barkhudarov [1961] L. Gor’kov and T. Melik-Barkhudarov, Contribution to the Theory of Superfluidity in an Imperfect Fermi Gas, Zh. Eskp. Theor. Fiz. 40, 1452 (1961).
  * Petrov  _et al._ [2004] D. S. Petrov, C. Salomon, and G. V. Shlyapnikov, Weakly bound dimers of fermionic atoms, Phys. Rev. Lett. 93, 090404 (2004).
  * Nozières and Schmitt-Rink [1985] P. Nozières and S. Schmitt-Rink, Bose condensation in an attractive fermion gas: From weak to strong coupling superconductivity, Journal of Low Temperature Physics 59, 195 (1985).
  * Veillette  _et al._ [2007] M. Y. Veillette, D. E. Sheehy, and L. Radzihovsky, Large-$n$ expansion for unitary superfluid fermi gases, Phys. Rev. A 75, 043614 (2007).
  * Chevy [2006] F. Chevy, Universal phase diagram of a strongly interacting Fermi gas with unbalanced spin populations, Phys. Rev. A 74, 063628 (2006).
  * Sa de Melo  _et al._ [1993] C. A. R. Sa de Melo, M. Randeria, and J. R. Engelbrecht, Crossover from BCS to Bose Superconductivity: Transition Temperature and Time-Dependent Ginzburg-Landau Theory, Phys. Rev. Lett. 71, 3202 (1993).
  * Diener  _et al._ [2008] R. B. Diener, R. Sensarma, and M. Randeria, Quantum fluctuations in the superfluid state of the bcs-bec crossover, Phys. Rev. A 77, 023626 (2008).
  * Braaten  _et al._ [2010] E. Braaten, D. Kang, and L. Platter, Short-time operator product expansion for rf spectroscopy of a strongly interacting fermi gas, Phys. Rev. Lett. 104, 223004 (2010).
  * Cooley and Tukey [1965] J. Cooley and J. Tukey, An Algorithm for the Machine Calculation of Complex Fourier Series, Mathematics of Computation 19, 297 (1965).
  * Haines and Jones [1988] G. V. Haines and A. G. Jones, Logarithmic fourier transformation, Geophysical Journal 92, 171 (1988).
  * Lang and Frank [2019] J. Lang and B. Frank, Fast logarithmic fourier-laplace transform of nonintegrable functions, Phys. Rev. E 100, 053302 (2019).
  * Atkinson [1989] E. K. Atkinson,  _An introduction to numerical analysis_ , 2nd ed. (Wiley, 1989).
  * Ömer Eǧecioǧlu  _et al._ [1989] Ömer Eǧecioǧlu, E. Gallopoulos, and Çetin K Koç, Fast computation of divided differences and parallel hermite interpolation, Journal of Complexity 5, 417 (1989).
  * Frank [2019] B. Frank,  _Thermodynamics and Transport in Fermi Gases near Unitarity_ , Ph.D. thesis, Technische Universität München (2019).
  * Bryan [1990] R. K. Bryan, Maximum entropy analysis of oversampled data problems, European Biophysics Journal 18, 165 (1990).
  * Sivia and Skilling [2006] D. Sivia and J. Skilling,  _Data Analysis: A Bayesian Tutorial_, Oxford science publications (OUP Oxford, 2006).
  * Castin [2008] Y. Castin, Basic theory tools for degenerate fermi gases, in  _Ultra-cold Fermi Gases, Proceedings of the International School of Physics ’Enrico Fermi’, Course 164, Varenna, 20-30 June 2006_ , edited by M. Inguscio, W. Ketterle, and C. Salomon (IOS Press, Amsterdam, 2008) pp. 289–349.
  * Norman  _et al._ [1998] M. R. Norman, M. Randeria, H. Ding, and J. C. Campuzano, Phenomenology of the low-energy spectral function in high-${T}_{c}$ superconductors, Phys. Rev. B 57, R11093 (1998).
  * Dagotto [1994] E. Dagotto, Correlated electrons in high-temperature superconductors, Rev. Mod. Phys. 66, 763 (1994).
  * Randeria [2010] M. Randeria, Pre-pairing for condensation, Nature Physics 6, 561 (2010).
  * Chien  _et al._ [2010] C.-C. Chien, H. Guo, Y. He, and K. Levin, Comparative study of bcs-bec crossover theories above ${T}_{c}$: The nature of the pseudogap in ultracold atomic fermi gases, Phys. Rev. A 81, 023622 (2010).
  * Löhneysen  _et al._ [2007] H. v. Löhneysen, A. Rosch, M. Vojta, and P. Wölfle, Fermi-liquid instabilities at magnetic quantum phase transitions, Rev. Mod. Phys. 79, 1015 (2007).
  * Lee [2018] S.-S. Lee, Recent developments in non-fermi liquid theory, Annual Review of Condensed Matter Physics 9, 227 (2018), https://doi.org/10.1146/annurev-conmatphys-031016-025531 .
  * Klein  _et al._ [2020] A. Klein, A. V. Chubukov, Y. Schattner, and E. Berg, Normal state properties of quantum critical metals at finite temperature, Phys. Rev. X 10, 031053 (2020).
  * Frank  _et al._ [2023] B. Frank, Z. H. Liu, F. F. Assaad, M. Vojta, and L. Janssen, Marginal fermi liquid at magnetic quantum criticality from dimensional confinement, Phys. Rev. B 108, L100405 (2023).
  * Frank and Piazza [2020] B. Frank and F. Piazza, $\omega/t$ scaling and ir/uv-mixing in ising-nematic quantum critical metals (2020), arXiv:2011.07076 [cond-mat.str-el] .
  * Prüfer  _et al._ [2018] M. Prüfer, P. Kunkel, H. Strobel, S. Lannig, D. Linnemann, C.-M. Schmied, J. Berges, T. Gasenzer, and M. K. Oberthaler, Observation of universal dynamics in a spinor bose gas far from equilibrium, Nature 563, 217 (2018).
  * Eigen  _et al._ [2018] C. Eigen, J. A. P. Glidden, R. Lopes, E. A. Cornell, R. P. Smith, and Z. Hadzibabic, Universal prethermal dynamics of bose gases quenched to unitarity, Nature 563, 221 (2018).
  * Erne  _et al._ [2018] S. Erne, R. Bücker, T. Gasenzer, J. Berges, and J. Schmiedmayer, Universal dynamics in an isolated one-dimensional bose gas far from equilibrium, Nature 563, 225 (2018).
  * Scazza  _et al._ [2017] F. Scazza, G. Valtolina, P. Massignan, A. Recati, A. Amico, A. Burchianti, C. Fort, M. Inguscio, M. Zaccanti, and G. Roati, Repulsive Fermi polarons in a resonant mixture of ultracold Li-6 atoms, Phys Rev. Lett. 118, 083602 (2017).
  * Knoll  _et al._ [2001] J. Knoll, Y. Ivanov, and D. Voskresensky, Exact conservation laws of the gradient expanded kadanoff–baym equations, Annals of Physics 293, 126 (2001).
  * Basov  _et al._ [2017] D. N. Basov, R. D. Averitt, and D. Hsieh, Towards properties on demand in quantum materials, Nature Materials 16, 1077 (2017).
  * Lang  _et al._ [2023] J. Lang, M. Buchhold, and S. Diehl, Field theory for the dynamics of the open $o(n)$ model (2023), arXiv:2310.06892 [cond-mat.stat-mech] .
  * Stojchevska  _et al._ [2014] L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Ultrafast switching to a stable hidden quantum state in an electronic crystal, Science 344, 177 (2014), https://www.science.org/doi/pdf/10.1126/science.1241591 .
  * Liu  _et al._ [2011] M. K. Liu, B. Pardo, J. Zhang, M. M. Qazilbash, S. J. Yun, Z. Fei, J.-H. Shin, H.-T. Kim, D. N. Basov, and R. D. Averitt, Photoinduced phase transitions by time-resolved far-infrared spectroscopy in ${\mathrm{v}}_{2}{\mathrm{o}}_{3}$, Phys. Rev. Lett. 107, 066403 (2011).
  * Zhang  _et al._ [2016] J. Zhang, X. Tan, M. Liu, S. W. Teitelbaum, K. W. Post, F. Jin, K. A. Nelson, D. N. Basov, W. Wu, and R. D. Averitt, Cooperative photoinduced metastable phase control in strained manganite films, Nature Materials 15, 956 (2016).
  * Vlietinck  _et al._ [2013] J. Vlietinck, J. Ryckebusch, and K. Van Houcke, Quasiparticle properties of an impurity in a fermi gas, Phys. Rev. B 87, 115133 (2013).
  * Van Houcke  _et al._ [2020] K. Van Houcke, F. Werner, and R. Rossi, High-precision numerical solution of the fermi polaron problem and large-order behavior of its diagrammatic series, Phys. Rev. B 101, 045134 (2020).
  * von Milczewski  _et al._ [2022] J. von Milczewski, F. Rose, and R. Schmidt, Functional-renormalization-group approach to strongly coupled bose-fermi mixtures in two dimensions, Phys. Rev. A 105, 013317 (2022).
  * Eugen [2023] D. Eugen, Spectral Functions of Ultracold Fermi Gases from Dyson-Schwinger Equations, MSc. thesis, Heidelberg University (2023).
  * Taylor [1972] J. R. Taylor,  _Scattering Theory_ (Wiley, 1972).
  * Braaten and Hammer [2006] E. Braaten and H.-W. Hammer, Universality in few-body systems with large scattering length, Physics Reports 428, 259 (2006).
  * Bruus and Flensberg [2001] H. Bruus and K. Flensberg,  _Many-Body quantum theory in condensed matter physics_ (Oxford University Press, 2001).
  * Phillips [2003] G. M. Phillips,  _Interpolation and Approximation by Polynomials_ , 1st ed. (Wiley, 2003).
  * Diessel  _et al._ [2022] O. K. Diessel, J. von Milczewski, A. Christianen, and R. Schmidt, Probing molecular spectral functions and unconventional pairing using raman spectroscopy (2022), arXiv:2209.11758 [cond-mat.quant-gas] .
  * Nogaki and Shinaoka [2023] K. Nogaki and H. Shinaoka, Bosonic nevanlinna analytic continuation, Journal of the Physical Society of Japan 92, 035001 (2023), https://doi.org/10.7566/JPSJ.92.035001 .