# Spectral properties and observables in ultracold Fermi gases

Eugen Dizer  Institut für Theoretische Physik, Universität Heidelberg,
Philosophenweg 16, 69120 Heidelberg, Germany  Jan Horak  Institut für
Theoretische Physik, Universität Heidelberg, Philosophenweg 16, 69120
Heidelberg, Germany  Jan M. Pawlowski  Institut für Theoretische Physik,
Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany  ExtreMe
Matter Institute EMMI, GSI, Planckstr. 1, 64291 Darmstadt, Germany

###### Abstract

We calculate non-perturbative self-consistent fermionic and bosonic spectral
functions of ultracold Fermi gases with the spectral functional approach. This
approach allows for a direct real-time computation of non-perturbative
correlation functions, and in the present work we use spectral Dyson-Schwinger
equations. We focus on the normal phase of the spin-balanced Fermi gas and
provide numerical results for the full fermionic and bosonic spectral
functions. The spectral functions are then used for the determination of the
equation of state, the Tan contact and ejection rf spectra at unitarity. These
results are compared to experimental data, the self-consistent T-matrix
approach and lattice results. Our approach offers a wide range of
applications, including the ab initio calculation of transport and spectral
properties of the superfluid phase in the BCS-BEC crossover.

##  I Introduction

Ultracold Fermi gases [1, 2, 3] represent a versatile tool for studying
strongly correlated quantum systems at low energies. The experimental control
over the system parameters allows for the investigation of a wide range of
phenomena, such as the BCS-BEC crossover [4, 5, 6, 7, 8], the formation of
polarons [9, 10, 11, 12], and the emergence of superfluidity [13, 14, 15].
Various important properties of the gas, such as transport and scattering
properties [16, 17, 18] or the excitation spectrum [19, 20], are encoded in
its spectral functions.

Strong correlations, e.g. close to unitarity [1], render the theoretical
description of interacting Fermi gases a non-perturbative problem. Therefore,
various non-perturbative approaches, such as Quantum Monte Carlo (QMC)
simulations [21, 22], self-consistent T-matrix theory [23, 24, 25, 26, 27] or
Dyson-Schwinger equation (DSE) [28, 29, 30] and functional renormalisation
group (fRG) [31, 32, 33] approaches, have been employed to describe the
system.

The computation of spectral functions requires access to the fermionic and
bosonic self-energies at real frequencies. While QMC simulations are
formulated at imaginary frequencies by construction, also functional methods
such as DSE and fRG are usually formulated in Euclidean spacetime due to
significantly reduced computational costs. However, this requires an analytic
continuation of the correlation functions from Matsubara to real frequencies,
which is an ill-conditioned numerical task accompanied by large systematic
uncertainties [34].

Recently, the spectral functional approach has been put forward [35]. Notably,
this approach allows for the direct computation of non-perturbative
correlation functions at real frequencies: with the spectral representations
of propagators and vertices in non-perturbative diagrams, the respective
frequency sums or integrals can be performed analytically, and the remaining
real spectral integrals are readily computed numerically. The spectral
functional approach has been developed within the framework of Dyson-Schwinger
equations, but it applies to all functional approaches. The spectral
functional renormalisation has been set up and used in [36, 37], and the
spectral approach is tailor made for application in 2PI resummation schemes
[38], where all correlation functions are constructed from non-perturbative
propagators. It has been successfully used for the non-perturbative
computation of spectral functions and bound state properties in scalar
theories [35, 39, 40], gauge theories [41, 42, 43, 44], and quantum gravity
[36].

In this work, we use the spectral functional approach for the computation of
non-perturbative self-consistent single-particle spectral functions in the
normal phase of a strongly interacting spin-balanced Fermi gas for all
scattering lengths $a$, and in particular close to a Feshbach resonance [45]
with $a\to\infty$. This constitutes the first application of the spectral
functional approach to non-relativistic systems. Specifically, we set up
spectral Dyson-Schwinger equations for the fermionic and bosonic self-energies
at real frequencies, yielding direct access to their spectral functions. With
the spectral functions we also determine the density equation of state [15],
the Tan contact [46] and ejection radio-frequency (rf) spectra [7] at
unitarity. Our results agree with results from the self-consistent T-matrix
approximation [23, 24], but the flexibility of the spectral functional
approaches gives access to more generic expansion schemes.

The present paper is organised as follows. In Section II, we introduce our
microscopic model for the spin-balanced Fermi gas. In Section III, we set up
the spectral Dyson-Schwinger equations for bosonic and fermionic propagators.
Our numerical results for the spectral functions as well as the density, rf
spectra and Tan contact are presented in Section IV. We conclude in Section V.

##  II Microscopic Model

We consider a non-relativistic two-component Fermi gas described by the
Euclidean classical action

$\displaystyle S[\psi]=\int_{0}^{\beta}d\tau\int
d^{3}x\,\Big{[}\sum_{\sigma=\uparrow,\downarrow}\psi_{\sigma}^{*}($
$\displaystyle\partial_{\tau}-\nabla^{2}-\mu_{\sigma})\psi_{\sigma}$
$\displaystyle+\lambda\,\psi_{\uparrow}^{*}\psi_{\downarrow}^{*}\psi_{\downarrow}\psi_{\uparrow}\Big{]}\,,$
(1)

where $\lambda$ is the coupling constant of the contact interaction. If used
on the mean field level, $\lambda$ would be related directly to the scattering
length $a$ via [33]

$\displaystyle\lambda=8\pi a\,.$ (2)

Instead of using the common split into bare coupling and cutoff contribution,
2 is readily converted into a cutoff-independent relation including quantum
effects, for more details see Section III and Appendix A, and in particular 33
and 34.

In 1, we also introduced the chemical potential of the fermion species,
$\mu_{\sigma}$ with $\sigma=(\uparrow,\downarrow)$. The fermionic fields
$\psi_{\sigma}(\tau,\bm{x})$ are Grassmann valued and depend on the Euclidean
time $\tau$, which is restricted to the circumference $\beta=1/T$, and the
spatial coordinates $\bm{x}$. In 1, we have also used natural units,
$\hbar=k_{B}=2\,m_{\psi}=1$, where $m_{\psi}$ is the mass of fermions.

The emergence of bosonic dimers and their condensation at low temperatures in
the fermionic ultracold gas can be efficiently accommodated by rewriting the
action 1 with a Hubbard–Stratonovich transformation. Thus, we will consider a
spin-balanced system ($\mu_{\uparrow}=\mu_{\downarrow}=\mu$), where the
contact interaction of fermions with opposite spin is substituted by the
exchange of bosonic dimers,

$\displaystyle S[\psi,\phi]=\int_{0}^{\beta}d\tau$ $\displaystyle\int
d^{3}x\,\Big{[}\sum_{\sigma=\uparrow,\downarrow}\psi^{*}_{\sigma}(\partial_{\tau}-\nabla^{2}-\mu)\psi_{\sigma}$
$\displaystyle+\nu\phi^{*}\phi-h\left(\phi^{*}\psi_{\uparrow}\psi_{\downarrow}-\phi\psi^{*}_{\uparrow}\psi^{*}_{\downarrow}\right)\,\Big{]}\,,$
(3)

where $h$ is the Feshbach coupling between the fermions and bosons, and $\nu$
is the detuning of the dimer. On the mean field level, similarly to 2, the
couplings and the detuning parameter would be related by

$\displaystyle\lambda=-\frac{h^{2}}{\nu}\,.$ (4)

For the full relation on the quantum level, see Section III and Appendix A,
and in particular 33 and 34.

It is readily shown that 3 reduces to 1 on the equation of motion (EoM) of
$\phi$, also using 4. Solving 4 for the detuning parameter on the full quantum
level, leads to 34 in Appendix A. Notably, while the four-fermion coupling
$\lambda$ gets strongly affected by fluctuations, the Feshbach coupling $h$ is
well-approximated by its classical value, as these fluctuation effects are
encoded in $\nu$. The latter is used to tune the (inverse) scattering length,
for a detailed discussion of all these properties see [33]. We also emphasise
that the absence of an explicit $s$-channel four-fermion scattering can be
sustained on the quantum level within the fRG approach with emergent
composites or rebosonisation, for the respective conceptual developments see
[47, 48, 49, 50], for an application to the present system see [51, 52].
Heuristically, this approach may be understood as a scale-dependent Hubbard-
Stratonovich transformation.

For a dilute ultracold Fermi gas, the effective range $r_{0}$ of the
interaction is of the order of the inverse van der Waals length
$l_{\mathrm{vdw}}$ and defines a physical momentum cutoff scale $\Lambda\sim
1/r_{0}$, which is much larger then all relevant physical scales ($\Lambda\gg
k_{F}$). Throughout this work, we will work in the limit of contact
interactions with zero effective range, such that $\Lambda$ can be sent to
infinity in the end. For a more detailed discussion of the renormalisation
scheme, we refer to Appendix A.

The complex bosonic field $\phi(\tau,\bm{x})$ has been introduced as an
efficient book-keeping device for the $s$-channel interaction between fermions
of opposite spin. For example, bubble resummations of $s$-channel diagrams,
but also beyond bubble resummations of $t$ and $u$-channel contributions are
encoded in two-point function diagrams or corrections to the Yukawa coupling.
More importantly, emergent dynamical $s$-channel dimer degrees of freedom are
already taken care of as is their condensation.

Hence, by tuning the dimensionless interaction strength $(k_{F}a)^{-1}$, the
system can be driven from a BCS-type superfluid to a Bose-Einstein condensed
state, allowing for the exploration of the BCS-BEC crossover [53, 54, 55, 56].
Here, the interaction strength is measured in terms of the Fermi momentum

$\displaystyle k_{F}=(3\pi^{2}n)^{1/3}\,,\qquad
T_{F}=\varepsilon_{F}=k_{F}^{2}\,,$ (5)

where $n$ is the total density, and we have also defined the Fermi
energy/temperature for later use. On the BCS side of the crossover,
$(k_{F}a)^{-1}<0$, the bosonic dimer $\phi$ describes weakly bound Cooper
pairs. On the BEC side, $(k_{F}a)^{-1}>0$, $\phi$ describes tightly bound
bosonic molecules. Bose-Einstein condensation of the bosonic pairs, i.e., a
non-vanishing expectation value of $\phi$, leads to superfluidity in the
system. In the following, we will consider the normal phase of the BCS-BEC
crossover with a special focus on the strongly correlated, unitary regime at
$(k_{F}a)^{-1}=0$.

##  III Spectral Dyson-Schwinger Equations

In this Section we briefly introduce the coupled spectral propagator DSEs
which will be used in obtaining all numerical and analytical results in this
work. We also specify the underlying approximation of the effective action.

###  III.1 Gap equations and the approximation of the effective action

Functional approaches, such as the DSEs, the fRG or two-particle irreducible
(2PI) resummation schemes, provide exact non-perturbative relations for
correlation functions in terms of full propagators and vertices (fRG) as well
as the classical vertices (DSE, 2PI). For reviews, see [57] (DSEs) and [58,
59, 60, 61, 62] (fRG) and [38] (2PI). The central object in these approaches
is the quantum effective action $\Gamma[\Phi]$, which is the generating
functional of one-particle irreducible (1PI) correlation functions, or the 2PI
effective action. Correlation functions of $n$ fields are obtained by $n$
functional derivatives with respect to the the mean fields,
$\Phi=(\psi_{\sigma},\psi_{\sigma}^{*},\phi,\phi^{*})$, where for the sake of
notational simplicity we have used the same expressions as for the classical
fields.

