Second sound in the crossover from the Bose-Einstein condensate to the Bardeen-Cooper-Schrieffer superfluid

Second sound is an entropy wave which propagates in the superfluid component of a quantum liquid. Because it is an entropy wave, it probes the thermodynamic properties of the quantum liquid. Here, we study second sound propagation for a large range of interaction strengths within the crossover between a Bose-Einstein condensate (BEC) and the Bardeen-Cooper-Schrieffer (BCS) superfluid, extending previous work at unitarity. In particular, we investigate the strongly-interacting regime where currently theoretical predictions only exist in terms of an interpolation in the crossover. Working with a quantum gas of ultracold fermionic 6Li atoms with tunable interactions, we show that the second sound speed varies only slightly in the crossover regime. By varying the excitation procedure, we gain deeper insight on sound propagation. We compare our measurement results with classical-field simulations, which help with the interpretation of our experiments.


Supplementary Note 1: Temperatures to the measurements in
and Fig. 2 In this section we present the temperatures to the measurements shown in Fig. 1 and  Table 1). We determine the temperatures by tting a second order virial expansion of the density distribution at the wings of the cloud 1 . To compare the absolute temperature with T c for various interaction strengths we use values for T c as shown in Supp.

Figure 1.
T c is not precisely known yet in the strongly interacting regime. In the limit of the BEC regime the BEC mean-eld model should give accurate values for critical temperature.
Closer towards the resonance we expect the diagrammatic t-matrix calculation to provide quite good values 2 . For the range in between (0.5 < (k F a) −1 < 3) we linearly interpolate between both T c curves. For the measurements on the resonance we compared our result with the temperature determined using the equation of state from Ref. 3 . We nd reasonable agreement between the temperatures obtained from the two approaches with deviations on the order of 5 − 10%.

Supplementary Note 2: C-eld simulation method
Here we present our simulation method that is used to simulate sound mode dynamics in a condensate of 6 Li molecules on the BEC side. The system is described by the Hamiltonian ψ andψ † are the bosonic annihilation and creation operator, respectively. The 3D interaction parameter is given by g = 4πa dd 2 /M , where a dd is the dimer-dimer scattering length and M the dimer mass. The external potential V (r) represents the cigar-shaped trap V trap (r) = M (ω 2 x x 2 + ω 2 r r 2 )/2. ω x and ω r are the axial and radial trapping frequencies, respectively. r = (y 2 + z 2 ) 1/2 is the radial coordinate.
To perform numerical simulations we discretize space with the lattice of 180×35×35 sites and the discretization length l = 0.5 µm, where l is chosen to be smaller than or comparable to the healing length ξ > 0.5 µm and the thermal de Broglie wavelength λ dB ≈ 1 − 1.5 µm.
Since λ dB determines the scale for thermal uctuations, the associated thermal energy should always be below the cuto energy introduced by the discretization length. We also note that in the opposite limit l > ξ, λ dB the simulation method would be inadequate to capture smalldistance excitations, such as vortices. In our c-eld representation we replace in Eq. 1 and in the equations of motion the operatorsψ by complex numbers ψ, see Ref. 4 . We sample the initial states in a grand-canonical ensemble of temperature T and chemical potential µ via a classical Metropolis algorithm. We obtain the time evolution of ψ(t) using the classical equations of motion. As our key observable, we calculate the density n(r, t) = |ψ(r, t)| 2 and average it over the thermal ensemble. For our simulations we use the trapping frequencies (ω x , ω r ) = 2π × (70 Hz, 780 Hz) that are higher than the experimental trap values. This is because the size of the simulation lattice is needed to be small in order to have a reasonable calculation time. We show below that this larger value of the trapping frequency does not aect our results of the sound velocity, which is determined from the sound propagation near the trap center. We choose the scattering length a dd and the trap central density To excite sound modes we add the perturbation H ex (t) = dr V (r, t)n(r), where n(r) is the density at the location r = (x, y, z). The excitation potential V (r, t) is given by where V 0 (t) is the time-dependent strength and σ is the width. The locations x 0 , z 0 are chosen to be the trap center. We excite sound modes following the scheme used in the experiment, where σ and V 0 are chosen such that the changes in the local density due to the excitation potential are consistent with the experiment. We calculate the density prolē n ex (x, t), which is integrated along the radial direction. For sound propagation we examine ∆n(x, t) = n ex (x, t)−n(x) /n(0), wheren(x) is the density prole of the unperturbed cloud integrated in the radial direction andn(0) is the maximum density. As we show in the main text, the time evolution of ∆n(x, t) displays excitation of second sound, which is identied by a vanishing sound velocity at R TF . We note that solitonic excitations are not expected as they involve a steep change of the local phase, whereas our excitation protocol modies the local density. Furthermore the size of solitonic wave features is on the order of the healing length of 0.5 µm whereas the rst and second sound features we observe have a width of a few tens of micrometers. We t ∆n(x, t) with a Gaussian to determine the second sound velocity u 2 at the trap center. We note that within the range of T /T c = 0.5−0.7 the velocity u 2 changes only negligibly with temperature compared to the experimental errorbars.
A. Low versus strong transverse trapping frequency To examine whether a stronger connement in the transverse direction aects the sim-ulation result of the sound velocity, we choose a lower transverse trapping frequency of ω r = 2π × 546 Hz and compare its result with that of ω r = 2π × 780 Hz, while we use the same axial trapping frequency ω x = 2π × 70 Hz and the same scattering length a dd = 840a 0 .
The simulated cloud consists of N ≈ 66, 000 and 55, 000 molecules for the systems of low and high trapping frequency, respectively. To excite sound modes we use the excitation frequency ω ex /ω r = 0.61 and a half-cycle modulation, which are the same as in the case of high frequency simulation. In Supp. Fig. 2  appears which diuses after a short propagation time. This happens due to the fact that this rst sound wave is created after the dark second sound wave but propagates with higher velocity within the superuid region. When it crosses the second sound wave it fades out in the simulation because mixing of rst and second sound leads to diusion. In addition, second sound signal seems to wash out for increasing temperatures. In fact, at a temperature near the transition temperature second sound becomes a diusive sound mode as discussed in Ref. 11 . This seems to also increase the diusion of rst sound, leading e.g. to the suppression of the additional bright rst sound in Supp. Fig. 3 d.
Although predicted in the simulations the diusion of rst sound is not observed in the experiment (see Fig. 4) as its wave is still clearly visible even beyond the Thomas-Fermi radius. The discrepancy between experiment and theory could be due to the inherent discrete nature of the simulation method, where the variation in the density and the phase is not smooth on the scale of the discretization length, causing an additional dispersion of the sound wave.
To understand the long time behavior we show in the lower row of Supp. Fig. 3  This reection is strongly visible for the lowest temperature of T = 120 nK in Supp. Fig. 3 a (upper row).

