#  Quasiparticle disintegration in fermionic superfluids

Senne Van Loon<EMAIL_ADDRESS>TQC, Universiteit Antwerpen,
Universiteitsplein 1, B-2610 Antwerpen, België  School of Physics, Georgia
Institute of Technology, Atlanta, Georgia 30332, USA  Jacques Tempere  TQC,
Universiteit Antwerpen, Universiteitsplein 1, B-2610 Antwerpen, België  Lyman
Laboratory of Physics, Harvard University, Cambridge, Massachusetts 02138, USA
Hadrien Kurkjian  TQC, Universiteit Antwerpen, Universiteitsplein 1, B-2610
Antwerpen, België  Laboratoire de Physique Théorique, Université de Toulouse,
CNRS, UPS, France

###### Abstract

We study the fermionic quasiparticle spectrum in a zero-temperature superfluid
Fermi gas, and in particular how it is modified by different disintegration
processes. On top of the disintegration by emission of a collective boson
($1\to 2$, subject of a previous study, PRL 124, 073404), we consider here
disintegration events where three quasiparticles are emitted ($1\to 3$). We
show that both disintegration processes are described by a $t$-matrix self-
energy (as well as some highly off-resonant vacuum processes), and we
characterize the associated disintegration continua. At strong coupling, we
show that the quasiparticle spectrum is heavily distorted near the $1\to 3$
disintegration threshold. Near the dispersion minimum, where the
quasiparticles remain well-defined, the main effect of the off-shell
disintegration processes is to shift the location of the minimum by a value
that corresponds to the Hartree shift in the BCS limit. With our approximation
of the self-energy, the correction to the energy gap with respect to the mean-
field result however remains small, in contrast with experimental
measurements.

##  I Introduction

Although very useful to characterize the state and dynamics of a many-body
system [1, 2] the concept of quasiparticles is quite often only an
approximation of its spectral properties. This is particularly so in three-
dimensional ergodic systems where quasiparticles are affected by several decay
channels whose energy-momentum conservation constraint is met by exploiting
the angular degrees of freedom [3]. The resonance of the Green’s function at
the eigenenergy of the quasiparticle is then broadened by a nonzero damping
rate and diluted by a smaller spectral weight. The proximity of a decay
threshold often distorts the quasiparticle dispersion relation and in some
special cases it can entirely wash away the resonance [4]. These universal
phenomena affect systems as diverse as normal [5, 6] or superfluid Fermi
liquids [7, 8, 9], superconductors [10], rotonic systems [11, 12], or nuclear
matter [13, 14, 15].

This article is devoted to fermionic pair condensates, which support two kinds
of elementary excitations: the fermionic quasiparticles, describing unpaired
fermions surrounded by the condensate of pairs, and the bosonic modes
describing the collective motion of the pairs, such as sound waves in neutral
gases [16, 17, 18], or plasma oscillations in charged systems [19, 10]. In
addition, a high-energy collective branch can exist slightly above the pair-
breaking threshold [20, 21, 22, 23, 24]. Throughout this article, we restrict
the word “quasiparticles” to the single-particle (fermionic) elementary
excitations, and “collective modes” for the pair excitations. The properties
of both kinds of elementary excitations are essential to understand transport
[25] and dissipation [26] phenomena, or the weakly-excited dynamics of these
fermionic condensates [27, 28, 29].

Neutral fermionic condensates are now routinely prepared using ultracold
fermions [30, 31, 32, 33, 34, 9]. On this highly versatile platform, one can
study the entire crossover from a BCS-like superfluid, to a Bose-Einstein
condensate (BEC) of closely bound pairs by tuning the $s$-wave scattering
length describing the interaction between fermionic atoms in different
hyperfine states. By changing the magnetic field over a Feshbach resonance, it
is even possible to reach the unitary regime where interactions are resonant.
Using the set of techniques available to ultracold atomic experiments, such as
Bragg [35, 36, 37, 38] or rf-spectroscopy [39, 40, 41, 42, 23], flat-bottom
potentials [43, 44] or interaction quenches, numerous experimental results
have been accumulated on the collective mode spectrum [36, 37, 38], single-
particle excitations [42, 45], and the energy gap [41, 38], some of them being
already beyond our state-of-the-art theoretical understanding.

At zero temperature, the mean-field BCS theory [46] describes non-interacting,
undamped quasiparticles whose gap coincides with the order parameter $\Delta$.
For collective modes, the RPA [19] predicts a similarly undamped phononic
collective branch [16, 17]. It is now established that beyond the mean-field
approximation the ideal quasiparticle picture only works at low energies. Away
from their dispersion minimum, the quasiparticles can disintegrate [47, 48,
49], resulting in a nonzero damping. The gapped energy spectrum of the
fermionic excitations makes sure that the main decay mechanism at low energies
is the disintegration by emission of a phonon [50, 51] (a sound wave of the
Fermi gas, no to be confused with a lattice phonon in solid-state physics).
This kind of interaction between fermionic quasiparticles and collective modes
was considered in several studies, either from the point of view of collective
modes [26, 50, 51, 52], or to study their effect on the quasiparticle
properties [53, 54, 49]. Although the order parameter cannot be immediately
related to the quasiparticle gap beyond BCS theory, we also note that there
exists several calculations of this quantity [55, 56, 41, 57] beyond the mean-
field BCS approximation.

In a previous study, we studied corrections to the quasiparticle dispersion
relation focusing on the effect of the disintegration by emission of a
collective mode ($1\to 2$). Here we add the effect of disintegration into
three quasiparticles ($1\to 3$). Both disintegration processes are naturally
described by the $t$-matrix self-energy [58, 54, 59]. In Sec. II, we show that
the $1\to 2$ and $1\to 3$ disintegrations are respectively associated with the
poles and the branch cuts of the pair propagator, we derive the corresponding
coupling amplitudes in the perturbative limit, and discuss the structure of
the disintegration continua. We then study the full quasiparticle Green’s
function (Sec. III), where the two processes interfere, and extract the
corrected quasiparticle energy and lifetime. At strong coupling, we observe a
strong distortion of the quasiparticle resonance around the $1\to 3$
disintegration threshold $3\Delta$. Finally, we characterize the undamped low-
energy region of the quasiparticle dispersion by measurable quantities such as
the gap, effective mass and location of the dispersion minimum.

##  II Single-particle propagator

We consider a two-component Fermi gas with neutral atoms of mass $m$,
interacting with a short-range potential. For a dilute and ultracold gas the
true interaction potential can be replaced by an effective attractive contact
potential $g_{0}\delta(\textbf{r}-\textbf{r}^{\prime})$. The coupling constant
$g_{0}$ should be renormalized [34]

$\frac{1}{g_{0}}=\frac{m}{4\pi\hbar^{2}a}-\frac{1}{V}\sum_{{\bf{k}}}\frac{m}{\hbar^{2}{\bf{k}}^{2}}$
(1)

to keep the pair propagator finite, where $a$ is the $s$-wave scattering
length that determines the interaction regime.

Quite generally, when a many-body system supports quasiparticles, they appear
as poles of the single-particle Green’s function $G$. In the superfluid phase
of a two-component Fermi gas, the single-particle Green’s function is a two-
by-two matrix with particle and hole Green’s functions on the diagonal and
anomalous Green’s function on the codiagonal. In this situation, the poles of
$G$ are given by 111We assume here that the matrix elements of $G^{-1}$ do not
vanish individually, which is generally the case.:

$\text{det}G^{-1}({\bf{k}},z_{\bf{k}})=0$ (2)

Here ${\bf{k}}$ is the wave vector of the quasiparticle (we consider a
homogeneous system for which ${\bf{k}}$ is a good quantum number) and
$z_{\bf{k}}$ its eigenenergy, possibly taking complex values. To fully
characterize the quasiparticle resonance in the Green’s function, one should
also introduce the matrix residue $Z_{\bf{k}}$ associated to the pole in
$z_{\bf{k}}$, such that

$G({\bf{k}},z)\underset{z\to z_{\bf{k}}}{\sim}\frac{Z_{\bf{k}}}{z-z_{\bf{k}}}$
(3)

Note that the full Green’s function $G({\bf{k}},z)$ contains in general more
information than summarized by the quasiparticle spectrum.

The mean-field BCS theory provides a zeroth order approximation of the Green’s
function:

$\displaystyle G^{(0)}({\bf{k}},z)=$
$\displaystyle\frac{1}{z+\epsilon_{{\bf{k}}}}\begin{pmatrix}-V_{{\bf{k}}}^{2}&U_{{\bf{k}}}V_{{\bf{k}}}\\\
U_{{\bf{k}}}V_{{\bf{k}}}&-U_{{\bf{k}}}^{2}\end{pmatrix}-\frac{1}{z-\epsilon_{{\bf{k}}}}\begin{pmatrix}U_{{\bf{k}}}^{2}&U_{{\bf{k}}}V_{{\bf{k}}}\\\
U_{{\bf{k}}}V_{{\bf{k}}}&V_{{\bf{k}}}^{2}\end{pmatrix},$ (4)
$\displaystyle(G^{(0)}({\bf{k}},z))^{-1}=$
$\displaystyle\begin{pmatrix}-z+\xi_{\bf{k}}&\Delta\\\
\Delta&-z-\xi_{\bf{k}}\end{pmatrix}.$ (5)

In this approximation, the quasiparticle energy is

$\epsilon_{{\bf{k}}}=\sqrt{\xi_{{\bf{k}}}^{2}+\Delta^{2}},$ (6)

with $\xi_{{\bf{k}}}=\hbar^{2}k^{2}/2m-\mu$ the energy of the free fermions,
$\mu$ the chemical potential, and $\Delta$ the order parameter, acting also as
a gap of the dispersion relation $\epsilon_{\bf{k}}$. Note that the minimum of
the dispersion relation is reached at the wave vector

$\hbar k_{m}^{(0)}=\sqrt{2m\mu}.$ (7)

The Bogoliubov coefficients

$U_{\bf{k}}=\sqrt{\frac{1}{2}\left(1+\frac{\xi_{{\bf{k}}}}{\epsilon_{{\bf{k}}}}\right)},\quad
V_{\bf{k}}=\sqrt{\frac{1}{2}\left(1-\frac{\xi_{{\bf{k}}}}{\epsilon_{{\bf{k}}}}\right)},$
(8)

can be interpreted as the weights of the quasiparticle on the particle and
hole channels: when $U_{\bf{k}}$ increases the quasiparticle resembles more a
particle excitation, whereas when $V_{\bf{k}}$ increases it resembles a hole
excitation. These coefficients allow us to diagonalize the Green’s function:

$B_{\bf{k}}\left(G^{(0)}\right)^{-1}B_{\bf{k}}^{\dagger}=\begin{pmatrix}-z+\epsilon_{\bf{k}}&0\\\
0&-z-\epsilon_{\bf{k}}\end{pmatrix},\qquad\text{ with
}B_{\bf{k}}=\begin{pmatrix}U_{\bf{k}}&V_{\bf{k}}\\\
-V_{\bf{k}}&U_{\bf{k}}\end{pmatrix}.$ (9)

The fact that mean-field theory describes quasiparticles with purely real
eigenenergies is a serious limitation, since the quasiparticle lifetime plays
an important role in many dissipation mechanisms. This is one of the
motivations to introduce a correction to the Green’s function in the form of a
self-energy:

$G^{-1}({\bf{k}},\mathrm{i}\hbar\omega_{n})=(G^{(0)}({\bf{k}},\mathrm{i}\hbar\omega_{n}))^{-1}-\Sigma({\bf{k}},\mathrm{i}\hbar\omega_{n}),$
(10)

where $\Sigma$ is defined a priori only for imaginary fermionic Matsubara
frequencies $\mathrm{i}\hbar\omega_{n}=\mathrm{i}(2n+1)\pi k_{\mathrm{B}}T$,
with $n\in\mathbb{Z}$. In this work, we consider the so-called $t$-matrix
self-energy [58, 54, 59]:

$\Sigma_{ss^{\prime}}({\bf{k}},\mathrm{i}\hbar\omega_{n})=-ss^{\prime}\frac{1}{\beta
V}\sum_{{\bf{q}},m}G^{(0)}_{s^{\prime},s}({\bf{q}}-{\bf{k}},\mathrm{i}\hbar\nu_{m}-\mathrm{i}\hbar\omega_{n})\Gamma_{ss^{\prime}}({\bf{q}},\mathrm{i}\hbar\nu_{m}),$
(11)

with $s,s^{\prime}\in\\{+,-\\}$.

The self-energy of Eq. (11) can be studied with various degrees of self-
consistency [59]. Here, we choose to use a completely non-self-consistent
approach that allows for an analytic study of the different processes
contributing to the quasiparticle spectrum. Then, the ladder-resummed pair
propagator $\Gamma({\bf{q}},\mathrm{i}\hbar\nu_{m})$, with $\nu_{m}=2m\pi
k_{\mathrm{B}}T/\hbar$ bosonic Matsubara frequencies, is given in terms of the
mean-field single-particle propagators

$\displaystyle\Gamma_{ss^{\prime}}({\bf{q}},\mathrm{i}\hbar\nu_{m})$
$\displaystyle=\left(-\frac{1}{g_{0}}\delta_{ss^{\prime}}+N_{ss^{\prime}}({\bf{q}},\mathrm{i}\hbar\nu_{m})\right)^{-1},$
(12) $\displaystyle N_{ss^{\prime}}({\bf{q}},\mathrm{i}\hbar\nu_{m})$
$\displaystyle=\frac{1}{\beta
V}\sum_{{\bf{k}},n}G^{(0)}_{ss^{\prime}}({\bf{k}}+\frac{{\bf{q}}}{2},\mathrm{i}\hbar\omega_{n}+\frac{\mathrm{i}\hbar\nu_{m}}{2})G^{(0)}_{-s^{\prime},-s}({\bf{k}}-\frac{{\bf{q}}}{2},\mathrm{i}\hbar\omega_{n}-\frac{\mathrm{i}\hbar\nu_{m}}{2}).$
(13)

where $N$ is the bare pair propagator. We note that $\Gamma$ is an even
function of its energy argument,
$\Gamma({\bf{q}},-z_{q})=\Gamma({\bf{q}},z_{q})$, which ensures that its
spectrum is symmetric about the imaginary axis.

To calculate the pair propagator, we make use of the mean-field equation for
the order parameter

$\Delta=-\frac{g_{0}}{\beta
V}\sum_{{\bf{k}},n}G^{(0)}_{+-}({\bf{k}},\mathrm{i}\hbar\omega_{n})$ (14)

to replace $g_{0}$ in Eq. (12). In this way, the equations explicitly only
depend on the natural parameters $\Delta$ and $\mu$, in favor of the $s$-wave
scattering length $a$. The interaction regime is then completely determined by
fixing the dimensionless coefficient $\mu/\Delta$. To relate the parameters of
the theory $\Delta$ and $\mu$ to the experimentally more relevant interaction
parameter $1/k_{\mathrm{F}}a$ and the Fermi wavevector $k_{\mathrm{F}}$
(effectively fixing the density through $n=k_{\mathrm{F}}^{3}/3\pi^{2}$), the
order parameter equation (14) has to be solved together with the number
equation, which in the mean-field approximation is given by

$n=\frac{1}{V}\sum_{\bf{k}}\left[1-\frac{\xi_{{\bf{k}}}}{\epsilon_{{\bf{k}}}}\right].$
(15)

One can replace this number equation to include beyond-mean-field corrections,
which can lead to more quantitatively correct results. To avoid making a
choice for the approximation of the number equation, we present our main
results in terms of the coefficient $\mu/\Delta$. Whenever we relate this to
the coupling parameter $1/k_{\mathrm{F}}a$, we use the mean-field number
equation, as it allows for an analytic solution [16].

At zero temperature, the Matsubara sum in the self-energy [Eq. (11)] becomes
an integral over the imaginary frequencies
$\sum_{m}\to\int\limits_{-\mathrm{i}\infty}^{\mathrm{i}\infty}\frac{\mathrm{d}{z_{q}}}{2\pi\mathrm{i}}$.
Replacing $G^{(0)}$ by its BCS form (4), we cast the self-energy in the form:

$\displaystyle\Sigma({\bf{k}},\mathrm{i}\hbar\omega_{n})=\frac{1}{V}\sum_{{\bf{q}}}\int\limits_{-\mathrm{i}\infty}^{\mathrm{i}\infty}\frac{\mathrm{d}{z_{q}}}{2\pi\mathrm{i}}\Bigg{[}$
$\displaystyle\frac{1}{z_{q}-\mathrm{i}\hbar\omega_{n}+\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}V_{{\bf{k}}-{\bf{q}}}^{2}\Gamma_{++}({\bf{q}},z_{q})&U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\Gamma_{+-}({\bf{q}},z_{q})\\\
U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\Gamma_{-+}({\bf{q}},z_{q})&U_{{\bf{k}}-{\bf{q}}}^{2}\Gamma_{--}({\bf{q}},z_{q})\end{pmatrix}$
(16)
$\displaystyle-\frac{1}{z_{q}-\mathrm{i}\hbar\omega_{n}-\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}-U_{{\bf{k}}-{\bf{q}}}^{2}\Gamma_{++}({\bf{q}},z_{q})&U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\Gamma_{+-}({\bf{q}},z_{q})\\\
U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\Gamma_{-+}({\bf{q}},z_{q})&-V_{{\bf{k}}-{\bf{q}}}^{2}\Gamma_{--}({\bf{q}},z_{q})\end{pmatrix}\Bigg{]}.$

