Linear response of a superfluid Fermi gas inside its pair-breaking continuum

We study the signatures of the collective modes of a superfluid Fermi gas in its linear response functions for the order-parameter and density fluctuations in the Random Phase Approximation (RPA). We show that a resonance associated to the Popov-Andrianov (or sometimes “Higgs”) mode is visible inside the pair-breaking continuum at all values of the wavevector q, not only in the (order-parameter) modulus-modulus response function but also in the modulus-density and density-density responses. At nonzero temperature, the resonance survives in the presence of thermally broken pairs even until the vicinity of the critical temperature Tc, and coexists with both the Anderson-Bogoliubov modes at temperatures comparable to the gap Δ and with the low-velocity phononic mode predicted by RPA near Tc. The existence of a Popov-Andrianov-“Higgs” resonance is thus a robust, generic feature of the high-energy phenomenology of pair-condensed Fermi gases, and should be accessible to state-of-the-art cold atom experiments.

www.nature.com/scientificreports www.nature.com/scientificreports/ functions, of the collective modes (called Popov-Andrianov modes hereafter) predicted by refs. 19,20,34 based on the analytic structure of the functions continuated to imaginary frequencies. There are two major obstacles 24,28,35 to the observation of the Popov-Andrianov-"Higgs" resonance in a conventional fermionic condensate. (i) So far the resonance has been clearly identified only in the modulus-modulus response function, whereas experiments (both in superconductors 24 and ultracold Fermi gases 13,22 ) usually excite or measure the density of the fermions. (ii) In a conventional fermionic condensate, where the resonance energy is above 2Δ and the resonance broadened by its coupling to the pair-breaking continuum, it is generally not known whether a quality factor and spectral weight large enough to allow for an observation can be reached. Most studies then look for situations where the damping by the continuum is absent, as in Charged-Density-Wave superconductors 25,[27][28][29][30] , inhomogeneous systems 36 or superfluids in unconventional lattice geometries 35 . Here, we show that the resonance is observable in the density-density and density-modulus response functions at strong coupling even in a conventional fermionic condensate. In those density responses, the spectral weights of the resonance tends to zero with the wavevector q while the quality factor decreases when q increases. Nevertheless we could identify an intermediary regime ( µ ≈ q m 2 at unitarity) where the resonance, and the characteristic quadratic dependence on q of its peak frequency, should be resolvable from the continuum background in an ultracold Fermi gas.
We study the response functions in Anderson's Random Phase Approximation (RPA) 37 . We use the formulation of ref. 38 in terms of bilinear quasiparticle operators that we generalize to nonzero temperature and to the presence of external drive fields. The RPA captures the coupling of the collective modes to the two fermionic continua (and the corresponding broadening of the resonances in the response functions) but neglects other couplings, in particular to the continua of two 39 , three 40,41 or more collective excitations. We show that in this approximation, the density fluctuations are sensitive to the fluctuations of Δ, so that both modulus and phase collective modes are visible in the density response, but that the converse is not true. We give explicit expressions of each element of the response function matrix 1,[42][43][44] , and show that they agree with path-integral based treatments 21,45 .
As the spectrum and response-function signatures of the low-energy collective modes is known in the RPA at zero 37,46,47 , nonzero temperature 38,48 and near the critical temperature T c 19,21,49 , we concentrate here on the high-energy ( 2 ω > Δ) modes. At zero temperature, we show that the resonance of the Popov-Andrianov-"Higgs" mode is visible not only in the modulus-modulus response 20 but also as a global extremum (in the region 2 ω > Δ) in the modulus-phase and modulus-density responses, and as a local extremum in the density-density response at strong coupling. As suggested by the analytic structure found in ref. 34 , we show that the branch remains observable at large q (in particular at µ ≈ q m 2 in the weak-coupling limit µ Δ ) with a quality factor below, but not much below unity.
At nonzero temperature, where the RPA captures the thermal population of the fermionic quasiparticle branches (and only of those branches) and describes the collective modes in the collisionless approximation, we show that the Popov-Andrianov resonance is not destroyed by the presence of thermally excited fermionic quasiparticles. On the contrary, the increase of temperature (which reduces Δ) favours the observability of the resonance in the density response functions by increasing the resonance spectral weight. The shape of the resonance is weakly affected by temperature, and for the order-parameter responses this shape is actually the same as at zero temperature for a slightly different interaction strength. Close to the critical temperature T c , we show that collective mode in the pair-breaking continuum branch is not hidden by the low-velocity phononic branch 21 as long as q m / / 2 2 2  µ Δ . This is in contrast with the Anderson-Bogoliubov branch, which disappears near T c according to the RPA.
Altogether our findings confirm the observability of the Popov-Andrianov-"Higgs" branch, which appears, after our study, as the strongest feature of the high-energy phenomenology of pair-condensed Fermi gases. It is observable in wide ranges of values of the interaction strength, exciting wavevector and temperature, and it is only weakly affected by the singularities caused in the response functions by the structure changes of the fermionic continuum. We are then optimistic about its observability, especially if the experiments can access one of the modulus response functions (modulus-modulus, modulus-phase or modulus-density). The modulus of the order-parameter can be excited by a Feshbach modulation of the scattering length 20,50 , after which the modulus-density response can be measured by absorption images as in refs. 13,22 . Alternatively the density can be excited by a Bragg pulse 13 or by shaking the confinement walls 22 , and the order-parameter modulus measured after a bosonization of the Cooper pairs. In the density-density response, it would be interesting to see if the peak observed in 13 above 2Δ has the characteristic behavior of the Popov-Andrianov-"Higgs" mode, that is, a quadratic dependence on q for both the peak frequency and its width.