Supplementary Note 3: Analytic description of the sound modes
In the following we present an analytic description of rst and second sound based on the two-uid hydrodynamic model for a uniform gas. The total density n of the gas is a sum of the superuid n s and normal uid density n n . The rst and second sound mode squared velocities are given by 5 where c 2 T = 1/M (∂p/∂n) T and c 2 2 = n s s 2 T /(n n c V ) representing the isothermal and entropic sound velocities, respectively. p is the pressure, s the entropy per unit mass, T the temperature, and c V = T (∂s/∂T ) n the heat capacity per unit mass. The quant- Here, rst and second sound can be described as a pressure and entropy wave, respectively.
To determine the second sound velocity u 2 , we calculate the entropy and the normal uid density dened as and respectively 5 is the excitation energy and k the wavevector. The upper and lower sign correspond to a Bose and Fermi gas, respectively.

A. BEC
We use the Bogoliubov theory, valid in the dilute limit, to analyze the regime k B T < gn, where gn is the mean-eld energy. The Bogoliubov spectrum is given by E k = k ( k + 2gn), where k = 2 k 2 /(2M ) is the free-particle spectrum. M is the molecular mass. To examine the decoupled modes in Eq. 4 we approximate E k by the linear spectrum E k ≈ ck, where c = gn/M is the Bogoliubov sound velocity. We obtain the entropy and the normal uid density, respectively, Here, This result is only valid at zero temperature, see Supp. Fig. 4a, where we show the full numerical solutions of Eq. 3 using the Bogoliubov description.
For k B T > gn instead we make use of a thermal gas description to determine s, c V , and n n , which are given by s = 2.568k B n n /(2M n), c V = 3s/2, and n n = n(T /T c ) 3/2 , respectively 5 .
In this regime, solving Eq. 3 the sound velocities read, u 2 is proportional to n s /n and can be approximated by u 2 = 1 − (T /T c ) 3/2 gn/M (see Supp. Fig. 4a).