The most important example is the full non-perturbative two-point function,
the full propagator $G_{\Phi}$, which is the pivotal building block of all
diagrammatic functional approaches. In terms of the 1PI effective action, it
is given is the inverse of the 1PI two-point function,
$G_{\Phi}=\left[\Gamma^{(2)}\right]^{-1}_{\Phi\Phi^{*}}$, which is a matrix
inverse in field space. In the normal phase, the expectation value of the
bosonic field is zero, i.e. $\phi=0$, and the propagator, if evaluated on the
EoMs, is diagonal in field space. Since we are dealing with a balanced system,
the propagators of the $\uparrow$ and $\downarrow$ species are the same and we
are left with only one fermion propagator
$G_{\psi_{\uparrow}}=G_{\psi_{\downarrow}}=G_{\psi}$ and one boson propagator
$G_{\phi}$. The gap equation or DSE for the inverse propagator
$\Gamma^{(2)}_{\Phi\Phi^{*}}$ involves the classical inverse propagator
$S^{(2)}_{\Phi\Phi^{*}}$, whose fermionic and bosonic components are obtained
from the microscopic action 3,

$\displaystyle S^{(2)}_{\psi\psi^{*}}(P)=-i\omega_{n}+\bm{p}^{2}-\mu\,,\qquad
S^{(2)}_{\phi\phi^{*}}(Q)=\nu\,.$ (6)

Here, $P=(\omega_{n},\bm{p})$ and $Q=(\epsilon_{n},\bm{q})$ with the fermionic
($\omega_{n}$) and bosonic ($\epsilon_{n}$) Matsubara frequencies
respectively,

$\displaystyle\omega_{n}=(2n+1)\pi T\,,\qquad\epsilon_{n}=2n\pi T\,.$ (7)

The coupled Dyson-Schwinger equations for the fermion and boson propagator are
depicted in Figure 1, while the diagrammatic notation is summarised in Figure
2. The renormalised fermion and boson gap equations are given by

$\displaystyle G^{-1}_{\psi}(P)$
$\displaystyle=-i\omega_{n}+\bm{p}^{2}-\mu+\Sigma_{\psi}(P)\,,$ $\displaystyle
G^{-1}_{\phi}(Q)$ $\displaystyle=-\frac{h^{2}}{8\pi a}-\Pi_{\phi}(Q)\,,$ (8)

with the renormalised fermionic and bosonic self-energy, $\Sigma_{\psi}(P)$
and $\Pi_{\phi}(Q)$, respectively. For more details on the spectral
renormalisation procedure see Appendix A.

Figure 1: Fermion (solid line) and boson (dotted line) propagator Dyson-
Schwinger equation for the balanced ultracold Fermi gas. Notation as defined
in Figure 2.

Now we specify the approximation of the effective action used for the
following computations. As already discussed below 4, the approximation of a
classical constant Yukawa coupling $h$ is working well as the fluctuations of
the four-fermi interaction are well-captured by that in the detuning
parameter. The latter is part of the full propagator and hence will be
accessed here. This approximation is readily lifted by also solving the
spectral DSE for the Yukawa coupling and the remnant four-fermion interaction
which will be considered elsewhere. This leads us to the following
approximation of the renormalised finite effective action,

$\displaystyle\Gamma[\psi,\phi]=$ $\displaystyle\,\int_{0}^{\beta}d\tau\int
d^{3}x\,\Biggl{[}\psi^{*}_{\sigma}\Bigl{(}\partial_{\tau}-\nabla^{2}-\mu+\Sigma_{\psi}\Bigr{)}\psi_{\sigma}$
$\displaystyle+\phi^{*}\Bigl{(}\nu-\Pi_{\phi}\Bigr{)}\phi-h\Bigl{(}\phi^{*}\psi_{\uparrow}\psi_{\downarrow}-\phi\psi^{*}_{\uparrow}\psi^{*}_{\downarrow}\Bigr{)}\,\Biggr{]}\,,$
(9)

where the spin sum over $\sigma=\uparrow,\downarrow$ is implied. Equation 9
includes the full frequency and momentum dependent fermionic and bosonic
propagator and keep the classical interaction part. This approximation
incorporates the $s$-channel bubble-resummation with full fermionic
propagators in $\Sigma_{\psi}$ as the bosonic dimer propagator simply
constitutes the bubble chain.

Figure 2: Diagrammatic notation used throughout this work: Lines stand for
full propagators, small black dots stand for classical vertices, and larger
blue dots stand for full vertices.

Using 3 and 9 on both sides of the gap equations 8 leads us to the
diagrammatic expression for the non-perturbative self-consistent and
renormalised self-energies at finite temperature,

$\displaystyle\Pi_{\phi}(Q)=$
$\displaystyle\,h^{2}\int_{\bm{p}}\left[T\sum_{\omega_{m}}G_{\psi}(P)G_{\psi}(Q-P)-\frac{1}{2\bm{p}^{2}}\right]\,,$
$\displaystyle\Sigma_{\psi}(P)=$
$\displaystyle\,h^{2}\int_{\bm{q}}\,T\sum_{\epsilon_{m}}G_{\phi}(Q)G_{\psi}(Q-P)\,.$
(10)

with

$\displaystyle\int_{\bm{q}}=\int\frac{d^{3}q}{(2\pi)^{3}}\,,$ (11)

and similarly for the integral over $\bm{p}$. The subtraction with
$-1/(2\bm{p}^{2})$ on the right hand side of the equation for $\Pi_{\phi}$ in
10 arises from the spectral renormalisation. Equation 10 is given by Figure 1
with two classical vertices in each diagram. Note also, that in the fermionic
self-energy $\Sigma_{\psi}(P)$, the sum is over bosonic Matsubara frequencies
$\epsilon_{m}=2m\pi T$, whereas in the bosonic self-energy $\Pi_{\phi}(Q)$,
the sum is over fermionic Matsubara frequencies $\omega_{m}=(2m+1)\pi T$.

Evidently, the present truncation with $\Gamma^{(3)}=S^{(3)}$ in the gap
equations, 10 derived from 9, is equivalent to an $s$-channel bubble-
resummation, also known as the self-consistent T-matrix approximation [27, 24,
63]. However, the present functional approach can be generalised
straightforwardly in a self-consistent way by taking the full three-point
function $\Gamma^{(3)}$ into account, as has been done in [35, 39] for a
scalar theory. Respective results will be presented elsewhere.

###  III.2 Spectral representation

The spectral representation of the propagator is at the core of the spectral
functional approach developed in [35]. In the present work, we employ the
Källén-Lehmann (KL) representation of the full propagators [24, 64, 65]

$\displaystyle
G(\omega_{n},\bm{p})=\int_{-\infty}^{\infty}d\lambda\,\frac{\rho(\lambda,\bm{p})}{-i\omega_{n}+\lambda}\,,$
(12)

where $\rho(\lambda,\bm{p})$ is the frequency and momentum-dependent spectral
function. In this way, the spectral function acts as a linear response
function of the propagator, encoding the energy spectrum of the theory.
Equation 12 leads to the following inverse relation between the spectral
function and the retarded propagator,

$\displaystyle\rho(\omega,\bm{p})=\frac{1}{\pi}\,\mathrm{Im}\,G^{R}(\omega,\bm{p})\,,$
(13)

where $G^{R}(\omega,\bm{p})=G(-i(\omega+i0^{+}),\bm{p})$ and $\omega$ is now a
real frequency. Note that 12 always defines the spectral function, related to
the statistical function by the fluctuation-dissipation theorem, whereas 12
only holds true for fields that define physical asymptotic states. While this
is true in the present ultracold gas, it is far from obvious for theories
defined in terms of unphysical degrees of freedom such as gauge fields, for
more details see [41, 42, 43]. In summary, the present ultracold gas is
optimally suited for application of spectral functional approaches.

The existence of a spectral representation restricts all non-analyticities of
the propagator to lie on the real frequency and momentum axis. The fermion
spectral function satisfies the sum rule

$\displaystyle\int_{-\infty}^{\infty}d\lambda\,\rho_{\psi}(\lambda,\bm{p})=1\,.$
(14)

The fermionic spectral functions satisfy $\rho_{\psi}(\omega,\bm{p})\geq 0$.
In our approach, the frequency and momentum dependence of the bosonic inverse
propagator is only given by the self-energy, see Equation 8. Accordingly, its
spectral function is not that of a physical state and is not normalised to
unity as the fermionic one. Moreover, the respective spectral function is not
positive semi-definite, but satisfies
$\mathrm{sgn}(\omega)\rho_{\phi}(\omega,\bm{p})\geq 0$ [65]. Note that the
negative sign of the boson spectral function for negative frequencies
guarantees the positivity of the boson momentum distribution function [66].

Using the spectral representation 12 for the fermion and boson propagator, the
self-energies can be written in terms of spectral loop integrals. We find for
the bosonic self-energy

$\displaystyle\Pi_{\phi}(\epsilon_{n},\bm{q})=h^{2}\int_{\bm{p}}\left[\int_{\lambda_{1},\lambda_{2}}\rho_{\psi}(\lambda_{1},\bm{p})\rho_{\psi}(\lambda_{2},\bm{q}-\bm{p})\right.$
$\displaystyle\left.\times\,T\sum_{\omega_{m}}\frac{1}{i\omega_{m}-\lambda_{1}}\frac{1}{i(\epsilon_{n}-\omega_{m})-\lambda_{2}}-\frac{1}{2\bm{p}^{2}}\right]\,.$
(15)

The fermionic self-energy takes the form

$\displaystyle\Sigma_{\psi}(\omega_{n},\bm{p})=$
$\displaystyle\,h^{2}\int_{\lambda_{1},\lambda_{2},\bm{q}}\rho_{\phi}(\lambda_{1},\bm{q})\rho_{\psi}(\lambda_{2},\bm{q}-\bm{p})$
$\displaystyle\times
T\sum_{\epsilon_{m}}\frac{1}{i\epsilon_{m}-\lambda_{1}}\frac{1}{i(\epsilon_{m}-\omega_{n})-\lambda_{2}}\,.$
(16)

In both, 15 and 16, we have used 11 and

$\displaystyle\int_{\lambda}=\int_{-\infty}^{\infty}d\lambda\,.$ (17)

Importantly, the Matsubara sum in 16 and 15 can be carried out analytically.
This leaves us with symbolic expressions in terms of the argument
$\epsilon_{n}$ and $\omega_{n}$ for both self-energies that can be evaluated
at any complex frequency. The explicit spectral integral expressions for
$\Sigma_{\psi}$ and $\Pi_{\phi}$ can be found in Appendix B and C.

###  III.3 Evaluation at real frequencies

The regularised and coupled DSEs in Equation 8 can be evaluated for arbitrary
complex frequencies. For the extraction of the spectral functions with 13, we
choose $\omega_{n}=-i(\omega+i\varepsilon)$ with $\varepsilon\rightarrow
0^{+}$. The limit $\varepsilon\rightarrow 0^{+}$ is performed analytically
using the relation

$\displaystyle\frac{1}{x\pm i0^{+}}=P\left(\frac{1}{x}\right)\mp
i\pi\delta(x)\,,$ (18)

where $P(1/x)$ denotes the principal value of $1/x$. This allows us to write
the imaginary part of the retarded self-energies as