BCS Theory At Nonzero Temperature
To derive the matrix of linear response functions, we use the formalism of ref. 38 , itself based on the RPA approach of Anderson 37 , and we generalize it to the presence of pairing and density exciting fields. We start by briefly recalling the formalism of BCS theory at nonzero temperature. In real and momentum space, the Hamiltonian of an isolated gas of fermions in two internal states σ = ↑ ↓ / with S-wave contact interactions is given by www.nature.com/scientificreports www.nature.com/scientificreports/ We use from now on the convention = = k 1 B . To introduce a momentum cutoff in a natural way, we discretize space into a cubic lattice of step l and impose periodic boundary conditions (in a volume = V L 3 ), which restrict the values of the wavevectors to D . The bare coupling constant g 0 is renormalized to reproduce the correct s-wave scattering length of the two-body problem: At the end of the calculation we take the lattice spacing l to 0, and thus g 0 tends to 0 to compensate the divergence of the integral on the right-hand-side of (3). BCS theory describes the equilibrium state at temperature T by the Gaussian state:Ẑ where Z is the partition function and the BCS Hamiltonian H BCS is obtained by treating the interactions in the mean-field approximation, i.e. by replacing the quartic interaction term g r r r r , through the introduction of the self-consistent pairing-field . This quadratic Hamiltonian can be diagonalized easily into a Hamiltonian describing fermionic elementary excitations on top of a ground state energy Here the eigenenergy of a fermionic excitation is The fermionic quasiparticle operators γ σ k, are obtain after a Bogoliubov rotation of the particle operators a k, σ as in the zero temperature case:ˆˆˆ † k kk k k with the Bogoliubov coefficients U k and V k : The difference with the zero temperature case lies in the average values of the bilinear operators: † which now depend on the Fermi-Dirac occupation number This thermal population of the quasiparticle modes also affect the gap equation: Thus, at nonzero temperature, BCS theory captures the effects due to the thermally excited fermionic quasiparticles (the broken pairs); it completely neglects that there are also thermally excited collective modes (some of which are gapless) which is a serious limitation, particularly at strong coupling.