While computing the self-energy by integrating over $z_{q}$ on the imaginary
axis is a frequent and perfectly conceivable strategy (including numerically),
here we prefer to deform the integration contour towards the real axis, as
shown in Figure 1. This strategy allows us to isolate the contribution of each
pole or branch cut of the pair propagator to the self-energy, and to interpret
them separately in terms of elementary decay processes affecting the
quasiparticle. The straightforward analytic continuation of $\Gamma$ to
complex frequencies ($\mathrm{i}\hbar\nu_{m}\rightarrow z_{q}$) reveals its
analytic structure on the real axis (in black on Figure 1). Firstly, there are
two real poles in $z_{q}=\pm\hbar\omega_{\bf{q}}$, representing the phononic
collective branch of a neutral Fermi gas [19]. Secondly, there are two gapped
branch cuts at energies where it is possible to break up pairs. Opposite in
energy, those branch cuts are bound by the pair-breaking threshold
$\pm\epsilon_{c}({\bf{q}})$ with

$\epsilon_{c}({\bf{q}})=\min_{{\bf{k}}}[\epsilon_{{\bf{k}}-{\bf{q}}/2}+\epsilon_{{\bf{k}}+{\bf{q}}/2}]=\begin{cases}2\Delta\text{
if }q\leq 2k_{m}^{(0)}\\\ 2\epsilon_{q/2}\text{ if }q>2k_{m}^{(0)}\end{cases}$
(17)

On top of the analytic structure of $\Gamma$, the free single-particle
propagator $G^{(0)}$ gives rise to two additional complex poles, in
$z_{q}=\mathrm{i}\hbar\omega_{n}-\epsilon_{{\bf{k}}-{\bf{q}}}$ [for the first
term in Eq. (16)] and $\mathrm{i}\hbar\omega_{n}+\epsilon_{{\bf{k}}-{\bf{q}}}$
(for the second term). To evade those poles and focus on the singularities of
the pair propagator, we choose different contours for these two terms: for the
first term of Eq. (16) we deform the contour to the positive real axis, (see
contour $\mathcal{C}_{1}$ in Fig. 1), while for the second term the contour is
deformed to the negative real axis (contour $\mathcal{C}_{2}$ in Fig. 1).

Figure 1: The analytic structure of the integrand of the self-energy
$\Sigma({\bf{k}},\mathrm{i}\hbar\omega_{n})$. The Matsubara integral over the
imaginary axis can be performed by analytically continuing the Matsubara
frequencies and computing a contour integral. The contours are chosen in such
a way to avoid the poles of the single particle propagator, resulting in the
two keyhole contours $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ used respectively
in the first and second term of Eq. (16).

This procedure allows us to identify two physically distinct contributions to
the self-energy. The residues of the poles of the pair propagator will give
rise to a coupling between the fermionic quasiparticles and the collective
modes, such as the emission of a phonon by the quasiparticle, a process which
has been studied in-depth in Ref. [48]. Conversely, the branch cut
contributions are related to four-fermion processes, such as the decay of a
quasiparticle into three, as we explain below.

We expect that these disintegration processes, described by the self-energy of
Eq. (11), are the main decay channels of the quasiparticles at low energy.
Unfortunately, there is no small parameter that can be used to estimate the
importance of higher-order processes. Such processes could include the
emission of two or more bosons by a quasiparticle, or the disintegration of a
quasiparticle into five or more. The former should lead to a small correction
at low temperature [61], while the latter only become resonant above
$5\Delta$, and should thus be suppressed at low energies.

###  II.1 Contribution of the pole

Figure 2: Analytic structure of the self-energy $\Sigma({\bf{k}},z)$ on the
real axis. The two branch cuts on the positive real axis represent the real
quasiparticle decay processes (orange stripes: the process $\gamma\to\gamma+b$
associated to the bosonic continuum $[\epsilon_{\mathrm{th}}^{1\rightarrow
2},+\infty[$, green stripes: the process $\gamma\to\gamma+\gamma+\gamma$
associated to the fermionic continuum $[\epsilon_{\mathrm{th}}^{1\rightarrow
3},+\infty[$, in terms the fermionic quasiparticle $\gamma$ and a bosonic
collective mode $b$). Below those thresholds, a real solution $z_{\bf{k}}$ of
(2) may exists depending on the value of ${\bf{k}}$. The two branches on the
negative real axis represent the far-off-shell virtual processes (orange
stripes: the process $\emptyset\to\gamma+\gamma+b$, green stripes: the process
$\emptyset\to\gamma+\gamma+\gamma+\gamma$).

The residue of the $z_{q}$-integral in the poles of $\Gamma$ is given by

$\displaystyle\Sigma^{\mathrm{p}}({\bf{k}},z)=\frac{1}{V}\\!\sum_{{\bf{q}}}$
$\displaystyle\frac{\hbar}{\partial_{\omega}\det\Gamma^{-1}({\bf{q}},\omega_{{\bf{q}}})}\Bigg{[}$
(18)
$\displaystyle\frac{1}{z-\hbar\omega_{{\bf{q}}}-\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}V_{{\bf{k}}-{\bf{q}}}^{2}\check{N}_{{--}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)&-U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}N_{{+-}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)\\\
-U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}N_{{+-}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)&U_{{\bf{k}}-{\bf{q}}}^{2}\check{N}_{{++}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)\end{pmatrix}$
$\displaystyle-\frac{1}{-z-\hbar\omega_{{\bf{q}}}-\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}U_{{\bf{k}}-{\bf{q}}}^{2}\check{N}_{{++}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)&U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}N_{{+-}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)\\\
U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}N_{{-+}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)&V_{{\bf{k}}-{\bf{q}}}^{2}\check{N}_{{--}}\left({{\bf{q}},\omega_{{\bf{q}}}}\right)\end{pmatrix}\Bigg{]},$

where we have defined $\check{N}_{\pm\pm}=N_{\pm\pm}-1/g_{0}$. The absence of
singularities of $\Sigma$ outside the real axis allowed us to analytically
continue it in a natural way from Matsubara to complex frequencies
$\mathrm{i}\hbar\omega_{n}\to z$. Looking at the poles of the integrand as a
function of the wave vector ${\bf{q}}$, we remark that the first term of
$\Sigma^{\rm p}$ in Eq. (18) has a branch cut on the positive real axis for
$z\in\\{\hbar\omega_{{\bf{q}}}+\epsilon_{{\bf{k}}-{\bf{q}}}\text{ for
}{\bf{q}}\in\mathbb{R}^{3}\\}$, while the second term has a symmetric branch
cut on the negative real axis (see the orange stripes on Fig. 2). Restricting
to $\text{Re}\,z>0$ without loss of generality, we relate these branch cuts to
two elementary processes contributing to the correction of the single-particle
Green’s function: the term with
$z-\hbar\omega_{{\bf{q}}}-\epsilon_{{\bf{k}}-{\bf{q}}}$ in the denominator
describes the emission of a collective excitation $b$ by the quasiparticle
($\gamma\to\gamma+b$, in short $1\to 2$), while the term with
$-z-\hbar\omega_{{\bf{q}}}-\epsilon_{{\bf{k}}-{\bf{q}}}$ depicts the
spontaneous emission out of vacuum of two quasiparticles and a collective
excitation ($\emptyset\to\gamma+\gamma+b$, in short $0\to 3$), as shown in
figures 3a and 3b. The emission process can be resonant when the quasiparticle
energy is larger than the threshold value

$\epsilon_{\mathrm{th}}^{1\rightarrow
2}({\bf{k}})=\min_{\bf{q}}[\epsilon_{{\bf{k}}-{\bf{q}}}+\hbar\omega_{{\bf{q}}}],$
(19)

such that an undamped quasiparticle can exist only at energies
$z_{\bf{k}}<\epsilon_{\mathrm{th}}^{1\rightarrow 2}$. Conversely, the
simultaneous emission out of vacuum ($0\to 3$) is a far-off-shell process
acting on the quasiparticle only through an energy shift.

Figure 3: Different processes contributing to the fermionic self-energy. The
full lines depict quasiparticles, while the dotted lines are collective modes.
Both (a) and (b) result from the residue of the poles of the self-energy
integrand of Eq. (16), and describe the emission of a collective mode by a
quasiparticle and the spontaneous emission of two quasiparticles and a
collective mode from the vacuum. The diagrams in (c) and (d) describe four-
fermion processes, coming from the contribution of the branch cut. Only in the
weak-coupling limit can an explicit coupling amplitude be written down for
these processes.

###  II.2 Contribution of the branch cut

For the part of the contour along the branch cut, one obtains the following
self-energy

$\Sigma^{\mathrm{bc}}({\bf{k}},z)=\frac{1}{V}\sum_{{\bf{q}}}\int\limits_{\epsilon_{c}({\bf{q}})}^{+\infty}\mathrm{d}{z_{q}}\Bigg{[}-\frac{1}{z_{q}-z+\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}V_{{\bf{k}}-{\bf{q}}}^{2}\rho_{++}({\bf{q}},z_{q})&U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\rho_{+-}({\bf{q}},z_{q})\\\
U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\rho_{-+}({\bf{q}},z_{q})&U_{{\bf{k}}-{\bf{q}}}^{2}\rho_{--}({\bf{q}},z_{q})\end{pmatrix}\\\
-\frac{1}{-z_{q}-z-\epsilon_{{\bf{k}}-{\bf{q}}}}\begin{pmatrix}U_{{\bf{k}}-{\bf{q}}}^{2}\rho_{--}({\bf{q}},z_{q})&-U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\rho_{+-}({\bf{q}},z_{q})\\\
-U_{{\bf{k}}-{\bf{q}}}V_{{\bf{k}}-{\bf{q}}}\rho_{-+}({\bf{q}},z_{q})&V_{{\bf{k}}-{\bf{q}}}^{2}\rho_{++}({\bf{q}},z_{q})\end{pmatrix}\Bigg{]},$
(20)

where we have introduced the spectral density of the pair propagator

$\rho_{ss^{\prime}}({\bf{q}},z_{q})=\frac{\Gamma_{ss^{\prime}}({\bf{q}},z_{q}-\mathrm{i}0^{+})-\Gamma_{ss^{\prime}}({\bf{q}},z_{q}+\mathrm{i}0^{+})}{2\pi\mathrm{i}}.$
(21)

This time, the first line of Eq. (20) has a branch cut for $z$ larger than the
fermionic disintegration threshold $\epsilon_{\mathrm{th}}^{1\rightarrow
3}=\min_{{\bf{q}}}\min_{\omega_{q}\geq\epsilon_{c}({\bf{q}})}[\omega_{q}+\epsilon_{{\bf{k}}-{\bf{q}}}]$.
This branch cut is shown as green stripes on the positive real axis of Fig. 2.
The symmetric branch cut on the negative real axis stems from the second line
of Eq. (20) and describes a far-off-shell process. Using the definition (17)
of the pair-breaking threshold $\epsilon_{c}({\bf{q}})$ one can rewrite
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ in a form which suggests its
physical origin:

$\epsilon_{\mathrm{th}}^{1\rightarrow
3}({\bf{k}})=\min_{{\bf{k}}_{1},{\bf{k}}_{2}}[\epsilon_{{\bf{k}}_{1}}+\epsilon_{{\bf{k}}_{2}}+\epsilon_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}]\geq
3\Delta.$ (22)

The associated resonance condition
$z_{\bf{k}}=\epsilon_{{\bf{k}}_{1}}+\epsilon_{{\bf{k}}_{2}}+\epsilon_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}$
indicates that the first line of $\Sigma^{\mathrm{bc}}$ describes the decay of
the quasiparticle of momentum ${\bf{k}}$ into three other quasiparticles of
momenta ${\bf{k}}_{1}$, ${\bf{k}}_{2}$ and
${\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}$ ($\gamma\to\gamma+\gamma+\gamma$, in
short $1\to 3$, see the diagram on Fig. 3c). The second line of
$\Sigma^{\mathrm{bc}}$ then corresponds to the far-off-shell process where 4
quasiparticles appear out of vaccum
($\emptyset\to\gamma+\gamma+\gamma+\gamma$, in short $0\to 4$, see diagram
3d). This term will not contribute to the imaginary part of the self-energy,
but will nevertheless result in an energy shift of the quasiparticles.

###  II.3 Coupling amplitudes

To confirm our intuition that $\Sigma^{\rm bc}$ describes 4-fermion processes,
we now express the corrected Green’s function in terms of the coupling
amplitudes associated to these processes. Due to the non-perturbative nature
of the ladder-resummed pair propagator, this is not possible in the general
case, but only in the perturbative limit. In this limit, we replace $g_{0}$ by
$g=4\pi\hbar^{2}a/m$ in the pair propagator [Eq. (12)] and expand for $g\to
0$:

$\Gamma_{ss^{\prime}}=-g^{2}\left(\frac{\delta_{ss^{\prime}}}{g}+\frac{\delta_{ss^{\prime}}}{V}\sum_{\bf{k}}\frac{m}{\hbar^{2}{\bf{k}}^{2}}+N_{ss^{\prime}}\right)+\mathcal{O}(g^{3}),$
(23)

such that the spectral density $\rho$ takes the following form

$\displaystyle\rho({\bf{q}},z_{q})\overset{g\to
0}{\longrightarrow}-\frac{g^{2}}{V}\sum_{{\bf{p}}}\Bigg{[}$
$\displaystyle\begin{pmatrix}U_{{\bf{p}}_{+}}^{2}U_{{\bf{p}}_{-}}^{2}&-U_{{\bf{p}}_{+}}V_{{\bf{p}}_{+}}U_{{\bf{p}}_{-}}V_{{\bf{p}}_{-}}\\\
-U_{{\bf{p}}_{+}}V_{{\bf{p}}_{+}}U_{{\bf{p}}_{-}}V_{{\bf{p}}_{-}}&V_{{\bf{p}}_{+}}^{2}V_{{\bf{p}}_{-}}^{2}\end{pmatrix}\delta(z_{q}-\epsilon_{{\bf{p}}_{+}}-\epsilon_{{\bf{p}}_{-}})$
(24)
$\displaystyle-\begin{pmatrix}V_{{\bf{p}}_{+}}^{2}V_{{\bf{p}}_{-}}^{2}&-U_{{\bf{p}}_{+}}V_{{\bf{p}}_{+}}U_{{\bf{p}}_{-}}V_{{\bf{p}}_{-}}\\\
-U_{{\bf{p}}_{+}}V_{{\bf{p}}_{+}}U_{{\bf{p}}_{-}}V_{{\bf{p}}_{-}}&U_{{\bf{p}}_{+}}^{2}U_{{\bf{p}}_{-}}^{2}\end{pmatrix}\delta(z_{q}+\epsilon_{{\bf{p}}_{+}}+\epsilon_{{\bf{p}}_{-}})\Bigg{]}$

where we have introduced the short-hand notation
${\bf{p}}_{\pm}={\bf{p}}\pm{\bf{q}}/2$. The Dirac delta functions in the
spectral density can be used to perform the integral over $z_{q}$ in the self-
energy $\Sigma^{\mathrm{bc}}$ of Eq. (20), which exposes the elementary
processes involving 4 fermionic quasiparticles:

$\displaystyle\Sigma^{\mathrm{bc}}({\bf{k}},z)\overset{g\to
0}{\longrightarrow}\frac{g^{2}}{V^{2}}$
$\displaystyle\sum_{{\bf{k}}_{1},{\bf{k}}_{2},{\bf{k}}_{3}}\Bigg{[}\frac{\delta_{{\bf{k}}_{1}+{\bf{k}}_{2}+{\bf{k}}_{3},{\bf{k}}}}{z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{{\bf{k}}_{3}}}\begin{pmatrix}-V_{{\bf{k}}_{1}}^{2}U_{{\bf{k}}_{2}}^{2}U_{{\bf{k}}_{3}}^{2}&U_{{\bf{k}}_{1}}V_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}V_{{\bf{k}}_{2}}U_{{\bf{k}}_{3}}V_{{\bf{k}}_{3}}\\\
U_{{\bf{k}}_{1}}V_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}V_{{\bf{k}}_{2}}U_{{\bf{k}}_{3}}V_{{\bf{k}}_{3}}&-U_{{\bf{k}}_{1}}^{2}V_{{\bf{k}}_{2}}^{2}V_{{\bf{k}}_{3}}^{2}\end{pmatrix}$
(25)
$\displaystyle+\frac{\delta_{{\bf{k}}_{1}+{\bf{k}}_{2}+{\bf{k}}_{3},-{\bf{k}}}}{-z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{{\bf{k}}_{3}}}\begin{pmatrix}U_{{\bf{k}}_{1}}^{2}V_{{\bf{k}}_{2}}^{2}V_{{\bf{k}}_{3}}^{2}&U_{{\bf{k}}_{1}}V_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}V_{{\bf{k}}_{2}}U_{{\bf{k}}_{3}}V_{{\bf{k}}_{3}}\\\
U_{{\bf{k}}_{1}}V_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}V_{{\bf{k}}_{2}}U_{{\bf{k}}_{3}}V_{{\bf{k}}_{3}}&V_{{\bf{k}}_{1}}^{2}U_{{\bf{k}}_{2}}^{2}U_{{\bf{k}}_{3}}^{2}\end{pmatrix}\Bigg{]}.$

Finally, we rewrite the Green’s function in the quasiparticle basis
$\tilde{G}^{-1}({\bf{k}},z)=B_{\bf{k}}G^{-1}({\bf{k}},z)B_{\bf{k}}^{\dagger}$.
In the quasiparticle-quasiparticle channel, we obtain

$\tilde{G}_{++}^{-1}({\bf{k}},z)\overset{g\to
0}{\longrightarrow}-z+\epsilon_{\bf{k}}+\delta\epsilon({\bf{k}},z),$ (26)

with a complex energy shift $\delta\epsilon({\bf{k}},z)$ that combines all the
processes of Fig. 3:

$\displaystyle\delta\epsilon({\bf{k}},z)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{V}\sum\limits_{\bf{q}}\left[\frac{\mathcal{A}^{2}_{{\bf{k}}-{\bf{q}},{\bf{q}}}}{z-\epsilon_{{\bf{k}}-{\bf{q}}}-\hbar\omega_{\bf{q}}}-\frac{\mathcal{B}^{2}_{{\bf{k}},{\bf{q}}}}{-z-\epsilon_{{\bf{k}}+{\bf{q}}}-\hbar\omega_{\bf{q}}}\right]$
(27)
$\displaystyle+\frac{1}{V^{2}}\sum_{{\bf{k}}_{1},{\bf{k}}_{2}}\left[\frac{(\mathcal{A}^{1\to
3}_{{\bf{k}}_{;}{\bf{k}}_{1},{\bf{k}}_{2}})^{2}}{z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}}-\frac{(\mathcal{A}^{4\to
0}_{{\bf{k}}_{,}{\bf{k}}_{1},{\bf{k}}_{2}})^{2}}{-z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{-{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}}\right].$

The coupling amplitudes $\mathcal{A}_{{\bf{k}}-{\bf{q}},{\bf{q}}}$ of the
emission process $1\to 2$ and $\mathcal{B}_{{\bf{k}},{\bf{q}}}$ of the vacuum
process $0\to 3$ are described by $\Sigma^{\rm p}$ and discussed in detail in
Ref. [48]. Here, we add the contribution of the fermionic disintegration
process $1\to 3$ and vacuum process $0\to 4$ from $\Sigma^{\rm bc}$. Their
respective coupling amplitudes are given by

$\displaystyle\mathcal{A}^{1\to 3}_{{\bf{k}};{\bf{k}}_{1},{\bf{k}}_{2}}$
$\displaystyle=$ $\displaystyle
g(U_{{\bf{k}}}U_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}V_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}-V_{{\bf{k}}}V_{{\bf{k}}_{1}}V_{{\bf{k}}_{2}}U_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}),$
(28) $\displaystyle\mathcal{A}^{4\to 0}_{{\bf{k}},{\bf{k}}_{1},{\bf{k}}_{2}}$
$\displaystyle=$ $\displaystyle
g(U_{{\bf{k}}}U_{{\bf{k}}_{1}}V_{{\bf{k}}_{2}}V_{-{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}+V_{{\bf{k}}}V_{{\bf{k}}_{1}}U_{{\bf{k}}_{2}}U_{-{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}).$
(29)

Note that the vacuum terms $0\to 3$ and $0\to 4$ appear with a minus sign in
Eq. (27) because they occur in the ground state of the system
$|\psi_{0}\rangle$ and not in the state containing one quasiparticle
$\hat{\gamma}_{{\bf{k}}}^{\dagger}|\psi_{0}\rangle$.

The amplitudes (28)–(29) were derived in Ref. [54] [see Eqs. (B7–9) therein]
by expressing the interaction Hamiltonian $\hat{H}_{\rm int}=g\int
d^{3}r\hat{\psi}_{\uparrow}^{\dagger}({\bf{r}})\hat{\psi}_{\downarrow}^{\dagger}({\bf{r}})\hat{\psi}_{\downarrow}({\bf{r}})\hat{\psi}_{\uparrow}({\bf{r}})$
in terms of the quasiparticle creation-annihilation operators. In this weak-
coupling limit the branch cut contribution to $\Sigma$ thus describes the
damping and energy-shift due to the 4-fermion processes treated to second
order in perturbation theory.

For completeness, we also give the Green’s function in the quasiparticle-
quasihole channel:

$\displaystyle\tilde{G}_{+-}^{-1}({\bf{k}},z)$ $\displaystyle\overset{g\to
0}{\longrightarrow}$
$\displaystyle-\frac{1}{V}\sum\limits_{\bf{q}}\left[\frac{\mathcal{A}_{{\bf{k}}-{\bf{q}},{\bf{q}}}\mathcal{B}_{{\bf{k}}-{\bf{q}},{\bf{q}}}}{z-\epsilon_{{\bf{k}}-{\bf{q}}}-\hbar\omega_{\bf{q}}}+\frac{\mathcal{A}_{{\bf{k}},{\bf{q}}}\mathcal{B}_{{\bf{k}},{\bf{q}}}}{-z-\epsilon_{{\bf{k}}+{\bf{q}}}-\hbar\omega_{\bf{q}}}\right]$
(30)
$\displaystyle-\frac{1}{V^{2}}\sum_{{\bf{k}}_{1},{\bf{k}}_{2}}\left[\frac{\mathcal{A}^{1\to
3}_{{\bf{k}}_{;}{\bf{k}}_{1},{\bf{k}}_{2}}\mathcal{A}^{4\to
0}_{{\bf{k}}_{,}{\bf{k}}_{1},{\bf{k}}_{2}}}{z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{{\bf{k}}_{3}}}+\frac{\mathcal{A}^{1\to
3}_{{\bf{k}}_{;}{\bf{k}}_{1},{\bf{k}}_{2}}\mathcal{A}^{4\to
0}_{{\bf{k}}_{,}{\bf{k}}_{1},{\bf{k}}_{2}}}{-z-\epsilon_{{\bf{k}}_{1}}-\epsilon_{{\bf{k}}_{2}}-\epsilon_{{\bf{k}}_{3}}}\right].$

This non-zero off-diagonal matrix element shows that the quasiparticle
described by $\tilde{G}$ is a mix of the original BCS quasiparticle and
quasihole (or in other word that its weight on the particle and hole channels
differ from the $U_{\bf{k}}$ and $V_{\bf{k}}$ prescribed by BCS theory). We
note that $\tilde{G}_{+-}^{-1}$ is no longer small outside the limit $g\to 0$,
such that $(i)$ the quasiparticle energy $z_{\bf{k}}$ is no longer given by
$\epsilon_{\bf{k}}+\delta\epsilon({\bf{k}},z)$, and $(ii)$ the contribution of
$\Sigma^{\rm p}$ and $\Sigma^{\rm bc}$ to $z_{\bf{k}}-\epsilon_{\bf{k}}$ can
no longer be disentangled. This is the case even in the BCS limit as shown by
our numerical results in the next section. As such, all the numerical results
we present in Sec. III make use of the full self-energy of Eq. (16), valid for
all couplings, and not the perturbative expressions given in this section
222All code used for our numerical calculations can be found on
https://github.com/hkurkjian/CodeArXiv_2111_04692.

###  II.4 Structure of the disintegration continua

We conclude this section by studying the two disintegration continua
$\mathcal{C}_{1\to
2}({\bf{k}})=\\{\epsilon_{{\bf{k}}-{\bf{q}}}+\hbar\omega_{{\bf{q}}}\text{, for
}{\bf{q}}\in\mathbb{R}^{3}\\}$ and $\mathcal{C}_{1\to
3}({\bf{k}})=\\{\epsilon_{{\bf{k}}_{1}}+\epsilon_{{\bf{k}}_{2}}+\epsilon_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}\text{,
for }{\bf{k}}_{1},{\bf{k}}_{2}\in\mathbb{R}^{3}\\}$ and their intersection as
a function of ${\bf{k}}$. The lower edge of $\mathcal{C}_{1\to 3}$ (red solid
curve in Fig. 4) has a simple analytic expression [49]:

$\epsilon_{\mathrm{th}}^{1\rightarrow 3}=\begin{cases}3\Delta\text{ if }k\leq
3k_{m}^{(0)}\\\ 3\epsilon_{k/3}\text{ if }k>3k_{m}^{(0)}\end{cases}$ (31)

This threshold is thus a three-body equivalent of the pair-breaking threshold
(17). Its expression is obtained by a similar reasoning: as long as $k$ is
below $3k_{m}^{(0)}$, the angular degrees of freedom can be used to
accommodate three wave vectors ${\bf{k}}_{1}$, ${\bf{k}}_{2}$,
${\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}$ of norm $k_{m}^{(0)}$. Beyond this point
the minimum of the continuum is reached by symmetry for equal wave vectors of
norm $k/3$. The internal structure of the three-quasiparticle continuum,
particularly complex due to the number of internal degrees of freedom, is
studied in detail in Appendix A.

The lower edge of $\mathcal{C}_{1\to 2}$ does not have in general an analytic
expression because the eigenenergy $\omega_{\bf{q}}$ of the bosonic collective
branch does not. For $k$ sufficiently close to the dispersion minimum
$k_{m}^{(0)}$ (more precisely when the group velocity
$\mathrm{d}\epsilon_{k}/\mathrm{d}k$ is below the speed of sound $c$ [63,
51]), there exists a region $[k_{r1},k_{r2}]$ where the continuum edge is
reached in $q=0$ such that $\epsilon_{\mathrm{th}}^{1\rightarrow
2}({\bf{k}})=\epsilon_{\bf{k}}$. Outside this region the edge is reached at
nonzero wave vector $q_{\rm th}$ (orange line in Fig. 4) and energy
$\omega_{q_{\rm th}}$ (black line). For $k$ sufficiently close to the
dispersion minimum $k_{m}^{(0)}$, this edge remains below
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}=3\Delta$ such that the resonant
energies are divided in two sectors: sector A where only the $1\to 2$
disintegration is energetically allowed and sector $B$ where both $1\to 2$ and
$1\to 3$ disintegrations are allowed.

In the BCS regime (Fig. 4 is for $\mu/\Delta=4$, $1/k_{\rm F}a\simeq-0.91$)
the bosonic branch has a termination point ($q_{\rm sup}\gtrsim 2k_{F}$ in the
BCS regime) where it touches the pair-breaking continuum edge and disappears.
This leads to a saturation of $q_{\rm th}$ at $k\approx k_{r2}+q_{\rm sup}$,
and correspondingly to a rapid increase of
$\epsilon_{\mathrm{th}}^{1\rightarrow 2}$. This opens a new sector C where the
$1\to 3$ disintegration is resonant while the $1\to 2$ is not. This sector
however never contains the unperturbed quasiparticle energy
$\epsilon_{\bf{k}}$.

Figure 4: Diagram of the disintegration continua in function of the wave
vector $k$ (in units of the mean-field dispersion minimum
$k_{m}^{(0)}=\sqrt{2m\mu}/\hbar$) and energy $\epsilon$. Blue and purple area
(regions A and B): the $1\to 2$ disintegration continuum
$\\{\epsilon_{{\bf{k}}-{\bf{q}}}+\omega_{\bf{q}}\text{ for
}{\bf{q}}\in\mathbb{R}^{3}\\}$ caused by the poles of the pair propagator and
bounded from below by $\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ (blue solid
curve). The orange and black solid lines are respectively the wave vector
$q_{\rm th}$ and energy $\omega_{q_{\rm th}}$ of the boson emitted at the
continuum threshold (which satisfy $\epsilon_{{\bf{k}}-{\bf{q}}_{\rm
th}}+\omega_{q_{\rm th}}=\epsilon_{\mathrm{th}}^{1\rightarrow 2}$). Red and
purple area (regions B and C): the $1\to 3$ disintegration continuum
$\\{\epsilon_{{\bf{k}}_{1}}+\epsilon_{{\bf{k}}_{2}}+\epsilon_{{\bf{k}}-{\bf{k}}_{1}-{\bf{k}}_{2}}\text{
for }{\bf{k}}_{1},{\bf{k}}_{2}\in\mathbb{R}^{3}\\}$ originating from the
branch cut of the pair propagator and bounded from below by
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (red solid curve). The two continua
overlap in sector B (purple area), while sectors A and C are defined by
$\epsilon_{\mathrm{th}}^{1\rightarrow
2}(k)<\epsilon<\epsilon_{\mathrm{th}}^{1\rightarrow 3}(k)$ and
$\epsilon_{\mathrm{th}}^{1\rightarrow
3}(k)<\epsilon<\epsilon_{\mathrm{th}}^{1\rightarrow 2}(k)$ respectively. The
dashed line shows the mean-field quasiparticle dispersion relation
$\epsilon_{\bf{k}}$. Here, we took an interaction strength $\mu/\Delta=4$
($1/k_{\rm F}a\simeq-0.91$).

We note that the quasiparticle spectral function
$\rho_{G}({\bf{k}},\varepsilon)=\text{Im}G({\bf{k}},\varepsilon+\mathrm{i}0^{+})$
is still not analytic within sector A, B, or C due to the internal structure
of the disintegration continua, which contain angular points not shown on Fig.
4 (see for instance Appendix A). Still, the boundaries of those regions
constitute major singularities of the spectral function, across which it
cannot be fitted.

##  III Quasiparticle spectrum

###  III.1 Quasiparticle Green’s function

Expressions (18) and (20) of the self-energies were all given in the particle-
hole basis, where the single-particle propagator has both a positive and
negative energy contribution, weighted by the Bogoliubov coefficients, as can
be seen at the mean-field level in Eq. (4). The particle-propagator $G_{++}$
then consists of a mixture of quasiparticles with positive energy and
quasiholes with negative energy. When studying the limit $g\to 0$ (Sec. II.3)
we explained that the Green’s function takes a more intuitive form in the
quasiparticle basis, after transformation with the mean-field Bogoliubov
matrix $B_{\bf{k}}$. Here, we generalize this quasiparticle basis as the basis
which diagonalizes the Green’s function. This diagonalisation takes the form
of a generalized Bogoliubov transformation:

$\mathcal{B}({\bf{k}},z)G^{-1}({\bf{k}},z)\mathcal{B}^{\dagger}({\bf{k}},z)=\begin{pmatrix}\mathcal{G}^{-1}({\bf{k}},z)&0\\\
0&-\mathcal{G}^{-1}({\bf{k}},-z)\end{pmatrix}.$ (32)

This defines the quasiparticle Green’s function

$\mathcal{G}^{-1}({\bf{k}},z)=-z-\mathcal{D}({\bf{k}},z)+\mathcal{E}({\bf{k}},z)$
(33)

and the energy functionals

$\displaystyle\mathcal{D}({\bf{k}},z)$ $\displaystyle=$
$\displaystyle\frac{\Sigma_{++}({\bf{k}},z)+\Sigma_{--}({\bf{k}},z)}{2},$ (34)
$\displaystyle\mathcal{E}({\bf{k}},z)$ $\displaystyle=$
$\displaystyle\sqrt{\left(\xi_{{\bf{k}}}+\frac{\Sigma_{--}({\bf{k}},z)-\Sigma_{++}({\bf{k}},z)}{2}\right)^{2}+(\Delta-\Sigma_{+-}({\bf{k}},z))^{2}}.$
(35)

The energy $\mathcal{E}({\bf{k}},z)$ plays the same role as the mean-field
quasiparticle energy in Eq. (8), whereas $\mathcal{D}({\bf{k}},z)$ acts as a
displacement of the reference energy (as long as its dependence on $z$ is
omitted). The transfer matrix

$\mathcal{B}({\bf{k}},z)=\begin{pmatrix}\mathcal{U}({\bf{k}},z)&\mathcal{V}({\bf{k}},z)\\\
-\mathcal{V}({\bf{k}},z)&\mathcal{U}({\bf{k}},z)\end{pmatrix}$ (36)

has the same structure as the mean-field transfer matrix, but contains altered
Bogoliubov coefficients

$\mathcal{U}({\bf{k}},z)=\sqrt{\frac{1}{2}\left(1+\frac{\mathcal{X}({\bf{k}},z)}{\mathcal{E}({\bf{k}},z)}\right)},\quad\mathcal{V}({\bf{k}},z)=\sqrt{\frac{1}{2}\left(1-\frac{\mathcal{X}({\bf{k}}.z)}{\mathcal{E}({\bf{k}},z)}\right)},$
(37)

with
$\mathcal{X}({\bf{k}},z)=\xi_{{\bf{k}}}+({\Sigma_{--}({\bf{k}},z)-\Sigma_{++}({\bf{k}},z)})/{2}$.
Clearly, in the quasiparticle basis, it is enough to study the quasiparticle
propagator $\mathcal{G}$, rather than the full matrix Green’s function. We
note that some usual experimental probes such as rf-spectroscopy [33] measure
the Green’s function in the particle-hole basis. However, key properties of
the system, in particular dissipative properties, are directly sensitive to
the quasiparticle spectrum.

###  III.2 Perturbative and self-consistent solutions

In terms of the quasiparticle Green’s function, the equation on the
quasiparticle energy takes the simple form

$\mathcal{G}^{-1}({\bf{k}},z_{\bf{k}})=0.$ (38)

Instead of solving this self-consistent equation on $z_{\bf{k}}$, it is
tempting, given the form (33) of the quasiparticle Green’s function, to
perform a perturbative approximation on the energy functionals $\mathcal{D}$
and $\mathcal{E}$, that is, evaluate them in
$z\rightarrow\epsilon_{{\bf{k}}}+\mathrm{i}0^{+}$. Within this approximation,
the solution of Eq. (38) reads

$z_{\bf{k}}\simeq\mathcal{E}({\bf{k}},\epsilon_{{\bf{k}}}+\mathrm{i}0^{+})-\mathcal{D}({\bf{k}},\epsilon_{{\bf{k}}}+\mathrm{i}0^{+})\equiv
E_{{\bf{k}}}-\frac{\mathrm{i}\hbar}{2}\Gamma_{{\bf{k}}},$ (39)

where we have separated the real and imaginary part to distinguish between the
corrected energy $E_{\bf{k}}$ and the damping rate $\Gamma_{\bf{k}}$ of the
quasiparticles. In this perturbative approximation, the damping rate
$\Gamma_{\bf{k}}$ is nonzero if the mean-field energy is inside one of the
disintegration continua, such that at least one decay process is resonant.
When the mean-field energy $\epsilon_{{\bf{k}}}$ is above
$\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ but below
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (sector A), the only contribution to
a finite lifetime comes from the boson-emission process $1\to 2$ of Fig. 3a,
while for $\epsilon_{{\bf{k}}}>\epsilon_{\mathrm{th}}^{1\rightarrow 3}$
(sector B), also the four-fermion process $1\to 3$ of Fig. 3c lowers the
lifetime. We note that $\epsilon_{\bf{k}}$ never enters sector C (the
inequality $\epsilon_{\mathrm{th}}^{1\rightarrow
2}>\epsilon_{{\bf{k}}}>\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ is never
fulfilled), such that the fermionic disintegration process $1\to 3$ never acts
as the sole damping channel [49]. In a neutral gas, sector C thus has little
practical importance, but the situation would be reversed in a charged
superfluid where the bosonic branch acquires a large gap corresponding to the
plasma frequency [19].

To look for the exact poles of the quasiparticle Green’s function, Eq. (38)
should be solved self-consistently (this is not to be confused with a self-
consistent treatment of the self-energy $G^{(0)}\to G$ in Eqs. (11)–(13),
which is beyond the scope of the present work). Close to the minimum of the
corrected fermionic branch, below the threshold energies
$\epsilon_{\mathrm{th}}^{1\rightarrow 2},\epsilon_{\mathrm{th}}^{1\rightarrow
3}$, a self-consistent solution of Eq. (38) can be found on the real axis,
indicating well-defined quasiparticles with an infinite lifetime in this case.
Once the solution enters the continuum, the self-consistent solution
$z_{\bf{k}}$ obtains an imaginary part, and the quasiparticles are damped.
Extracting this complex solution requires in principle an analytic
continuation of $\mathcal{G}(z)$ through its branch cuts on the real axis (see
Fig. 2). Unfortunately, in this problem, we do not have access to an analytic,
or partially analytic [24, 52, 64], expression of the spectral density
$\mathcal{G}({\bf{k}},\varepsilon+\mathrm{i}0^{+})-\mathcal{G}({\bf{k}},\varepsilon-\mathrm{i}0^{+})$
on which we could rely to extend the function to the lower-half complex plane.
For this reason, we estimate the complex solution from its residuals just
above the real axis in $\varepsilon+\mathrm{i}0^{+}$. To do so, we fit the
quasiparticle Green’s function to a Lorentzian resonance

$\mathcal{G}({\bf{k}},\varepsilon+\mathrm{i}0^{+})\underset{\textrm{fit in
sector }X\textrm{ to
}}{\to}\frac{\Upsilon_{\bf{k}}^{(X)}}{z_{\bf{k}}^{(X)}-\varepsilon},$ (40)

and extract the fitted complex eigenenergy $z_{\bf{k}}^{(X)}$ and associated
residue $\Upsilon_{\bf{k}}^{(X)}$. Since the Green’s function shows a sharp
singularity at the continuum thresholds $\epsilon_{\mathrm{th}}^{1\rightarrow
2}$ and $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$, the fitting domain is
always comprised in one of the analyticity sectors $X=A,B$ from Fig. 4. We
note that the existence of several fitting sectors implies that two solutions
$z_{\bf{k}}^{(A)}$ and $z_{\bf{k}}^{(B)}$ coexist for the same value of $k$.
This reflects the structure of the analytic continuation, which has (at least)
two separate Riemann sheets corresponding to the continuation through sector
$A$ or $B$ [52, 65, 64] separated by branching points. This is thus not an
artifact of our fitting strategy. We shall see in the discussion of Fig. 7
that the coexistence of the two solutions reflects a physical phenomenon,
particularly near the point where $\text{Re}z_{\bf{k}}$ passes the fermionic
threshold $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$.

###  III.3 Numerical results

Figure 5: Quasiparticle spectrum in the BCS regime ($\mu/\Delta=4$,
$1/k_{\mathrm{F}}a\simeq-0.9$). The quasiparticle spectral density
$\mathrm{Im}[\mathcal{G}({\bf{k}},\varepsilon+\mathrm{i}0^{+})]$ is shown as a
function of the wavenumber $k$ and energy $\varepsilon$ in units of the order
parameter $\Delta$. This is compared to the self-consistent energy (red line
found by solving Eq. (38) in the undamped region, and by fitting the Green’s
function to Eq. (40) elsewhere), the perturbative energy given by Eq. (39)
(blue dashed line), and the mean-field energy (black dotted line). The
threshold energies $\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ (dark grey) and
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (light grey) are also shown to
distinguish the different regions of possible resonances. In the lower panel,
the imaginary part of the eigenenergy is plotted, showing the damping rate of
the quasiparticles.

In Fig. 5 we show results of the quasiparticle spectrum at $\mu/\Delta=4$
(thus in the shallow BCS regime, corresponding to
$1/k_{\mathrm{F}}a\simeq-0.9$). The imaginary part of the Green’s function is
exactly zero below the threshold energy $\epsilon_{\mathrm{th}}^{1\rightarrow
2}$. There, a real solution of $\mathcal{G}^{-1}({\bf{k}},z_{\bf{k}})=0$ can
be found, represented as the part of the red curve below
$\epsilon_{\mathrm{th}}^{1\rightarrow 2}$. It is clear that the minimum of the
energy is shifted toward higher values of the wavenumber $k$ with respect to
the mean-field energy, while the correction to the energy gap remains small. A
more in-depth study of the low-energy spectrum can be found in the next
section. Once the self-consistent solution reaches the first threshold
$\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ (on both sides of the minimum), the
eigenenergy $z_{\bf{k}}$ becomes complex, which translates into a broadened
peak in the spectral density. We keep track of the complex pole using the
fitted energy $z_{\bf{k}}^{(A)}$ from Eq. (40). The real part perfectly
continues the undamped solution (such that we represent it by the same red
curve on the top panel), illustrating the efficiency of our fitting strategy.
Further away from the dispersion minimum, the resonance reaches the second
threshold $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ and the red curve switches
from $z_{\bf{k}}^{(A)}$ to $z_{\bf{k}}^{(B)}$ to describe the resonance in
region $B$. The transition is again very smooth, and only translates into a
small kink in the damping rate (bottom panel). This is in stark contrast to
previous results of Ref. [48] taking into account only the $1\to 2$
disintegration, where the damping rate is sharply peaked when the eigenenergy
reaches the threshold $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$. This peaked
behavior is washed out by including the now resonant $1\to 3$ process in the
self-energy, which sharply reduces the quasiparticle lifetime. At larger $k$,
the damping rate peaks around $\hbar k\simeq 5.5\sqrt{2m\Delta}$, after which
it exhibits a large $1/k$ tail. In the limit $k\rightarrow\infty$, where the
interaction between the quasiparticle and the rest of the superfluid becomes
negligible, the quasiparticle energy tends to the kinetic energy of a free
fermion: $z_{\bf{k}}\sim\hbar^{2}k^{2}/2m$.

The self-consistent eigenenergy can be compared to the perturbative result of
Eq. (39), which is shown as a blue dashed line in Fig. 5. Both results remain
close to each other, although the energy gap is lowered in the perturbative
case with respect to the self-consistent solution. Most notably, the
perturbative spectrum incorrectly predicts a finite lifetime of the
quasiparticles at their energy minimum. This is due to the fact that the
resonance condition in this case is controlled by the mean-field energy
$\epsilon_{{\bf{k}}}$, which does not account for the shift in the energy
minimum. This is resolved with the self-consistent solution, which accurately
predicts well-defined quasiparticles at the energy minimum.

Figure 6: Bogoliubov coefficients in the BCS regime ($\mu/\Delta=4$,
$1/k_{\mathrm{F}}a\simeq-0.9$). The altered Bogoliubov weights defined in Eq.
(37) are shown in function of the wavenumber using the eigenenergy from
different approximations used in our study, namely the self-consistent (red
line) and perturbative (blue dashed line) solutions. The corrected
coefficients are shifted with respect to the BCS result (black dotted line),
in the same way as the minimum of the dispersion (see Fig. 5).

In order to reveal the underlying mixing of the particle and hole channels
contributing to the quasiparticle spectrum, we show in Fig. 6 the altered
Bogoliubov coefficients of Eq. (37). Similar to the BCS case, before the
dispersion minimum $\mathcal{V}$ is close to unity and $\mathcal{U}$ is small,
indicating that the quasiparticles resemble hole excitations, while at higher
${\bf{k}}$ values the situation is reversed, such that the quasiparticles are
well approximated by particle excitations. The region where the mixing is
largest corresponds to the dispersion minimum and is consequently shifted to
higher $k$ values compared to the BCS prediction. Note that when $z_{\bf{k}}$
is complex the self-consistent values of $\mathcal{U}$ and $\mathcal{V}$ are
evaluated in $\text{Re}\,z_{\bf{k}}$, rather than in $z_{\bf{k}}$ itself, as
we don’t have access to the analytical continuation of the self-energy. This
approximation works well as long as the imaginary part of $z_{\bf{k}}$ remains
small.

Figure 7: Quasiparticle spectrum at unitarity ($\mu/\Delta\simeq 0.86$). To
correctly account for different possible branches, the energy domain for which
we fit the quasiparticle Green’s function is split up for each value of the
wavenumber $k$. Three different regions can be recognized, separated by the
threshold energies $\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ (dark grey) and
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (light gray). In the first region
($\varepsilon<\epsilon_{\mathrm{th}}^{1\rightarrow 2}$) the quasiparticle
propagator is real, and a real self-consistent pole $z_{\bf{k}}$ can be found
(green line). Above the threshold $\epsilon_{\mathrm{th}}^{1\rightarrow 2}$,
the complex eigenenergy is computed from a fit as in Eq. (41), using only
energy values $\epsilon_{\mathrm{th}}^{1\rightarrow
2}<\varepsilon<\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (green line, sector
A), or $\epsilon_{\mathrm{th}}^{1\rightarrow 3}<\varepsilon$ (red line, sector
B). We furthermore compare the self-consistent solutions with the perturbative
(blue dashed line) and mean-field energy (black dotted line). In the lower
panel, we show the imaginary part of the eigenenergy, with the same color
code.

We perform a similar analysis at unitarity in Fig. 7, where the quasiparticle
spectral density is shown and compared with the eigenenergies from the self-
consistent and perturbative methods. Here, however, we notice a more abrupt
transition when the resonance reaches the second threshold
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$. To make this clear, we show on Fig.
8 a cross section of the Green’s function at $\hbar k=2.16\sqrt{2m\Delta}$
where we noticed a marked angular point in
$\epsilon=\epsilon_{\mathrm{th}}^{1\rightarrow 3}$. At the same time, the
resonance in both sectors $A$ and $B$ loses its Lorentzian behavior (in other
words, the fit in Eq. (40) loses accuracy), as visible in the large shoulders
on the side of the resonance peaks in Fig. 8. Incidentally, the perturbative
method predicts a large damping rate in this regime (blue dashed curve on the
bottom panel of Fig. 7). Overall, this behavior indicates that the
quasiparticle is no longer a well defined object when its energy approaches
$3\Delta$ at strong-coupling. This is in qualitative agreement with the
results of Ref. [48].

One possibility to account for the increasing complexity of the quasiparticle
Green’s function, and in particular to model the shoulders that the resonance
develops in this regime, is to use a fit function describing two interfering
resonances:

$\mathcal{G}({\bf{k}},\varepsilon+\mathrm{i}0^{+})\underset{\textrm{fit in
sector }X\textrm{ to
}}{\to}\frac{\Upsilon_{\bf{k}}^{(X,1)}}{z_{\bf{k}}^{(X,1)}-\varepsilon}+\frac{\Upsilon_{\bf{k}}^{(X,2)}}{z_{\bf{k}}^{(X,2)}-\varepsilon}.$
(41)

where $X=A,B$ and the second poles (not explicitly shown in Fig. 7) have a
lower spectral weight than the first ones
$|\Upsilon_{\bf{k}}^{(X,2)}|<|\Upsilon_{\bf{k}}^{(X,1)}|$ by convention. This
second pole in the energy sector $A$ ($\epsilon_{\mathrm{th}}^{1\rightarrow
2}<\varepsilon<\epsilon_{\mathrm{th}}^{1\rightarrow 3}$) can be seen as a
remnant of the first pole from sector $B$
($\epsilon_{\mathrm{th}}^{1\rightarrow 3}<\varepsilon$) and vice-versa. Note
that the main fitted solutions $z_{\bf{k}}^{(A,1)}$ (in green) and
$z_{\bf{k}}^{(B,1)}$ (in red) continue to formally exist outside their
respective energy sector, with however a decreasing fitting accuracy, and less
physical significance [52, 64]. For this reason, we display the
eigenfrequencies as dotted lines outside their respective energy sector on
Fig. 7. We note an important mismatch between $z_{\bf{k}}^{(A,1)}$ and
$z_{\bf{k}}^{(B,1)}$, in the real and especially the imaginary part, which
reflects the repulsion exerted on the resonance peak by the $3\Delta$
threshold.

This non trivial behavior near $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ is
missed by the perturbative solution which incorrectly predicts a smooth entry
into the fermionic $1\to 3$ continuum (although with a degraded quality factor
$E_{\bf{k}}/\Gamma_{\bf{k}}\approx 1$). Another discrepancy between the
perturbative and self-consistent energies is visible at low $k$, in the
decreasing part of the dispersion: there, the mean-field energy is at the
disintegration threshold
$\epsilon_{\bf{k}}=\epsilon_{\mathrm{th}}^{1\rightarrow 2}$, such that
perturbative approach predicts undamped quasiparticles. Conversely, the self-
consistent solution $z_{\bf{k}}^{(A)}$ is pushed up in energy and thus
acquires a (small) nonzero damping rate.

Figure 8: Cross section of the real (top panel) and imaginary (bottom panel)
part of the quasiparticle Green’s function at unitarity and $\hbar
k=2.16\sqrt{2m\Delta}$. In order to provide a good fit of the propagator, we
split up the energy domain in resonance domains delimited by the threshold
energies $\epsilon_{\mathrm{th}}^{1\rightarrow 2}$ and
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ shown as black dash-dotted lines.
The fitted function according to Eq. (41) is shown as a blue dashed line below
$\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (sector A) and as a green dotted
line above (sector B). The real parts of the eigenenergies found from the fit
are indicated on the axis with $a_{j}\equiv\mathrm{Re}\,z_{\bf{k}}^{(A,j)}$
and $b_{j}\equiv\mathrm{Re}\,z_{\bf{k}}^{(B,j)}$

###  III.4 Low-energy properties

The previous section has shown that the quasiparticles corresponding to the
exact zeros of Eq. (38) remain well-defined around their energy minimum
everywhere in the BCS-BEC crossover. It is thus useful to explicitly examine
the low-energy properties of the quasiparticle branch, which we do in this
last section.

In order to extract experimentally relevant properties from the quasiparticle
branch, we fit a quadratic dispersion around the minimum of the eigenenergy

$\epsilon_{k}^{\mathrm{fit}}=\epsilon^{\ast}+\frac{\hbar^{2}(k-k_{\mathrm{m}}^{\ast})^{2}}{2m^{\ast}},$
(42)

and study the effective energy gap $\epsilon^{\ast}$, location of the minimum
$k_{\mathrm{m}}^{\ast}$, and effective mass $m^{\ast}$ for different values of
the interaction. The results are shown in Figs. 9 and 10.

Figure 9: The quadratic shift of the energy minimum compared to the mean-field
value. Our self-consistent calculation of the shift (red circles) is compared
to the Hartree [56, 66, 67, 57] (black dotted line) and Galitskii [59, 68]
(blue dashed line) results in the BCS limit.

In Fig. 9 we show the (squared) shift of the dispersion minimum with respect
to the mean-field value $k_{\mathrm{m}}^{(0)}$. This shift is strictly
positive throughout the BEC-BCS crossover (as long as
$k_{\mathrm{m}}^{\ast}\neq 0$), such that the dispersion minimum is reached at
shorter wavelengths. In the BCS limit, our numerical results recover the
Hartree shift
$(k_{\mathrm{m}}^{\ast})^{2}-(k_{\mathrm{m}}^{(0)})^{2}=-4k_{\mathrm{F}}^{3}a/3\pi$,
as shown by the dotted asymptote. Our inconsistent $t$-matrix self-energy
however does not capture the next-to-leading order correction to the Hartree
shift obtained by Galitskii [59, 68] (dashed line in Fig. 9). The appearance
of the Hartree shift can be understood heuristically by retaining only the
leading order term in the expansion Eq. (23) of the pair propagator for $g\to
0$: $\Gamma_{ss^{\prime}}=-g\delta_{ss^{\prime}}+\mathcal{O}(g^{2})$. Then,
the off-diagonal elements of the self-energy become zero, while
$\Sigma_{++}(k)=g/\beta V\sum_{{\bf{k}},n}G^{(0)}_{++}(k)$. This sum can be
performed analytically, to find $\Sigma_{++}=-gn/2$, with $n$ the total
fermion density. This is translated in the fermionic energy by a shift of the
chemical potential $\mu\rightarrow\mu-
gn/2=\mu-4\epsilon_{\mathrm{F}}k_{\mathrm{F}}a/3\pi$, better known as the
Hartree shift [56, 66, 67, 57]. Note that although this Hartree shift tends to
0 in the BCS limit, it becomes much larger than the gap (which vanishes
exponentially with $k_{F}a$). Near the dispersion minimum, this forbids a
linearisation of Eq. (38) where $\Sigma$ would be treated as an infinitesimal
[48]. Although in general it is not possible to separate the contributions of
the different disintegration processes, the appearance of the Hartree shift in
the BCS limit can be understood to be mainly coming from
$\Sigma^{\mathrm{bc}}$, as it is not present when only including the $1\to 2$
process [48].

Figure 10: The quasiparticle gap (top panel) and effective mass (bottom panel)
in function of the interaction $1/k_{\mathrm{F}}a$ (bottom axis) or
$\mu/\Delta$ (top axis). The red circles show the self-consistent solutions,
which are compared to the perturbative results (blue crosses) and the mean-
field values (black dotted line). It is clear that the perturbative
calculation leads to wrong results in the BCS limit, as it is not possible to
treat the self-energy $\Sigma$ as a small correction.

Looking then at the energy gap in Fig. 10, we observe a fairly small
correction throughout the BCS-BEC crossover. We underline that this figure
should be interpreted by keeping in mind the difference between the energy gap
$\epsilon^{\ast}$ and the order parameter $\Delta$. While most theoretical
works (see e.g. Refs. [55, 56, 57]) focus on the order parameter, which is a
variational parameter determining the superfluid phase transition,
experimental results usually measure (twice) the energy gap (see e.g. Refs.
[41, 36]), the energy needed to break a condensed pair. While both are the
same on the mean-field level, this is no longer obviously true when
fluctuations are taken into account.

In the BCS limit, the corrected gap tends to the mean-field gap,
$\epsilon^{*}/\Delta\to 1$. This appears as a serious limitation of our
inconsistent $t$-matrix self-energy, since Gor’kov and Melik-Barkhudarov [55]
have predicted a lowering of the order parameter of a factor $\approx 2.2$,
while experimental results have reported a correction of the same magnitude on
the energy gap [38]. To account for this reduction of the gap, it is necessary
to add diagrams to the quasiparticle self-energy to include a similar
correction as in Ref. [55]. For a consistent theory, the equation for the
order parameter should then also be modified. In this case, it is still an
open question if the corrected gap and order parameter are equal, even in the
weak-coupling limit. In the BCS limit, we also note that the perturbative
result for $\epsilon^{\ast}$ differs much from the self-consistent one. This
can be seen as a consequence of the Hartree-shift: the energy displacement
$\text{Re}(z_{\bf{k}})-\epsilon_{\bf{k}}$ near the dispersion minimum is large
compared to $\Delta$ (as visible on Fig. 5), such that the replacement
$z_{\bf{k}}\to\epsilon_{\bf{k}}+\mathrm{i}0^{+}$ in Eq. (39) has a large
effect on the energy correction. At unitarity, we find the gap to be larger
than the mean-field prediction $\epsilon^{\ast}\simeq 1.06\Delta$. This
contrasts with our previous study [48] limited to the bosonic decay $1\to 2$,
and with experimental results observing a reduction of the gap [41, 38].
Possible explanations of this quantitative mismatch are the non-self-
consistent definition of the self-energy and the omission of the Gorkov-Melik
Barkhudarov contribution to the self-energy.

Finally, it is useful to examine the low-energy spectrum in detail in the BCS
limit, which can be seen in Fig. 11 for $\mu/\Delta=100$
($1/k_{\mathrm{F}}a\simeq-3$). There, as explained above, the self-energy
mainly produces a shift of the energy minimum, while the correction to the gap
remains small. We can align the computed eigenenergies by shifting the
wavenumber with respect to the location of the minimum
$k_{\mathrm{m}}^{\ast}$, to reveal that the structure of $z_{\bf{k}}$ remains
close to the BCS energy (this can also be seen in Fig. 10, as both the gap and
effective mass tend to the mean-field result in the BCS limit). Due to the
large shift, the self-consistent energy enters the continuum at
${\bf{k}}$-values where the threshold $\epsilon_{\mathrm{th}}^{1\rightarrow
2}$ is very close to $\epsilon_{\mathrm{th}}^{1\rightarrow 3}$, such that the
energy interval in sector A is small (see Fig. 4). Therefore, the damping rate
rapidly increases as both the $1\rightarrow 2$ and $1\rightarrow 3$ processes
are resonant. The perturbative approach, however, breaks down in the BCS limit
as it fails to capture the Hartree shift, as explained above. Most notably,
the perturbative approximation predicts a finite lifetime of the quasiparticle
excitations at the energy minimum, as the resonance condition is controlled by
the mean-field energy $\epsilon_{{\bf{k}}}$. For this reason, we have to shift
the perturbative damping rate by $k_{\mathrm{m}}^{(0)}$ instead of
$k_{\mathrm{m}}^{\ast}$ used for the eigenfrequency to be able to compare with
the self-consistent solution in Fig. 11. Doing this, we see that the damping
rate remains small when the BCS energy enters sector A and starts to grow
rapidly when it reaches the second threshold energy
$\epsilon_{{\bf{k}}}>\epsilon_{\mathrm{th}}^{1\rightarrow 3}$ (sector B).

Figure 11: Quasiparticle spectrum in the BCS regime ($\mu/\Delta=100$,
$1/k_{\mathrm{F}}a\simeq-3$). The eigenenergies and spectral density are
shifted to facilitate the comparison between the different curves. The self-
consistent solution is shown in red, and both its real (top panel) and
imaginary part (bottom panel) are shifted with respect to
$k_{\mathrm{m}}^{\ast}$ computed from $z_{\bf{k}}$. The perturbative result
(blue dashed line) incorrectly predicts a finite quasiparticle lifetime at the
dispersion minimum, but to be able to compare with the self-consistent result,
the real part $E_{\bf{k}}$ is shifted with respect to $k_{\mathrm{m}}^{\ast}$,
while the damping rate $\Gamma_{\bf{k}}$ is shifted to $k_{\mathrm{m}}^{(0)}$.
The black dotted line is the mean-field energy, and the vertical dash-dotted
lines depict the shifted $k$-values where the mean-field energy
$\epsilon_{{\bf{k}}}$ and self-consistent energy $z_{\bf{k}}$ cross the
threshold energies respectively in black and red.

##  IV Conclusion

We have performed an in-depth study of the quasiparticle spectrum in a
superfluid Fermi gas and identified the disintegration processes described by
the $t$-matrix self-energy. We find evidence of multiple excitation branches
at strong coupling close to the threshold energy where the four-fermion
processes of disintegration into 3 quasiparticles become resonant. Using
momentum-resolved rf-spectroscopy [42], the quasiparticle spectrum can
directly be probed, such that our results could be experimentally verified. At
weak-coupling, our approach captures the Hartree shift of the dispersion
minimum, but not the Gor’kov-Melik Barkhudarov correction to the gap.
Generalizing the correction of Ref. [55] to the quasiparticle spectrum thus
appears as a necessary continuation of the present study.

Although the calculation was done within the framework of ultracold fermionic
gases, similar decay processes occur in other quantum many-body systems where
both phononic and rotonic excitations are present [5, 6, 10, 11, 12, 13, 14,
15]. In superconductors, besides the intrinsic processes considered here, the
quasiparticle lifetime is limited also by extraneous processes such as
impurities scattering [26] and emission of lattice phonons [10]. However, the
$3\Delta$ threshold should still be visible, and the $1\to 3$ disintegration
process should be the main intrinsic disintegration process of the electron
gas, owing to the gapped nature of the collective plasma branch.

##  Appendix A Integration over the internal structure of the $1\to 3$
disintegration continuum

In this appendix we study the internal structure of the fermionic
disintegration process $1\to 3$ (Fig. 3c), in relation to the numerical
evaluation of $\Sigma^{\rm bc}$ in Eq. (20). Rather than using the symmetric
formulation as in Eq. (22), the form of the integral incites us to pair up two
of the emitted quasiparticles. The process is then resonant for a given value
of the fermion momentum $k$ and energy $z_{k}$ if one can find ${\bf{q}}$ and
$\omega_{q}$ such that:

$z_{k}-\omega_{q}=\epsilon_{{\bf{k}}-{\bf{q}}},$ (43)

where $\omega_{q}$ has to be inside the pair-breaking continuum
$[\epsilon_{c}(q),+\infty[$. (Everywhere in this appendix, energies and
momenta are in units of $\Delta$ and $k_{\Delta}=\sqrt{2m\Delta}/\hbar$, and
$k_{0}=\sqrt{2m\mu}/\hbar$ denotes the mean-field dispersion minimum). Whether
this resonance condition is met or not has a huge impact on the integrand of
Eq. (20), in particular on its imaginary part. This motivates a precise
identification of the resonance domains.

Numerically, we integrate successively over the angle
$u={\bf{k}}\cdot{\bf{q}}/kq$ (the integral is done analytically as explained
in Appendix B), then over $\omega_{q}$ and finally over $q$. As we integrate
over $u$, the energy $\epsilon_{{\bf{k}}-{\bf{q}}}$ of the unpaired emitted
quasiparticle goes from $\epsilon_{k-q}$ to $\epsilon_{k+q}$ (black solid
curves in Fig. 12), either monotonically (for $q<|k-k_{0}|$ or $q>k+k_{0}$) or
non-monotonically (for $|k-k_{0}|<q<k+k_{0}$) if it passes through the energy
minimum (horizontal dotted line in Fig. 12). Depending on the value of
$z_{k}-\omega_{q}$, there can then be either 0, 1 or 2 resonance angles. We
then integrate over $\omega_{q}$, which corresponds to a descending vertical
path on Fig. 12, as shown by the red arrow. Along the integration, the number
of resonance angles changes if $\omega_{q}$ passes through
$z_{k}-\epsilon_{k-q}$ or $z_{k}-\epsilon_{k+q}$ and if $\omega_{q}$ passes
$z_{k}-1$ while in the blue (2R) region. To reach a good precision on the
integral, we split it at those values. Depending on the values of $k,z_{k}$
and the last integration variable $q$, there are five possible configurations
in which the integral over $\omega_{q}$ may be split; they are denoted by
greek letters ($\alpha$, $\beta$, $\gamma$, $\delta$ and $\epsilon$),
described in Tab. 1, and shown in colors on the right panel of Fig. 12 in
function of $q$ and $z_{k}-2$. In the example shown by the red arrow of Fig.
12, there are 0 resonance angles from $\omega_{q}=\epsilon_{c}(q)$ to
$z_{k}-\epsilon_{k+q}$, 1 from $z_{k}-\epsilon_{k+q}$ to
$z_{k}-\epsilon_{k-q}$, 2 from $z_{k}-\epsilon_{k-q}$ to $z_{k}-1$, and again
0 from $z_{k}-1$ to $+\infty$, which corresponds to the configuration
$\epsilon$ in Tab. 1.

$k_{4}\leq k\leq k_{3}$, type IV ($k_{0}=2$, $k=0.92)$

Figure 12: Resonance diagrams and integration over $\omega_{q}$ and $q$. The
left panel shows the number of resonance angles in function of $q$ and
$\omega_{q}$ (the $y$-axis is $z_{k}-\omega_{q}$): either 0 (white region), 1
(yellow region) or 2 (blue region) resonance angles, depending on whether the
function $u\mapsto\epsilon_{{\bf{k}}-{\bf{q}}}$ is monotonous or not, and
whether $z_{k}-\omega_{q}$ lies within its energy-range or not (black solid
line: value $\epsilon_{k\mp q}$ of the function in $u=\pm 1$). As we then
integrate over $\omega_{q}$ in the pair-breaking continuum ($z_{k}-\omega_{q}$
varies from $z_{k}-\epsilon_{c}(q)$ to $-\infty$, see the red vertical arrow
on the left panel), the number of resonance angles changes according to a
configuration $\alpha$, $\beta$, $\gamma$, $\delta$ or $\epsilon$ (see Tab. 1)
shown in color on the right panel. As we finally integrate over $q$ (red solid
curve on the right panel), the succession of greek-letter configurations gives
rise to one of the superconfiguration of Tabs. 2 and 3. This
superconfiguration itself changes when $z_{k}$ passes one of the energy lines
shown as thin dashed horizontal curves. This figure is drawn for
$k_{4}<k<k_{3}$ (vertical sector number IV in Fig. 13) and $k_{0}=2$.
Configurations of the $\omega_{q}$ integral
---
Name | Energy boundaries | 0,1 or 2 resonances
$0$ | $2$ |  |  |  | 0R |  |  |  |
$\alpha$ | $2$ | $z_{k}-\epsilon_{k-q}$ |  |  | 1R | 0R |  |  |
$\beta$ | $2$ | $z_{k}-\epsilon_{k-q}$ | $z_{k}-1$ |  | 1R | 2R | 0R |  |
$\gamma$ | $2$ | $z_{k}-1$ |  |  | 2R | 0R |  |  |
$\delta$ | $2$ | $z_{k}-\epsilon_{k+q}$ | $z_{k}-\epsilon_{k-q}$ |  | 0R | 1R | 0R |  |
$\epsilon$ | $2$ | $z_{k}-\epsilon_{k+q}$ | $z_{k}-\epsilon_{k-q}$ | $z_{k}-1$ | 0R | 1R | 2R | 0R |
Table 1:  Resonance configurations of the integral over $\omega_{q}$ (at fixed
$k$, $z_{k}$ and $q$ but after integration over $u$). The table reads as
follows: the configuration $\beta$ has 1 resonance angle (1R) for
$2<\omega_{q}<z_{k}-\epsilon_{k-q}$, 2 angles (2R) for
$z_{k}-\epsilon_{k-q}<\omega_{q}<z_{k}-1$ and 0 angle (0R) for
$\omega_{q}>z_{k}-1$. The configurations symbolized by the prime letters
$\alpha^{\prime}$,…$\epsilon^{\prime}$ are deduced from $\alpha$,…$\epsilon$
by changing the value of the continuum threshold:
$2\to\epsilon_{c}(q>2k_{0})$. The tilde configuration $\tilde{\alpha}$,
$\tilde{\beta}$, $\tilde{\delta}$ and $\tilde{\epsilon}$ are obtained by
swapping $\epsilon_{k-q}$ and $\epsilon_{k+q}$ in the energy boundaries of
$\alpha$, $\beta$, $\delta$, $\epsilon$.

The succession of the $\alpha$, $\beta$, $\gamma$, $\delta$ or $\epsilon$
configurations as we finally integrate from $q=0$ to $+\infty$ gives rise to
one of the superconfigurations denoted by capital letters from A to J and
described in Tabs. 2 and 3 (for $k<k_{0}$ and $k>k_{0}$ respectively). The
boundary values of $q$ where the (greek-letter) configuration changes are then
given in Tabs. 4 and 5 (again for $k<k_{0}$ and $k>k_{0}$ respectively). Note
that in all the zoology of superconfigurations, the set of resonant wave
numbers $q$ is always either the empty set or a (connected) interval (which
can easily be read on Tabs. 2 and 3 as the interval between the lower bound of
the first configuration and the upper bound of the last one, excluding the 0
configuration). Numerically, the splitting of the integral over $q$ at those
boundaries is thus less crucial than for the integral over $\omega_{q}$,
except at the boundaries where the resonance totally disappears.

Superconfigurations of the $q$-integral for $k<k_{0}$
---
Name | Boundaries | Configuration
A | $0$ | $q_{1}$ | $k_{-}$ | $q_{2}$ |  | $q_{2}^{\prime}$ | $k_{+}$ | $q_{3}^{\prime}$ | $0$ | $\tilde{\alpha}$ |  | $\tilde{\beta}$ | $\gamma$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
B | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $q_{2}$ |  | $q_{2}^{\prime}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\alpha}$ |  | $\tilde{\beta}$ | $\gamma$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{B}^{\prime}$ | $0$ | $q_{1}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{c}$ | $q_{2}$ | $k_{+}$ | $q_{3}^{\prime}$ | $0$ | $\tilde{\alpha}$ | $\tilde{\beta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
C | $0$ | $k_{-}$ | $q_{1}^{\prime}$ | $q_{2}$ |  | $q_{2}^{\prime}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ |  | $\tilde{\beta}$ | $\gamma$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{C}^{\prime}$ | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{c}$ | $q_{2}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\alpha}$ | $\tilde{\beta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{C}^{\prime\prime}$ | $0$ | $q_{1}$ | $q_{2}^{\prime}$ | $k_{-}$ | $q_{c}$ | $q_{2}$ | $k_{+}$ | $q_{3}^{\prime}$ | $0$ | $\tilde{\alpha}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
D | $0$ | $k_{-}$ | $q_{c}$ | $q_{2}$ | $q_{1}^{\prime}$ | $q_{2}^{\prime}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\gamma$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{D}^{\prime}$ | $0$ | $k_{-}$ | $q_{1}^{\prime}$ | $q_{2}^{\prime}$ | $q_{c}$ | $q_{2}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\tilde{\beta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{D}^{\prime\prime}$ | $0$ | $q_{1}^{\prime}$ | $q_{2}^{\prime}$ | $k_{-}$ | $q_{c}$ | $q_{2}$ | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\alpha}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{D}^{\prime\prime\prime}$ | $0$ | $q_{1}$ | $q_{2}^{\prime}$ | $k_{-}$ | $q_{c}$ | $k_{+}$ | $q_{2}$ | $q_{3}^{\prime}$ | $0$ | $\tilde{\alpha}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\delta$ | $\alpha\text{ or }\alpha^{\prime}$
E | $0$ | $k_{-}$ | $q_{c}$ | $q_{2}$ |  |  | $k_{+}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ |  |  | $\beta$ | $\alpha\text{ or }\alpha^{\prime}$
$\text{E}^{\prime}$ | $0$ | $q_{1}^{\prime}$ | $q_{2}^{\prime}$ | $k_{-}$ | $q_{c}$ | $k_{+}$ | $q_{2}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\alpha}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ | $\delta$ | $\alpha\text{ or }\alpha^{\prime}$
F | $0$ | $k_{-}$ | $q_{c}$ | $k_{+}$ |  |  | $q_{2}$ | $q_{3}^{\prime}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\epsilon$ |  |  | $\delta\text{ or }\delta^{\prime}$ | $\alpha\text{ or }\alpha^{\prime}$
Table 2:  For $k<k_{0}$, table of the resonance superconfigurations for the
integration over $q$ (at fixed $k$ and $z_{k}$ and after integration over $u$
and $\omega_{q}$). The first column is the name of the configuration in Roman
letters, the second column is the list of $p+1$ boundaries values of $q$ (see
Table 4 for their meaning) splitting the integration interval $[0,+\infty[$
into $p$ domains, and the third column is the configuration itself: a list of
$p$ greek letters (plus the final configuration $0$ not written explicitly),
each one denoting the shape of the resonance in the subintegral over
$\omega_{q}$ (see Table 1). Superconfigurations of the $q$-integral for
$k>k_{0}$
---
Name | Boundaries | Configuration
G | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{\rm m}$ |  |  | $0$ | $\alpha$ | $\beta$ | $\gamma\text{ or }\gamma^{\prime}$ |  |
J | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{3}^{\prime}$ | $q_{3}^{\prime\prime}$ | $q_{\rm m}$ | $0$ | $\alpha$ | $\beta$ | $\gamma\text{ or }\gamma^{\prime}$ | $\beta\text{ or }\beta^{\prime}$ | $\gamma^{\prime}$
A | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{3}^{\prime}$ | $k_{+}$ | $q_{4}^{\prime}$ | $0$ | $\alpha$ | $\beta$ | $\gamma\text{ or }\gamma^{\prime}$ | $\beta\text{ or }\beta^{\prime}$ | $\alpha^{\prime}$
B | $0$ | $q_{1}$ | $k_{-}$ | $q_{2}^{\prime}$ | $q_{3}^{\prime}$ | $k_{+}$ | $q_{4}^{\prime}$ | $\delta$ | $\alpha$ | $\beta$ | $\gamma\text{ or }\gamma^{\prime}$ | $\beta\text{ or }\beta^{\prime}$ | $\alpha^{\prime}$
D | $0$ | $k_{-}$ | $q_{1}$ | $q_{2}^{\prime}$ | $q_{3}^{\prime}$ | $k_{+}$ | $q_{4}^{\prime}$ | $\delta$ | $\epsilon\text{ or }\epsilon^{\prime}$ | $\beta$ | $\gamma\text{ or }\gamma^{\prime}$ | $\beta\text{ or }\beta^{\prime}$ | $\alpha^{\prime}$
E | $0$ | $k_{-}$ | $q_{1}$ | $k_{+}$ |  |  | $q_{4}^{\prime}$ | $\delta$ | $\epsilon\text{ or }\epsilon^{\prime}$ |  |  | $\beta\text{ or }\beta^{\prime}$ | $\alpha^{\prime}$
F | $0$ | $k_{-}$ | $k_{+}$ | $q_{1}$ |  |  | $q_{4}^{\prime}$ | $\delta$ | $\epsilon\text{ or }\epsilon^{\prime}$ |  |  | $\delta^{\prime}$ | $\alpha^{\prime}$
$\text{B}^{\prime}$ | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $k_{+}$ | $q_{4}^{\prime}$ |  |  | $0$ | $\alpha$ | $\beta\text{ or }\beta^{\prime}$ |  |  | $\alpha^{\prime}$
$\text{D}^{\prime}$ | $0$ | $q_{1}^{\prime}$ | $k_{-}$ | $k_{+}$ | $q_{4}^{\prime}$ |  |  | $\delta$ | $\alpha$ | $\beta\text{ or }\beta^{\prime}$ |  |  | $\alpha^{\prime}$
H | $0$ | $q_{1}^{\prime}$ | $q_{1}^{\prime\prime}$ |  |  |  |  | $0$ | $\alpha\text{ or }\alpha^{\prime}$ |  |  |  |
Table 3:  Same table as Table 2 but for $k>k_{0}$. See Table 5 for the meaning of the boundary values of $q$. Boundary momentum | Solution of | in the interval | condition of existence
---|---|---|---
$q_{1}$ | $\tilde{z}_{k}=\epsilon_{k+q}$ | $[0,k_{0}-k]$ | $\epsilon_{k}>\tilde{z}_{k}>1$
$q_{2}$ | $\tilde{z}_{k}=\epsilon_{k+q}$ | $[k_{0}-k,+\infty[$ | $\tilde{z}_{k}>1$
$q_{1}^{\prime}$ | $\tilde{z}_{k}=\epsilon_{k-q}$ | $[0,k]$ | $\epsilon_{0}>\tilde{z}_{k}>\epsilon_{k}$
$q_{2}^{\prime}$ | $\tilde{z}_{k}=\epsilon_{k-q}$ | $[k,k+k_{0}]$ | $\epsilon_{0}>\tilde{z}_{k}>1$
$q_{3}^{\prime}$ | $\tilde{z}_{k}=\epsilon_{k-q}$ | $[k+k_{0},+\infty[$ | $\tilde{z}_{k}>1$
$q_{c}$ | $\epsilon_{k+q}=\epsilon_{k-q}$ | $[k-k_{0},k+k_{0}]$ | $k<k_{0}$
Table 4:  Table of the boundary momenta (used in Table 2) in the region $k<k_{0}$. We use here the notation $\tilde{z}_{k}=z_{k}-2$. Boundary momentum | Solution of | in the interval | condition of existence
---|---|---|---
For all $k$
$q_{1}$ | $z_{k}=\epsilon_{k+q}+\epsilon_{c}(q)$ | $[0,+\infty[$ | $\tilde{z}_{k}>\epsilon_{k}$
$q_{\rm m}$ | $z_{k}=1+\epsilon_{c}(q)$ | $[2k_{0},+\infty[$ | $\tilde{z}_{k}>1$
$q_{4}^{\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $[k+k_{0},+\infty[$ | $\tilde{z}_{k}>\omega_{3}$
specific to $2k_{0}>k>k_{0}$
$q_{1}^{\prime}$ | $z_{k}=\epsilon_{k-q}+2$ | $[0,k-k_{0}]$ | $\epsilon_{k}>\tilde{z}_{k}>1$
$q_{2}^{\prime}$ | $z_{k}=\epsilon_{k-q}+2$ | $[k-k_{0},k]$ | $\epsilon_{0}>\tilde{z}_{k}>1$
$q_{3}^{\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $[k,k+k_{0}]$ | $\epsilon_{0}>\tilde{z}_{k}>\omega_{2}$
$q_{3}^{\prime\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $]q_{3}^{\prime},k+k_{0}]$ | $\omega_{3}>\tilde{z}_{k}>\omega_{2}$
specific to $3k_{0}>k>2k_{0}$
$q_{1}^{\prime}$ | $z_{k}=\epsilon_{k-q}+2$ | $[0,k-k_{0}]$ | $\epsilon_{k}>\tilde{z}_{k}>1$
$q_{2}^{\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $[k-k_{0},k+k_{0}]$ | $\omega_{3}>\tilde{z}_{k}>1$
specific to $k>3k_{0}$
$q_{1}^{\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $[0,k-k_{0}]$ | $\epsilon_{k}>\tilde{z}_{k}>\omega_{2}^{\prime}$
$q_{1}^{\prime\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $]q_{1}^{\prime},k-k_{0}]$ | $\omega_{3}^{\prime}>\tilde{z}_{k}>\omega_{2}^{\prime}$
$q_{2}^{\prime}$ | $z_{k}=\epsilon_{k-q}+\epsilon_{c}(q)$ | $[k-k_{0},k+k_{0}]$ | $\omega_{3}>\tilde{z}_{k}>\omega_{3}^{\prime}$
Table 5:  Table of the boundary momenta (used in Table 3) in the region
$k>k_{0}$.

Fig. 13 shows the superconfigurations in colors in function of $k$ and
$z_{k}$. The energy lines which separate them in the $k,z_{k}$ plane (among
which the lower edge of the continuum $\epsilon_{\mathrm{th}}^{1\rightarrow
3}$) are given in Tab. 6.

Figure 13: Diagram of the resonance superconfigurations for $k\leq k_{0}$
(left panel) and $k\geq k_{0}$ (right panel). The configurations denoted by
Roman letters (see Tab. 2 for their signification) are shown in colors and
delimited by energy lines shown in solid curves: in black $\epsilon_{k}$,
$\epsilon_{2k-k_{0}}$ and $\epsilon_{2k+k_{0}}$, in blue
$\epsilon_{0}=\sqrt{k_{0}^{2}+1}$, in red the energies $\omega_{i}$ listed in
Tab. 6 (note that $\omega_{4}$ and $\omega_{4}^{\prime}$ continue
$\epsilon_{2k+k_{0}}$ and $\epsilon_{2k-k_{0}}$ respectively at $k>k_{0}$ and
$k>k_{12}$). The vertical lines delimit the $k$-intervals where the same
succession of configurations is found when increasing $z_{k}$ from 3 to
$+\infty$. To each sector between two vertical lines corresponds a diagram
similar to Fig. 12.

The only energy line posing a slight difficulty is $\omega_{2}^{\prime}$
(which separates superconfigurations J and G, see Fig. 13). At this energy,
the displaced pair-breaking continuum is tangent to the function
$q\mapsto\epsilon_{k-q}$, such that

$\displaystyle\omega_{2}^{\prime}-2\epsilon_{q_{\rm tg}/2}$ $\displaystyle=$
$\displaystyle\epsilon_{k-q_{\rm tg}},$ (44) $\displaystyle v_{q_{\rm tg}/2}$
$\displaystyle=$ $\displaystyle v_{k-q_{\rm tg}}\qquad\text{with
}v_{k}=\partial\epsilon_{k}/\partial k$ (45)

and where $q_{\rm tg}$ belongs to the interval $[k_{0},2k_{0}]$. This leads to
the polynomial equation on $q_{\rm tg}$:

$16q_{\rm tg}^{8}-32kq_{\rm tg}^{7}+8\left(3k^{2}-5k_{0}^{2}\right)q_{\rm
tg}^{6}-8k\left(k^{2}-9k_{0}^{2}\right)q_{\rm
tg}^{5}+\left(k^{4}-50k^{2}k_{0}^{2}+33k_{0}^{4}+21\right)q_{\rm tg}^{4}\\\
+4k\left(4k^{2}k_{0}^{2}-12k_{0}^{4}-9\right)q_{\rm
tg}^{3}+\left[-2k^{4}k_{0}^{2}+k^{2}\left(28k_{0}^{4}+25\right)-10\left(k_{0}^{6}+k_{0}^{2}\right)\right]q_{\rm
tg}^{2}\\\ -8k\left(k_{0}^{4}+1\right)(k-k_{0})(k+k_{0})q_{\rm
tg}+\left(k_{0}^{4}+1\right)\left(k^{2}-k_{0}^{2}\right)^{2}=0.$ (46)

Energy line | Expression
---|---
$\omega_{1}$ | $\text{ solution of }\epsilon_{k+q}=\epsilon_{k-q}\text{ for }q\in[0,k-k_{0}]$
$\omega_{2}$ | ${\epsilon_{k-q_{\rm tg}+2\epsilon_{q_{\rm tg}/2}}}\text{ with $q_{\rm tg}$ solution of }P_{8}(q_{\rm tg})=0\text{ in }[k_{0},2k_{0}]\text{(see Eq.~{}\eqref{P8q})}$
$\omega_{2}^{\prime}=\epsilon_{\mathrm{th}}^{1\rightarrow 3}-2$ | $3\epsilon_{k/3}-2\text{ (lower edge of the continuum)}$
$\omega_{3}$ | $\epsilon_{c}(k+k_{0})-1$
$\omega_{3}^{\prime}$ | $\epsilon_{c}(k-k_{0})-1$
$\omega_{4}$ | $\epsilon_{2k+k_{0}}+\epsilon_{c}(k+k_{0})-2$
$\omega_{4}^{\prime}$ | $\epsilon_{2k-k_{0}}+\epsilon_{c}(k-k_{0})-2$
$\omega_{4}$ | $\epsilon_{2k+k_{0}}+\epsilon_{c}(k+k_{0})-2$
Table 6:  Energy lines separating the resonance superconfigurations in Fig.
13.

Finally the remarkable values of $k$ at which the energy lines cross or
anticross (the vertical lines on Fig. 13) are gathered in Tab. A. To each
interval between successive $k_{i}$ values corresponds a type of diagram
(denoted by Roman numbers) similar to Fig. 12.

Name | $\text{Value in function of }k_{0}$ | Energy lines crossing
---|---|---
$k_{0}^{-}$ |  | $\omega_{1}=\epsilon_{k}=\epsilon_{2k-k_{0}}=1$
$k_{1}$ | $k_{0}/\sqrt{2}$ | $\omega_{1}=\epsilon_{0}\text{ (avoided crossing, }\omega_{1}\leq\epsilon_{0}\text{)}$
$k_{2}$ | $3k_{0}/5$ | $\epsilon_{2k-k_{0}}=\omega_{1}$
$k_{3}$ | $k_{0}/2$ | $\epsilon_{2k-k_{0}}=\epsilon_{0}\text{ (avoided crossing, }\epsilon_{2k-k_{0}}\leq\epsilon_{0}\text{)}$
$k_{4}$ | $k_{0}/\sqrt{5}$ | $\epsilon_{k}=\omega_{1}$
$k_{5}$ | $k_{0}/3$ | $\epsilon_{2k-k_{0}}=\epsilon_{k}$
$k_{6}$ | $(\sqrt{2}-1)k_{0}/2$ | $\epsilon_{2k+k_{0}}=\epsilon_{0}$
$k_{7}$ | $k_{0}/5$ | $\epsilon_{2k+k_{0}}=\epsilon_{k}$
$k_{0}^{+}$ |  | $\omega_{2}=\omega_{3}=\epsilon_{k}=\epsilon_{2k+k_{0}}=1$
$k_{8}$ | $(\sqrt{2}+1)k_{0}/2$ | $\epsilon_{2k+k_{0}}=\epsilon_{0}$
$k_{9}$ | $\sqrt{2}k_{0}$ | $\epsilon_{k}=\epsilon_{0}$
$k_{10}$ | $\sqrt{4k_{0}^{2}+2\sqrt{k_{0}^{4}+2\sqrt{1+k_{0}^{4}}-2}}-k_{0}$ | $\omega_{3}=\epsilon_{0}$
$k_{11}$ | $2k_{0}$ | $\epsilon_{0}=\omega_{2}\text{($\omega_{2}$ disappear beyond this point)}$
$k_{12}$ | $3k_{0}$ |
Table 7:  Values of $k$ where energy lines cross.

##  Appendix B Angular integration in the self-energy

In this appendix, we derive analytic formulas of the angular integrals
appearing in the self-eigenenergies $\Sigma^{\rm p}$ and $\Sigma^{\rm bc}$. In
Eq. (18) the pair propagator is independent of the direction of ${\bf{q}}$ and
goes out of the integral over $u={\bf{k}}\cdot{\bf{q}}/kq$. The same is true
for the spectral density of the pair propagator in Eq. (20). This leaves only
elementary functions to integrate over $u$. We thus introduce the integrals

$\displaystyle I_{U^{2}|V^{2}}(z)$ $\displaystyle=$
$\displaystyle\int_{-1}^{1}\frac{\mathrm{d}u}{2}\left(1\pm\frac{\xi_{{\bf{k}}-{\bf{q}}}}{\epsilon_{{\bf{k}}-{\bf{q}}}}\right)\frac{\Delta}{z+\epsilon_{{\bf{k}}-{\bf{q}}}}$
(47) $\displaystyle I_{UV}(z)$ $\displaystyle=$
$\displaystyle\int_{-1}^{1}\frac{\mathrm{d}u}{2}{\frac{\Delta}{\epsilon_{{\bf{k}}-{\bf{q}}}}}\frac{\Delta}{z+\epsilon_{{\bf{k}}-{\bf{q}}}},$
(48)

and write the matrix elements of the self-energy in the form

$\displaystyle\Sigma_{++}^{\mathrm{p}}({\bf{k}},z_{k})\\!$
$\displaystyle=\\!-\frac{1}{4\pi^{2}}\int\limits_{0}^{+\infty}\frac{q^{2}dq}{\partial_{\omega}\det\Gamma^{-1}(q,\omega_{{\bf{q}}})}\left[I_{V^{2}}(\omega_{{\bf{q}}}-z_{k})\check{N}_{{--}}\left({q,\omega_{{\bf{q}}}}\right)-I_{U^{2}}(\omega_{{\bf{q}}}+z_{k})\check{N}_{{++}}\left({q,\omega_{{\bf{q}}}}\right)\right]$
(49) $\displaystyle\Sigma_{+-}^{\mathrm{p}}({\bf{k}},z_{k})\\!$
$\displaystyle=\\!\frac{1}{4\pi^{2}}\int\limits_{0}^{+\infty}\frac{q^{2}dq}{\partial_{\omega}\det\Gamma^{-1}(q,\omega_{{\bf{q}}})}\left[I_{UV}(\omega_{{\bf{q}}}-z_{k})+I_{UV}(\omega_{{\bf{q}}}+z_{k})\right]N_{{+-}}\left({q,\omega_{{\bf{q}}}}\right)$
(50) $\displaystyle\Sigma_{++}^{\mathrm{bc}}({\bf{k}},z_{k})\\!$
$\displaystyle=\\!-\frac{1}{4\pi^{2}}\int\limits_{0}^{+\infty}q^{2}dq\int\limits_{\epsilon_{c}(q)}^{+\infty}d\omega_{q}\left[I_{V^{2}}(\omega_{q}-z_{k})\rho_{++}(q,\omega_{q})-I_{U^{2}}(\omega_{q}+z_{k})\rho_{--}(q,\omega_{q})\right]$
(51) $\displaystyle\Sigma_{+-}^{\mathrm{bc}}({\bf{k}},z_{k})\\!$
$\displaystyle=\\!-\frac{1}{4\pi^{2}}\int\limits_{0}^{+\infty}q^{2}dq\int\limits_{\epsilon_{c}(q)}^{+\infty}d\omega_{q}\left[I_{UV}(\omega_{q}-z_{k})+I_{UV}(\omega_{q}+z_{k})\right]\rho_{+-}(q,\omega_{q}).$
(52)

As before, we work here in units of $\Delta$, setting
$\hbar^{2}k^{2}/2m\Delta\to k^{2}$, $\hbar^{2}q^{2}/2m\Delta\to q^{2}$ and

$\frac{{\xi}_{{\bf{k}}-{\bf{q}}}}{\Delta}=\xi_{0}-2kqu\rightarrow{\xi}_{{\bf{k}}-{\bf{q}}}\quad\text{with}\quad\xi_{0}=k^{2}+q^{2}-\mu.$
(53)

We use a Euler substitution to rationalize the integrand, setting

${\epsilon}_{{\bf{k}}-{\bf{q}}}-{\xi}_{{\bf{k}}-{\bf{q}}}=x\quad\Longleftrightarrow\quad
u=\frac{x^{2}+2\xi_{0}x-1}{4kqx}$ (54)
$\mathrm{d}u=\frac{x^{2}+1}{4kqx^{2}}\mathrm{d}x.$ (55)

Note in passing that this change of variable is monotonous (and increasing).
In the new variable $x$, the integrals become

$\displaystyle I_{U^{2}}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{2kq}\int_{x_{\rm min}}^{x_{\rm
max}}\frac{\mathrm{d}x}{xP(x)}$ (56) $\displaystyle I_{V^{2}}(z)$
$\displaystyle=$ $\displaystyle\frac{1}{2kq}\int_{x_{\rm min}}^{x_{\rm
max}}\mathrm{d}x\frac{x}{P(x)}$ (57) $\displaystyle I_{UV}(z)$
$\displaystyle=$ $\displaystyle\frac{1}{2kq}\int_{x_{\rm min}}^{x_{\rm
max}}\frac{\mathrm{d}x}{P(x)}$ (58)

with the new integration boundaries

$\displaystyle x_{\rm min}$ $\displaystyle=$
$\displaystyle\epsilon_{k+q}-\xi_{k+q}$ (59) $\displaystyle x_{\rm max}$
$\displaystyle=$ $\displaystyle\epsilon_{k-q}-\xi_{k-q}$ (60)

and the polynomial

$P(x)=x^{2}+2zx+1.$ (61)

The roots of $P$ are

$\displaystyle x_{1}$ $\displaystyle=$
$\displaystyle-z+\sqrt{z^{2}-1}\underset{z=\omega+\mathrm{i}0^{+}}{=}-\omega+\sqrt{\omega^{2}-1}+\mathrm{i}\text{sg}(\omega)0^{+}$
(62) $\displaystyle x_{2}$ $\displaystyle=$
$\displaystyle-z-\sqrt{z^{2}-1}\underset{z=\omega+\mathrm{i}0^{+}}{=}-\omega-\sqrt{\omega^{2}-1}-\mathrm{i}\text{sg}(\omega)0^{+}.$
(63)

Effecting the partial fraction decomposition of the integrals, we finally get

$\displaystyle I_{U^{2}}(\omega)$ $\displaystyle=$
$\displaystyle\frac{-I_{1}-\omega I_{2}+2I_{3}}{4kq}$ (64) $\displaystyle
I_{V^{2}}(\omega)$ $\displaystyle=$ $\displaystyle\frac{I_{1}-\omega
I_{2}}{4kq}$ (65) $\displaystyle I_{UV}(\omega)$ $\displaystyle=$
$\displaystyle\frac{I_{2}}{4kq},$ (66)

with

$I_{1}=\int_{x_{\rm min}}^{x_{\rm
max}}\mathrm{d}x\left(\frac{1}{x-x_{1}}+\frac{1}{x-x_{2}}\right)=\text{ln}\left|\frac{(x_{\rm
max}-x_{1})(x_{\rm max}-x_{2})}{(x_{\rm min}-x_{1})(x_{\rm
min}-x_{2})}\right|\\\ +\mathrm{i}\pi\text{sg}(\omega)\left(\Theta(x_{\rm
max}-x_{1})\Theta(x_{1}-x_{\rm min})-\Theta(x_{\rm
max}-x_{2})\Theta(x_{2}-x_{\rm min})\right)$ (67) $I_{2}=\int_{x_{\rm
min}}^{x_{\rm
max}}\frac{\mathrm{d}x}{\sqrt{\omega^{2}-1}}\left(\frac{1}{x-x_{1}}-\frac{1}{x-x_{2}}\right)\\\
={\small\begin{cases}\frac{1}{\mathrm{i}\sqrt{1-\omega^{2}}}\text{ln}{\frac{(x_{\rm
max}-x_{1})(x_{\rm min}-x_{2})}{(x_{\rm max}-x_{2})(x_{\rm
min}-x_{1})}}=\frac{2}{\sqrt{1-\omega^{2}}}\left(\text{arg}(x_{\rm
max}-x_{1})-\text{arg}(x_{\rm min}-x_{1})\right)\quad\text{if}\quad
1>|\omega|\\\
\frac{1}{\sqrt{\omega^{2}-1}}\Big{[}\text{ln}\left|{\frac{(x_{\rm
max}-x_{1})(x_{\rm min}-x_{2})}{(x_{\rm max}-x_{2})(x_{\rm
min}-x_{1})}}\right|\\\ +\mathrm{i}\pi\text{sg}(\omega)\left(\Theta(x_{\rm
max}-x_{1})\Theta(x_{1}-x_{\rm min})+\Theta(x_{\rm
max}-x_{2})\Theta(x_{2}-x_{\rm min})\right)\Big{]}\text{if}\quad
1<|\omega|\end{cases}}$ (68) $I_{3}=\text{ln}\frac{x_{\rm max}}{x_{\rm min}}.$
(69)

##  Appendix C UV contributions to $\Sigma^{\rm bc}$

In this appendix, we compute analytically the contributions of the high
energies ($\omega_{q}\gg\epsilon_{k},\Delta$) and short wavelengths ($q\gg
k,k_{F},\sqrt{2m\Delta}$) to the self-energy $\Sigma^{\rm bc}$. We consider
two limits: either $\omega_{q}$ tends to $+\infty$ at fixed $q$, or
$\omega_{q}$ and $q$ both tend to $+\infty$ with a fixed ratio
$\omega_{q}/q^{2}$. Again, we use straightforward units of $\Delta$ (for
instance $k_{\Delta}a\to a$ with $k_{\Delta}=\sqrt{2m\Delta}$), except for the
matrix elements of the pair propagator, where we include a factor
$(2\pi)^{3}$: $(2\pi)^{3}\Delta
N_{s,s^{\prime}}/(k_{\Delta}^{3}V)\to\check{N}_{s,s^{\prime}}$ and
correspondingly
$(k_{\Delta}^{3}V)\Gamma_{s,s^{\prime}}/(2\pi)^{3}\Delta\to\check{\Gamma}_{s,s^{\prime}}$.

We first give the equivalent of the angular integrals computed in App. B. In
the first limit, $\omega_{q}\gg q^{2},k^{2},|\mu|,1$, one has

$\displaystyle I_{U^{2}}$ $\displaystyle\sim$
$\displaystyle\frac{{\frac{1}{x_{\rm min}}-\frac{1}{x_{\rm
max}}}}{4kq\omega_{q}}$ (70) $\displaystyle I_{V^{2}}$ $\displaystyle\sim$
$\displaystyle\frac{{x_{\rm max}-x_{\rm min}}}{4kq\omega_{q}}$ (71)
$\displaystyle I_{UV}$ $\displaystyle\sim$
$\displaystyle\frac{\log\frac{x_{\rm max}}{x_{\rm min}}}{4kq\omega_{q}}$ (72)

and in the second limit, $\omega_{q}\approx q^{2}\gg k^{2},|\mu|,1$:

$\displaystyle I_{U^{2}}$ $\displaystyle\sim$
$\displaystyle\frac{2}{\omega_{q}+q^{2}}$ (73) $\displaystyle I_{V^{2}}$
$\displaystyle\sim$ $\displaystyle\frac{1}{2(\omega_{q}+q^{2})q^{4}}$ (74)
$\displaystyle I_{UV}$ $\displaystyle\sim$
$\displaystyle\frac{1}{(\omega_{q}+q^{2})q^{2}}.$ (75)

Next, we expand the pair propagator in the UV limit. Explicitly, the matrix
elements of the bare propagator in the cartesian basis (as opposed to the
phase-modulus basis [18, 51]) take the form:

$\displaystyle\check{N}_{++}(z,{\bf{q}})$ $\displaystyle=$
$\displaystyle\check{M}_{11}(z,{\bf{q}})=-\frac{\pi^{2}}{a}+\int
d^{3}k\left[\frac{U_{+}^{2}U_{-}^{2}}{z-\epsilon_{+}-\epsilon_{-}}-\frac{V_{+}^{2}V_{-}^{2}}{z+\epsilon_{+}+\epsilon_{-}}+\frac{1}{2k^{2}}\right]$
(76) $\displaystyle\check{N}_{--}(z,{\bf{q}})$ $\displaystyle=$
$\displaystyle\check{M}_{22}(z,{\bf{q}})=-\frac{\pi^{2}}{a}+\int
d^{3}k\left[\frac{V_{+}^{2}V_{-}^{2}}{z-\epsilon_{+}-\epsilon_{-}}-\frac{U_{+}^{2}U_{-}^{2}}{z+\epsilon_{+}+\epsilon_{-}}+\frac{1}{2k^{2}}\right]$
(77) $\displaystyle\check{N}_{+-}(z,{\bf{q}})$ $\displaystyle=$
$\displaystyle\check{M}_{12}(z,{\bf{q}})=\int
d^{3}k\left[\frac{U_{+}U_{-}V_{+}V_{-}}{z+\epsilon_{+}+\epsilon_{-}}-\frac{U_{+}U_{-}V_{+}V_{-}}{z-\epsilon_{+}-\epsilon_{-}}\right]$
(78)

with the notation $\epsilon_{\pm}=\epsilon_{{\bf{k}}\pm{\bf{q}}/2}$ and
similarly for $U_{\pm}$ and $V_{\pm}$. Note that we gave in passing the
equivalence between our notation $N_{s,s^{\prime}}$ and the notation $M_{ij}$
of Refs.[18, 51]. At large $\omega_{q}$, those integrals are dominated by the
large-$k$ region. We thus expand at $k^{2}\gg 1,|\mu|/\Delta$ (but a priori
$q\approx k$):

$\displaystyle U_{{\bf{k}}}$ $\displaystyle=$ $\displaystyle 1+O(k^{-4})$ (79)
$\displaystyle V_{{\bf{k}}}$ $\displaystyle=$
$\displaystyle\frac{1}{2k^{2}}+O(k^{-4})$ (80)
$\displaystyle\epsilon_{{\bf{k}}+{\bf{q}}/2}+\epsilon_{{\bf{k}}-{\bf{q}}/2}$
$\displaystyle=$ $\displaystyle 2k^{2}+q^{2}/2+O(1).$ (81)

With this, we obtain the expansions of the matrix elements in the limit
$\omega_{q},q^{2}\gg,k^{2},|\mu|,1$:

$\displaystyle\text{Re}\check{N}_{--}(\omega_{q},{\bf{q}})$ $\displaystyle=$
$\displaystyle\pi^{2}\left[\frac{\sqrt{2\omega_{q}+q^{2}}}{2}-\frac{1}{a}\right]+O(1)$
(82)
$\displaystyle\text{Im}\check{N}_{--}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle\frac{\pi^{2}}{q^{7}}n_{--}\left(\frac{2\omega_{q}}{q^{2}}\right)+O(q^{-8})\quad\text{with}\quad
n_{--}(\alpha)=-\frac{4\sqrt{\alpha-1}}{(\alpha-2)^{2}\alpha^{2}}-\frac{\text{ln}\left[\frac{\sqrt{\alpha-1}+1}{\sqrt{\alpha-1}-1}\right]^{2}}{\alpha^{3}}$
$\displaystyle\check{N}_{++}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle=-\pi^{2}\left[\mathrm{i}\frac{\sqrt{2\omega_{q}-q^{2}}}{2}+\frac{1}{a}\right]+O(1)$
(83) $\displaystyle\check{N}_{+-}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle\frac{\pi^{2}}{q^{3}}n_{+-}\left(\frac{2\omega_{q}}{q^{2}}\right)+O(q^{-4})$
$\displaystyle\text{with}\quad
n_{+-}(\alpha)=\frac{4}{\pi}\mathcal{P}\int_{0}^{+\infty}\frac{t\mathrm{d}t\text{Argth}\left(\frac{2t}{t^{2}+1}\right)}{(1+t^{2})^{2}-\alpha^{2}}+\frac{\mathrm{i}}{\alpha}\text{Argth}\left(\frac{2\sqrt{\alpha-1}}{\alpha}\right).$

For the pair propagator $\Gamma$, this leads to

$\displaystyle\text{Re}\check{\Gamma}_{--}(\omega_{q},{\bf{q}})$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}\left[\frac{\sqrt{2\omega_{q}+q^{2}}}{2}-\frac{1}{a}\right]}+O(q^{-3})$
(85)
$\displaystyle\text{Im}\check{\Gamma}_{--}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}q^{9}}\gamma_{22}\left({2\omega_{q}}/{q^{2}}\right)+O(q^{-10})$
$\displaystyle\text{with}\quad\gamma_{22}(\alpha)=-\frac{4}{\alpha+1}\left(n_{--}(\alpha)-\frac{2\text{Re}n_{+-}^{2}(\alpha)}{\sqrt{\alpha-1}}\right)$
$\displaystyle\check{\Gamma}_{++}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle-\frac{1}{\pi^{2}\left[\mathrm{i}\frac{\sqrt{2\omega_{q}-q^{2}}}{2}+\frac{1}{a}\right]}+O(q^{-2})$
(87) $\displaystyle\check{\Gamma}_{+-}(\omega_{q}+\mathrm{i}0^{+},{\bf{q}})$
$\displaystyle=$
$\displaystyle\frac{4}{\pi^{2}q^{5}}\frac{n_{+-}(2\omega_{q}/{q^{2}})}{\sqrt{4\omega_{q}^{2}/{q^{4}}-1}}+O(q^{-6}).$
(88)

In the more stringent limiting case $\omega_{q}\gg q^{2}$, those asymptotic
behaviors simplify to

$\displaystyle\text{Re}\check{N}_{--}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}\sqrt{\omega_{q}}}{\sqrt{2}}+O(1)$ (89)
$\displaystyle\text{Im}\check{N}_{--}$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}}{\sqrt{2}{\omega_{q}}^{7/2}}+O(\omega_{q}^{-4})$
(90) $\displaystyle\text{Re}\check{N}_{++}$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}}{a}+O({\omega_{q}}^{-1/2})$ (91)
$\displaystyle\text{Im}\check{N}_{++}$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}\sqrt{\omega_{q}}}{\sqrt{2}}+O(1)$ (92)
$\displaystyle\check{N}_{+-}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{\sqrt{2}{\omega_{q}}^{3/2}}(\mathrm{i}-1)+O(\omega_{q}^{-2})$
(93) $\displaystyle\text{Re}\check{\Gamma}_{--}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\pi^{2}\sqrt{\omega_{q}}}+O({\omega_{q}}^{-1})$
(94) $\displaystyle\text{Im}\check{\Gamma}_{--}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\pi^{2}\omega_{q}^{9/2}}+O({\omega_{q}}^{-5})$
(95) $\displaystyle\text{Re}\check{\Gamma}_{++}$ $\displaystyle=$
$\displaystyle-\frac{2}{\pi^{2}a\omega_{q}}+O(\omega_{q}^{-3/2})$ (96)
$\displaystyle\text{Im}\check{\Gamma}_{++}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\pi^{2}\sqrt{\omega_{q}}}+O(\omega_{q}^{-1})$
(97) $\displaystyle\Gamma_{12}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\pi^{2}\omega_{q}^{5/2}}(1+\mathrm{i})+O(\omega_{q}^{-3}).$
(98)

We then extract the UV behavior of the integrand of Eqs. (51)–(52), first at
finite $q$:

$\displaystyle-\check{\rho}_{++}(q,\omega_{q})I_{V^{2}}(\omega_{q}\pm z_{k})$
$\displaystyle\sim$ $\displaystyle\frac{\sqrt{2}\left(x_{\rm max}-x_{\rm
min}\right)}{4\pi^{3}kq\omega_{q}^{3/2}}$ (99)
$\displaystyle-\check{\rho}_{--}(q,\omega_{q})I_{U^{2}}(\omega_{q}\pm z_{k})$
$\displaystyle\sim$ $\displaystyle\frac{\sqrt{2}\left(x_{\rm min}^{-1}-x_{\rm
max}^{-1}\right)}{4\pi^{3}kq\omega_{q}^{11/2}}$ (100)
$\displaystyle-\check{\rho}_{+-}(q,\omega_{q})I_{UV}(\omega_{q}\pm z_{k})$
$\displaystyle\sim$ $\displaystyle\frac{\sqrt{2}\left(\text{log}(x_{\rm
max}/x_{\rm min})\right)}{4\pi^{3}kq\omega_{q}^{7/2}}.$ (101)

With this, we compute analytically the integral over $\omega_{q}$ from a UV
cutoff $\Omega_{0}\gg 1$ to $+\infty$:

$\displaystyle J_{++}(q)$ $\displaystyle\equiv$
$\displaystyle-\int_{\Omega_{0}}^{+\infty}\check{\rho}_{++}(q,\omega_{q})I_{V^{2}}(\omega_{q}\pm
z_{k})\underset{\Omega_{0}\to+\infty}{\sim}\frac{\sqrt{2}\left(x_{\rm
max}-x_{\rm min}\right)}{2\pi^{3}kq\Omega_{0}^{1/2}}$ (102) $\displaystyle
J_{--}(q)$ $\displaystyle\equiv$
$\displaystyle-\int_{\Omega_{0}}^{+\infty}\check{\rho}_{--}(q,\omega_{q})I_{U^{2}}(\omega_{q}\pm
z_{k})\underset{\Omega_{0}\to+\infty}{\sim}\frac{\sqrt{2}\left(x_{\rm
min}^{-1}-x_{\rm max}^{-1}\right)}{18\pi^{3}kq\Omega_{0}^{9/2}}$ (103)
$\displaystyle J_{+-}(q)$ $\displaystyle\equiv$
$\displaystyle-\int_{\Omega_{0}}^{+\infty}\check{\rho}_{+-}(q,\omega_{q})I_{UV}(\omega_{q}\pm
z_{k})\underset{\Omega_{0}\to+\infty}{\sim}\frac{\sqrt{2}\left(\text{log}(x_{\rm
max}/x_{\rm min})\right)}{10\pi^{3}kq\Omega_{0}^{5/2}}.$ (104)

Finally we compute analytically the integral over $q$ from $Q\gg 1$ to
$+\infty$. At such large $q$, the subintegral over $\omega_{q}$ runs from
$q^{2}/2$ ($\alpha=1$) to infinity. Setting $\Omega_{0}=q^{2}/2$ in Eqs.
(102–104), and expanding for $q\to+\infty$, we obtain:

$\displaystyle
q^{2}J_{++}(q)\underset{q\to+\infty}{\sim}\frac{2}{\sqrt{3}\pi^{3}q^{5}}\quad\text{such
that}\quad\int_{Q}^{+\infty}q^{2}J_{++}(q)$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{3}\pi^{3}Q^{4}}$ (105)
$\displaystyle\int_{Q}^{+\infty}q^{2}J_{--}(q)$ $\displaystyle=$
$\displaystyle O(1/Q^{8})$ (106)
$\displaystyle\int_{Q}^{+\infty}q^{2}J_{+-}(q)$ $\displaystyle=$
$\displaystyle O(1/Q^{6}).$ (107)

###### Acknowledgements.

Discussions with C. A. R. Sá de Melo are gratefully acknowledged. SVL was
supported by a Fellowship of the Belgian American Educational Foundation. The
computational resources and services used in this work were provided by the
HPC core facility CalcUA of the Universiteit Antwerpen, and VSC (Flemish
Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and
the Flemish Government.

## References

  * Nozières and Pines [1966] P. Nozières and D. Pines,  _The theory of quantum liquids_ (W.A. Benjamin, New York, 1966).
  * Fetter and Walecka [1971] A. L. Fetter and J. D. Walecka,  _Quantum theory of many-particle systems_ (McGraw-Hill, San Francisco, 1971).
  * Landau and Khalatnikov [1949] L. Landau and I. M. Khalatnikov, Teoriya vyazkosti Geliya-II, Zh. Eksp. Teor. Fiz. 19, 637 (1949).
  * Pitaevskii [1959] L. Pitaevskii, Properties of the spectrum of elementary excitations near the disintegration threshold of the excitations, Zh. Eksp. Teor. Fiz. 36, 1168 (1959), [Sov. Phys. JETP, Vol. 9, No. 4, p. 830 (1959)].
  * Nozières [1964] P. Nozières,  _Theory of interacting Fermi systems_ (W.A. Benjamin, New York, 1964).
  * Godfrin  _et al._ [2012] H. Godfrin, M. Meschke, H.-J. Lauter, A. Sultan, H. M. Böhm, E. Krotscheck, and M. Panholzer, Observation of a roton collective mode in a two-dimensional Fermi liquid, Nature 483, 576 (2012).
  * Tucker and Wyatt [1992] M. A. H. Tucker and A. F. G. Wyatt, Four-phonon scattering in superfluid 4 He, Journal of Physics: Condensed Matter 4, 7745 (1992).
  * Vollhardt and Woelfle [1990] D. Vollhardt and P. Woelfle,  _The Superfluid Phases Of Helium 3_ (Taylor & Francis, 1990).
  * Zwerger [2011] W. Zwerger, ed.,  _The BCS-BEC crossover and the unitary Fermi gas_ (Springer Science & Business Media, 2011).
  * Kaplan  _et al._ [1976] S. B. Kaplan, C. C. Chi, D. N. Langenberg, J. J. Chang, S. Jafarey, and D. J. Scalapino, Quasiparticle and phonon lifetimes in superconductors, Phys. Rev. B 14, 4854 (1976).
  * Fåk  _et al._ [2012] B. Fåk, T. Keller, M. E. Zhitomirsky, and A. L. Chernyshev, Roton-Phonon Interactions in Superfluid ${}^{4}\mathrm{He}$, Phys. Rev. Lett. 109, 155305 (2012).
  * Chomaz  _et al._ [2018] L. Chomaz, R. M. W. van Bijnen, D. Petter, G. Faraoni, S. Baier, J. H. Becher, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Observation of roton mode population in a dipolar quantum gas, Nature Physics 14, 442 (2018).
  * Schulze  _et al._ [2001] H.-J. Schulze, A. Polls, and A. Ramos, Pairing with polarization effects in low-density neutron matter, Phys. Rev. C 63, 044310 (2001).
  * Cao  _et al._ [2006] L. G. Cao, U. Lombardo, and P. Schuck, Screening effects in superfluid nuclear and neutron matter within Brueckner theory, Phys. Rev. C 74, 064301 (2006).
  * 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, Physics Reports 738, 1 (2018), the BCS–BEC crossover: From ultra-cold Fermi gases to nuclear systems.
  * Marini  _et al._ [1998] M. Marini, F. Pistolesi, and G. Strinati, Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions, The European Physical Journal B - Condensed Matter and Complex Systems 1, 151 (1998).
  * Combescot  _et al._ [2006] R. Combescot, M. Y. Kagan, and S. Stringari, Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover, Physical Review A 74, 042717 (2006).
  * 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).
  * Anderson [1958] P. W. Anderson, Random-Phase Approximation in the Theory of Superconductivity, Phys. Rev. 112, 1900 (1958).
  * Andrianov and Popov [1976] V. A. Andrianov and V. N. Popov, Gidrodinamičeskoe dejstvie i Boze-spektr sverhtekučih Fermi-sistem, Teoreticheskaya i Matematicheskaya Fizika 28, 341 (1976), [English translation: Theoretical and Mathematical Physics, 1976, 28:3, 829–837].
  * Sooryakumar and Klein [1980] R. Sooryakumar and M. V. Klein, Raman Scattering by Superconducting-Gap Excitations and Their Coupling to Charge-Density Waves, Phys. Rev. Lett. 45, 660 (1980).
  * Méasson  _et al._ [2014] M.-A. Méasson, Y. Gallais, M. Cazayous, B. Clair, P. Rodière, L. Cario, and A. Sacuto, Amplitude Higgs mode in the $2H-{\text{NbSe}}_{2}$ superconductor, Phys. Rev. B 89, 060503 (2014).
  * Behrle  _et al._ [2018] A. Behrle, T. Harrison, J. Kombe, K. Gao, M. Link, J. S. Bernier, C. Kollath, and M. Köhl, Higgs mode in a strongly interacting fermionic superfluid, Nature Physics 10.1038/s41567-018-0128-6 (2018).
  * Kurkjian  _et al._ [2019] H. Kurkjian, S. N. Klimin, J. Tempere, and Y. Castin, Pair-Breaking Collective Branch in BCS Superconductors and Superfluid Fermi Gases, Phys. Rev. Lett. 122, 093403 (2019).
  * Zhou and Ma [2021] H. Zhou and Y. Ma, Thermal conductivity of an ultracold fermi gas in the bcs-bec crossover, Scientific Reports 11, 1228 (2021).
  * Artemenko and Volkov [1975] S. Artemenko and A. Volkov, Collective excitations with a sound spectrum in superconductors, Zh. Eksp. Teor. Fiz. 69, 1764 (1975), [Sov. Phys. JETP, Vol. 42, No. 5, p. 896, 1975].
  * Volkov and Kogan [1973] A. Volkov and C. M. Kogan, Collisionless relaxation of the energy gap in superconductors, Zh. Eksp. Teor. Fiz. 65, 2038 (1973), [Sov. Phys. JETP, Vol. 38, No. 5, p. 1018].
  * Gurarie [2009] V. Gurarie, Nonequilibrium Dynamics of Weakly and Strongly Paired Superconductors, Phys. Rev. Lett. 103, 075301 (2009).
  * Klimin  _et al._ [2015] S. N. Klimin, J. Tempere, G. Lombardi, and J. T. Devreese, Finite temperature effective field theory and two-band superfluidity in Fermi gases, The European Physical Journal B 88, 122 (2015).
  * 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).
  * 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).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Ketterle and Zwierlein [2008] W. Ketterle and M. W. Zwierlein, Making, probing and understanding ultracold Fermi gases, in  _Ultra-cold Fermi Gases_, Proceedings of the International School of Physics “Enrico Fermi”, Course CLXIV, edited by M. Inguscio, W. Ketterle, and C. Salomon (IOS Press, Amsterdam, 2008).
  * Stoof  _et al._ [2009] H. T. Stoof, K. B. Gubbels, and D. Dickerscheid,  _Ultracold quantum fields_ (Springer, 2009).
  * Veeravalli  _et al._ [2008] G. Veeravalli, E. Kuhnle, P. Dyke, and C. J. Vale, Bragg Spectroscopy of a Strongly Interacting Fermi Gas, Phys. Rev. Lett. 101, 250403 (2008).
  * Hoinka  _et al._ [2017] S. Hoinka, P. Dyke, M. G. Lingham, J. J. Kinnunen, G. M. Bruun, and C. J. Vale, Goldstone mode and pair-breaking excitations in atomic Fermi superfluids, Nature Physics 13, 943 (2017).
  * Kuhn  _et al._ [2020] C. C. N. Kuhn, S. Hoinka, I. Herrera, P. Dyke, J. J. Kinnunen, G. M. Bruun, and C. J. Vale, High-Frequency Sound in a Unitary Fermi Gas, Phys. Rev. Lett. 124, 150401 (2020).
  * Biss  _et al._ [2021] H. Biss, L. Sobirey, N. Luick, M. Bohlen, J. J. Kinnunen, G. M. Bruun, T. Lompe, and H. Moritz, Excitation Spectrum and Superfluid Gap of an Ultracold Fermi Gas, arXiv:2105.09820 (2021).
  * Gupta  _et al._ [2003] S. Gupta, Z. Hadzibabic, M. W. Zwierlein, C. A. Stan, K. Dieckmann, C. H. Schunck, E. G. M. van Kempen, B. J. Verhaar, and W. Ketterle, Radio-Frequency Spectroscopy of Ultracold Fermions, Science 300, 1723 (2003).
  * Shin  _et al._ [2007] Y. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, Tomographic rf Spectroscopy of a Trapped Fermi Gas at Unitarity, Phys. Rev. Lett. 99, 090403 (2007).
  * 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).
  * 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).
  * Mukherjee  _et al._ [2017] B. Mukherjee, Z. Yan, P. B. Patel, Z. Hadzibabic, T. Yefsah, J. Struck, and M. W. Zwierlein, Homogeneous Atomic Fermi Gases, Phys. Rev. Lett. 118, 123401 (2017).
  * 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).
  * Sagi  _et al._ [2015] Y. Sagi, T. E. Drake, R. Paudel, R. Chapurin, and D. S. Jin, Breakdown of the Fermi Liquid Description for Strongly Interacting Fermions, Phys. Rev. Lett. 114, 075301 (2015).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Microscopic Theory of Superconductivity, Phys. Rev. 106, 162 (1957).
  * 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).
  * Van Loon  _et al._ [2020] S. Van Loon, J. Tempere, and H. Kurkjian, Beyond Mean-Field Corrections to the Quasiparticle Spectrum of Superfluid Fermi Gases, Phys. Rev. Lett. 124, 073404 (2020).
  * Castin [2020] Y. Castin, Marche au hasard d’une quasi-particule massive dans le gaz de phonons d’un superfluide à très basse température, Comptes Rendus. Physique 21, 571 (2020).
  * Kulik  _et al._ [1981] I. O. Kulik, O. Entin-Wohlman, and R. Orbach, Pair susceptibility and mode propagation in superconductors: A microscopic approach, Journal of Low Temperature Physics 43, 591 (1981).
  * Kurkjian and Tempere [2017] H. Kurkjian and J. Tempere, Absorption and emission of a collective excitation by a fermionic quasiparticle in a Fermi superfluid, New Journal of Physics 19, 113045 (2017).
  * Klimin  _et al._ [2019] S. N. Klimin, J. Tempere, and H. Kurkjian, Phononic collective excitations in superfluid Fermi gases at nonzero temperatures, Phys. Rev. A 100, 063634 (2019).
  * Lerch  _et al._ [2008] N. Lerch, L. Bartosch, and P. Kopietz, Absence of Fermionic Quasiparticles in the Superfluid State of the Attractive Fermi Gas, Phys. Rev. Lett. 100, 050403 (2008).
  * Haussmann  _et al._ [2009] R. Haussmann, M. Punk, and W. Zwerger, Spectral functions and rf response of ultracold fermionic atoms, Physical Review A 80, 063612 (2009).
  * Gor’kov and Melik-Barkhudarov [1958] L. Gor’kov and T. Melik-Barkhudarov, Contribution to the theory of superfluidity in an imperfect Fermi gas, Zh. Eksp. Teor. Fiz. 40, 1452 (1958), [Sov. Phys. JETP, 13, 1018 (1958)].
  * 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).
  * Pisani  _et al._ [2018] L. Pisani, P. Pieri, and G. C. Strinati, Gap equation with pairing correlations beyond the mean-field approximation and its equivalence to a Hugenholtz-Pines condition for fermion pairs, Physical Review B 98, 104507 (2018).
  * Haussmann [1993] R. Haussmann, Crossover from BCS superconductivity to Bose-Einstein condensation: A self-consistent theory, Zeitschrift für Physik B Condensed Matter 91, 291 (1993).
  * 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).
  * Note [1] We assume here that the matrix elements of $G^{-1}$ do not vanish individually, which is generally the case.
  * Castin  _et al._ [2017] Y. Castin, A. Sinatra, and H. Kurkjian, Landau Phonon-Roton Theory Revisited for Superfluid He 4 and Fermi Gases, Physical Review Letters 119, 260402 (2017).
  * Note [2] All code used for our numerical calculations can be found on https://github.com/hkurkjian/CodeArXiv_2111_04692.
  * Zhang and Liu [2011] Z. Zhang and W. V. Liu, Finite-temperature damping of collective modes of a BCS-BEC crossover superfluid, Phys. Rev. A 83, 023617 (2011).
  * Klimin  _et al._ [2021] S. N. Klimin, J. Tempere, and H. Kurkjian, Collective excitations of superfluid Fermi gases near the transition temperature, Phys. Rev. A 103, 043336 (2021).
  * Castin and Kurkjian [2020] Y. Castin and H. Kurkjian, Branche d’excitation collective du continuum dans les gaz de fermions condensés par paires : étude analytique et lois d’échelle, Comptes Rendus. Physique 21, 253 (2020).
  * 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 CLXIV, edited by M. Inguscio, W. Ketterle, and C. Salomon (IOS Press, Amsterdam, 2008).
  * Kinnunen [2012] J. J. Kinnunen, Hartree shift in unitary Fermi gases, Phys. Rev. A 85, 012701 (2012).
  * Galitskii [1958] V. M. Galitskii, The energy spectrum of a non-ideal Fermi gas, Zh. Eksp. Teor. Fiz. 34, 151 (1958), [Sov. Phys. JETP, 7, 104 (1958)].