Sound amplitudes
Besides the sound velocity, our analytic description can be used to determine the amplitudes of the propagating sound modes, described as 6 where δñ(x, t) is the density variation created by the excitation potential. δñ(x ± u 1/2 t) represent wave packets of rst and second sound with weights W 1/2 . The relative weight is given by 6 We determine W 2 /W 1 by numerically solving Eq. 3 for the regimes k B T < gn and k B T > gn using the Bogoliubov and thermal gas description, respectively. We show these results in Supp. Fig. 4b. The Bogoliubov description of the weight works only for k B T gn. We note that at higher temperatures terms beyond Bogoliubov are needed to account for the thermal damping of the modes. The Bogoliubov description thus leads to an overestimation of the weight at high temperatures. For temperatures above the mean-eld energy the weight is described by the thermal gas description, which we use to estimate the relative weight of the two modes in the main text. Please note that the thermal description gives unphysical solutions for k B T /gn → 1. In the experiment presented in g. 3 of the main text k B T /gn ≈ 0.8 in the central region and therefore the Bogoliubov description should be valid.

B. BCS
A condensate of an interacting Fermi gas is described by the BCS spectrum E k = ξ 2 k + ∆ 2 , with ξ k = 2 k 2 /(2m) − µ, where µ is the chemical potential and ∆(T ) the gap. At low k B T ∆, we use µ ≈ E F and expand ξ k near the Fermi surface, i.e.
. The entropy in Eq. 5 results in with which is the gap at zero temperature 8 . With Eq. 12 we determine s = S/(mN tot ) and c V .
The normal uid density in Eq. 6 gives Using s, c V , and n n in Eq. 4 we obtain the second sound velocity which is valid for T < T c . The BCS critical temperature is given by k B T c = (γ/π)∆ 0 = 0.567∆ 0 , which depends on the interaction parameter (k F a) −1 . We show in the main text the result u 2 at various interactions on the BCS side (see Fig. 2). u 2 vanishes at zero temperature contrary to the BEC superuids. We note that this result is consistent with Ref. 9 .

Supplementary Note 4: BEC mean-eld model
To estimate the density distribution of a partially Bose condensed cloud in the BEC regime we carry out a self-consistent calculation where the condensate phase is treated within the Thomas-Fermi approximation and for the normal phase we use a standard thermodynamical approach. Specically, we solve the following set of coupled equations 10 Here, λ dB is the thermal de Broglie wavelength, g = 4π 2 a dd /M is the coupling constant, T is the temperature and V ext (r) is the external potential consisting of the harmonic trapping potential and the repulsive potential of the excitation beam, µ s and µ n are the chemical potentials of the superuid and the normal uid part, respectively. For the calculation we set µ n = min[V ext (r) + 2gn s (r) + 2gn n (r)] which ensures that the normal gas reaches the critical density n n,crit = Li 3/2 (1)/λ 3 dB at the Thomas-Fermi radius. This way, the number of normal uid atoms is xed. µ s is chosen such that the total atom number matches the experimental value.
Equation 16 represents the Thomas-Fermi approximation where we take into account the repulsive mean-eld potential of the normal uid part. Equation 17 is the density distribution of a thermal bosonic cloud, again including the additional mean-eld potential produced by the atoms. By self-consistently solving the coupled equations we obtain the density distributions of the superuid and the normal uid gas as shown in Supp. Fig. 5.
The repulsive excitation beam pushes the atoms away from the trap center which creates a density prole with two peaks of the same height. Interestingly in our self-consistent calculations we nd that the peak density with and without the excitation beam is almost the same. This holds both for the line density and the 3D density. This allows for extracing the peak density in our experiments, when the excitation laser is present, from reference absorption images of an unperturbed cloud when no excitation laser is present. To do this, we apply the inverse Abel transformation to the reference absorption images to reconstruct the 3D density prole. Note that this transformation is in generally valid only for rotationally symmetric clouds, which is the case when no exciation beam is present. We have veried that we obtain the same density using our self-consistent calculations where we input the temperature, the total number of atoms, the trapping frequencies and the scattering length. This is possible since the interaction parameter of (k F a) −1 = 1.61 is still close to the BEC regime. Supplementary Figure 5 a-c shows there is good agreement for the calculated and measured density distributions, which validates this approach.