RPA Equations of Motion in Presence of Drive Fields
To study the linear response of the gas, we introduce, on top of the Hamiltonian (2) of the isolated gas, a quadratic Hamiltonian describing the experimental driving of the system: Here the fields σ u r ( ), coupled to the density of spin σ fermions, describe for instance a Bragg excitation of the gas 13 . The complex field φ r ( ) coupled to the quantum pairing field r r can be imposed for instance by a Feshbach-modulation of the interaction strength 50 . An excitation coupled to the phase of ˆ † † ψ ψ ↑ ↓ r r ( ) ( ) can be achieved using a time-and space-dependent Josephson junction as proposed in ref. 21 . This drive Hamiltonian decomposes into a sum of Fourier components of the momentum q transferred to system: In the framework of linear response theory, we seek the response of the system to first order in the fields φ and u σ . We thus neglect the quantum fluctuations in the terms of the equations of motion deriving from Ĥ drive :ˆˆˆˆâ The resulting equations of motion of the bilinear particle operators are given in Appendix A. We give here the equations of motion in their simplest form, which is in the quasiparticle basis. At the level of the bilinear operators, the Bogoliubov rotation (8) and (9) . Performing this change of basis on the equations of motion, we get: (the subtraction of the mean-field average value matters only when q = 0). At the linear order, the sole effect of the drive fields is thus to shift the collective quantities which enter in the equations of motion: Note that one recovers the zero temperature system Eqs. (14)(15)(16) of ref. 38 by setting f 0 k = (in which case the equations of motion of the ˆ † γ γ operators become trival).

Linear Response to a Periodic Drive
Matrix of response functions of a driven system. We now assume that the system is driven at a fixed We can then replace the time derivatives ∂ i t in Eqs. (22)(23)(24)(25) by i0 ω + + . Rederiving with respect to time and resumming the system to form the collective quantities (26)(27)(28)(29) yields the 4-dimensional linear system where  is the identity matrix. We have introduced the density and polarisation fluctuations and reparametrized the fluctuations of the order-parameter as:ˆê We treat the phase θ q of the order-parameter as an infinitesimal and therefore linearize the exponential in (31) and (32), which is consistent with our symmetry-breaking approach where the expansion is done around the mean-field state with a real Δ. The linear response matrix 1,42-44 , which relates the fluctuations of the density and order-parameter to the infinitesimal drive fields, is then V g 0 www.nature.com/scientificreports www.nature.com/scientificreports/ To describe the experimental behavior of a driven system, one usually concentrates on the imaginary part of χ, which describes the energy absorbed by the system 1 (whereas the real part describes the energy refracted or reflected by the system).
The matrix χ is expressed in terms of the 4 × 4 matrix Π of the pair correlation functions 52 , computed in the BCS thermal state (4): where we generalize the notations of refs. 34,51 : Here a and b are one of the functions + W , W − , w + or − w of k and q, and ε ± kq is short-hand for ε ε k q k q /2 /2 ± + − . The first and second matrix in the right-hand side of (36) are the contribution of respectively the quasiparticle-quasiparticle and quasiparticle-quasihole continua to Π. In our unpolarized system, the polarization fluctuations ˆ− ↑ ↓ n n q q are entirely decoupled from the other collective fields. Note that the response functions computed here for a driven system also give access, through a Laplace transform 34 , to the time response following a perturbation localized in time.
Remark that, up to some signs, the matrix Π has a tensor-product structure when expressed in terms of the vector = :

Eigenenergies of the collective modes.
The response of the system should diverge when the drive frequency coincides with the eigenfrequency ω q of a collective mode; to find those eigenfrequencies, one should thus search for the poles of χ, in other words the zero of its denominator: When Π has a branch cut on the real axis (which occurs for all  ω ∈ at nonzero temperature and for min ( ) at zero temperature), this equation cannot have a real solution. Its analytic continuation to the lower-half complex plane may however have solutions describing damped collective modes. Numerical and analytic methods to continue the matrix Π through its branch cuts have been described in refs. 20,21,34 . Note that the bare density-density response function Π 33 may have poles in the complex plane, which remain poles of the dressed response χ 33 .
Explicit expressions of the response functions in the limit g 0 → 0. In the limit of zero lattice spacing ( → l 0), g 0 tends to 0, Π 11 and Π 22 are equivalent to V g / 0 , while 33 Π and 44 Π have a finite limit (to interpret physically the elements of Π the reader can use the correspondance θ ρ → Δ p 1, 2, 3, 4 , , , ). We thus have the equivalences www.nature.com/scientificreports www.nature.com/scientificreports/   (41) to compute the matrix product in χ we obtain: We have used the notation Π = Π − ∼ ). For the density response functions we have in particular:  which coincides with the explicit expressions obtained by refs. 21,45 in the path-integral formalism. At weak coupling ( µ Δ → / 0) and q O( ) = Δ , the modulus-phase and modulus-density matrix elements 12 Π and 23 Π vanish such that the collective modes are either pure modulus modes (if their eigenenergy solves ).

