# Supplementary Material: Higgs oscillations in a unitary Fermi superfluid

###### Contents

  1. 1 Experimental protocol
    1. 1.1 Sample Preparation
    2. 1.2 Bragg Spectroscopy Calibration
    3. 1.3 Determination of Density
  2. 2 Order parameter dynamics within time-dependent BCS theory
    1. 2.1 Analytical solution in the small-amplitude regime
    2. 2.2 Numerical simulations
  3. 3 Application of time-dependent BCS theory to the experimental setup
    1. 3.1 Ramp-rate dependence
    2. 3.2 Inhomogeneous broadening of the signal
      1. 3.2.1 Qualitative discussion in the limit of a tightly focused beam
      2. 3.2.2 Numerical solutions in the general case
    3. 3.3 Averaging over the Bragg time
  4. 4 Comparing theoretical predictions and experimental measurements
    1. 4.1 Magnitude averaged over different time windows
    2. 4.2 Comparisons including the sensitivity of the Bragg pulse

##  1 Experimental protocol

###  1.1 Sample Preparation

In our experiments, we cool a balanced mixture of fermionic 6Li atoms in the
$|F=1/2,m_{\mathrm{F}}=\pm 1/2\rangle(\equiv|1\rangle-|2\rangle$, blue circles
in Fig. 2, main text) or
$|F=1/2,m_{\mathrm{F}}=1/2\rangle-|F=3/2,m_{\mathrm{F}}=-3/2\rangle(\equiv|1\rangle-|3\rangle$,
green squares) hyperfine states in a 100 W, 1075 nm single beam optical dipole
trap. Degeneracy is reached through evaporative cooling by smoothly lowering
the trap laser power at the magnetic field where the $\it{s}$-wave scattering
length diverges, $\it{a}$ $\rightarrow$ $\infty$. Subsequently, the atom cloud
is transferred to an oblate harmonic potential, formed by a combination of
optical and magnetic fields. The oblate trap is formed between two anti-nodes
of a cylindrically focused, 532-nm (blue-detuned), TEM01 mode laser beam [1,
2, 3], where the two anti-nodes are separated by $\approx$ 90 $\mu$m
propagating along the y-direction and the 1/e2 radius in the x-direction is
$\approx$ 1.0 mm. The optical potential produces the confinement in the $z$
direction and a very weak anti-confinement in the $x$-$y$ plane. The residual
magnetic field curvature from the Feshbach coils provides highly harmonic and
cylindrically symmetric confinement in the $x$-$y$ plane, which dominates the
anti-trapping of the optical potential. The measured trapping frequencies are
$\omega_{z}/2\pi$ = 105 Hz(1) and ($\omega_{x},\omega_{y}$) = 2$\pi$ $\times$
(24.5, 24.5) Hz at a magnetic field of $B$ = 832.2 G (corresponding to the
$|1\rangle-|2\rangle$ Feshbach resonance [4]). The asymmetry in the trapping
potential is $|\omega_{x}-\omega_{y}|/\omega_{r}\lesssim 0.01$. Note that
$\omega_{r}\propto\sqrt{B}$ so the radial confinement also changes when we
tune interactions. Typically, we produce clouds with temperatures of 0.09
$T/T_{\mathrm{F}}$, where $T_{\mathrm{F}}$ is the Fermi temperature and
$\it{N}$ = 3$\times$105 atoms per spin state.
To determine the temperature of the atom clouds (horizontal axis of Figs. 2
and 3 in the main text) we follow the quench procedure as described in the
main text, however, instead of probing the clouds with Bragg lasers following
the quench the cloud is held for 500 ms and an $in-situ$ absorption image is
taken (after this hold time, no dynamics in the density distribution are
detected [5]). The temperature can be determined by fitting the known equation
of state (EoS) for the pressure of a unitary Fermi gas [6] to the line
densities of a trapped atom cloud at equilibrium. Throughout this study we
probe temperatures ranging from 0.1 to 0.18 $T/T_{\mathrm{F}}$. The higher
temperatures are obtained by varying the endpoint of the initial evaporative
cooling process, which loads a higher temperature cloud into the oblate
trapping potential and the temperature is extracted following the same
procedure described above.

###  1.2 Bragg Spectroscopy Calibration

The post-quench dynamics are probed using two tightly-focused Bragg laser
beams, as shown in Fig. 1(a), that intersect in the centre of the cloud where
the density is near-uniform [7, 8]. Two-photon Bragg spectroscopy is both
energy and momentum selective thus, by setting the frequency difference to
half of the atomic recoil frequency ($\omega_{r}/2=\hbar k^{2}/4m$), pairs
with zero center-of-mass momentum are resonantly excited [9, 8]. These pairs
then begin moving with a velocity $\hbar\mathbf{k}/m$ where,
$\mathbf{k}=\mathbf{k}_{a}-\mathbf{k}_{b}$ is the difference of the wave-
vectors of the two Bragg lasers. We use $k\simeq 4k_{F}$ to ensure that
$\omega$ is large compared to $E_{\mathrm{F}}$, and the pulse duration (50
$\mu$s) is relatively short with respect to the dynamics. The center of mass
displacement $\Delta X_{\mathrm{CoM}}$ following time-of-flight expansion,
reveals the momentum imparted to the atoms by the Bragg lasers. Fig. 1(b)
shows an absorption image of an atom cloud following the experimental sequence
described in the main text. A large number of atoms are seen to be displaced
from the centre of the cloud towards the right of the image. The crescent
shape in Fig. 1(b) is a result of collisions between scattered atoms and the
unperturbed atoms in the cloud which occur following the Bragg pulse [10].

We have measured the sensitivity of the centre of mass displacement, $\Delta
X_{\mathrm{CoM}}$ at unitarity as the pair condensate fraction varies. To
determine the centre of mass displacement, $\Delta X_{\mathrm{CoM}}$, we
prepare clouds at a range of initial temperatures by varying the endpoint of
the evaporation and then applying a Bragg pulse (to the clouds at equilibrium)
and measure the response as described above. To determine the pair condensate
fraction a degenerate cloud at 832.2 G is produced at a given temperature and
in equilibrium. Subsequently, we jump the magnetic field far onto the BEC side
of the Feshbach resonance in $\approx$ 50 $\mu$s and simultaneously turn off
the optical ($z$) confinement. The weakly bound pairs are converted to tightly
bound molecules, whose centre of mass momentum is preserved [11]. The weakly
interacting molecules then expand ballistically along $z$ before the magnetic
field is ramped back to 832.2 G in 2.5 ms to dissociate the molecules. The
total expansion time is approximately one quarter of the radial trapping
period, to best reveal the momentum distribution. Finally, an absorption image
is taken and fitted with a bi-modal distribution to determine the condensate
fraction.

Figure 1: (Color online) (a) Schematic of the experimental setup for focused
beam Bragg spectroscopy. Two laser beams with wave vectors $\mathbf{k}_{a}$
and $\mathbf{k}_{b}$ and frequencies $\omega_{a}$ and $\omega_{b}$ are focused
into the center of a trapped atom cloud. The beams have a 1/e2 radius of 15
$\mu$m and intersect at an angle of 2$\theta$ = 72.6∘ (b) Absorption image of
an atom cloud following an excitation with a Bragg frequency
$\omega/(2\pi)=+50$ kHz. The dashed rectangle indicates the region used for
determination of the centre of mass displacement. (c) Measured center of mass
displacement, $\Delta X_{\mathrm{CoM}}$ as a function of measured pair
condensate fraction.

Fig. 1(c) shows the response of the $\Delta X_{\mathrm{CoM}}$ of the cloud as
the pair condensate fraction, and therefore the temperature, is varied. As the
condensate fraction is reduced, the centre of mass displacement, $\Delta
X_{\mathrm{CoM}}$ also decreases. The Bragg signal shows a clear dependence on
the condensate fraction, which reflects the change in the number of condensed
pairs available to take part in the scattering process.