$\displaystyle\mathrm{Im}\,\Pi^{R}_{\phi}(\omega,\bm{q})$ $\displaystyle=\pi
h^{2}\int_{\lambda,\bm{p}}\rho_{\psi}(\omega-\lambda,\bm{p})\rho_{\psi}(\lambda,\bm{q}-\bm{p})$
$\displaystyle\times\left[1-n_{F}(\omega-\lambda)-n_{F}(\lambda)\right]\,,$
(19) $\displaystyle\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega,\bm{p})$
$\displaystyle=\pi
h^{2}\int_{\lambda,\bm{q}}\rho_{\phi}(\omega+\lambda,\bm{q})\rho_{\psi}(\lambda,\bm{q}-\bm{p})$
$\displaystyle\times\left[-n_{B}(\omega+\lambda)-n_{F}(\lambda)\right]\,.$
(20)

Here we have introduced the Fermi-Dirac distribution $n_{F}$ and Bose-Einstein
distribution $n_{B}$ with

$\displaystyle n_{F}(x)=\frac{1}{e^{x/T}+1}\,,\qquad
n_{B}(x)=\frac{1}{e^{x/T}-1}\,.$ (21)

A detailed derivation of the above formulae is deferred to Appendix B and
Appendix C. Note that the pole of $n_{B}(\omega+\lambda)$ in the fermion self-
energy is exactly cancelled by the zero crossing of the boson spectral
function $\rho_{\phi}$, see Section III.2.

From the imaginary part of the fermion self-energy, e.g., the real part can be
obtained via the Kramers-Kronig relation

$\displaystyle\mathrm{Re}\,\Sigma^{R}_{\psi}(\omega,\bm{p})$
$\displaystyle=\frac{1}{\pi}\,P\int_{\lambda}\frac{\mathrm{Im}\,\Sigma^{R}_{\psi}(\lambda,\bm{p})}{\lambda-\omega}\,.$
(22)

For the real part of the boson self-energy, the renormalisation of $\nu$ has
to be taken into account. The problematic vacuum part can be treated
analytically, see Appendix D. This set of equations allow for an iterative
calculation of the retarded self-energies and therefore, a direct
determination of the spectral functions which are presented in the next
section. For a detailed discussion of the numerical implementation, see
Appendix D.

##  IV Results

In this Section, we discuss the numerical results for the fully self-
consistent spectral functions as well as its use for the computation of the
density equation of state, the Tan contact and ejection rf spectra at
unitarity. In Section IV.1 we focus on the novel results obtained in real
frequencies. In Section IV.2 we use the spectral functions from Section IV.1
for computing radio-frequency (rf) spectra in comparison to experimental data.
Finally, in Section IV.3 we compute the equation of state and use the large
momentum dependence of the density for extracting the Tan contact.

###  IV.1 Spectral functions

We present numerical results for the self-consistent spectral functions of the
spin-balanced Fermi gas in the normal phase at different scattering lengths
and, in particular, at unitarity. Physical properties of the single-particle
fermion spectral function have been discussed extensively in the literature
[67, 20, 24]. In this work, we also discuss the bosonic two-particle (dimer)
spectral function which has not received much attention in the past. For
details on the numerical implementation, we refer the reader to Appendix D.

Figure 3 shows results for the fermionic spectral function $\rho_{\psi}$ at
$\bm{p}=0$ for the unitary Fermi gas at different temperatures, which are used
for the computation of the rf spectra in Section IV.2. Furthermore, Figure 4
shows results for the full frequency and momentum-dependent fermion spectral
function $\rho_{\psi}$, and Figure 5 shows the corresponding bosonic dimer
spectral functions $\rho_{\phi}$ for different interaction strengths
$1/(k_{F}a)$, with the Fermi momentum $k_{F}$ defined in 5. In Figures 3, 4
and 5, the frequency $\omega$ and momentum $\bm{p}$ are measured in the Fermi
energy $\varepsilon_{F}$ and $k_{F}$, respectively, and the fermionic
(bosonic) spectral functions have units of $\varepsilon_{F}^{-1}$
($\varepsilon_{F}^{-1/2}$). We state all results in dimensionless form. Note
that the bosonic dimer spectral function is not normalised and depends on the
choice of the Feshbach coupling $h$. This is connected to it being a pure
interaction exchange boson and not a real particle. To eliminate the
$h$-dependence, we multiply $\rho_{\phi}$ by $h^{2}/(8\pi)$ since
$h^{2}G_{\phi}/(8\pi)$ is the relevant quantity, which is related to the
scattering amplitude, see Appendix A.

Figure 3: Results for the fermionic spectral function
$\rho_{\psi}(\omega,\bm{0})\,\varepsilon_{F}$ for the unitary Fermi gas at
different temperatures, which are used for the computation of the rf spectra.

Our results are computed directly in real-time and agree qualitatively with
the reconstruction results from previous works [24] for the fermion spectral
function. Thus, our results corroborate the use of reconstruction methods for
spectral functions. It is well-known that numerical reconstruction problems
are ill-conditioned or at least have a high condition number. Moreover, the
respective ambiguities are specifically large in the low frequency regime at
finite temperature. The computation of dynamical properties such as the shear
viscosity requires a precise determination of the low frequency regime and its
transport peaks, which can be achieved with our spectral functional approach.
We envisage a two-step process with a direct use of Euclidean benchmark
results for a quality check of the direct spectral computation as well as
using its structural low frequency properties as an inductive bias for
spectral reconstructions.

(a) $(k_{F}a)^{-1}=-0.5$

(b) $(k_{F}a)^{-1}=0$

(c) $(k_{F}a)^{-1}=0.5$

Figure 4: Results for the fermionic spectral function
$\rho_{\psi}(\omega,\bm{p})\,\varepsilon_{F}$ for $T/T_{F}=0.56$ at different
interaction strengths $(k_{F}a)^{-1}$.
(a) BCS regime: $\beta\mu=0.5$ (b) unitarity: $\beta\mu=0.13$ (c) BEC regime:
$\beta\mu=-0.54$.

(a) $(k_{F}a)^{-1}=-0.5$

(b) $(k_{F}a)^{-1}=0$

(c) $(k_{F}a)^{-1}=0.5$

Figure 5: Results for the bosonic dimer spectral function
$h^{2}\rho_{\phi}(\omega,\bm{p})\,\sqrt{\varepsilon_{F}}/(8\pi)$ for
$T/T_{F}=0.56$ at different $(k_{F}a)^{-1}$.
(a) BCS regime: $\beta\mu=0.5$ (b) unitarity: $\beta\mu=0.13$ (c) BEC regime:
$\beta\mu=-0.54$.

The results for the bosonic dimer spectral function encode important
information about the physics of the ultracold Fermi gas and are consistent
with previous studies. First, we observe a very broad peak structure for weak
attractive interactions in the BCS regime, and a very sharp peak structure
with thermal broadening for strong attractive interactions in the BEC regime.
The latter property signals the tightly bound molecules on the BEC side. On
this side of the BCS-BEC phase diagram at $1/(k_{F}a)=0.5$, the boson spectral
function is very sharp and the system can be described as normal Bose liquid.
On the BCS side of the crossover at $1/(k_{F}a)=-0.5$, the fermion spectral
function is sharper and the system is described by a normal Fermi liquid.

Figure 6: Calculated ejection rf spectra $I(\omega)$ for the spin-balanced
unitary Fermi gas as a function of the reduced temperature $T/T_{F}$. Results
of this work (solid lines) are compared to experimental data from MIT
(circles) [7]. A Fourier broadening of $0.1\varepsilon_{F}$ to account for the
finite experimental resolution, and a right-shift by $0.09\varepsilon_{F}$ to
account for the final state interaction were applied.

Additionally, the boson spectral function reveals crucial information about
the critical region of the phase transition to the superfluid state. The onset
of superfluidity is marked by the divergence of the boson propagator at zero
frequency and momentum, $G^{-1}_{\phi}(0,\bm{0})=0$. This property is known as
the Thouless criterion [68]. Thus, the closer and sharper the peak gets at
zero frequency and momentum, the closer the system is to the phase transition,
until the spectral function eventually diverges at the critical temperature.
We refrain from a respective analysis and only present results in the normal
phase.

Computations at and below the phase transition are beyond the scope of the
present paper and shall be considered elsewhere. Note, however, that
respective work in the spectral functional approach has already been done in
relativistic systems both with the DSE and the fRG [35, 39, 40]. The
approximations there go beyond that used in the present work. In [35, 39, 40],
the full $s$-channel resummation of the four-point function has been
considered, as well as additional diagrams that arise in the presence of a
condensate $\phi_{0}$. In contrast to the Monte Carlo (MC) integration routine
used here, we used efficient integration routines in [35, 39, 40], aiming at
optimal convergence for low-dimensional integrals. The use of adaptive MC
integrators in the present work has to be understood as a test case for
higher-dimensional computations as will arise for non-trivial vertex
approximations beyond mere sums of momentum channels. This may be relevant in
the presence of competing order effects.

In summary we emphasise, that contrary to respective statements in the very
recent work [69], the numerical implementation of the present spectral
functional approach certainly is not less efficient than that used there.
Indeed, one of the reasons for setting up the spectral functional approach was
its numerical efficiency, which by now has been proven in abundance in [36,
35, 39, 40, 41, 42, 43, 44]. These successful tests include its convergence in
coupled sets of integral equations, as well as the evaluation of bound state
properties in the broken phase, as well as through the phase transition from
the symmetric to the broken phase. Specifically, self-consistent real-time
computations in the broken phase have so far solely been done in the present
approach and standard Bethe-Salpeter type computations. We add that, given the
structural similarities of the approaches used in [69], as well as [70], to
the spectral functional approach, it should be possible to extend the former
approaches towards the broken regime along the lines in [35, 39, 40].

###  IV.2 Radio-frequency spectroscopy

Figure 7: Peak positions ($E_{p}=-\omega_{p}$) and full width at half maximum
$\Gamma$ as a function of the reduced temperature $T/T_{F}$. The red vertical
dashed lines mark the superfluid phase transition.

Now we use our numerical results for the spectral functions to compute
experimentally measurable radio-frequency (rf) spectra [71, 72, 73]. We
compare the results from our approach, obtained for the spin-balanced Fermi
gas at unitarity with recent experimental data from MIT [7]. Specifically, we
apply the relation

$\displaystyle
I(\omega)=\int_{\bm{q}}\rho_{\psi}(\varepsilon_{\bm{q}}-\omega-\mu,\bm{q})\,n_{F}(\varepsilon_{\bm{q}}-\omega-\mu)\,,$
(23)

for the computation of rf spectra $I(\omega)$ from the fermion spectral
functions $\rho_{\psi}$, see [71, 74]. In Equation 23,
$\varepsilon_{\bm{q}}=\bm{q}^{2}$ is the classical fermion dispersion.

Note that the chemical potential $\mu(T)$ is temperature-dependent and has to
be determined self-consistently from the density equation 25 below, see also
[75, 23, 76, 77, 78]. More explicitly, the number density $n=1/(3\pi^{2})$ is
fixed by the choice $k_{F}=1$, and temperature is measured in units of
$T_{F}$. Consequently, the chemical potential $\mu(T)$ has to be chosen such
that the density remains constant. One readily observes that the rf spectrum
is normalised to the total density $n$ [76],

$\displaystyle n=2\int_{\lambda}I(\lambda)\,.$ (24)

Figure 6 shows our results in comparison with the experimental data from MIT
[7]. Apart from adjusting the peak heights, no fitting parameter have been
used. In order to account for the finite rectangular rf pulse duration and,
thus, a finite experimental resolution, a Fourier broadening of
$0.1\varepsilon_{F}$ has been applied. Additionally, the curves have been
right-shifted by an amount of $0.09\varepsilon_{F}$ to compensate residual
final state interactions [79].

Even after taking into account all these possible factors, the calculated rf
spectra do not fit the experimental data for higher temperatures very well.
This is also apparent in the comparison of the peak position and full width as
shown in Figure 7.

This discrepancy was already observed in [79] for the case of a highly spin-
imbalanced unitary Fermi gas with the non-self-consistent T-matrix approach,
and persists in the very recent update [80]. It has been argued in [79], the
misalignment may be caused by the missing vertex corrections, that may also be
relevant for pseudogap effects above $T_{c}$ [81]. Another possible
explanation might be the absence of a trap average, see e.g. [82].

Figure 8: Large $p$ behaviour of the momentum density $n(p)$ of the balanced
unitary Fermi gas at $T/T_{F}=0.9$ ($\beta\mu=-0.5$). The $1/p^{4}$ tail for
$p\gg k_{F}$ is clearly visible and the contact $C$ can be determined from the
fit. Lattice data from [83], see also [84].

###  IV.3 Equation of state and Tan contact

Finally, we use the present results for the spectral functions to compute the
equation of state for the unitary Fermi gas. Thermodynamic quantities, such as
the total particle density, can be calculated precisely in imaginary
frequencies without the need of analytic continuation. For this reason, it is
a good way to validate the new spectral approach against well tested and
robust imaginary-time calculations [23, 85, 84, 86, 87, 88, 89].

Figure 9: Universal contact $C\,\lambda^{4}_{T}$, where
$\lambda_{T}=\sqrt{4\pi\beta}$ is the thermal wavelength, of the spin-balanced
unitary Fermi gas as a function of dimensionless chemical potential
$\beta\mu$. Our results in real frequencies compare well with Luttinger-Ward
[23, 85] and Bold Diagrammatic Monte Carlo [46].

The total density $n$ of fermions at finite chemical potential $\mu$ and
temperature $T$ can be calculated from the spectral function via [78]

$\displaystyle
n=2\int_{\bm{p}}T\sum_{\omega_{n}}G_{\psi}(\omega_{n},\bm{p})e^{i\omega_{n}0^{+}}=2\int_{\bm{p}}n(\bm{p})\,,$
(25)

where the factor two accounts for both fermion species and $n(\bm{p})$ is the
momentum distribution function,

$\displaystyle
n(\bm{p})=\int_{\lambda}\rho_{\psi}(\lambda,\bm{p})\,n_{F}(\lambda)\,.$ (26)

For recent precision computations with functional approaches in the full phase
structure of ultracold gases, see [90]. The large-momentum tail of $n(\bm{p})$
is related to the Tan contact [91, 92]

$\displaystyle\lim_{p\rightarrow\infty}n(\bm{p})=\frac{C}{p^{4}}\,,$ (27)

where $p=\|\bm{p}\|$ and $C$ is the contact parameter that can be determined
from the fit. For a different computation in the DSE approach, see [30]. We
have compared the present results for the momentum density $n(\bm{p})$ and the
total density $n$ with lattice data generated for this purpose with a
Metropolis algorithm [83]. For results in the literature see also e.g. [84,
86, 88], using simulations with complex Langevin equations. The results for
$n(\bm{p})$ in the present spectral DSE approach compare well to the benchmark
results from the lattice simulation [83]. Moreover, in Figure 8 we also show a
fit to the large momentum asymptotics and we deduce that the asymptotic
$1/p^{4}$ tail is approached for

$\displaystyle\|\bm{p}\|\gtrsim p_{\textrm{as}}\,,\qquad
p_{\textrm{as}}=8k_{F}\,.$ (28)

The Tan contact $C$ is the slope of the asymptotics and, within the present
approach, its value for $T/T_{F}=0.9$ ($\beta\mu=-0.5$) is given by

$\displaystyle C/k^{4}_{F}\approx 0.084\,.$ (29)

Equation 29 and further values for the Tan contact $C(\beta\mu)$ for other
$\beta\mu$ are shown in Figure 9 and will be used for the following
determination of the density equation of state which is an additional
benchmark result: Figure 10 shows the results for the normalised density
$n/n_{0}$ as a function of dimensionless chemical potential $\beta\mu$ in
comparison with other approaches. The density $n_{0}$ of the non-interacting
Fermi gas is given by

$\displaystyle
n_{0}=2\int_{\bm{p}}n_{F}(p^{2}-\mu)=-2\mathrm{Li}_{3/2}(-e^{\beta\mu})/\lambda_{T}^{3}\,,$
(30)

where $\lambda_{T}=\sqrt{4\pi/T}$ is the thermal wavelength and
$\mathrm{Li}_{3/2}$ is a poly-logarithm function [93].

The computation of the total density 25 through the spectral function requires
its large momentum tail. For an efficient computation we split the momentum
integral in 25 into parts with radial momenta smaller and larger than
$p_{\textrm{as}}$ defined in 28,

$\displaystyle n=\sum_{\pm}\Delta n^{\pm}_{p_{\textrm{as}}}\,,\qquad\Delta
n^{\pm}_{p_{\textrm{as}}}=2\int_{\|\bm{p}\|\gtrless
p_{\textrm{as}}}n(\bm{p})\,.$ (31)

The low momentum part $\Delta n^{-}_{p_{\textrm{as}}}$ is computed directly
from the spectral function, while we shall use the results for the asymptotic
behaviour of $n(\bm{p})$ for the large momentum part $\Delta
n^{+}_{p_{\textrm{as}}}$: Equation 27 provides us with an analytic form of the
large momentum distribution and we obtain

$\displaystyle\Delta
n^{+}_{p_{\textrm{as}}}=2\int_{\|\bm{p}\|>p_{\textrm{as}}}\frac{C}{p^{4}}=\frac{C}{\pi^{2}p_{\textrm{as}}}\,.$
(32)

Here, we have used $p_{\textrm{as}}=10k_{F}$ in the explicit computation. Our
real-time method reproduces the results of the well tested Luttinger-Ward
approach [23, 85].

Figure 10: Normalised density $n/n_{0}$ of the spin-balanced unitary Fermi gas
as a function of dimensionless chemical potential $\beta\mu$. Results directly
obtained in real frequencies (this work) in comparison with experimental data
from MIT [15], lattice [83], see also [86, 88], and Luttinger-Ward [23, 85].

##  V Conclusion & Outlook

In the present work, we have extended the spectral functional approach [35]
for the non-perturbative and self-consistent computation of spectral
properties in non-relativistic systems. The approach was used for the direct
real-time computation of fermionic single-particle and bosonic dimer spectral
functions in the normal phase of a three-dimensional ultracold gas for a wide
range of scattering lengths. The present results are benchmarked with
observables in the unitary limit of a spin-balanced Fermi gas: the equation of
state, the Tan contact and ejection rf spectra in comparison to other
theoretical results and the experiment.

The present approach opens the path towards a wide range of applications,
including transport properties and the ab initio calculation of spectral
functions in the superfluid phase of ultracold Fermi gases [94, 95], and
through the phase transition, including the evaluation of critical properties
of the system. Its potential for computations in the broken phase of quantum
systems as well as in the critical regime has already been demonstrated in
relativistic systems and the respective computational advances can directly be
applied here. Furthermore, the present approach can also be readily applied to
general spin and mass-imbalance to study polaron physics or Fermi-Bose
mixtures [96, 97], where even the bosonic dimer spectral function can be
measured using molecular injection spectroscopy [98]. Finally, we envisage the
inclusion of vertex corrections or other classes of diagrams in systematic
extensions towards full quantitative precision. We hope to report on these
matters soon, but in particular on results in the superfluid phase.

Note added: After this work was completed, two works with complementary
approaches for real-time computations on the Keldysh contour were published,
[69, 70], that have been triggered partially by the present work. The
respective results compare well with our results, and in combination this
provides a respective non-trivial systematic error control in particular
within improvements. Furthermore, another independent work for the computation
of self-consistent spectral functions of Fermi polarons was published very
recently [80].

### Acknowledgements

We thank M. Bauer for providing us with lattice data for the equation of state
[83]. We thank F. Attanasio, M. Bauer, J. Braun, R. Kapust, J. Lang, and T.
Enss and R. Schmidt for discussions, and M. W. Zwierlein, B. Mukherjee and B.
Frank for providing their data. E. Dizer thanks J. Lang also for comparison of
numerical data of the present approach and Keldysh contour results. This work
is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the
Heidelberg STRUCTURES Excellence Cluster) and the Collaborative Research
Centre SFB 1225 (ISOQUANT).

##  Appendix A Renormalisation

The four-fermi model with the action 1 or equivalently its bosonised version
with the action 3 requires regularisation and renormalisation. While the
present model anyway has a physical cutoff given by the van der Waals length,
see also the discussion below 4, the respective standard renormalisation
procedure is accommodated in the present non-perturbative spectral approach
with spectral renormalisation, see [35, 37].

In the present approximation the renormalisation procedure is easily
implemented, following the standard renormalisation procedure, see e.g. [99,
100, 3]. To begin with, the physical renormalisation condition can be
extracted by the quantum analogue of the classical relation between the two-
body scattering length $a$ in the vacuum and the full two-to-two scattering
coupling at large distances, $\lambda$. In the present approximation this is
simply the bubble-resummed four-fermi interaction and 2 turns into

$\displaystyle-8\pi a=h^{2}G_{\phi}(0,\bm{0})\,.$ (33)

This can be written in terms of the full physical detuning parameter

$\displaystyle\nu=G^{-1}_{\phi}(0,\bm{0})=\Gamma^{(2)}_{\phi\phi^{*}}(0,\bm{0})\,.$
(34)

Equation 34 is the quantum analogue of the classical relation 4 and
constitutes a renormalisation condition of the dimer gap equation. A
renormalisation at vanishing momentum and frequency in vacuum yields

$\displaystyle\Pi_{\phi}(0,\bm{0})=0\,,$ (35)

which marks the detuning parameter in the effective action 9 as the physical
one in 34. The renormalisation condition 35 yields the renormalised boson DSE,

$\displaystyle G_{\phi}^{-1}(\epsilon_{n},\bm{q})$
$\displaystyle=-\frac{h^{2}}{8\pi a}$
$\displaystyle-h^{2}\int_{\bm{p}}\left[T\sum_{\omega_{m}}G_{\psi}(P)G_{\psi}(Q-P)-\frac{1}{2\bm{p}^{2}}\right]\,.$
(36)

We have used the $T$-independence of $h$ in the present approximation together
with the fact, that the fermion propagator is classical in the present
approximation, see [29, 100]. In this approximation the renormalised vacuum
result for $\Pi_{\phi}(Q)$ is readily obtained from 36 by evaluating the self-
energy loop integral with the classical fermion propagators,

$\displaystyle\Pi_{\phi}(\epsilon_{n},\bm{p})$
$\displaystyle=h^{2}\int_{\bm{q}}\left[\frac{1}{-i\epsilon_{n}+\varepsilon_{\bm{q}}+\varepsilon_{\bm{p-q}}}-\frac{1}{2\varepsilon_{\bm{q}}}\right]$
$\displaystyle=-\frac{h^{2}}{8\pi}\sqrt{-\frac{i\epsilon_{n}}{2}+\frac{\bm{p}^{2}}{4}}\,,$
(37)

for general complex frequencies $\epsilon_{n}\in\mathbbm{C}$. This concludes
our discussion of the renormalisation of the bosonised theory within the
present approximation.

##  Appendix B Boson self-energy calculation

In this Appendix, we discuss the explicit computations and analytic results
for the boson self-energy. Starting from 15, we define

$\displaystyle\Pi_{\phi}(\epsilon_{n},\bm{p})=h^{2}\int_{\bm{q}}\left[\int_{\lambda_{1},\lambda_{2}}\rho_{\psi}(\lambda_{1},\bm{q})\rho_{\psi}(\lambda_{2},\bm{p}-\bm{q})\right.$
$\displaystyle\left.\times\,I(\epsilon_{n},\lambda_{1},\lambda_{2})-\frac{1}{2\bm{q}^{2}}\right]\,,$
(38)

with the analytic sum over fermionic Matsubara frequencies
$\omega_{m}=(2m+1)\pi T$,

$\displaystyle I(\epsilon_{n},\lambda_{1},\lambda_{2})$
$\displaystyle=T\sum_{\omega_{m}}\frac{1}{i\omega_{m}-\lambda_{1}}\frac{1}{i(\epsilon_{n}-\omega_{m})-\lambda_{2}}$
$\displaystyle=\frac{1-n_{F}(\lambda_{1})-n_{F}(\lambda_{2})}{-i\epsilon_{n}+\lambda_{1}+\lambda_{2}}\,,$
(39)

Performing the analytic continuation $i\epsilon_{n}\rightarrow\omega+i0^{+}$
and taking the imaginary part, we end up with Section III.3.

For the first iteration, analytic expressions for the self-energy at finite
and zero temperature can be derived. Inserting the classical fermion spectral
function $\rho_{\psi}(\lambda,\bm{p})=\delta(\lambda-\bm{p}^{2}+\mu)$, we
obtain the well-known expression for the retarded boson self-energy [53, 101]

$\displaystyle\Pi^{R}_{\phi}(\omega,\bm{p})=$ $\displaystyle
h^{2}\int_{\bm{q}}\Biggl{[}\frac{1-n_{F}(\varepsilon_{\bm{q}}-\mu)-n_{F}(\varepsilon_{\bm{p-q}}-\mu)}{-\omega+\varepsilon_{\bm{q}}+\varepsilon_{\bm{p-q}}-2\mu-i0^{+}}$
$\displaystyle\hskip 28.45274pt-\frac{1}{2\varepsilon_{\bm{q}}}\Biggl{]}\,,$
(40)

where $\varepsilon_{\bm{p}}=\bm{p}^{2}$ is the classical momentum dispersion.
The boson self-energy can be separated in a temperature-independent and a
temperature-dependent part,

$\displaystyle\Pi^{R}_{\phi}=\Pi^{R,0}_{\phi}+\Pi^{R,T}_{\phi}\,.$ (41)

After the shift $\bm{q}\rightarrow\bm{q}+\bm{p}/2$, the angular integration is
trivial, and we obtain the finite expression for the vacuum part,

$\displaystyle\Pi^{R,0}_{\phi}(\omega,\bm{p})$
$\displaystyle=\int_{\bm{q}}\left[\frac{h^{2}}{-\omega+\varepsilon_{\bm{q}}+\varepsilon_{\bm{p-q}}-2\mu-i0^{+}}-\frac{h^{2}}{2\varepsilon_{\bm{q}}}\right]$
$\displaystyle=-\frac{h^{2}}{8\pi}\sqrt{-y-i0^{+}}\,,$ (42)

where we have defined $y=\omega/2-p^{2}/4+\mu$. The temperature-dependent part
can be obtained analogously. We note that the contribution from both Fermi
distributions is identical in the spin-balanced case and find

$\displaystyle\Pi^{R,T}_{\phi}(\omega,\bm{p})$
$\displaystyle=-h^{2}\int_{\bm{q}}\frac{n_{F}(\varepsilon_{\bm{q}}-\mu)+n_{F}(\varepsilon_{\bm{p-q}}-\mu)}{-\omega+\varepsilon_{\bm{q}}+\varepsilon_{\bm{p-q}}-2\mu-i0^{+}}$
$\displaystyle=-\frac{h^{2}}{4\pi^{2}}\int_{0}^{\infty}\frac{2\chi(q)\,q^{2}\,dq}{q^{2}-y-i0^{+}}\,,$
(43)

with the angle-integrated function

$\displaystyle\chi(q)=\int_{-1}^{1}\frac{dx}{2}\,n_{F}\left(\left[\bm{q}\pm\bm{p}/2\right]^{2}-\mu\right)=\begin{cases}n_{F}\left(q^{2}-\mu\right)&\quad\textrm{for}\quad
p=0\\\\[8.61108pt]
\frac{T}{2pq}\ln\left[\frac{n_{F}\left(\mu-q^{2}_{+}\right)}{n_{F}\left(\mu-q^{2}_{-}\right)}\right]&\quad\textrm{for}\quad
p\neq 0\end{cases}\,,$ (44)

with $q_{\pm}=q\pm p/2$. Here, $x=\cos(\theta_{\bm{pq}})$ and
$\theta_{\bm{pq}}$ is the angle between the vectors $\bm{p}$ and $\bm{q}$.
Note that this function yields a non-zero contribution in the vacuum for
$\mu>0$. In the limit of $T\rightarrow 0$, the Fermi functions
$n_{F}(x)\rightarrow\theta(-x)$, and the expression simplifies to

$\displaystyle\chi^{T=0}(q)=\begin{cases}\theta\left(\mu-q^{2}\right)&\quad\textrm{for}\quad
p=0\\\\[8.61108pt]
\frac{\theta\left(\mu-q^{2}_{-}\right)}{2pq}\left(\mu-q^{2}_{-}-(\mu-q^{2}_{+})\theta\left(\mu-q^{2}_{+}\right)\right)&\quad\textrm{for}\quad
p\neq 0\end{cases}\,.$ (45)

The real and imaginary part can be obtained analytically by using

$\displaystyle\frac{1}{x-i0^{+}}=P\left(\frac{1}{x}\right)+i\pi\delta(x)\,.$
(46)

There are the following two cases: If $y<0$, the integrals for the real part
are well-defined and $\mathrm{Im}\,\Pi^{R,T}_{\phi}=0$. If $y\geq 0$, the
imaginary part is given analytically by

$\displaystyle\mathrm{Im}\,\Pi^{R,T}_{\phi}(\omega,\bm{p})=-\frac{h^{2}}{4\pi}\int_{0}^{\infty}dq\,2\chi(q)q^{2}\delta(q^{2}-y)=-\frac{1}{4\pi}\sqrt{y}\chi(\sqrt{y})\,,$
(47)

and the real part $\mathrm{Re}\,\Pi^{R,T}_{\phi}$ can be obtained numerically
via a one-dimensional principal value integral, see [79].

For completeness, we also give the analytic expressions for the real part at
$T=0$. Since the contribution from both Fermi functions is identical, we
obtain for $\text{Re}\,\Pi^{R,T=0}_{\phi}$ at zero temperature ($\mu>0$),

$\displaystyle\text{Re}\,\Pi^{R,T=0}_{\phi}(\omega,\bm{p})=\,\frac{h^{2}}{4\pi^{2}}\left[\sqrt{\mu}-\frac{y-\mu+\frac{p^{2}}{4}}{2p}\log\left(\frac{y-\xi^{2}_{+}}{y-\xi^{2}_{-}}\right)-\sqrt{|y|}\begin{cases}\text{arctanh}\left(\frac{\xi_{-}}{\sqrt{y}}\right)+\text{arctanh}\left(\frac{\xi_{+}}{\sqrt{y}}\right)&\quad\textrm{for}\quad
y\geq 0\\\\[8.61108pt]
\text{arctan}\left(\frac{\xi_{-}}{\sqrt{|y|}}\right)+\text{arctan}\left(\frac{\xi_{+}}{\sqrt{|y|}}\right)&\quad\textrm{for}\quad
y<0\end{cases}\right]\,,$ (48)

with $\xi_{\pm}=\sqrt{\mu}\pm p/2$. The $p=0$ limit is given by

$\displaystyle\text{Re}\,\Pi^{R,T=0}_{\phi}(\omega,\bm{0})$
$\displaystyle=\frac{h^{2}}{2\pi^{2}}\begin{cases}\sqrt{\mu}-\sqrt{y}\,\text{arctanh}\left(\sqrt{\frac{\mu}{y}}\right)&\quad\textrm{for}\quad
y\geq 0\\\\[8.61108pt]
\sqrt{\mu}-\sqrt{|y|}\arctan\left(\sqrt{\frac{\mu}{|y|}}\right)&\quad\textrm{for}\quad
y<0\end{cases}\,.$ (49)

(a) $\mathrm{Re}\,\Pi_{\phi}^{R}(\omega,\bm{p})/\sqrt{\varepsilon_{F}}$

(b) $\mathrm{Im}\,\Pi_{\phi}^{R}(\omega,\bm{p})/\sqrt{\varepsilon_{F}}$

Figure 11: Real and imaginary part of the retarded boson self-energy
$\Pi_{\phi}^{R}(\omega,\bm{p})$ for the balanced Fermi gas at unitarity and
$T/T_{F}=0.56$ ($\beta\mu=0.13$).

The contribution from the next iterations has to be computed numerically and
gives rise to a finite correction on top of the analytic result. For more
details on the numerical treatment, see Appendix D. Finally, an exemplary plot
of the fully converged boson self-energy $\Pi_{\phi}^{R}(\omega,\bm{p})$ is
shown in Figure 11.

##  Appendix C Fermion self-energy calculation

For the fermion self-energy, the Matsubara sum can also be calculated
analytically. However, there is no analytic result for the first iteration
since the classical boson spectral function is not defined. Therefore the
fermion self-energy has to be computed numerically. Starting from Equation 16,
we define

$\displaystyle\Sigma_{\psi}(\omega_{n},\bm{p})=\int\limits_{\lambda_{1},\lambda_{2},\bm{q}}\rho_{\phi}(\lambda_{1},\bm{q})\rho_{\psi}(\lambda_{2},\bm{p}-\bm{q})I(\omega_{n},\lambda_{1},\lambda_{2})\,,$
(50)

with the analytic Matsubara sum over bosonic frequencies $\epsilon_{m}=2m\pi
T$,

$\displaystyle
I(\omega_{n},\lambda_{1},\lambda_{2})=T\sum_{\epsilon_{m}}\frac{1}{i\epsilon_{m}-\lambda_{1}}\frac{1}{i(\epsilon_{m}-\omega_{n})-\lambda_{2}}=\frac{-n_{B}(\lambda_{1})-n_{F}(\lambda_{2})}{-i\omega_{n}+\lambda_{1}-\lambda_{2}}\,.$
(51)

Performing the analytic continuation $i\omega_{n}\rightarrow\omega+i0^{+}$ and
taking the imaginary part, we recover Section III.3. Note that the pole in
$n_{B}(\lambda_{1})$ at $\lambda_{1}=0$ is exactly cancelled by the zero-
transition of $\rho_{\phi}(\lambda_{1})$.

For the numerical computation of the fermion self-energy, the high-frequency
tails contribute significantly. Therefore, we adopt a semi-analytic treatment
which is discussed in the next section and improves the determination of the
real part via Kramers-Kronig substantially. An exemplary plot of the fully
converged fermion self-energy $\Sigma_{\psi}^{R}(\omega,\bm{p})$ is shown in
Figure 12.

(a) $\mathrm{Re}\,\Sigma_{\psi}^{R}(\omega,\bm{p})/\varepsilon_{F}$

(b) $\mathrm{Im}\,\Sigma_{\psi}^{R}(\omega,\bm{p})/\varepsilon_{F}$

Figure 12: Real and imaginary part of the retarded fermion self-energy
$\Sigma_{\psi}^{R}(\omega,\bm{p})$ for the balanced Fermi gas at unitarity and
$T/T_{F}=0.56$ ($\beta\mu=0.13$).

##  Appendix D Numerical implementation

The numerical iteration is done with a standard iteration procedure for
integral equations. The first iteration step is done for the bosonic
propagator and only involves the classical fermion propagators. The following
one for the fermion propagator includes the results of the first step for the
bosonic propagator and again uses the classical fermion propagator. All the
remaining ones use the fully numerical results from the previous steps. The
respective workflow can be summarised as follows,

  1. 1.

First iteration of the boson propagator:  Start with the analytic expression
for $\text{Im}\,\Pi^{R}_{\phi}$ and calculate $\text{Re}\,\Pi^{R}_{\phi}$ on a
finite grid. To simplify the calculation, treat the divergent vacuum part
analytically, see below. Choose grid size large enough such that the numerical
corrections are small. This can be quantified in dependence of the
temperature.

  2. 2.

First iteration of the fermion propagator:  Take the classical fermion
spectral function (delta peak) and the semi-analytic boson spectral function
from the first iteration and calculate $\text{Im}\,\Sigma^{R}_{\psi}$ with
III.3 on a finite grid. From this, obtain $\text{Re}\,\Sigma^{R}_{\psi}$ as
described below.

  3. 3.

Further iterations of the fermion propagator:  Take the numerical spectral
functions for the fermion and boson and calculate the self-energy as in step
2.

  4. 4.

Further iterations of the boson propagator:  Take numerical spectral functions
for the fermions and calculate $\text{Im}\,\Pi^{R}_{\phi}$ on a finite grid.
Compute the difference to the non-self-consistent result and obtain the real
part. Outside the grid, glue smoothly to the semi-analytical non-self-
consistent result.

For the numerical calculation of the two-dimensional functions an adaptive
method is used. The sampling points are chosen automatically based on the
functional form. Less samples are taken in slow varying regions and more
samples are taken in faster varying regions. Additionally, the calculation of
these sampling points can be parallelised over multiple cores. The resulting
sampling points are then linearly interpolated. It is useful to simplify the
interpolation by transforming the functions onto the quadratic dispersion
relation. For the bosonic self-energy, this step can significantly improve the
resolution of sharp edges in the function. For the fermionic self-energy, it
can additionally help to ensure a large enough distance of the grid boundary
to the main peak, such that asymptotic behaviour is guaranteed for all
momenta.

The numerical integration of the three-dimensional self-energy loop integrals
is performed mainly with an adaptive Monte Carlo method. This allows for
maximal flexibility when dealing with highly peaked integrands in a
multidimensional space and can be generalised easily. However, it comes with
the downside of a long runtime in comparison to standard adaptive routines in
low dimensions like here. For this specific case, a different adaptive
integration routine using sparse grids is far more efficient. The one-
dimensional principal value integral for the real part of the self-energy is
computed efficiently using an adaptive quadrature integration.

For the representation of the numerical self-energy on the finite grid, an
adaptive grid is chosen for the imaginary and real part, separately. Typical
boundaries in frequency are $\omega=[-200,200]\,\varepsilon_{F}$ and in
momentum $p=[0,10]\,k_{F}$. Approximately 10.000-20.000 grid points are needed
to obtain stable numerical results.

###  D.1 Boson spectral function

In this Section, we detail the numerical calculation and representation of the
boson spectral function $\rho_{\phi}$. Since the vacuum part of the bosonic
self-energy is problematic, the numerical procedure requires suitable
subtraction schemes and analytic treatment.

As seen above, the vacuum contribution to the boson self-energy is given by

$\displaystyle\Pi^{0}_{\phi}(\epsilon_{n},\bm{p})=-\frac{h^{2}}{8\pi}\sqrt{-\frac{i\epsilon_{n}}{2}+\frac{\bm{p}^{2}}{4}-\mu}\,,$
(52)

With 52, the (retarded) imaginary part follows as

$\displaystyle\mathrm{Im}\,\Pi^{R,0}_{\phi}(\omega,\bm{p})=\frac{h^{2}}{8\pi}\sqrt{\frac{\omega}{2}-\frac{\bm{p}^{2}}{4}+\mu}\,.$
(53)

In order to obtain the full real part $\mathrm{Re}\,\Pi^{R}_{\phi}$ of the
self-energy via Kramers-Kronig numerically, we can subtract Equation 53 from
the total numerical imaginary part $\mathrm{Im}\,\Pi^{R}_{\phi}$ and find

$\displaystyle\mathrm{Re}\,\Pi^{R}_{\phi}(\omega,\bm{p})$
$\displaystyle=\mathrm{Re}\,\Pi^{R,0}_{\phi}(\omega,\bm{p})$
$\displaystyle+\frac{1}{\pi}\int_{\lambda}\frac{\mathrm{Im}\,\Pi^{R}_{\phi}(\lambda,\bm{p})-\mathrm{Im}\,\Pi^{R,0}_{\phi}(\lambda,\bm{p})}{\lambda-\omega}\,,$
(54)

where the real part of the vacuum solution is

$\displaystyle\mathrm{Re}\,\Pi^{R,0}_{\phi}(\omega,\bm{p})=-\frac{h^{2}}{8\pi}\sqrt{-\frac{\omega}{2}+\frac{\bm{p}^{2}}{4}-\mu}\,.$
(55)

With these formulas, we can obtain $\mathrm{Im}\,\Pi^{R}_{\phi}$ and
$\mathrm{Re}\,\Pi^{R}_{\phi}$ on a finite grid numerically. However, the
calculation of the real part requires also information about the high
frequency tails outside the numerical grid. In this case, the imaginary part
outside the grid is approximated by the analytic formula of the non-self-
consistent self-energy discussed in the previous Section,

$\displaystyle\mathrm{Im}\,\Pi^{R}_{\phi}(\omega,\bm{p})-\mathrm{Im}\,\Pi^{R,0}_{\phi}(\omega,\bm{p})\approx\mathrm{Im}\,\Pi^{R,T}_{\phi}(\omega,\bm{p})\,,$
(56)

for large $\omega$. Thus, the bosonic spectral function outside the grid is
approximated by the non-self-consistent (first iteration) spectral function
$\rho^{(1)}_{\phi}$. A different subtraction scheme would be to subtract the
whole non-self-consistent imaginary part straight away and only deal with
differences to the first iteration. In principle, both methods are equivalent
and work similar. For practical reasons, we choose the latter subtraction
scheme.

Another trick to improve the numerical calculation was already mentioned
above. Since the boson self-energy follows a quadratic dispersion relation, it
is practical to transform the functions onto the dispersion relation before
interpolation. This way, a better sampling and interpolation of the important
regions can be achieved. Afterwards, the interpolated function is shifted back
to the correct dispersion relation.

Since the fermion spectral functions in the self-energy integrals are very
peaked for lower temperatures and larger momenta, it might be useful to
subtract a broadened classical spectral function, or the first iteration, from
the further iterations in order to improve the integration for larger momenta.
The contribution from the subtracted spectral function has then to be taken
into account. Since the fermion spectral functions are represented on a finite
numerical grid, the bare delta peak contributions from outside the grid has to
be taken into account too. This is similar to the peak-tail splits from [35].

###  D.2 Calculation of the real part via Kramers-Kronig

It turns out that the calculation of the real part via Kramers-Kronig relation
is very sensitive to the high frequency tails of the imaginary part. We
accommodate for this fact with an extrapolation of the high frequency tail,
which allows us to estimate its contribution. For constant momentum $p$, the
large-frequency asymptotics of the fermion self-energy is captured by the
power law

$\displaystyle\lim_{\omega\rightarrow\infty}\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega)=a\,\omega^{-1/2}\,,$
(57)

where the constant $a$ is determined by a fitting routine. From this fit one
can calculate the missing contribution for the real part at a frequency
$\omega$ via

$\displaystyle\delta\mathrm{Re}\,\Sigma^{R}_{\psi}(\omega)=\frac{2a}{\pi\sqrt{\Delta}}\,{}_{2}F_{1}\left(\frac{1}{2},\frac{1}{2};\frac{3}{2};-\frac{\omega}{\Delta}\right)\,,$
(58)

where $\Delta=\omega_{\mathrm{max}}-\omega$ and $\omega_{\mathrm{max}}$ is the
largest frequency of the grid, and ${}_{2}F_{1}(a,b;c;z)$ being the
hypergeometric function [93]. Using this large-frequency contribution from
outside the grid, improves the determination of the fermion real part
significantly.

This asymptotic behaviour can be shown by similar arguments as in Ref. [73].
We consider the large frequency behaviour and rewrite the imaginary part of
the fermion self-energy III.3 as

$\displaystyle\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega)\sim\int_{\lambda_{1},\lambda_{2},\bm{q}}$
$\displaystyle\rho_{\phi}(\lambda_{1})\rho_{\psi}(\lambda_{2})\delta(\omega-\lambda_{1}+\lambda_{2})$
$\displaystyle\times\left[n_{B}(\lambda_{1})+n_{F}(\lambda_{2})\right]\,.$
(59)

In the limit $\omega\rightarrow\infty$, the delta function contributes only if
(a) $\lambda_{1}$ large and $\lambda_{2}$ small or (b) $\lambda_{2}$ large
negative and $\lambda_{1}$ small. In case (b), however, the fermion spectral
function $\rho_{\psi}$ vanishes for small momenta and
$\lambda_{2}\rightarrow-\infty$. Thus, only case (a) contributes and we are
left with

$\displaystyle\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega)\sim\int_{\lambda_{2},\bm{q}}\rho_{\phi}(\omega+\lambda_{2})\rho_{\psi}(\lambda_{2})n_{F}(\lambda_{2})\,,$
(60)

where we have used $n_{B}(\lambda_{1})\rightarrow 0$ for
$\lambda_{1}\rightarrow\infty$. As seen above, the large frequency behaviour
of the boson spectral function is dominated by $\rho_{\phi}(\omega)\sim
1/\sqrt{\omega}$. Additionally, the integrand vanishes for large negative
values of $\lambda_{2}$, since $\rho_{\psi}(\lambda_{2})$ vanishes, and for
large positive values of $\lambda_{2}$, since $n_{F}(\lambda_{2})$ vanishes.
Thus, the range of $\lambda_{2}$ is limited and we can write

$\displaystyle\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega)\sim\int_{\lambda_{2},\bm{q}}\rho_{\psi}(\lambda_{2})n_{F}(\lambda_{2})/\sqrt{\omega}\sim
1/\sqrt{\omega}\,,$ (61)

where we used that
$\int_{\lambda_{2},\bm{q}}\rho_{\psi}(\lambda_{2})n_{F}(\lambda_{2})=n$ is
finite in the last step. With the same argument one can show that
$\mathrm{Im}\,\Sigma^{R}_{\psi}(\omega)$ is mostly suppressed for large
negative frequencies.

###  D.3 Iterative procedure

Using III.3 and III.3, and the Kramers-Kronig relation, the retarded self-
energies can be computed in every iteration step and fed back into the next,
after extraction of the spectral function via Equation 13. The iterative
procedure is initialised with the classical fermion spectral function,

$\displaystyle\rho^{(0)}_{\psi}(\omega,\bm{p})=\delta(\omega-\bm{p}^{2}+\mu)\,.$
(62)

This initial guess is inserted into the spectral form of $\Pi_{\phi}(P)$ to
obtain the first iteration of the boson spectral function $\rho_{\phi}^{(1)}$.
Then, $\rho_{\phi}^{(1)}$ together with $\rho^{(0)}_{\psi}$ are inserted into
the spectral integral of $\Sigma_{\psi}(P)$ to obtain the first iteration of
the fermion spectral function $\rho_{\psi}^{(1)}$. Now, $\rho_{\psi}^{(1)}$
can be used to obtain the next iteration for the boson spectral function
$\rho_{\phi}^{(2)}$, and so on. In general, $\rho_{\psi}^{(i)}$ is used to
obtain $\rho_{\phi}^{(i+1)}$, and $\rho_{\phi}^{(i+1)}$ and
$\rho_{\psi}^{(i)}$ are used to obtain $\rho_{\psi}^{(i+1)}$. This iteration
is repeated until convergence is reached. We observe that around 5-20
iterations are needed to obtain a converged result with

$\displaystyle\int_{\lambda}\left\|\rho^{(i)}(\lambda,\bm{p})-\rho^{(i-1)}(\lambda,\bm{p})\right\|\lesssim
0.005\,,\quad\forall\bm{p}\,.$ (63)

As mentioned in the main text, the convergence is worse closer to the critical
temperature and the spectral functions may oscillate between intermediate
solutions. One can improve the convergence by iterating twice over the
fermions or updating the spectral functions only partially.

## References

  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008), arXiv:0704.3011 [cond-mat.other] .
  * Giorgini  _et al._ [2008] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of ultracold atomic Fermi gases, Rev. Mod. Phys. 80, 1215 (2008), arXiv:0706.3360 [cond-mat.other] .
  * Zwerger [2016] W. Zwerger, Strongly Interacting Fermi Gases, (2016), arXiv:1608.00457 [cond-mat.quant-gas] .
  * Greiner  _et al._ [2003] M. Greiner, C. A. Regal, and D. S. Jin, Emergence of a molecular Bose–Einstein condensate from a Fermi gas, Nature 426, 537 (2003).
  * Chin  _et al._ [2004] C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, J. H. Denschlag, and R. Grimm, Observation of the Pairing Gap in a Strongly Interacting Fermi Gas, Science 305, 1128 (2004).
  * Perali  _et al._ [2004] A. Perali, P. Pieri, L. Pisani, and G. C. Strinati, BCS-BEC crossover at finite temperature for superfluid trapped Fermi atoms, Phys. Rev. Lett. 92, 220404 (2004), arXiv:cond-mat/0311309 .
  * 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).
  * Strinati  _et al._ [2018] G. C. Strinati, P. Pieri, G. Röpke, P. Schuck, and M. Urban, The BCS–BEC crossover: From ultra-cold Fermi gases to nuclear systems, Phys. Rept. 738, 1 (2018), arXiv:1802.05997 [cond-mat.quant-gas] .
  * Schirotzek  _et al._ [2009] A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein, Observation of Fermi Polarons in a Tunable Fermi Liquid of Ultracold Atoms, Phys. Rev. Lett. 102, 230402 (2009).
  * Nascimbène  _et al._ [2009] S. Nascimbène, N. Navon, K. J. Jiang, L. Tarruell, M. Teichmann, J. McKeever, F. Chevy, and C. Salomon, Collective Oscillations of an Imbalanced Fermi Gas: Axial Compression Modes and Polaron Effective Mass, Phys. Rev. Lett. 103, 170402 (2009).
  * Massignan  _et al._ [2014] P. Massignan, M. Zaccanti, and G. M. Bruun, Polarons, dressed molecules and itinerant ferromagnetism in ultracold Fermi gases, Reports on Progress in Physics 77, 034401 (2014).
  * Scazza  _et al._ [2022] F. Scazza, M. Zaccanti, P. Massignan, M. M. Parish, and J. Levinsen, Repulsive Fermi and Bose Polarons in Quantum Gases, Atoms 10, 55 (2022), arXiv:2204.06984 [cond-mat.quant-gas] .
  * Zwierlein  _et al._ [2005] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H. Schunck, and W. Ketterle, Vortices and superfluidity in a strongly interacting Fermi gas, Nature 435, 1047 (2005).
  * Zwierlein  _et al._ [2006] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ketterle, Fermionic Superfluidity with Imbalanced Spin Populations, Science 311, 492 (2006), arXiv:cond-mat/0511197 .
  * Ku  _et al._ [2012] M. J. H. Ku, A. T. Sommer, W. C. Lawrence, and M. W. Zwierlein, Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas, Science 335, 563 (2012).
  * Taylor and Randeria [2010] E. Taylor and M. Randeria, Viscosity of strongly interacting quantum fluids: spectral functions and sum rules, Phys. Rev. A 81, 053610 (2010), arXiv:1002.0869 [cond-mat.quant-gas] .
  * Enss  _et al._ [2011] T. Enss, R. Haussmann, and W. Zwerger, Viscosity and scale invariance in the unitary Fermi gas, Annals Phys. 326, 770 (2011), arXiv:1008.0007 [cond-mat.quant-gas] .
  * 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), arXiv:2003.10338 [cond-mat.quant-gas] .
  * Greiner  _et al._ [2005] M. Greiner, C. A. Regal, and D. S. Jin, Probing the Excitation Spectrum of a Fermi Gas in the BCS-BEC Crossover Regime, Phys. Rev. Lett. 94, 070403 (2005), arXiv:cond-mat/0407381 .
  * Palestini  _et al._ [2012] F. Palestini, A. Perali, P. Pieri, and G. C. Strinati, Dispersions, weights, and widths of the single-particle spectral function in the normal phase of a Fermi gas, Phys. Rev. B 85, 024517 (2012).
  * Burovski  _et al._ [2006] E. Burovski, N. Prokof’ev, B. Svistunov, and M. Troyer, Critical Temperature and Thermodynamics of Attractive Fermions at Unitarity, Phys. Rev. Lett. 96, 160402 (2006), arXiv:cond-mat/0602224 .
  * Magierski  _et al._ [2009] P. Magierski, G. Wlazlowski, A. Bulgac, and J. E. Drut, The Finite Temperature Pairing Gap of a Unitary Fermi Gas by Quantum Monte Carlo, Phys. Rev. Lett. 103, 210403 (2009), arXiv:0801.1504 [cond-mat.stat-mech] .
  * 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), arXiv:cond-mat/0608282 .
  * 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).
  * Hu  _et al._ [2008] H. Hu, X.-J. Liu, and P. D. Drummond, Comparative study of strong-coupling theories of a trapped Fermi gas at unitarity, Phys. Rev. A 77, 061605 (2008).
  * 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).
  * Hanai and Ohashi [2014] R. Hanai and Y. Ohashi, Self-Consistent T-Matrix Approach to an Interacting Ultracold Fermi Gas with Mass Imbalance, Journal of Low Temperature Physics 175, 272 (2014).
  * 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).
  * Diehl and Wetterich [2006] S. Diehl and C. Wetterich, Universality in phase transitions for ultracold fermionic atoms, Phys. Rev. A 73, 033615 (2006), arXiv:cond-mat/0502534 .
  * Boettcher  _et al._ [2013] I. Boettcher, S. Diehl, J. M. Pawlowski, and C. Wetterich, Tan contact and universal high momentum behavior of the fermion propagator in the BCS-BEC crossover, Phys. Rev. A 87, 023606 (2013), arXiv:1209.5641 [cond-mat.quant-gas] .
  * Diehl  _et al._ [2007a] S. Diehl, H. Gies, J. M. Pawlowski, and C. Wetterich, Flow equations for the BCS-BEC crossover, Phys. Rev. A 76, 021602 (2007a), arXiv:cond-mat/0701198 .
  * Diehl  _et al._ [2007b] S. Diehl, H. Gies, J. M. Pawlowski, and C. Wetterich, Renormalisation flow and universality for ultracold fermionic atoms, Phys. Rev. A 76, 053627 (2007b), arXiv:cond-mat/0703366 .
  * Diehl  _et al._ [2010] S. Diehl, S. Floerchinger, H. Gies, J. M. Pawlowski, and C. Wetterich, Functional renormalization group approach to the BCS-BEC crossover, Annalen Phys. 522, 615 (2010), arXiv:0907.2193 [cond-mat.quant-gas] .
  * Jarrell and Gubernatis [1996] M. Jarrell and J. E. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo data, Physics Reports 269, 133 (1996).
  * Horak  _et al._ [2020] J. Horak, J. M. Pawlowski, and N. Wink, Spectral functions in the ${\phi}^{4}$-theory from the spectral Dyson-Schwinger equations, Phys. Rev. D 102, 125016 (2020).
  * Fehre  _et al._ [2023] J. Fehre, D. F. Litim, J. M. Pawlowski, and M. Reichert, Lorentzian Quantum Gravity and the Graviton Spectral Function, Phys. Rev. Lett. 130, 081501 (2023), arXiv:2111.13232 [hep-th] .
  * Braun  _et al._ [2022] J. Braun _et al._ , Renormalised spectral flows, (2022), arXiv:2206.10232 [hep-th] .
  * Berges [2004] J. Berges, Introduction to nonequilibrium quantum field theory, AIP Conf. Proc. 739, 3 (2004), arXiv:hep-ph/0409233 .
  * Horak  _et al._ [2023] J. Horak, F. Ihssen, J. M. Pawlowski, J. Wessely, and N. Wink, Scalar spectral functions from the spectral fRG, (2023), arXiv:2303.16719 [hep-th] .
  * Eichmann  _et al._ [2023] G. Eichmann, A. Gómez, J. Horak, J. M. Pawlowski, J. Wessely, and N. Wink, Bound states from the spectral Bethe-Salpeter equation, (2023), arXiv:2310.16353 [hep-ph] .
  * Horak  _et al._ [2022a] J. Horak, J. M. Pawlowski, and N. Wink, On the quark spectral function in QCD, (2022a), arXiv:2210.07597 [hep-ph] .
  * Horak  _et al._ [2022b] J. Horak, J. M. Pawlowski, and N. Wink, On the complex structure of Yang-Mills theory, (2022b), arXiv:2202.09333 [hep-th] .
  * Horak  _et al._ [2021] J. Horak, J. Papavassiliou, J. M. Pawlowski, and N. Wink, Ghost spectral function from the spectral Dyson-Schwinger equation, Phys. Rev. D 104, 074017 (2021), arXiv:2103.16175 [hep-th] .
  * Horak [2023] J. Horak,  _Realtime properties of QCD_ , Ph.D. thesis, U. Heidelberg (2023).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Rossi  _et al._ [2018] R. Rossi, T. Ohgoe, E. Kozik, N. Prokof’ev, B. Svistunov, K. Van Houcke, and F. Werner, Contact and Momentum Distribution of the Unitary Fermi Gas, Phys. Rev. Lett. 121, 130406 (2018), arXiv:1303.6245 [cond-mat.quant-gas] .
  * Gies and Wetterich [2002] H. Gies and C. Wetterich, Renormalization flow of bound states, Phys. Rev. D 65, 065001 (2002), arXiv:hep-th/0107221 .
  * Pawlowski [2007] J. M. Pawlowski, Aspects of the functional renormalisation group, Annals Phys. 322, 2831 (2007), arXiv:hep-th/0512261 .
  * Floerchinger and Wetterich [2009] S. Floerchinger and C. Wetterich, Exact flow equation for composite operators, Phys. Lett. B 680, 371 (2009), arXiv:0905.0915 [hep-th] .
  * Fu  _et al._ [2020] W.-j. Fu, J. M. Pawlowski, and F. Rennecke, QCD phase structure at finite temperature and density, Phys. Rev. D 101, 054032 (2020), arXiv:1909.02991 [hep-ph] .
  * Floerchinger  _et al._ [2008] S. Floerchinger, M. Scherer, S. Diehl, and C. Wetterich, Particle-hole fluctuations in the BCS-BEC Crossover, Phys. Rev. B 78, 174528 (2008), arXiv:0808.0150 [cond-mat.supr-con] .
  * Floerchinger  _et al._ [2010] S. Floerchinger, M. M. Scherer, and C. Wetterich, Modified Fermi-sphere, pairing gap and critical temperature for the BCS-BEC crossover, Phys. Rev. A 81, 063619 (2010), arXiv:0912.4050 [cond-mat.quant-gas] .
  * Nozieres and Schmitt-Rink [1985] P. Nozieres and S. Schmitt-Rink, Bose condensation in an attractive fermion gas: From weak to strong coupling superconductivity, J. Low Temp. Phys. 59, 195 (1985).
  * Sá de Melo  _et al._ [1993] C. A. R. Sá 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).
  * Chen  _et al._ [2005] Q. Chen, J. Stajic, S. Tan, and K. Levin, BCS–BEC crossover: From high temperature superconductors to ultracold superfluids, Physics Reports 412, 1 (2005).
  * Perali  _et al._ [2002] A. Perali, P. Pieri, G. C. Strinati, and C. Castellani, Pseudogap and spectral function from superconducting fluctuations to the bosonic limit, Phys. Rev. B 66, 024510 (2002).
  * Swanson [2010] E. S. Swanson, A Primer on Functional Methods and the Schwinger-Dyson Equations, AIP Conf. Proc. 1296, 75 (2010), arXiv:1008.4337 [hep-ph] .
  * Kopietz  _et al._ [2010] P. Kopietz, L. Bartosch, and F. Schütz,  _Introduction to the Functional Renormalization Group_, Vol. 798 (2010).
  * Metzner  _et al._ [2012] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schonhammer, Functional renormalization group approach to correlated fermion systems, Rev. Mod. Phys. 84, 299 (2012), arXiv:1105.5289 [cond-mat.str-el] .
  * Scherer  _et al._ [2011] M. M. Scherer, S. Floerchinger, and H. Gies, Functional renormalization for the BCS-BEC crossover, Phil. Trans. Roy. Soc. Lond. A 368, 2779 (2011), arXiv:1010.2890 [cond-mat.quant-gas] .
  * Boettcher  _et al._ [2012] I. Boettcher, J. M. Pawlowski, and S. Diehl, Ultracold atoms and the Functional Renormalization Group, Nucl. Phys. B Proc. Suppl. 228, 63 (2012), arXiv:1204.4394 [cond-mat.quant-gas] .
  * Dupuis  _et al._ [2021] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, The nonperturbative functional renormalization group and its applications, Phys. Rept. 910, 1 (2021), arXiv:2006.04853 [cond-mat.stat-mech] .
  * 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).
  * Abrikosov  _et al._ [1975] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski,  _Methods of Quantum Field Theory in Statistical Physics_ (Dover, 1975).
  * Fetter and Walecka [1971] A. L. Fetter and J. D. Walecka,  _Quantum Theory of Many-Particle Systems_ (McGraw-Hill, 1971).
  * Fratini and Pieri [2013] E. Fratini and P. Pieri, Single-particle spectral functions in the normal phase of a strongly attractive Bose-Fermi mixture, Phys. Rev. A 88, 013627 (2013).
  * Reichl and Mueller [2015] M. D. Reichl and E. J. Mueller, Quasiparticle dispersions and lifetimes in the normal state of the BCS-BEC crossover, Phys. Rev. A 91, 043627 (2015).
  * Thouless [1960] D. J. Thouless, Perturbation theory in statistical mechanics and the theory of superconductivity, Ann. Phys. 10, 553 (1960).
  * Johansen  _et al._ [2023] C. H. Johansen, B. Frank, and J. Lang, Spectral functions of the strongly interacting 3D Fermi gas, (2023), arXiv:2311.03953 [cond-mat.quant-gas] .
  * Enss [2023] T. Enss, Particle and pair spectra for strongly correlated Fermi gases: a real-frequency solver, (2023), arXiv:2311.05443 [cond-mat.quant-gas] .
  * Punk and Zwerger [2007] M. Punk and W. Zwerger, Theory of rf-Spectroscopy of Strongly Interacting Fermions, Phys. Rev. Lett. 99, 170404 (2007).
  * Stewart  _et al._ [2008] J. Stewart, J. Gaebler, and D. Jin, Using photoemission spectroscopy to probe a strongly interacting Fermi gas, Nature 454, 744 (2008).
  * Schneider  _et al._ [2009] W. Schneider, V. B. Shenoy, and M. Randeria, Theory of Radio Frequency Spectroscopy of Polarized Fermi Gases, (2009), arXiv:0903.3006 [cond-mat.other] .
  * Tsuchiya  _et al._ [2010] S. Tsuchiya, R. Watanabe, and Y. Ohashi, Photoemission spectrum and effect of inhomogeneous pairing fluctuations in the BCS-BEC crossover regime of an ultracold Fermi gas, Phys. Rev. A 82, 033629 (2010).
  * Veillette  _et al._ [2008] M. Veillette, E. G. Moon, A. Lamacraft, L. Radzihovsky, S. Sachdev, and D. E. Sheehy, Radio-frequency spectroscopy of a strongly imbalanced Feshbach-resonant Fermi gas, Phys. Rev. A 78, 033614 (2008).
  * Tajima and Uchino [2019] H. Tajima and S. Uchino, Thermal crossover, transition, and coexistence in Fermi polaronic spectroscopies, Phys. Rev. A 99, 063606 (2019).
  * Kagamihara  _et al._ [2019] D. Kagamihara, D. Inotani, and Y. Ohashi, Shear Viscosity and Strong-Coupling Corrections in the BCS–BEC Crossover Regime of an Ultracold Fermi Gas, J. Phys. Soc. Jap. 88, 114001 (2019), arXiv:1907.12414 [cond-mat.quant-gas] .
  * Pisani  _et al._ [2023] L. Pisani, M. Pini, P. Pieri, and G. C. Strinati, Peaks and widths of radio-frequency spectra: An analysis of the phase diagram of ultra-cold Fermi gases, (2023), arXiv:2311.00479 [cond-mat.quant-gas] .
  * Hu and Liu [2022] H. Hu and X.-J. Liu, Fermi polarons at finite temperature: Spectral function and rf spectroscopy, Phys. Rev. A 105, 043303 (2022), arXiv:2201.07872 [cond-mat.quant-gas] .
  * Hu and Liu [2023] H. Hu and X.-J. Liu, Spectral function of Fermi polarons at finite temperature from a self-consistent many-body $T$-matrix approach in real frequency, (2023), arXiv:2311.11554 [cond-mat.quant-gas] .
  * Li  _et al._ [2023] X. Li _et al._ , Observation and quantification of pseudogap in unitary Fermi gases, (2023), arXiv:2310.14024 [cond-mat.quant-gas] .
  * Schmidt  _et al._ [2012] R. Schmidt, T. Enss, V. Pietilä, and E. Demler, Fermi polarons in two dimensions, Phys. Rev. A 85, 021602 (2012).
  * Bauer [2023] M. Bauer, Lattice data from a Metropolis benchmark code, (2023).
  * Drut  _et al._ [2011] J. E. Drut, T. A. Lahde, and T. Ten, Momentum Distribution and Contact of the Unitary Fermi gas, Phys. Rev. Lett. 106, 205302 (2011), arXiv:1012.5474 [cond-mat.stat-mech] .
  * Frank [2018] B. Frank,  _Thermodynamics and Transport in Fermi Gases near Unitarity_ , Ph.D. thesis, TU München, T34 (2018).
  * Drut  _et al._ [2012] J. E. Drut, T. A. Lahde, G. Wlazlowski, and P. Magierski, The Equation of State of the Unitary Fermi Gas: An Update on Lattice Calculations, Phys. Rev. A 85, 051601 (2012), arXiv:1111.5079 [cond-mat.quant-gas] .
  * Drut [2012] J. E. Drut, Improved lattice operators for non-relativistic fermions, Phys. Rev. A 86, 013604 (2012), arXiv:1203.2565 [cond-mat.stat-mech] .
  * Rammelmüller  _et al._ [2018] L. Rammelmüller, A. C. Loheac, J. E. Drut, and J. Braun, Finite-temperature equation of state of polarized fermions at unitarity, Phys. Rev. Lett. 121, 173001 (2018), arXiv:1807.04664 [cond-mat.quant-gas] .
  * Jensen  _et al._ [2020] S. Jensen, C. N. Gilbreth, and Y. Alhassid, Contact in the Unitary Fermi Gas across the Superfluid Phase Transition, Phys. Rev. Lett. 125, 043402 (2020), arXiv:1906.10117 [cond-mat.quant-gas] .
  * Faigle-Cedzich  _et al._ [2023] B. M. Faigle-Cedzich, J. M. Pawlowski, and C. Wetterich, Towards quantitative precision in ultracold atoms with functional renormalisation, (2023), arXiv:2307.14787 [cond-mat.quant-gas] .
  * Tan [2008] S. Tan, Large momentum part of a strongly correlated Fermi gas, Annals of Physics 323, 2971 (2008).
  * Braaten and Platter [2008] E. Braaten and L. Platter, Exact Relations for a Strongly Interacting Fermi Gas from the Operator Product Expansion, Phys. Rev. Lett. 100, 205301 (2008).
  * Abramowitz and Stegun [1972] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions_ (National Bureau of Standards, 1972).
  * Engelbrecht  _et al._ [1997] J. R. Engelbrecht, M. Randeria, and C. A. R. Sá de Melo, BCS to Bose crossover: Broken-symmetry state, Phys. Rev. B 55, 15153 (1997).
  * Pieri  _et al._ [2004] P. Pieri, L. Pisani, and G. C. Strinati, BCS-BEC crossover at finite temperature in the broken-symmetry phase, Phys. Rev. B 70, 094508 (2004).
  * 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), arXiv:2104.14017 [cond-mat.quant-gas] .
  * Guo  _et al._ [2023] Y. Guo, H. Tajima, T. Hatsuda, and H. Liang, BCS-BCS crossover between atomic and molecular superfluids in a Bose-Fermi mixture, Phys. Rev. A 108, 023304 (2023), arXiv:2302.04617 [cond-mat.quant-gas] .
  * 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] .
  * Schonenberg  _et al._ [2017] L. M. Schonenberg, P. C. Verpoort, and G. J. Conduit, Effective-range dependence of two-dimensional Fermi gases, Phys. Rev. A 96, 023619 (2017).
  * 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), arXiv:1104.1379 [cond-mat.quant-gas] .
  * Punk [2010] M. Punk,  _Many-particle physics with ultracold gases_ , Ph.D. thesis, TU München, T34 (2010).