Angular points of the response functions.
We conclude this section by remarking that the response functions have the same angular points as the spectral density The quasiparticle-quasiparticle part of the spectral density (which originates from the first matrix in (36) with ( ) in the denominator) is nonzero when the "pair-breaking" resonance condition is satisfied Physically, when this resonance condition is met, the drive field can break the pair of total momentum q into unpaired fermions of momenta k q/2 + and − k q/2. As a function of the increasing drive frequency ω, this www.nature.com/scientificreports www.nature.com/scientificreports/ resonance condition is (i) satisfied by no wavevector when ω ω < 1 (with 2 1 ω = Δ at low-q), (ii) satisfied for a connected set of wavevectors around the dispersion minimum of the BCS branch for 1 2 ω ω ω < < , (iii) satisfied by two connected sets of wavevectors, one in the increasing and one in the decreasing part of the BCS branch for 2 3 ω ω ω < < and (iv) satisfied for a connected set of wavevectors in the increasing part of the branch for ω ω > 3 . These three boundary frequencies 1 ω , 2 ω and ω 3 will appear as angular points in the spectral function and therefore in the response functions.
The quasiparticle-quasihole part of the spectral density (which originates from the second matrix in (36) in the denominator) is nonzero when the "absorption-emission" resonance condition is satisfied In this case, the drive field is not breaking the Cooper pairs but simply transferring more energy to the unpaired fermions already created by thermal agitation. This requires much less energy, which is why the quasiparticle-quasihole continuum is gapless. As a function of ω, this resonance condition can be met on i ( ) two disconnected sets of wavevectors, one in the increasing and one in the decreasing part of the BCS branch for ω ω < ph , ii ( ) a single connected set of wavevectors in the increasing part of the branch for ph ω ω > . With ph ω , we have found the four angular point of the response function ( ) ].