###  1.3 Determination of Density

To observe the Higgs oscillation most clearly the Bragg scheme addresses atoms
at the centre of the trap, where the density is near homogeneous, as in
previous studies [7, 8]. We define the mean density $\bar{n}$ in the Bragg
volume as

$\bar{n}=\frac{\int\Omega_{\mathrm{Br}}(\mathbf{r})n(\mathbf{r})\mathrm{d}^{3}\mathbf{r}}{\int\Omega_{\mathrm{Br}}(\mathbf{r})\mathrm{d}^{3}\mathbf{r}},$
(1)

where $\Omega_{\mathrm{Br}}(\mathbf{r})$ is the spatially dependent two-photon
Rabi frequency, which is proportional to the geometric mean of the intensities
of the two laser beams $\sqrt{I_{a}(\mathbf{r})I_{b}(\mathbf{r})}$ and
$n(\mathbf{r})$ is the 3D density distribution. The mean density sets the
relevant energy scale and the Higgs oscillation frequency via Fermi energy
$E\mathrm{{}_{F}}=\frac{\hbar}{2m}(3\pi^{2}\bar{n})^{2/3}$, and the normalized
frequency $\hbar\omega_{\rm H}/2E\mathrm{{}_{F}}$, where $m$ is the atomic
mass and $\omega_{\rm H}$ is the Higgs oscillation frequency.

Figure 2: (Color online) A selection of measured centre of mass displacement
($\Delta X_{\mathrm{CoM}}$) vs hold time, $t_{h}$ in absolute units of $\mu$s
in (a) and relative to the local Fermi time $\tau_{F}=\hbar/E_{F}$ in (b) for
atom clouds prepared at different temperatures. Also plotted are the fitted
sinusoidal functions used to determine the frequency and damping of the
amplitude oscillation. Different temperature clouds can have different
densities leading to the different absolute oscillation frequencies, (a). When
the time axis is taken relative to the Fermi time, (b) the relative
frequencies agree more closely, as seen in the main text.

To determine the mean density, the true (3D) density distribution
$n$($\mathbf{r}$) of the trapped cloud and the intensity product of the Bragg
laser beams with their intersecting Gaussian spatial profiles are combined.
The 3D density $n$($\mathbf{r}$) is extracted through applying an inverse Abel
transform to absorption images of unperturbed trapped clouds. Images of
trapped clouds are taken using a high intensity ($I/I_{\mathrm{sat}}>10$) and
short imaging pulse length of 1 $\mu s$. The inverse Abel transform method
uses a Fourier decomposition where the radial density distribution is expanded
in a Fourier series. This method requires no direct differentiation and allows
reconstruction of the central density without singularities. From this, we are
able to determine $\bar{n}$ with an error of approximately 7.5 %, which leads
to a 5 % error in the Fermi energy.

In Fig. 2, we plot the direct measured centre of mass displacement ($\Delta
X_{\mathrm{CoM}}$) vs. time along with the fitted sinusoid, for the same data
as in Fig. 2 of the main text. These traces show how the total Bragg signal
decreases with increasing temperature, the magnitude of the oscillations falls
off with increasing temperatures and the frequency changes for clouds prepared
at different temperatures. In Fig 2(a) the hold time is in ms and in 2(b) in
units of the local Fermi time $\tau_{F}=\hbar/E_{F}$. The absolute oscillation
frequency is set by the local Fermi energy, which in turn is set by the local
density. Clouds at different temperatures may have different densities due to
changes in the preparation sequence. This leads to the variations in the
absolute oscillation frequency in (a), which become more consistent when
plotted as a function of the scaled time in (b).

##  2 Order parameter dynamics within time-dependent BCS theory

We recall here what is theoretically expected for the time evolution of a
superfluid Fermi gas after a rapid change of the interatomic interactions at
zero temperatures. We consider a two-component ($\uparrow$ and $\downarrow$)
Fermi gas interacting through a pairwise single-channel $s$-wave interaction,
which captures broad, entrance-channel dominated Feshbach resonances [12].
Assuming the gas is homogeneous in a cubic volume $V$ (see Sec. 3 for
consideration of inhomogeneities in the trapped gas) the momentum
representation of the Hamiltonian for this system is

$\hat{H}=\sum_{\mathbf{k}\sigma}\frac{k^{2}}{2m}\hat{a}^{\dagger}_{\mathbf{k}\sigma}\hat{a}_{\mathbf{k}\sigma}+\frac{g_{0}}{V}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}\hat{a}^{\dagger}_{\mathbf{k}+\mathbf{q},\uparrow}\hat{a}^{\dagger}_{-\mathbf{k},\downarrow}\hat{a}_{-\mathbf{k}^{\prime},\downarrow}\hat{a}_{\mathbf{k}^{\prime}+\mathbf{q},\uparrow},$
(2)

where the interaction strength $g_{0}$ is renormalized to reproduce the
correct $s$-wave scattering length $a$ of the two-body problem. In the
remaining sections, we use the convention $\hbar=1$ throughout.

The equations of motion for the momentum distribution $n_{\mathbf{k}}$ and
pairing function $c_{\mathbf{k}}$ (Eqs. (1) and (2) in the main text) are
derived from the Hamiltonian in Eq. (2), using the BCS mean-field
approximation. At $t=0$ the system is at equilibrium at scattering length
$a_{i}$ and temperature $T<T_{c,i}$, we can therefore fix the initial
conditions for $n_{\mathbf{k}}$ and $c_{\mathbf{k}}$ using the BCS ground-
state solutions, as discussed in the main text

$\displaystyle n_{\mathbf{k}}(t=0)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\frac{\xi_{\mathbf{k},i}}{\epsilon_{\mathbf{k},i}}F_{\beta}(\epsilon_{\mathbf{k},i})\right),$
(3) $\displaystyle c_{\mathbf{k}}(t=0)$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{i}}{2\epsilon_{\mathbf{k},i}}F_{\beta}(\epsilon_{\mathbf{k},i}),$
(4)

where $\xi_{\mathbf{k},i}=k^{2}/2m-\mu_{i}$,
$\epsilon_{\mathbf{k},i}=\sqrt{\xi_{\mathbf{k},i}^{2}+\Delta_{i}^{2}}$,
$\mu_{i}$ and $\Delta_{i}$ are calculated at the initial scattering length
$a_{i}$, and where
$F_{\beta}(\epsilon)=\tanh(\beta\epsilon/2)=1-2(\mathrm{exp}(\beta\epsilon)+1)^{-1}$
is the thermal distribution with $\beta=1/k_{\mathrm{B}}T$. The limit of zero
temperature corresponds to $F_{\beta}(\epsilon)\to 1$ for $\epsilon>0$.
Conversely, the regime of temperatures close to the critical temperature
$T_{c}$ is found by taking the limit $\Delta\to 0$ [13].

###  2.1 Analytical solution in the small-amplitude regime

For shallow quenches, that is, small deviations of $n_{\mathbf{k}}$ and
$c_{\mathbf{k}}$ from their equilibrium values, the time-dependent BCS
equations (Eqs. (1) and (2) of the main text) can be solved analytically. At
zero-temperature, this is a well-studied problem [14, 15, 16]. Here we only
detail the expression of the normalized amplitude $A_{\rm th}$ of the
oscillations (see Eq. (4) of the main text).

In the final equilibrium state (at $a=a_{f}$ and $T=0$) and in the modulus-
phase basis [17], the linear-response matrix of the order-parameter is given
by

$M(\omega)=\begin{pmatrix}\omega^{2}f(\omega)&\omega g(\omega)\\\ \omega
g(\omega)&(\omega^{2}-4\Delta_{f}^{2})f(\omega)\end{pmatrix},$ (5)

where $\Delta_{f}$ and $\mu_{f}$ are the gap and chemical potential,
respectively, of the final equilibrium state, and we have introduced the
functions:

$\displaystyle f(\omega)$ $\displaystyle=$
$\displaystyle\frac{\Delta_{f}}{V}\sum_{\mathbf{k}}\frac{1}{2\epsilon_{\mathbf{k}}(\omega^{2}-4\epsilon_{\mathbf{k}}^{2})},$
(6) $\displaystyle g(\omega)$ $\displaystyle=$
$\displaystyle\frac{\Delta_{f}}{V}\sum_{\mathbf{k}}\frac{E_{\mathbf{k}}-\mu_{f}}{\epsilon_{\mathbf{k}}(\omega^{2}-4\epsilon_{\mathbf{k}}^{2})}.$
(7)

The amplitude-amplitude response function introduced in Eq. (3) of the main
text follows directly from $M$ through the relation

$\chi_{|\Delta||\Delta|}(\omega)=(M^{-1}(\omega))_{22}=\frac{f(\omega)}{(\omega^{2}-4\Delta_{f}^{2})f^{2}(\omega)-g^{2}(\omega)}.$
(8)

As mentioned in the main text, the imaginary part
$\chi_{|\Delta||\Delta|}^{\prime\prime}$ of this function has a squareroot
divergence in $\omega=2\Delta$ for $\mu_{f}>0$ and a square root cancellation
for $\mu_{f}<0$. We denote by $f_{l}$ the spectral weight of this squareroot
singularity:

$\chi_{|\Delta||\Delta|}^{\prime\prime}(\omega)\underset{\omega\to\omega_{\rm
th}}{\sim}\begin{cases}f_{l}\sqrt{\frac{\omega_{\rm th}}{\omega-\omega_{\rm
th}}}\text{ when }\mu_{f}>0\text{ and }\omega_{\rm th}=2\Delta_{f}\\\
f_{l}\sqrt{\frac{\omega-\omega_{\rm th}}{\omega_{\rm th}}}\text{ when
}\mu_{f}<0\text{ and }\omega_{\rm
th}=2\sqrt{\Delta_{f}^{2}+\mu_{f}^{2}}\end{cases}.$ (9)

In general, linear response theory expresses time-dependent quantities as
Laplace transforms of frequency responses. For the evolution of the amplitude
of the order parameter, with initial conditions Eqs. (3) and (4), this gives

$\frac{\Delta(t)-\Delta_{f}}{\Delta_{i}-\Delta_{f}}=1-\frac{1}{\chi_{|\Delta||\Delta|}(0)}\int_{+\infty+\mathrm{i}\eta}^{-\infty+\mathrm{i}\eta}\frac{\mathrm{e}^{-\mathrm{i}zt}}{2\mathrm{i}\pi
z}\chi_{|\Delta||\Delta|}(z).$ (10)

This integral is deformed to a contour enclosing the branch cut
$[2\Delta,+\infty)$ of $\chi_{|\Delta||\Delta|}$, as well as the pole of the
integrand in $z=0$. Since only the contribution of the pole survives at long
times, one sees that $\lim_{t\to+\infty}\Delta(t)-\Delta_{f}=0$, in other
words

$\Delta_{\infty}=\Delta_{f}$ (11)

in the small-amplitude regime and at $T=0$. From the contribution of the
branch cut we then derive the version of Eq. (4) of the main text that is
valid at all times:

$\frac{\Delta(t)-\Delta_{f}}{\Delta_{i}-\Delta_{f}}=\frac{2}{\pi}\int_{2\Delta_{f}}^{+\infty}\frac{\cos\omega
t}{\omega}\frac{\chi^{\prime\prime}_{|\Delta||\Delta|}(\omega)}{\chi_{|\Delta||\Delta|}(0)}\mathrm{d}\omega.$
(12)

At long times, this frequency integral is evaluated by approximating
$\chi^{\prime\prime}_{|\Delta||\Delta|}$ by its behavior (Eq. (9)) near the
pair-breaking threshold. This yields the expression of $A_{\rm th}$ (see Eq.
(4) of the main text):

$A_{\rm
th}\underset{T=0}{=}\frac{2f_{l}}{\sqrt{\pi}\chi_{|\Delta||\Delta|}(0)}$ (13)

with $\chi_{|\Delta||\Delta|}(0)$ the static modulus response and $f_{l}$ the
spectral weight of the pair-breaking threshold.

###  2.2 Numerical simulations

Outside the small amplitude regime, we evolve the equations of motion (Eqs.
(1) and (2) of the main text) using a Runge-Kutta method [18]. The momentum-
space integrals in the gap ($\Delta=(g_{0}/V)\sum_{\mathbf{k}}c_{\mathbf{k}}$)
and number ($N=2\sum_{\mathbf{k}}n_{\mathbf{k}}$) equations are evaluated in
spherical coordinates using a fine grid in momentum space111Numerically, we
evaluate the momentum-space integral of a generic function $f(k)$ using $\int
d^{3}k/(2\pi)^{3}f(k)\rightarrow\sum_{i}w_{i}k_{i}^{2}f(k_{i})/2\pi^{2}$,
where $k_{i}$ and $w_{i}$ are the abscissas and weights of the Gauss-Legendre
quadrature method [18]. Because $n_{\mathbf{k}}$ and $c_{\mathbf{k}}$ become
small after a few $k_{F}$, we split the integration in two intervals:
$[0,10k_{F}]$ and $[10k_{F},\Lambda]$, where $\Lambda=350k_{F}$ is a momentum
cut-off (large enough such that numerical results are insensitive to it). In
every interval, we take $1500$ points, so that the first interval is much more
dense than the second one, where $n_{\mathbf{k}}$ and $c_{\mathbf{k}}$ are
$\approx 0$.. The conservation of the total number of atoms $N$ provides a
check on the integrity of the numerics. Similar dynamical problems near
resonance have been studied in Refs. [19, 20].

In this section, we consider an infinitely fast quench from
$1/(k_{F}a_{i})=-0.18$ to $1/(k_{F}a_{f})=0$ (see Sec. 3 for a discussion of
the nonzero ramping time). As shown in Fig. 3, the relative variation of the
equilibrium order parameter is small at $T=0$
($\delta\Delta\equiv|\Delta_{f}-\Delta_{i}|\approx 0.16\Delta_{f}$) but
increases with temperature, tending to 1 at $T_{c}$. Therefore, this limits
the applicability of the analytical results presented in Sec. 2.1 to the small
temperature regime.

Figure 3:  Temperature dependence of the difference in the equilibrium order
parameter from the initial and final state
$\delta\Delta\equiv|\Delta_{f}-\Delta_{i}|$. The initial and final state
correspond to $1/(k_{F}a_{i})=-0.18$ and $1/(k_{F}a_{f})=0$, respectively. The
temperature is rescaled by the critical temperature at $a_{i}$, that is
$T_{c,i}=0.399T_{F}$

.

Using numerical simulations, we can instead explore the entire temperature
range $0\leq T\lesssim T_{c,i}$. For example, Fig. 4(a) shows the time
evolution of the amplitude of the order parameter $\Delta$ at temperature
$T=0.25T_{c,i}$. As expected from literature (see e.g. Ref. [16] for $T=0$)
and discussed in the main text, the order parameter asymptotes to a value
$\Delta_{\infty}$ a bit less than the expected final state equilibrium value
$\Delta_{f}$ and oscillates at a frequency $\omega_{\rm H}=2\Delta_{\infty}$.
Moreover, we expect a power-law damping coefficient $\gamma=0.5$, which is
found if we neglect the very short time dynamics, as shown in Fig. 4(b) for
different temperatures. Finally, in Fig. 4(c), we show the temperature
dependence of the Higgs oscillation frequency (orange dash-dotted line) and
the final order parameter (pink solid), and compare it to the experimental
data. The temperature is scaled to the critical temperature of the initial
scattering length (for the experimental data we have used
$T_{c,i}=0.15T_{F}$). The comparison is not meant to be quantitative, as BCS
theory is known to overestimate both the superfluid gap and critical
temperature at unitarity [13], but rather to confront the trend in the
vicinity of $T_{c}$.

Figure 4: (a) Numerical result for the time evolution of the amplitude of the
order parameter $\Delta(t)$ at $T=0.25T_{c,i}$ after a fast quench from
$1/(k_{F}a_{i})=-0.18$ to $1/(k_{F}a_{f})=0$. The frequency of oscillation is
$\omega_{\rm H}=2\Delta_{\infty}\lesssim 2\Delta_{f}$. (b) Damping coefficient
$\gamma$ for $\Delta(t)$ as a function of the starting point of the fit
$t_{\mathrm{st}}$ for different temperatures. Here, we have used as a fitting
function $\Delta(t)=\Delta_{\infty}+A\cos(\omega_{\rm H}t+\phi)/t^{\gamma}$ in
the range $[t_{\mathrm{st}},t_{\mathrm{end}}]$, where the final point
$t_{\mathrm{end}}$ is fixed. (c) Oscillation frequency (dashed dot line) as a
function of temperature $T/T_{c,i}$. The dashed (purple) line corresponds to
$\Delta_{f}$ in Fermi units. The data points are from the experimental results
in Fig. 3(a) of the main text.

##  3 Application of time-dependent BCS theory to the experimental setup

In this Section, we consider the impact of various experimental effects on the
amplitude oscillations and their detection by Bragg spectroscopy. This
includes the following: the presence of a nonzero-duration initial ramp of the
magnetic field (Sec. 3.1), density inhomogeneities due to the trapping
potential (Sec. 3.2), and time averaging over a Bragg pulse of nonzero
duration (Sec. 3.3).

###  3.1 Ramp-rate dependence

To mimic experimental protocols for Feshbach resonances, we include in our
numerics an initial linear magnetic field ramp: $B(t)=B_{i}+Rt$, where
$R=dB/dt$ is the ramp rate. The presence of the ramp delays the start of the
Higgs oscillation. This delay can be quantified using a timescale
$t_{q}=|B_{f}-B_{i}|/R$, which must be compared with the Fermi scale
$\tau_{F}$ and the characteristic timescale of the order parameter dynamics
$t_{\Delta}=1/\Delta$ [21, 22]. When $t_{q}$ becomes much longer than
$t_{\Delta}$, $|\Delta(t)|$ adiabatically evolves to $\Delta_{f}$.

We vary the ramp rate $R$ to study its impact on the early-time dynamics of
$|\Delta(t)|$. Fig. 5 shows how for slower changes of $B$ (smaller $R$) the
amplitude of the oscillations becomes smaller leading to a decrease in the
visibility of the Higgs mode, and the early-time behaviour of $|\Delta|$ has a
more parabolic shape that cannot be fitted using a function $\sim\cos(2\Delta
t)/\sqrt{t}$ (see Eq. (4) in the main text). In the inset of Fig. 5, we plot
the same quantities but we include a delay equal to half ramp time $t_{q}/2$.
Accounting for this delay, acts to align the oscillation phases.

Figure 5: Early-time evolution of the amplitude of the order parameter
$\Delta(t)$ after a quench from $1/(k_{F}a_{i})=-0.18$ to $1/(k_{F}a_{f})=0$
for an infinitely fast ramp (black solid line) and for different ramp time
$t_{q}$ (coloured lines) at zero temperature. We indicate with the bottom
coloured arrows, the values of the ramp time $t_{q}$ corresponding to ramp
rate $R$ (see the text). Inset: Same as (a) but the curves are delayed by half
the ramp time $t_{q}/2$.

###  3.2 Inhomogeneous broadening of the signal

As discussed in Sec. 1.1, during the sample preparation, the atoms are loaded
into an oblate harmonic potential. In this section, we show how this causes an
inhomogeneous broadening of the signal, which can lead to a sharp increase of
the damping coefficient $\gamma$ from $1/2$ to 2 depending on the width of the
Bragg beam.

We begin by writing the external trapping potential as

$V_{\mathrm{ho}}(\mathbf{r})=\frac{1}{2}m\left(\omega_{r}^{2}r^{2}+\omega_{z}^{2}z^{2}\right),$
(14)

where $\omega_{r}$ and $\omega_{z}$ are the trapping angular frequencies in
the radial and axial directions, respectively [23]. In the local density
approximation (LDA), the trapped system can be treated as a collection of
locally uniform subsystems222The local density approximation remains valid
while the harmonic oscillator length, which is the typical length at which
$n(\mathbf{r})$ varies, is much larger than both the healing length
$1/\sqrt{2m\mu}$ and pair correlation length [24] with the latter diverging
near $T_{c}$ as $k_{F}/m\Delta$, marking an obvious breakdown of the
approximation.. This means that at equilibrium the equation of state for the
superfluid can be solved after replacing $\mu$ by a local
$\mu(\mathbf{r})=\mu-V_{\mathrm{ho}}(\mathbf{r})$ and the number equation
becomes

$N=\int d^{3}rn(\mathbf{r}),$ (15)

where $n(\mathbf{r})=2\sum_{\mathbf{k}}n_{\mathbf{k}}(\mathbf{r})$ is the
local density. Consequently, $\Delta(\mathbf{r})$, $1/(k_{F}(\mathbf{r})a)$,
and $k_{\mathrm{B}}T/E_{\mathrm{F}}(\mathbf{r})$ also acquire a dependence on
$\mathbf{r}$ [25, 26, 27, 28]. At unitarity, $1/(k_{F}(\mathbf{r})a)=+\infty$
becomes uniform, but $\Delta$ still depends on $\mathbf{r}$ because
$T/T_{c}(\mathbf{r})$ does.

At zero temperature, for a weakly interacting gas ($|a|\ll
k_{\mathrm{F}}^{-1}$), $n(\mathbf{r})$ is well-approximated by the Thomas-
Fermi distribution:

$n^{(0)}(\mathbf{r})=\frac{8N}{\pi^{2}R_{r}^{2}R_{z}}\left[1-\left(\frac{r}{R_{r}}\right)^{2}-\left(\frac{z}{R_{z}}\right)^{2}\right]^{3/2},$
(16)

where
$R_{i}=\sqrt{2E_{F}^{\mathrm{ho}}/(m\omega_{i}^{2})}=a_{\mathrm{ho}}(24N)^{1/6}\omega_{\mathrm{ho}}/\omega_{i}$
($i=\\{r,z\\}$) are the Thomas Fermi radii,
$\omega_{\mathrm{ho}}=(\omega_{r}^{2}\omega_{z})^{1/3}$ is the geometric mean
of the angular frequencies,
$a_{\mathrm{ho}}=\sqrt{1/(m\omega_{\mathrm{ho}})}=(24N)^{1/6}/k_{F}^{\mathrm{ho}}$
is the oscillator length, and
$k_{F}^{\mathrm{ho}}=\sqrt{2mE_{F}^{\mathrm{ho}}}$ is the Fermi momentum in a
trapped system.

Equation 16 is altered by increasing the scattering length with respect to
$k_{\mathrm{F}}^{\mathrm{ho}}$ as shown in Fig. 6(a), where we compare
$n^{(0)}({\bf r})$ with the density profiles for
$1/(k_{F}^{\mathrm{ho}}a)=-0.18$ and $0$ using the trapping parameters of Sec.
1.1333When $a\to\infty$, the dimensionless parameter $1/(k_{F}(\mathbf{r})a)$
loses spatial dependence, and the equation of state for the Fermi gas is
expected to be universal in terms of density scales and is written in terms of
the Bertsch parameter $\xi_{\mathrm{B}}$ [29, 13]. Therefore, using
$\mu(\mathbf{r})/E_{F}(\mathbf{r})=\xi_{\mathrm{B}}$ and Eq. 16, the density
profile can be approximated by
$n(\mathbf{r})=\xi_{\mathrm{B}}^{-3/2}n^{(0)}(\mathbf{r})$, which well
describes experimental data in the unitary regime when a beyond-mean field
$\xi_{\mathrm{B}}$ is used [30].. In Fig. 6, we have utilized the weighted
radial coordinate

$\rho=\frac{\sqrt{(\omega_{r}r)^{2}+(\omega_{z}z)^{2}}}{\omega_{\mathrm{ho}}},$
(17)

with corresponding (rescaled) Thomas-Fermi readius
$\rho_{\mathrm{TF}}=\sqrt{2E_{\mathrm{F}}^{\mathrm{ho}}/(m\omega_{\mathrm{ho}}^{2})}$.
From Fig. 6(a), it is clear that interactions lead to a narrowing of the
density distribution. It is noteable that Eq. (17) also shows how the case of
a cigar-shape trap can be reduced to the spherical case by homothety.

Figure 6: Particle density $n(\rho)$ in a harmonic trap as a function of the
weighted radial coordinate $\rho$ in Fermi units. (a) At zero temperature for
a non-interacting system (dashed grey), for
$1/k_{\mathrm{F}}^{\mathrm{ho}}a=-0.18$ (green dash dotted), and at unitarity
(black solid). (b) For fixed $1/(k_{\mathrm{F}}^{\mathrm{ho}}a)=-0.18$ and
different temperatures (different colors).

The influence of temperature on the density profile is shown in Fig. 6(b).
Deviations from the zero-temperature profile (black solid curve) begin on the
edges of the cloud and progress towards the center as temperature is
increased.

Within the local density approximation, we use the density distribution
$n(\mathbf{r})$ to calculate the global order parameter by averaging over the
cloud profile. The spatial weight is the product of the local density,
$n(\mathbf{r})$, and local intensity of the probe beam, $I(\mathbf{r})$,

${\bar{\Delta}(t)}=\frac{1}{\bar{n}}\int
d^{3}rn(\mathbf{r})I(\mathbf{r})\Delta(\mathbf{r}),$ (18)

such that

$\frac{1}{\bar{n}}\int d^{3}rn(\mathbf{r})I(\mathbf{r})=1,$ (19)

where $\bar{n}$ is the mean density. The finite width of the probe beam allows
us to consider smaller volumes of the cloud, mimicking the experimental Bragg
scheme described in Sec. 1.3. We approximate the probe beam with a Gaussian
profile, such that

$I(\mathbf{r})=I(0)\,\exp\left(-\frac{|\mathbf{r}|^{2}}{2\lambda^{2}}\right),$
(20)

where $I(0)$ is a normalization constant and $\lambda$ is the effective width
of the beam which is fixed to match the experimental mean density (Eq. 1).

Using Eq. (4) of the main text for the time-evolution of $\Delta$, Eq. (18)
becomes

${\bar{\Delta}(t)}=\frac{1}{\bar{n}}\int
d^{3}rn(\mathbf{r})I(\mathbf{r})\left(\Delta_{\infty}(\mathbf{r})+\frac{\cos(2\Delta(\mathbf{r})t+\phi)}{\sqrt{\nu(\Delta(\mathbf{r}))t}}\right).$
(21)

with the rewriting of the oscillation amplitude $\nu=2\Delta/(A_{\rm
th}(\Delta_{i}-\Delta_{\infty}))^{2}$. The density average of functions
oscillating at frequencies $2\Delta(\mathbf{r})$ leads to a blurring of the
Higgs signal and consequently to a reduced contrast of the oscillations.

####  3.2.1 Qualitative discussion in the limit of a tightly focused beam

When the beam is well focused around the center of the trap, one can
approximate the density using

$n(\rho)=n_{0}-\alpha\rho^{2},$ (22)

with $\alpha=-(1/2)d^{2}n/d\rho^{2}$ (see Eq. (16)). The small variation of
the density around trap center causes a spatial variation of the order
parameter $\Delta(\rho)=\Delta_{0}-\Delta^{\prime}\alpha\rho^{2}$, and of the
magnitude of the Higgs oscillations
$\nu(\rho)=\nu_{0}-\nu^{\prime}\alpha\rho^{2}$ (where
$\Delta^{\prime}=d\Delta/dn$ and $\nu^{\prime}=d\nu/dn$). To simplify the
discussion, we consider here only a probe which respects the cylindrical
symmetry of the trapping potential $I(\mathbf{r})=I(\rho)$, as defined in Eq.
(20). This leads to

${\bar{\Delta}(t)}={{\bar{\Delta}_{\infty}}}+\frac{4\pi
n_{0}}{\bar{n}\sqrt{\nu_{0}t}}\int_{0}^{+\infty}\rho^{2}d\rho
I(\rho)(1-c\rho^{2})\cos(2\Delta_{0}t-\Delta^{\prime}\alpha\rho^{2}t+\phi),$
(23)

where the term proportional to $c=\alpha(1/n_{0}-\nu^{\prime}/2\nu_{0})$ is
included for completeness but becomes negligible when the beam waist tends to
0. In this approximation, the averages of slowly-varying quantities simply
coincide with their value at the center of the trap ($\bar{n}=n_{0}$,
$\bar{\nu}=\nu_{0}$ etc.), and only the average of the oscillatory part
remains.

Performing the radial integral in Eq. (23) for the Gaussian profile Eq. (20)
gives

${\bar{\Delta}(t)}={{\bar{\Delta}_{\infty}}}+\frac{1}{\sqrt{\nu_{0}t}}\text{Re}\left[\mathrm{e}^{2\mathrm{i}\Delta_{0}t+\mathrm{i}\phi}f(t)\right],$
(24)

where the complex-valued blurring function

$f(t)=\left(\frac{1-\mathrm{i}t/t_{\rm
d}}{1+t^{2}/t_{\mathrm{d}}^{2}}\right)^{3/2}$ (25)

depends on the characteristic timescale

$t_{\rm
d}=\frac{1}{4\lambda^{2}\alpha\Delta^{\prime}}\approx\frac{1}{\Delta(\lambda)-\Delta_{0}}$
(26)

Assuming that this timescale is large ($\Delta_{0}t_{\rm d}\gg 1$), the Higgs
oscillations will still display the expected $1/\sqrt{t}$ attenuation law at
short times ($t\ll t_{\rm d}$). However, beginning at $t\approx t_{\rm d}$ the
decay is quickened and the phase shifts from the original $\phi=\pi/4$ phase.
At later times ($t\gg t_{\rm d}$), one finds instead a $1/t^{2}$ attenuation
law:

${\bar{\Delta}(t)}\underset{t\to+\infty}{\simeq}{{\bar{\Delta}_{\infty}}}-\frac{t_{\rm
d}^{3/2}}{\sqrt{\nu_{0}}t^{2}}\sin 2\Delta_{0}t.$ (27)

This shows that inhomogeneous broadening can seriously impair the observation
of the power-law damping exponent $\gamma=1/2$ or $3/2$.

####  3.2.2 Numerical solutions in the general case

More generally, we solve numerically444 To perform the integration in Eq. 18,
we sample $n(\mathbf{r})$ using a “accept-reject method” [31], which generates
a set of random numbers distributed according to the initial distribution:
$n(\mathbf{r})\rightarrow\\{n_{i},$ with $i=1,\cdots,N_{\mathrm{p}}\\}$. The
average on Eq. 18 can be then performed using
$\bar{\Delta}(t)=\frac{1}{\sum_{i}^{N_{\mathrm{p}}}I(n_{i})}\sum_{i=1}^{N_{\mathrm{p}}}\,I(n_{i})\Delta(n_{i},t),$
(28) where $N_{\mathrm{p}}$ is the total of number of extracted values
$n_{i}$.  Eq. (18) for different values of $\lambda$ and, consequently, for
different $\bar{n}$. Fig. 7 shows the evolution of the superfluid order
parameter at zero temperature after a quench from
$1/k_{F}^{\mathrm{ho}}a=-0.18$ to unitarity after averaging over a probe beam
of width $\lambda=0.06\rho_{\rm TF}$ such that $\bar{n}=0.91n_{0}$, compared
to the case where the probe beam covers the entire volume. We observe that
this latter case yields strongly-damped oscillations, which would be difficult
to observe. Instead, the experimentally-relevant average over a smaller volume
leads to distinct oscillations. Moreover, because $\lambda=0.06\rho_{\rm TF}$
is small compared to the effective Thomas-Fermi radius $\rho_{\mathrm{TF}}$,
the characteristic time $t_{d}$ is much larger than the oscillation period
$2\pi/\omega_{\mathrm{H}}$, and therefore, for times $t<t_{\rm d}$, we do not
expect the attenuation law $1/t^{2}$ (as described by Eq. (27)) to play a role
within the experimentally relevant time window.

Figure 7: Time evolution of the order parameter after averaging over the total
density profile of the trap (green solid line) and over a reduced volume with
corresponding $\bar{n}=0.91n_{0}$ (blue dashed line) at zero temperature.
Black dotted and yellow dashed dotted lines indicate the power law curves
$1/\sqrt{t}$, which is found at short times, and $1/t^{2}$ from Eq. 27, which
characterizes longer times than the maximum time in the figure ($t\gg
t_{\mathrm{d}}$).

We address also the impact of inhomogeneous broadening on the damping law in
the time region $0<t<t_{\mathrm{d}}$. Fig. 8 shows trends in the local
oscillation maxima $\bar{\Delta}(t_{n})$ at discrete times $t_{n}$ measured in
units of the oscillation frequency $\omega_{\rm H}$ and at temperature
$T=0.66T_{\mathrm{c,i}}$. We fit the first few oscillations to a power-law
damping, thereby mimicking the analysis performed on the experimental data.
For smaller values of $\lambda$ compared to $\rho_{\mathrm{TF}}$, the power
law decay $1/\sqrt{t}$ is still a good fit for the short-time behavior of
$\bar{\Delta}$, as shown by the case with $\bar{n}=0.98n_{0}$ in Fig. 8. For
larger values of $\lambda$, we find an intermediate behavior, characterized by
a power-law decay $1/t^{\beta}$ with $0.5<\beta<2.0$.

Figure 8: Data points: Maximum values of $(\bar{\Delta}(t)/\Delta_{\infty}-1)$
for different $\bar{n}$ or equivalently $\lambda$ at $T=0.66T_{\mathrm{c,i}}$.
The solid black line represents the homogeneous power-law behavior
$1/\sqrt{t}$ (see Eq. (4) in the main text). Dashed lines: Fitting functions
$1/t^{\beta}$ for times $t<t_{d}$. Gray solid line: power law $1/t^{2}$
recovered only in the case with lower $\bar{n}$.

###  3.3 Averaging over the Bragg time

Experimental measurements are taken using two-photon Bragg scattering as
illustrated in Fig. 1 of the main text. To avoid over-broadening the sharp
features of the response function, the duration $t_{\mathrm{B}}=50\mu s$ of
the Bragg pulse should be as long as possible. However, in an out-of-
equilibrium system, using pulses that are too long allows the system to evolve
significantly during the probe time violating the timescale separation
assumptions underlying Eq. (5) of the main text and possibly blurring the
oscillatory signal. Explicitly, the observable $\Delta$ reported at the time
$t$ is actually the result of a time-average:

$\braket{\Delta(t)}_{\mathrm{Bragg}}=\frac{1}{t_{\mathrm{B}}}\int_{t-t_{\mathrm{B}}/2}^{t+t_{\mathrm{B}}/2}\Delta(t^{\prime})dt^{\prime}.$
(29)

Fig. 9 shows the comparison between the time-evolution of (spatially-averaged)
$\Delta$ for $t_{B}=0$ and $t_{B}=50\mu s$. Although the time-averaging
reduces the contrast of the oscillations by reducing their the magnitude, it
is less critical than the spatial averaging as it does not affect the
attenuation law.

Figure 9: Comparison between the time and space averaged order parameter
($\langle\bar{\Delta}(t)\rangle$, yellow solid line) and the same quantity
only averaged over space ($\bar{\Delta}(t)$, blue dashed line). The time
average is performed over the Bragg time $t_{\mathrm{B}}$ (see Eq. (29)) and
the density average corresponds to the case $\bar{n}=0.91n_{0}$ (see Eq.
(18)). The temperature is $T=0.66T_{\mathrm{c,i}}$.

##  4 Comparing theoretical predictions and experimental measurements

The comparisons with experiment shown in Figs. 2(b) and 3(b) of the main text
use theoretical predictions for $\braket{\bar{\Delta}(t)}$ that combine both
the spatial and temporal averages discussed in Secs. 3.2 and 3.3 and the
dephasing introduced by the finite ramp time $t_{q}$ discussed in Sec. 3.1.
From Eq. (5) of the main text, we write the experimental signal as
${S(t)}\approx\alpha+\beta{\Delta(t)}$, where $\alpha$ and $\beta$ depend on
the sensitivity with an explicit expression of those quantities given below.
This dependency is removed by taking the ratio

$\frac{{S(t)}-S_{\infty}}{S_{f}-S_{i}}\sim\frac{\braket{\bar{\Delta}(t)}-\Delta_{\infty}}{\Delta_{f}-\Delta_{i}}.$
(30)

To avoid the well-known overestimation of BCS theory of the critical
temperature $T_{c}$ and superfluid gap $\Delta$, we compare theoretical
predictions and experiment signal as a function of $\omega_{\rm H}t$ and
$T/T_{c,i}$ using experimental values of $\omega_{\rm H}/\epsilon_{F}$, fitted
from the bare experimental signal, and $T_{c,i}/T_{F}$, estimated as the
temperature at which Higgs oscillations are no longer observed.

###  4.1 Magnitude averaged over different time windows

To quantify the agreement between experiment and theory, we compute at each
temperature the root-mean-square of the signal (30):

$A_{\rm
rms}=\sqrt{\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}\mathrm{d}t\left(\frac{S(t)-S_{\infty}}{S_{f}-S_{i}}\right)^{2}},$
(31)

As discussed in the main text, the choice of the averaging time window
$[t_{1},t_{2}]$ significantly affects the gap between measurements and
theoretical predictions, the best agreement being found at early times. In
Fig. 10, we show that the agreement is better from $[t_{1},t_{2}]=[\tau_{\rm
H},10\tau_{\rm H}]$ than in the window $[3\tau_{\rm H},20\tau_{\rm H}]$ used
in the main text.

Figure 10: Comparison between the theoretical (diamonds and x’s) and
experimental (squares and plus) averaged magnitudes of the oscillations,
$A_{\mathrm{rms}}$, in two different time ranges in units of $\tau_{\rm H}$.
Shown as dashed-dotted lines is the magnitude of the oscillations obtain by
integrating Eq. (12), that is, in the linear regime and without any of the
experimental effects discussed in Sec. 3.

###  4.2 Comparisons including the sensitivity of the Bragg pulse

In this subsection, we predict the experimental Bragg signal from our
calculation of the density-density response within the Random Phase
Approximation (RPA), that is, without scaling out the sensitivity. As
discussed in the main text (see Eq. (5)), the Bragg signal is sensitive to the
variations of $\Delta(t)$ through the dimensionless sensitivity

$\sigma=\frac{\Delta}{\chi_{nn}}\left(\frac{\partial\chi_{nn}}{\partial\Delta}\right)_{n,T},$
(32)

This thermodynamic quantity is computed by comparing equilibrium states having
the same equilibrium density and temperature but slightly different values of
the order parameter $\Delta$ (or equivalently, slightly different values of
the scattering length $a$). The calculation of $\chi_{nn}$ within RPA is a
standard problem [32, 33]. In practice we have used the explicit expressions
of Ref. [34] (see in particular Eqs. (36) and (46) therein) to evaluate
numerically $\chi_{nn}$ as a function of the excitation frequency $\omega$ and
wavevector $q$.

To illustrate how the density-density response is sensitive to $\Delta$, Fig.
11 shows $\chi_{nn}(\omega)$ in the regime of large excitation wavevector
($q=4k_{F}$) used in the experiment. The feature which characterizes the
superfluid phase is a sharp edge at the pair-breaking threshold:

$\omega_{\rm
th}=2\sqrt{\Delta^{2}+\left(\frac{q^{2}}{8m}-\mu\right)^{2}}\quad\text{ (for
}{q^{2}}/{8m}>\mu).$ (33)

As this threshold location changes with $\Delta$, the density response
$\chi_{nn}$ for fixed $\omega/\epsilon_{F}$ varies sharply. However, as
temperature increases towards $T_{c}$, the spectral weight of the edge drops
in favor of the broad atomic scattering continuum, characteristic of the ideal
(normal) Fermi gas. This fading of the pair-breaking edge consequently reduces
the sensitivity $\sigma$.

Figure 11:  The imaginary part of the density-density response function
$\chi_{nn}$ in density units as a function of the excitation frequency
$\omega$ at fixed excitation wave number $q=4k_{F}$ and temperature varying
from $T=0$ (blue curves), $T=0.9T_{c}$ (purple curves) to $T=0.97T_{c}$ (red
curves). The threshold of the pair-breaking continuum (black dotted vertical
lines, given by Eq. (33)) is sensitive to the value of $\Delta$ and thus
depends slightly on temperature, tending to $2({q^{2}}/{8m}-\mu)$ in the limit
$T\to T_{c}$ where $\Delta\to 0$. For comparison, we show the kinetic energy
of atoms (mass $m$) and pairs (mass $2m$) in vacuum (blue vertical dotted
lines, independent of $T$). At the pair-breaking threshold, $\chi_{nn}$ has a
sharp peak indicating resonance with the condensed pairs. As temperature
increases, the number of pairs decreases in favor of unpaired atoms. This
results in a decreased spectral weight of the peak at $\omega_{\rm th}$ and
correspondingly in a growth of a broad single-atom scattering continuum
centered around $\omega=q^{2}/2m$. This broad peak is due to the imaginary
part of the normal density-density propagator $\Pi_{nn}$ (see Eq. (46) in
[34]) which is shown by the red crosses at the highest temperature
$T=0.97T_{c}$. These curves are drawn at unitarity ($1/|a|=0$) which
corresponds at $T=0$ and with the mean-field equation-of-state to
$\mu/\Delta\simeq 0.86$ and $\epsilon_{F}/\Delta\simeq 1.46$. At $T=T_{c}$,
this corresponds to $\mu(T_{c})/T_{c}=1.50448$,
$\epsilon_{F}(T_{c})/T_{c}\simeq 0.89$.

In theory, the vertical tangent of $\omega\mapsto\chi_{nn}(\omega,\Delta)$ in
$\omega_{\rm th}^{+}$, causes a divergence555To understand this divergence and
also to obtain rigorous numerical results for the blue curve in Fig. 12, we
write $\chi_{nn}(\omega,\Delta)=\Theta(\omega-\omega_{\rm
th}(\Delta))f(\omega-\omega_{\rm th}(\Delta),\Delta)$ where the function
$\delta\omega\mapsto f(\delta\omega,\Delta)$ is defined on $[0,+\infty[$.
Taking the derivative with respect to $\Delta$ (and using $f(0,\Delta)=0$), we
have $\partial\chi_{nn}/\partial\Delta=\partial
f/\partial\Delta-(\mathrm{d}\omega_{\rm th}/\mathrm{d}\Delta)\times(\partial
f/\partial\delta\omega)$. The vertical tangent of $\chi_{nn}$ in $\omega_{\rm
th}^{+}$ implies a divergence of $\partial f/\partial\delta\omega$ in
$\delta\omega=0^{+}$, and therefore a divergence of the sensitivity. of the
sensitivity at the edge (see the blue curve in Fig. 12), which should favor
the detection of Higgs oscillations. In practice, the finite duration of the
Bragg pulse limits the spectral resolution on $\chi_{nn}$, thereby limiting
the maximal accessible value of the sensitivity to roughly $\sigma\lesssim 2$
(red curve in Fig. 12). This is obtained by broadening $\chi_{nn}$ according
to the convolution formula

$\chi_{nn}^{\rm
br}(\omega_{0})=\int_{-\infty}^{+\infty}\mathrm{d}\omega\chi_{nn}(\omega)f_{\delta\omega}(\omega_{0}-\omega)$
(34)

where $f_{\delta\omega}(\omega)$ is a broadening function parametrized by its
energy-width $\delta\omega$. In the comparison performed in Fig. 13 we have
taken $\delta\omega/\epsilon_{F}=\tau_{\rm F}/t_{\rm B}\simeq 0.55$, using the
Fermi time $\tau_{\rm F}\simeq 27.4402\mu s$ associated to the density in the
trap center at unitarity, and a Gaussian broadening profile:

$f_{\delta\omega}(\omega)=\frac{1}{\sqrt{2\pi}\delta\omega}\mathrm{e}^{-\omega^{2}/2(\delta\omega)^{2}}.$
(35)

However the choice of a Gaussian broadening profile is not crucial for this
discussion. Compared to our previous experimental schemes [35], the spectral
width is relatively larger here, due to the necessity of maintaining a Bragg
duration much smaller than $2\pi/\omega_{\rm H}$. As visible in Fig. 12, the
selected excitation frequency $\omega=\hbar^{2}q^{2}/4m$ is fairly optimal for
the broadened sensitivity whereas $\omega=\omega_{\rm th}$ is preferable in
the absence of broadening.

Figure 12:  The bare (blue curve) and broadened (red curve) sensitivities as a
function of the excitation frequency for $q=4k_{F}$, $T=0$, and $1/k_{F}a=0$.

Using the sensitivity (32), we directly predict the scaled Bragg response
$S(t)/S_{\infty}$ from

$\frac{\braket{\bar{\chi}_{nn}^{\prime\prime}(t)}}{\bar{\chi}_{nn}^{\prime\prime}(+\infty)}=1+\sigma_{f}\frac{\Delta_{\infty}}{\Delta_{f}}\left(\frac{\braket{\bar{\Delta}(t)}}{\Delta_{\infty}}-1\right),$
(36)

where we have assumed that the width of the Bragg laser $\lambda$ is small
enough such that the dependence of the sensitivity
$\sigma_{f}\equiv\sigma(\Delta_{f})$ on $\mathbf{r}$ can be neglected. Note
that the linearization of $\chi_{nn}$ performed in Eq. (36) is valid only in
the small amplitude approximation, hence at low $T/T_{c}$.

In Fig. 13, we compare the theoretical curve Eq. (36) to the experimental
measurements of $S(t)/S_{\infty}$. We consider the lowest temperature case,
$T=0.10T_{F}$, for which $\delta\Delta$ is small enough with respect to
$\Delta_{f}$ to justify the use of the small-amplitude regime (we have
$|\Delta_{\infty}-\Delta_{f}|\simeq 0.11\Delta_{f}$). At the latest times
considered in the experiments, we find a small deviation of our theoretical
prediction from the homogeneous $1/t^{1/2}$ power-law damping, with $\gamma$
increasing to $0.645\pm 0.005$. This deviation could be seen if we exclude the
first peak of the oscillation from the theoretical fitting window,
consistently with our analysis in Sec. 3.2 (compare with Fig. 8 and the inset
of Fig. 13). Instead, using the experimental fitting window, which also
includes the first peak, this deviation is nearly compensated, as discussed in
the main text.

Figure 13: Comparison between the experimental Bragg signal $S(t)/S_{\infty}$
and its theoretical prediction Eq. (36). The theoretical curve has been
obtained at fixed $T=0.66T_{c,i}$ by including a linear ramp in $50\mu s$ with
$t_{h}\ =t-t_{q}/2$ corresponding to the theoretical delay found in Sec. 3.1,
the inhomogeneous broadening, the Bragg time averaging, and the Fourier
broadening. The experimental points correspond to $T=0.10T_{F}$ and the mean
density is $\bar{n}=(0.955\pm 0.018)n_{0}$ (see main text). We show the
deviation from the $1/t^{0.5}$ power-law damping (grey dashed line). The
power-law $1/t^{0.645\pm 0.005}$ (dashed-dotted purple curve) has been
obtained by fitting the theoretical curve in the time range
$9\leq\omega_{H}t_{h}\leq 23$ (from the second peak to the last experimental
point). Inset: Damping coefficient as a function of the starting time
$t_{\mathrm{st}}$ of the fitting window used for the theoretical curve.

## References

  * [1] Smith, N. L., Heathcote, W. H., Hechenblaikner, G., Nugent, E. & Foot, C. J.  Quasi-2d confinement of a BEC in a combined optical and magnetic potential.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 38, 223–235 (2005).
  * [2] Rath, S. P. _et al._ Equilibrium state of a trapped two-dimensional bose gas.  _Phys. Rev. A_ 82, 013609 (2010).
  * [3] Dyke, P. _et al._ Criteria for two-dimensional kinematics in an interacting fermi gas.  _Phys. Rev. A_ 93, 011603(R) (2016).
  * [4] Zürn, G. _et al._ Precise characterization of ${}^{6}\mathrm{Li}$ feshbach resonances using trap-sideband-resolved rf spectroscopy of weakly bound molecules.  _Phys. Rev. Lett._ 110, 135301 (2013).
  * [5] Dyke, P. _et al._ Dynamics of a fermi gas quenched to unitarity.  _Phys. Rev. Lett._ 127, 100405 (2021).
  * [6] Ku, M. J. H., Sommer, A. T., Cheuk, L. W. & Zwierlein, M. W.  Revealing the superfluid lambda transition in the universal thermodynamics of a unitary fermi gas.  _Science_ 335, 563–567 (2012).
  * [7] Hoinka, S. _et al._ Goldstone mode and pair-breaking excitations in atomic fermi superfluids.  _Nature Physics_ 13, 943–946 (2017).
  * [8] Carcy, C. _et al._ Contact and sum rules in a near-uniform fermi gas at unitarity.  _Phys. Rev. Lett._ 122, 203401 (2019).
  * [9] Lingham, M. G., Fenech, K., Hoinka, S. & Vale, C. J.  Local observation of pair condensation in a fermi gas at unitarity.  _Phys. Rev. Lett._ 112, 100404 (2014).
  * [10] Veeravalli, G., Kuhnle, E., Dyke, P. & Vale, C. J.  Bragg spectroscopy of a strongly interacting fermi gas.  _Phys. Rev. Lett._ 101, 250403 (2008).
  * [11] Regal, C. A., Greiner, M. & Jin, D. S.  Observation of resonance condensation of fermionic atom pairs.  _Phys. Rev. Lett._ 92, 040403 (2004).
  * [12] Chin, C., Grimm, R., Julienne, P. & Tiesinga, E.  Feshbach resonances in ultracold gases.  _Rev. Mod. Phys._ 82, 1225–1286 (2010).
  * [13] Strinati, G. C., Pieri, P., Röpke, G., Schuck, P. & Urban, M.  The bcs–bec crossover: From ultra-cold fermi gases to nuclear systems.  _Physics Reports_ 738, 1 – 76 (2018).
  * [14] Volkov, A. & Kogan, S. M.  Collisionless relaxation of the energy gap in superconductors.  _Zh. Eksp. Teor. Fiz_ 65, 2038 (1973).
  * [15] Gurarie, V.  Nonequilibrium dynamics of weakly and strongly paired superconductors.  _Phys. Rev. Lett._ 103, 075301 (2009).
  * [16] Yuzbashyan, E. A., Dzero, M., Gurarie, V. & Foster, M. S.  Quantum quench phase diagrams of an $s$-wave BCS-BEC condensate.  _Phys. Rev. A_ 91, 033628 (2015).
  * [17] Kurkjian, H., Klimin, S. N., Tempere, J. & Castin, Y.  Pair-Breaking Collective Branch in BCS Superconductors and Superfluid Fermi Gases.  _Phys. Rev. Lett._ 122, 093403 (2019).
  * [18] Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T. _et al._ _Numerical Recipes_ , vol. 3 (Cambridge University Press, Cambridge, 1989).
  * [19] Holland, M., Kokkelmans, S. J. J. M. F., Chiofalo, M. L. & Walser, R.  Resonance superfluidity in a quantum degenerate fermi gas.  _Phys. Rev. Lett._ 87, 120406 (2001).
  * [20] Kokkelmans, S. J. J. M. F., Milstein, J. N., Chiofalo, M. L., Walser, R. & Holland, M. J.  Resonance superfluidity: Renormalization of resonance scattering theory.  _Phys. Rev. A_ 65, 053617 (2002).
  * [21] Barankov, R. A., Levitov, L. S. & Spivak, B. Z.  Collective rabi oscillations and solitons in a time-dependent bcs pairing problem.  _Phys. Rev. Lett._ 93, 160401 (2004).
  * [22] Scott, R. G., Dalfovo, F., Pitaevskii, L. P. & Stringari, S.  Rapid ramps across the bec-bcs crossover: A route to measuring the superfluid gap.  _Phys. Rev. A_ 86, 053604 (2012).
  * [23] Pitaevskii, L. & Stringari, S.  _Bose-Einstein Condensation and Superfluidity_ , vol. 164 (Oxford University Press, 2016).
  * [24] Marini, M., Pistolesi, F. & Strinati, G.  Evolution from bcs superconductivity to bose condensation: analytic results for the crossover in three dimensions.  _The European Physical Journal B_ 1, 151–159 (1998).
  * [25] Perali, A., Pieri, P. & Strinati, G. C.  Shrinking of a condensed fermionic cloud in a trap approaching the bose-einstein condensation limit.  _Phys. Rev. A_ 68, 031601 (2003).
  * [26] Salasnich, L., Manini, N. & Parola, A.  Condensate fraction of a fermi gas in the bcs-bec crossover.  _Phys. Rev. A_ 72, 023621 (2005).
  * [27] Chiofalo, M. L., Kokkelmans, S. J. J. M. F., Milstein, J. N. & Holland, M. J.  Signatures of resonance superfluidity in a quantum fermi gas.  _Phys. Rev. Lett._ 88, 090402 (2002).
  * [28] Haussmann, R. & Zwerger, W.  Thermodynamics of a trapped unitary fermi gas.  _Phys. Rev. A_ 78, 063602 (2008).
  * [29] Ho, T.-L.  Universal thermodynamics of degenerate quantum gases in the unitarity limit.  _Phys. Rev. Lett._ 92, 090402 (2004).
  * [30] Perali, A., Pieri, P. & Strinati, G. C.  Quantitative comparison between theoretical predictions and experimental results for the bcs-bec crossover.  _Phys. Rev. Lett._ 93, 100404 (2004).
  * [31] Krauth, W.  _Statistical Mechanics Algorithms and Computations_ (Oxford University Press, Oxford, 2006).
  * [32] Bruun, G. M. & Mottelson, B. R.  Low energy collective modes of a superfluid trapped atomic fermi gas.  _Phys. Rev. Lett._ 87, 270403 (2001).
  * [33] Minguzzi, A., Ferrari, G. & Castin, Y.  Dynamic structure factor of a superfluid fermi gas.  _The European Physical Journal D_ 17, 49–55 (2001).
  * [34] Kurkjian, H., Tempere, J. & Klimin, S. N.  Linear response of a superfluid fermi gas inside its pair-breaking continuum.  _Scientific Reports_ 10(1), 11591 (2020).
  * [35] Hoinka, S. _et al._ Goldstone mode and pair-breaking excitations in atomic Fermi superfluids.  _Nature Physics_ 13, 943–946 (2017).