Long Wavelength Limit
In the long wavelength limit ( Δ q m 2 ) the solutions of the collective mode Eq. (40) and the behavior of the response functions can be studied analytically. Below the pair-breaking continuum (at energies lower than ε ε min ( ) the problem has been studied in-depth at zero and nonzero temperature. At zero temperature a real solution of (40) corresponding to the Anderson-Bogoliubov sound branch can be found 37,46,47 . This branch appears as a pole in Reχ and as a Dirac peak in Imχ because the zero-temperature response functions are free of branch cuts below the pair-breaking continuum. This resonance was observed experimentally in the density-density response function by Bragg spectroscopy 13 and in the low-q limit it can be identified with hydrodynamic first sound 8,12,22 . When the dispersion is supersonic (i.e. the function ω q q is convex) the resonance is broadened by intrinsic effects not captured by RPA 39,41 .
By contrast, the behavior of the response functions at high energy ( ε ε min ( ) ) is known in literature only at zero temperature. There, the collective mode Eq. (40) has a complex root departing quadratically with q from 2Δ 19 , and the modulus-modulus response shows correspondingly a resonance at energies above 2Δ 20 . Here, we show analytically that (within the RPA) the same resonance exists at nonzero temperature, even until the vicinity of T c . This result is in agreement with Popov-Andrianov's claim that the collective mode has the same energy at zero temperature and when T T c → .
General case. We first compute the matrix elements ij Π on the upper-half complex plane, that is for z Im 0 > . The long wavelength expansion can be perform using the method exposed in ref. 20 : since the energy is expected to depart quadratically from 2Δ, one parametrizes it as The window [2 , ] 2 ω Δ between the first two angular points of the branch cut corresponds in the limit → q 0 to ζ ∈ [0, 1]. With respect to the zero temperature case 20 , the phase-phase and modulus-modulus matrix elements are simply multiplied by a factor T tanh( /2 )  (1,1) or (2,2) and α Π = Π ij iǰ for the other matrix elements, the nondimensionalization factor α being π Δ Δ L m (2 / 2 ) 3 (the dimensionless response functions χ ij will accordingly be multiplied by the appropriate power of α: α −1 for 11 χ , χ 22 and χ 12 , α 0 for χ 13 and 23 χ and α 1 for χ 33 ). We have written the complex functions f 11 and f 22 of ζ in such a way that their spectral density (their imaginary part in the limit Im 0 ζ → + ) is directly given by their last term. The modulus-phase matrix element is independent of ζ (it can be approximated by its value in q 0 = , = Δ z 2 ) but, unlike 11 Π and 22 Π , depends on temperature through the two ratios T / Δ and T / µ : www.nature.com/scientificreports www.nature.com/scientificreports/ˇˇπ where we have introduced the improper integral (P denotes the Cauchy principal part): Note that the quasiparticle-quasihole integrals (38) have been negligible in deriving expressions (51)(52)(53).
To find a root of the collective mode Eq. (40), one analytically continues Π 11 and Π 22 (and hence the determinant of Π) from upper to lower half-complex plane (the forms given in Eqs. (51) and (52) At zero temperature, this equation was derived in ref. 19 in the weak-coupling limit and ref. 20 in the general case. In this low-q limit, the only dependence on temperature is through the ζ-independent second term of (55).
Close to the phase transition. Unlike in the phononic regime (q 0 → with z cq = ) 21 , no dramatic phenomenon occurs for the collective mode of the pair-breaking continuum when T tends to the critical temperature T c . We recall that in the RPA, the limit of the phase transition from the superfluid phase corresponds to The RPA also assumes an infinite fermionic quasiparticle lifetime and thus describes the collective modes and their damping by the fermionic continua in the collisionless approximation. The function g 12 tends to a finite nonzero constant in the limit → T T c : . In fact, the resonance near T c for a given value of k a 1/ F has exactly the same shape as the T 0 = resonance for a lower value (corresponding to weaker-coupling) of k a 1/ F . Using an equation-of-state to relate k a 1/ F to both µ T / c c and T T ( 0)/ ( 0) µ = Δ = 21 , the corresponding values a 0 and a c of the scattering length at = T 0 and T c are found by solving:   Thus, the response functions have exactly the same shape (they coincide up to a proportionality factor) near T c as they have at = T 0 for the slightly different value of the interaction strength given by (59). In Fig. 1, we show how the shape of the order-parameter response functions ( 11 χ , χ 22 and 12 χ ) change when going from the BCS limit ( k a to the threshold of the BEC regime where µ vanishes. Exploiting the equivalence (59), the figure describes together the crossover at T 0 = and T T c → . Irrespectively of the interaction regime, the phase-phase response is a monotonously increasing function of the drive frequency and only reflects the incoherent response of the pair-breaking continuum, without collective effects. Conversely, both the modulus-modulus and modulus-phase response functions display a maximum that can be interpreted as a collective mode in the BCS limit (black curves) and up until unitarity (blue curves). As explained in ref. 20 , this maximum can be fitted to extract the frequency and damping rate of the collective mode to a good precision. The fit function to use is , where the complex parameters z q , Z q and C q represent respectively the complex energy of the collective mode, its residue, and an incoherent flat background. A remarkable effect of this background C q is to displace the location of the www.nature.com/scientificreports www.nature.com/scientificreports/ maximum of 22 χ and χ 12 to respectively ζ . 0 4 and 0 1 ζ . in a very broad interaction range. The variations of ζ s (which decreases when increasing the coupling strength) are thus not visible by simply looking at the maximum location. Soon after unitarity, the resonance in χ 22 and χ 12 disappears and only a sharp feature near 2 ω = Δ remains. This abrupt lower edge of the continuum is in ζ = 0 so it is not departing quadratically with q from 2Δ (see also the color Fig. 6 in the BEC regime) as the Popov-Andrianov resonance does in the BCS regime, and it can no longer be interpreted as a collective mode. As understood in ref. 20 , this is because the complex root z q of the collective mode Eq. (40) has a real part below 2Δ (i.e. ζ < Re 0 s ) and does no longer trigger a resonance inside the pair-breaking continuum.
Density matrix elements in the long wavelength limit. We now study the density responses of the system χ i3 , = i 1, 2, 3 in the long wavelength limit → q 0 at energies above but close to 2Δ. For → q 0 at fixed ζ, the three matrix elements needed to compute the density response functions are given by:ˇˇδ     /2 ), located in this limit around the minimum k 0 of the BCS branch. For those wavevectors, we set 0 and expand the integrand in (37) at fixed K. This yields the leading order contribution to Π 11 , 22 Π and 13 δΠ . For Π 12 , 23 δΠ and δΠ 33 the leading order is dominated by the wavevectors away from k 0 and is obtained by expanding directly in powers of q at fixed k (with a contribution of the quasiparticle-quasihole integrals from Eq. (38)). For δΠ 33 specifically, the subleading order O q ( ) 3 (which matters for the imaginary part of the response function Im 33 χ ), is obtained by subtracting the leading one and then using the reparametrisation of the wavevectors, Eq. (69).
Using the expansions Eqs. (63-65), we obtain the expressions of the density response functions:χ  In Fig. 2, we show the density response functions on the BCS side of the crossover. Unlike for the order-parameter response functions, no exact correspondance between zero temperature and the transition temperature can be found by changing the interaction strength (this is due to the temperature dependence of g 23 ), so we show separately the functions at T 0 = (in solid curves) and at T T c → (in dashed curves). The difference between the = T 0 and T T c → curves (after the appropriate rescaling) remains however fairly small, and tends to 0 in the BCS limit (black curves). Remarkably, a minimum characteristic of the Popov-Andrianov collective mode is visible in all three density responses. In 23 χ , this minimum is a global minimum (for ζ ∈ [0, 1]) which exist (as www.nature.com/scientificreports www.nature.com/scientificreports/ in 22 χ and 12 χ ) from the BCS limit up until unitarity. For the density-density and density-phase responses 33 χ and 13 χ , this minimum is a local minimum, which exists close to unitarity (blue curves) around ζ = . 0 1. Because of the decoupling between the phase-density fluctuations and the modulus fluctuations in the weak-coupling limit, this minimum disappears from 13 χ and 33 χ when → −∞ k a 1/ F (black curves). After unitary, when approaching the BEC regime (orange curves), the resonances in all three density responses are replaced by a sharp edge in 2 ω → Δ + ( 0 ζ → + ). This is the same phenomenon as in the order-parameter response functions.

Coexistence with the phononic collective modes near T c .
To compare the Popov-Andrianov resonance to the other collective effects of a superfluid Fermi gases near T c , we show in Fig. 3 the response functions from ω = 0 up until ω > Δ 3 in the strong coupling regime and temperature close to T c . The sharpest feature in both the order-parameter and density responses is the resonance, at very low energy (that is at uq ω = with a velocity u T T c ∝ − ), of the collisionless phononic collective mode found in ref. 21 . Still at phononic energies ω ∝ q, the density-density response function shows a broad peak caused entirely by 33 Π (shown as a black dashed line) and also noticed in ref. 21 . This is simply the peak of the Lindhard function, which exists also in the normal phase. Finally, inside the first window [2 , ] 2 ω Δ of analyticity of the pair-breaking continuum, all response functions show the peak characteristic of the Popov-Andrianov resonance, whose shape matches the one shown on Figs. 1 and 2. Due to the absence of rescaling with the wavevector q in Fig. 3, the peak is much more intense in the modulus-modulus response, and to a lesser extent in the modulus-density response, than in the density-density response.
Experimental protocol. Our results suggest a very simple experimental protocol to observe the resonance: using a Bragg spectroscopic measurement as in ref. 13 , one should observe that the first extremum above 2Δ varies quadratically (both in location and width) with q, a behavior which can be viewed as the fingerprint of the Popov-Andrianov-Higgs mode. The optimal interaction regime is around unitarity and the optimal wavevector q 0 1). The dimensionless response functions are shown this time without rescaling and from ω = 0 until far inside the pair-breaking continuum. The vertical dashed lines show the angular points of the fermionic continua, from left to right ph ω , ω 1 (= Δ 2 ) and 2 ω . On the top panel, the black dashed line is the pure density contribution Im 33 Π to the density-density response (see Eq. (46)) and the inset is a zoom on the behaviour near 2 ω = Δ.
www.nature.com/scientificreports www.nature.com/scientificreports/ where the shape of the response functions changes dramatically. In the bottom left panel, we also show the contribution of the pure density-density fluctuations 33 Π to the total density response (black and blue dotted line). Bottom right panel: ˇχ Π − Im( ) 33 33 is shown in colors as a function of ω and q (the color scale is logarithmic) after division by q 3 . The angular points 2 ω and ω ω ≥  www.nature.com/scientificreports www.nature.com/scientificreports/ is around . × Δ m 0 5 2 (q should not be too small to avoid the q 3 cancellation of χ 33 near 2Δ but not too large either otherwise the minimum is reabsorbed by the continuum edge, see the lower panels of Fig. 5).
Alternatively, the resonance could be observed through the modulus-density response function 23 χ by (i) exciting the order-parameter modulus δ Δ through a modulation of the scattering length at frequency ω and wavelength π q 2 / and (ii) measuring the intensity of the density modulation δρ at wavelength q 2 / π . This should be easier than the scheme of ref. 20 which proposed to measure χ 22 by interferometry. Using the symmetry of the response matrix χ, one can also excite the density δρ (by a Bragg pulse 13 or using the trapping potential 22 ) and measure the order-parameter modulus δ Δ either by interferometry or by bosonizing the Cooper pairs through a fast sweep of the scattering length, as was done in ref. 23 .

At Shorter Wavelengths
Outside the long wavelength limit, that is µ ≈ Δ q m m 2 , 2 when µ and Δ are comparable (see Eq. (90) in ref. 34 for a more detailed discussion of the limit of validity of the long-wavelength limit), we study the response functions by performing numerically the integral over internal wavevectors k in Eqs. (37) and (38) (see Appendix B for more details on the numerical implementation).
At zero temperature. Weak-coupling regime. On the left panel of Fig. 4, we show the modulus-modulus response at relatively weak-coupling ( / 10 µ Δ = ) and zero temperature as a function of ω (rescaled as in the long wavelength section) for increasing values of the wavevector q. On the right panel, we show the same dispersion relation but in colors, with q on the x-axis and ω on the y-axis. The Popov-Andrianov resonance we have characterized at low q remains as a broader and shallower maximum as q increases (see the rescaling of the x and y-axis on the left panel of Fig. 4) that travels roughly quadratically through the continuum. In the modulus-modulus response function, the augmentation of wavevector is thus unfavourable for the observation of the resonance in the pair-breaking continuum. Note that the location of the maximum is discontinuous when crossing 2 ω and 3 ω (which both decrease with q), but remains a monotonously increasing function of q. The non-monotonic behavior of the collective mode eigenfrequency z q found in the analytic continuation through the interval ω Δ [2 , ] 2 of the real axis 20 is thus not reflected on the response function. In fact the angular points ω 2 and 3 ω only slightly affect the shape of the resonance when they cross it (see in particular the black curve on the left panel of Fig. 4). This is consistent with the finding of ref. 34  . The same robustness towards the choice of the real axis interval through which the analytic continuation is made was noticed by ref. 21 for the phononic modes. It is a sign that the Popov-Andrianov collective mode is a fundamental physical phenomenon, which does not depend on a specific configuration of the fermionic continuum.
Strong-coupling regime. Conversely, the increase of q favours the observability of the resonance in both the modulus-density and density-density response functions at strong coupling. On Fig. 5, we show 23 χ and χ 33 (as well as χ 22 ) at unitarity (µ Δ . / 086 ) and still at zero temperature. As long as it does not encounter the singularity in 3 ω , a smooth extremum (in 23 χ and χ 33 it's a minimum) whose location increases quadratically with q remains visible. The resonance broadens with q, but this is compensated by a deepening of the resonance peak roughly as q in 23 χ and as q 3 in χ 33 . The resonance in 33 χ is caused by the order-parameter contribution χ − Π 33 33 to the density-density fluctuations (compare the blue dotted and the blue solid line on the bottom left panel of Fig. 5), in which it is a global minimum as a function of ω (rather than a local minimum in 33 χ ). To emphasize the dispersion of the resonance, we thus plot on the bottom right panel of Fig. 5, Im Im 33 33 χ Π − divided by q 3 in www.nature.com/scientificreports www.nature.com/scientificreports/ colors as a function of ω and q. The global extremum of ω χ Π − Im( ) 33 33 is shown as a function of q in white solid line. As long as it stays in the window [2 , ] 2 ω Δ , it varies approximatively quadratically with q. Contrarily to what happens at weak-coupling, the resonance shape at strong coupling is much distorted when going through the singularity ω 3 . This effect is particularly visible on the modulus-modulus and modulus-density responses (upper panels of Fig. 5) where the resonance seems broken in ω 3 such that the smooth extremum has disappeared in favour of a sharp extremum in ω 3 . On the color plot of Fig. 5, the quadratic growth of the resonance frequency is also visibly halted when it encounters the angular point in 3 ω . This is not surprising since the poles found in the analytic continuation through windows [ In the BEC regime. In the BEC regime (that is for us when µ < 0), the lower-edge of the pair-breaking continuum is no longer flat at low q, but increases quadratically with q. Although a pole can be found in the analytic continuation through the interval ω + ∞ [ , [ 3 (the only one available when µ < 0), its real part always stays below ω 3 , such that no smooth peak appears in the response function. Instead there is only a sharp feature pinned at the lower-edge of the continuum. Figure 6, shows the example of the modulus-modulus response function (the other responses have a similar behavior) at / 1 µ Δ = − ( k a 1/ 13 F . ). This sharp feature can hardly be interpreted as a collective mode and only reflects the incoherent response of the fermionic continuum when the pairs are tightly bound.
near T c . At nonzero temperature and even near T c , we have shown in section V that the Popov-Andrianov resonance exists in the limit → q 0 and is almost insensitive to the quasiparticle-quasihole contributions (38) to the fluctuation matrix Π. This is no longer the case at higher q. The angular point ω ph of the quasiparticle-quasihole continuum in particular destroys the resonance as it increases (initially linearly) with q. This effect is illustrated on Fig. 7 showing the modulus-modulus response function near T c : for q m / 2 0 12 Δ = . (orange dashed curve on Fig. 7) the lower tail of the resonance is trimmed by the angular point at ph ω , and for q m / 2 0 3 Δ = . (long-dashed green curve) it is completely hidden. This can be understood by a simple reasoning: near T c , ph ω varies as q m 2 / µ at low q 21 , such that it reaches 2Δ for 1/4 . The long wavelength limit near T c is thus limited to q m /2 / 2 2 µ Δ (as in the weak-coupling case at = T 0 see Eq. (90) in ref. 34 ).

conclusion
We have computed the response function matrix of a superfluid Fermi gas in the Random Phase Approximation at nonzero temperature, and used it to study the observability of the order-parameter collective modes. We have shown that the appearance of a resonance inside the pair-breaking continuum associated to the Popov-Andrianov-"Higgs" mode is a very robust phenomenon which concerns not only the modulus-modulus response function but also the modulus-density and density-density responses, which are easier to measure. At weak-coupling the resonance is observable at all values of the wavevector q and is only weakly sensitive to the angular points created in the response functions by the changes of structure of the fermionic continuum. At nonzero temperature, we www.nature.com/scientificreports www.nature.com/scientificreports/ have shown analytically that the resonance is not destroyed by the presence of excited fermionic quasiparticles, and retains approximatively the same shape as when T 0 = . It also coexists with the low-velocity phononic collective mode which RPA predicts near T c . The spectral weight of the resonance is enhanced in the modulus-density and density-density responses when T increases, which should favour its observability.

Appendix A: Derivation of the equations of motion
We give here a few additional steps leading to the equations of motion (22)(23)(24)(25). In the particle basis, the equations of motion take the form: For ω ω ω < < 1 2 , this quantity is comprised in [0,1] (such that the resonance in (B2) is reached) for ∈ k k k [ , ] 1 2 with k 1 and k 2 solutions of ε ε ω + = instead of the wavenumber k, and ω = Δ t argch( /2 ) instead of the drive frequency, then using the Dirac delta to integrate analytically over the scattering angle u, we have  In the quasiparticle-quasihole spectral density (B4), we give the resonance angle expressed in terms of k m q m /2 /8 , with k k , 3 2 the two solutions of ω − =− + − k q k q /2 /2 ε ε . Using the variable ω ξ = y /2 instead of the wavenumber k, and ω = Δ t arccos( /2 ) instead of the drive frequency, then using the Dirac delta to integrate analytically over the scattering angle u, we have: