Collective excitations of a charged Fermi superfluid in the BCS-BEC crossover

We consider collective excitations in the superfluid state of Fermi condensed charged gases. The dispersion and damping of collective excitations at nonzero temperatures are examined, and the coexistence and interaction of different branches of collective excitations: plasma oscillations, pair-breaking Higgs modes, and Carlson-Goldman phonon-like excitations are taken into account. The path integral methods for superfluid Fermi gases and for Coulomb gas are combined into a unified formalism that extends the Gaussian fluctuation approximation to account for plasmonic modes. This approximation of Gaussian pair and density fluctuations is able to describe all branches of collective excitations existing in a charged superfluid. The spectra of collective excitations are determined in two ways: from the spectral functions and from the complex poles of the fluctuation propagator. A resonant avoided crossing of different modes is shown. It is accompanied by resonant enhancement of the response provided by the pair-breaking modes due to their interaction with plasma oscillations. This may facilitate the experimental observation of the pair-breaking modes.


I. INTRODUCTION
For decades, collective excitations in neutral and charged superfluids have been the subject of great interest in condensed matter physics. Their manifestations are found in a wide range of phenomena, from superconductors and quantum gases to nuclear systems and neutron stars [1][2][3][4][5]. Collective excitations are important for both experiment and theory because they determine the response spectra of condensed systems. Interest in collective excitations in superfluid and superconducting systems has been reinforced by experiments to study their response properties.
The present work focuses on collective excitations in charged superfluid Fermi fluids and superconductors. Several branches of collective excitations are a subject of the present study. The gapless soundlike mode [6][7][8][9], called Anderson-Bogoliubov mode, is well specified in neutral superfluid systems such as cold atomic gases. In superconductors, this gapless mode is affected by the Coulomb interaction and pushed up to the plasma mode. In the long-wavelength limit, the plasma mode is gapped, and the size of the gap is the same both in the superfluid/superconducting and in the normal state [6,10]. Near the transition temperature T c , the other gapless mode can appear in BCS superconductors, discovered by Carlson and Goldman [11,12].
The plasma mode is associated with oscillations of the particle density and is therefore well resolved in the density response. For the pair field response, the belonging of different branches of collective excitations to the pure amplitude and phase responses is only asymptotically exact in the far BCS limit. Plasma and Anderson-Bogoliubov collective excitations are revealed in the pair field response through oscillations of the superfluid phase. There exists also an amplitude mode due to oscillations of the pair field modulus, attributed in many papers to the Higgs mechanism [3,4] and revisited recently in Ref. [13] as the pairbreaking mode. In Ref. [14], spectra of both Anderson-Bogoliubov and pair-breaking modes have been experimentally studied for neutral atomic Fermi superfluids. The pair-breaking collective excitations in superconductors were theoretically predicted long ago [15], but have only recently been discovered experimentally [16].
In this paper special attention is paid to the interaction between different branches of collective excitations in charged Fermi superfluids. In the theory of collective excitations in superconductors, plasma frequency is usually assumed to be very large with respect to the superconducting gap. Here, we focus on the other case, when they are comparable to each other. It may be realized in strong-coupling superconductors where the BCS-BEC crossover regime can exist, particularly in iron-based superconductors. Consequently, the treatment is performed using methods suitable for the crossover. Also, collective excitations in high-T c superconductors [17] can reveal an interaction of plasma and pair-field branches.
The energy spectrum of charge carriers in high-T c superconductors is substantially different from the single-band 3D picture exploited in the present study, but the extension of the formalism to a two-dimensional and multiband system is straightforward. It is a subject of the future investigation.
In the BCS-BEC crossover, neither the plasma, nor the Anderson-Bogoliubov, nor the pair-breaking mode can be attributed exactly to modulus or phase responses, because of amplitude-phase mixing, which is also a subject of attention in this treatment. Another important point of the study is to investigate collective excitations at nonzero temperatures, where the spectrum of excitations and the picture of an interplay of different branches is richer than at zero temperature.
In the present work we apply the Gaussian fluctuation approach, which is well established for a description of collective excitations in neutral superfluids of cold atomic gases [7-9, 18, 19]. It can be straightforwardly extended to charged superfluid Fermi gases by addition of a Hubbard-Stratonovich field which describes oscillations of the particle density [20,21] and is promising also for application to superconductors. The Gaussian fluctuation method is equivalent to the extraction of the excitation spectra from the linear response within the random phase approximation (RPA) which is also frequently used in the theory of collective excitations [6,10,12,13].

A. Effective action and Gaussian fluctuation approximation
We consider a charged Fermi gas using the path-integral formalism. The thermodynamic properties of an interacting Fermi gas are determined by the partition function represented through the path integral over Grassmann field variables ψ σ , ψ σ with the action functional S, where the model attractive interaction for the pairing channel is expressed by a contact potential with the coupling constant g < 0. Fermions are assumed to have the spin 1/2, σ = (↑, ↓) are the spin projections. The part of the action describing the Coulomb interaction is with the particle density and the Coulomb interaction potential where ǫ 0 is the permittivity of free space, ε is the high-frequency dielectric constant of a medium. In the absence of the Coulomb interaction, (2) turns into the widely used fermionic action [24].
We choose the units: = 1, m = 1/2 and k F ≡ (3π 2 n) 1/3 = 1, which leads to the equality E F = 1 for free-fermion Fermi energy E F ≡ 2 k 2 F 2m . For a charged gas, one more input parameter appears: the effective charge e/ √ 4πǫ 0 ε. It can be expressed through the dimensionless parameter, the Coulomb α 0 having some analogy with the Fröhlich electronphonon interaction constant α: In these units, the bare plasma frequency ω p = e 2 n/ǫ 0 εm is expressed as It is worth discussing how can the chosen model potential be relevant to experiments on collective excitations in charged Fermi superfluids. In order to clarify the subject of the study, we note that this term means superconductors within the present work. Another class of systems which might represent an interest for the application of the used theoretical method and its future development, are ultracold ionic or atom-ionic gases [25]. However, they do not yet realize a superfluid state for a charged component, while the present treatment involves pair field as an essential element. Therefore at the present state-of-art of experiment, other charged fermion systems are not considered in the present investigation.
In rather early works [7,24], the effective action approach is used exploiting only the contact model pairing interaction, assuming that observable effects of a complete true interaction potential can be summarized through a single parameter, the effective scattering length. In these works, the discussion was therefore performed in the context of both superconductors and atomic superfluids. This approach however was not able to describe the plasma branch of collective excitations, because it needs to account for the Coulomb interaction explicitly.
For a Coulomb gas without pairing, the path-integral approach using the effective action with the Coulomb interaction described by (3) and (5) [23] appears to be equivalent to the random phase approximation. It effectively describes the plasma collective excitations in the normal state of a gas of electrically charged fermions.
The straightforward way to consider collective excitations taking into account both plasma oscillations and excitations of the pair field is to include the Coulomb repulsion potential as a separate term in the model interaction potential. The same model potential as in the present treatment, which is a sum of pairing interaction and Coulomb potentials, has been already used in a series of preceding publications, e. g., using the contact [10,22] or finite-width [20,21] pairing interaction potential, or an effective pairing potential provided by the electron-phonon interaction [6].
The repulsion between electrons is screened in superconductors, so that the effective repulsion can differ from the Coulomb potential (5). It should be noted however that the random phase approximation (which is equivalent to the Gaussian fluctuation approximation applied below) leads, in particular, to dynamic screening of the Coulomb interaction. As RPA leads to the account of dynamic screening, it is logical to write down the Coulomb potential in the non-screened form, to avoid double counting. See also the discussion of Eq. (9) in Ref. [6], where screening should be taken into account for exchange terms. The exchange terms are neglected in the present work as well as in the preceding paper [26]. They can be non-negligible at strong-coupling. Here, we suggest that the calculated spectra of collective excitations remain qualitatively correct in the BCS-BEC crossover except maybe the BEC regime, which is beyond the scope of the present work.
The other part of the interaction potential, which is responsible for the pairing channel, is applied here in the form of a contact potential. As long as calculated results are expressed in terms of the scattering length, a true shape of the potential has no significance, because we have no aim to derive this interaction from the first principles. The absence of such derivation of course makes the theory more phenomenological than first-principle theories.
However its experimental relevance can be kept if we match the scattering length with, for example, the ratio ∆| T =0 /E F , which can be an independent input parameter (what is done in figures with numeric results of the present work). Also, the same method can be used to compare different theoretical approaches to each other.
The effective bosonic action for a Fermi superfluid is obtained after introducing two auxiliary bosonic fields: the pair field Ψ , Ψ and the density field Φ, by adding them to the fermionic action as follows: The next step is the Hubbard-Stratonovich (HS) transformation, which shifts the bosonic fields in order to remove the fermionic interaction in S. After this, we use the Nambu representation of fermionic spinors, determined as The extended action (7) after the HS shift is then where the quadratic form with the inverse Nambu matrix is given by: It differs from the analogous matrix for a neutral Fermi superfluid by the presence of the terms provided by the Coulomb interaction, which are proportional to the density field Φ.
The integration over fermionic variables in the partition function with the action S ′ ext is performed formally exactly and leads to the partition function expressed as the path integral over bosonic pair and density fields, with the effective bosonic action For the subsequent derivation, we apply the Fourier representation for fermionic and bosonic fields, with the fermion and boson Matsubara frequencies The effective bosonic action then takes the form In order to consider thermodynamics and response of the superfluid fermionic system with the effective bosonic action (15), the lowest-order approximation is the saddle-point one, which determines macroscopic values of the field variables (Ψ, Φ) from the least action We apply trial saddle-point values of the pair and density fields to be uniform in space, to consider collective excitations on top of a uniform background. Also coordinate-dependent saddle-point solutions are possible [6,9], but they are beyond the scope of the present work.
The uniform saddle-point value for the density field appears to be equal to zero. Therefore we arrive at the gap equation for the saddle-point value of the pair field, which takes the same form as for the neutral superfluid, where E k = ξ 2 k + ∆ 2 is the BCS excitation energy and ξ k = k 2 − µ the free-fermion energy. We note that in this formulation, our saddle-point approximation misses the exchange scattering contributions which typically add a term of the form Eq. (17) where V (q) = q 2 /8π is the Fourier transform of the Coulomb potential. The coupling constant of the contact interaction is renormalized through the scattering length a s [24]: and X (E) is the function The next approximation takes into account the fluctuations about the saddle point, where ϕ is the pair fluctuation field. Introducing the amplitude-phase representation with, respectively, amplitude and phase fluctuations λ q,m and θ q,m . These Fourier components correspond to real amplitude and phase fields in time and space, so thatλ q,m = λ −q,−m andθ q,m = θ −q,−m . The resulting Gaussian pair and density fluctuation (GPDF) action is given by: The matrix elements of the pair-field (GPF) part of the GPDF action are equivalent to those obtained in preceding works, e. g. Ref. [18], and read: The other matrix elements have been derived in Refs. [26,27]. Explicitly, they are: which satisfy the symmetry properties and the diagonal matrix element is B. Analytic determination of collective excitation spectra Energies and damping factors of collective excitations can be determined through complex poles of the GPDF propagator, following to the procedure proposed by Nozières [28]. The formalism remains the same as in Refs. [9,13,19]. For completeness, we reproduce here its main steps. Formally, the complex poles of the GPDF propagator (det K) −1 are determined by the equation in the complex z plane, Due to the branch cut at the real axis, these poles can become visible only after the analytic continuation of the propagator through the branch cut to the next sheet of the Riemann surface. Otherwise, they are hidden by the branch cut as behind a mirror wall.
The analytic continuation is performed using the spectral density function. At the real axis, it is determined by: The spectral density is analytic on the real axis except maybe a finite number of points. In any chosen interval between these points it can be straightforwardly continued analytically to complex z plane. Let us label these intervals by the index n, and denote the spectral density in each interval as ρ (n) K (ω), so that the analytic continuation of the spectral density from each interval is ρ (n) K (z). The analytic continuation det K ↓ (z) of det K (z) through the branch cut in the n-th interval is then The set of eigenfrequencies and damping factors for collective modes is therefore determined by roots of the equations det K (n) for all intervals n. The bounds between different intervals in which ρ K (ω) is analytic are the angular points. They indicate a change in the configuration of the resonant wave vectors for one of the resonance conditions, The case ω > 0 is considered, without loss of generality. at which different angular-point frequencies coincide. Fig. 1 shows an example of angular points and, correspondingly, intervals for the analytic continuation. A detailed description of them can be found in Ref. [19]. Briefly, the angularpoint frequency ω 1 is the pair-breaking continuum edge. The angular points ω (pp) 2 (particleparticle) and ω (ph) 2 (particle-hole) correspond to, respectively, local extrema of the energies E k+ q 2 + E k− q 2 and E k+ q 2 − E k− q 2 as functions of k at k q. The frequency ω 3 is the energy of the BCS pair E k− q 2 + E k+ q 2 at zero momentum k. The angular point ω Because of multiplicity of intervals, some roots coming from the analytic continuation through different intervals may be physically equivalent, i. e., correspond physically to the same modes. In this case, a selection of appropriate roots must be performed using physical reasoning, as discussed below in Sec. III.

C. Response functions
The spectra of collective excitations can be rather qualitatively but reliably extracted from spectral functions for several types of the response. Also, spectral functions give an information on relative magnitudes of peaks provided by different branches of collective excitations. Here, we begin with two spectral functions for the pair field response: the modulus-modulus and phase-phase spectral function, determined using the bosonic Green's functions in the Matsubara representation, The spectral functions are obtained after the analytic continuation from a sequence of bosonic Matsubara frequencies to the complex plane near the real axis, Because the pair-density fluctuation action is quadratic, the bosonic averages are explicitly calculated in a standard way. The resulting spectral weight functions are then given by: .
The density-density response function χ ρρ is determined as through the Green's function G ρ (q, z) which is the analytic continuation to the complex z plane of the Matsubara Green's function for Fourier components of the fermion density The averages in (43) are obtained using the generating functional with the auxiliary source field υ, using the relation The terms in (44) containing the source field, are included to the fermionic action (2), only resulting in addition of υ to the chemical potential. Therefore the generating functional (44) is calculated analytically in the same way as in Subsec. II A: performing the HS transformation and integrating over fermionic fields. As a result, we arrive at the extended effective fluctuation action, Here, S GP F is the part of the GPDF action, which only describes pair fluctuations, and used in the theory of neutral Fermi superfluids, The relation (46) allows us to express the Green's function (43) defined originally through the fermionic average (43), in terms of averages of bosonic fields, As a result, after the calculation of bosonic averages, we arrive at a remarkably compact where det K GP F is the determinant of the GPF 2 × 2 matrix, det K GP F = K 1,1 K 2,2 + K 2 1,2 , and det K is the determinant of the whole fluctuation matrix K: Consequently, according to (42), the density-density spectral weight function is: In the limit α 0 → 0, this density-density spectral weight function continuously turns to that for a neutral Fermi superfluid.

III. COLLECTIVE EXCITATIONS AT ZERO TEMPERATURE
We now present our results on collective excitations at T = 0, in which case frequencies below 2∆ are free from damping. The high temperature case (T → T c ), which exhibits which is particularly useful when the parameter space (in our case the 4D space spanned by q/k F , 1/k F a, T /T c and ω p /∆) is large. In the past [9,13], it was shown that the poles in the analytic continuation (together with their residues) are usually a good summary of the behavior of the response functions, even in the unconventional case where these pole have a large imaginary part, comparable to their eigenfrequency. However, such highly damped and distorted resonances are not elementary excitations strictly speaking [29].
In Fig. 2, we compare the roots of det K to the contour plots of the spectral functions (restricting ourselves to the diagonal modulus-modulus, phase-phase and density-density responses). We fix the plasma frequency to ω p = 1.5∆, the interaction strength to the BCS regime (1/k F a = −1 or equivalently ∆ ≈ 0.2084E F ), and we vary the pair momentum q to explore the dispersion of the modes. The dispersions are overall similar to what was found in the BCS limit [26] (reminded by the dashed-dotted curves in Fig. 2a).
a. Plasma branch Below the pair-breaking continuum, detK ↓ has an undamped root ω I (black solid curve in Fig. 2a) which produces a Dirac peak in the response functions (overlayed yellow curves in Figs. 2 b, c, d). This mode departs from ω p at q = 0, such that at low q it can be attributed without ambiguity to the plasma mode. We remind that the existence of a plasma mode in q = 0 and ω = ω p is guaranteed by a sum rule [12]. In accordance with the calculation by Anderson [6], this demonstrates that the gapless Goldstone mode disappears (at zero temperature at least) due to long-range Coulomb interactions in a charged superfluid. For the chosen value of ω p , the plasma mode has a positive dispersion but we remind that when ω p approaches 2∆ from below, the dispersion becomes nontrivial, particularly showing a downward dispersion with a minimum at some nonzero q [26]. At larger q, the plasma branch does not enters smoothly into the pair-breaking continuum but splits into multiple peaks both above and below 2∆. In particular, it generates a sharp peak at frequency 2∆ + (fairly constant in function of q), well visible in Fig. 2b, and partially explained by the presence of a root z (2) II of det K ↓ in window II (blue dotted curve in Fig. 2a). In the BCS limit, this root belongs to the density-phase sector (it solves 1,3 = 0), contrarily to the pair-breaking mode z II (red solid curve in Fig. 2a) which belongs to the modulus sector. We note that a multiple resonance also appears in the phase-phase response (Fig. 2c), with a secondary peak visible around ω (pp) 2 at low q. Absent in the neutral case, this peak is a signature of long-range interactions. Remarkably, due to non-vanishing modulus-phase and modulus-density couplings (K 1,2 and K 2,3 respectively), the plasma mode also appears in the modulus-modulus channel, which was not the case in the BCS limit. In the interaction regime considered here, it has a small spectral weight at low q.
b. Pair-breaking branch A pair-breaking ("Higgs") mode (z (1) II or in short z pb ) is found in the analytic continuation through window II. At low-q its dispersion is qualitatively the same as in a neutral Fermi superfluid [9,13], departing quadratically from 2∆. The branch (overlayed in purple in Fig. 2d ) dominates the modulus-modulus response inside the pairbreaking continuum. The most remarkable difference between Fig. 2 and the neutral case [13,26] is the notable rise in the frequency of the modulus mode for q/ √ 2m∆ > 1 (compare the red and blue dashed curves), which we interpret as due to a repulsive interaction with the plasma mode. Since this interaction is carried by the matrix elements K 2,3 (modulusdensity) and K 1,2 (modulus-phase), it vanishes in the BCS limit, such that the modulus mode in the neutral and charged cases coincide in this limit.
To further illustrate the mixing of pair-breaking and plasma mode when their frequencies are in resonance, we show in Fig. 3 the poles and the spectral functions in function of the plasma frequency at fixed q = 0.4 √ 2m∆. At this wavevector, the bare Higgs eigenfrequency (calculated at ω p = 0) is ω pb = 2.1∆. We overlay to the contour plots the "bare" frequency of the plasma branch ω 2 p + c 2 q 2 (c is the sound velocity of the neutral gas) which accurately predicts the location of the resonance away from the interval [2∆, ω (pp) 2 ]. The interaction between plasma and modulus modes reveals itself through the intensity increase of the resonance in the [2∆, ω (pp) 2 ] sector of the modulus-modulus response, clearly visible in Fig. 3 d. Rather than a change in eigenenergy or damping rate (note the relative flatness of ω pb and γ pb , red curves in Figs. 3a and 3e), this increase is caused by the transfer of the spectral weight of the plasma branch (which is small but nonzero in the modulus-modulus channel) to the former modulus mode, causing the emergence of a more intense mixed modulus-plasma mode (note the increase of the residue modulus Z (1) II in Fig. 3 f ) . This enhancement is promising for an experimental observation of modulus "Higgs" collective excitations in charged superfluids and superconductor. ]. Note that this effect is not due to an interaction between the modulus and plasma branches, as it is also observed in the BCS limit [26] where the two branches are decoupled.

IV. VICINITY OF THE PHASE TRANSITION
We now turn to the high-temperature regime |T − T c | ≪ T c which (as in the neutral case [19]) differs much from the zero-temperature case. Fig. 4 shows the dispersion and contour plots at T = 0.99T c and 1/k F a = −1 and ω p = 1.5∆ = 1.75 × 10 −2 E F . We note that having ω p comparable to ∆ close to the phase transition is not the typical experimental situation of superconductors (where ω p is rather fixed in units of E F so large compared to ∆ near T c ), this may be approximately the case in highly-layered materials such as cuprates where at zero-temperature, the plasma frequency in the transverse c-direction is much less than the gap.
To understand the spectrum near T c , one should first notice that density fluctuations decouple from the pair-field fluctuations, in both modulus and phase (in other words, K 1,3 and K 2,3 tends to 0 at T c ). This reflects the situation in the normal phase where the pair susceptibility is decoupled from the density-density response function. Thus, in χ ρρ (Fig. 4 b), the plasma branch appears very close to its normal limit (shown by the overlaid white solid curve in Fig. 4 b). Unlike at T = 0, it is barely sensitive to the structure of the pairbreaking continuum: no resonance splitting is visible around 2∆, ω (pp) 2 or ω 3 . Several roots of det K are associated to plasma modes in the analytic continuation but they reconnect very well at the angular points (compare the green curve ω Ib and the blue curve ω IIb in Fig. 4 a) indicating that they are physically equivalent. Their remains however an important distortion of the plasma branch near ω nothing else than the distortion due to the upper-edge of the particle-hole continuum of the normal phase. This edge is shown by the overlaid white dashed curve in Fig. 4 b and is approximated by ω As it turns into its normal limit, the plasma resonance also looses its spectral weight in the pair field channels (modulus and phase), which is visible in Fig. 4, c,d. Within our approximation, Coulomb interactions do not intrinsically enter into the pair-field propagator K GPF , such that above T c (when the ρ-∆ matrix elements K 1,3 and K 2,3 vanish) this propagator coincides with the pair susceptibility of the neutral and normal Fermi gas with short range interactions. We exclude here the exotic limiting case in which ω p would tend to 0 with |T −T c | faster than ∆ such that the plasma frequency would remain comparable to that of the phononic branches. This susceptibility (studied in [19]) displays a pairing mode with a quadratic dispersion in the regime 1/ξ ≪ q ≪ k F (with ξ ≈ 2m∆/k F the pair coherence length), evolving into a double phononic branch (corresponding to ω (1) Ia and ω (2) Ia in Fig. 4 a) when q ≈ 1/ξ. A modulus mode near 2∆ also survives at wavelengths comparable to the pair coherence length. Due to the decoupling from Coulomb interactions, this low-energy collective physics emerges near T c as it would in the neutral case: except for the residual plasma branch still appearing in the phase response, the contour plots Fig. 4, c and d nearly coincide with their neutral equivalent. At nonzero temperature, the conjecture of Anderson concerning the disappearance of phononic branches in charged systems is thus limited to the density channel: in pair-field channels, a "collisionless second sound" [19] exists even in presence of Coulomb interactions. As can be seen from Fig. 4, more than one gapless collective excitations are predicted in the BCS-BEC crossover regime contrary to the far BCS limit studied in preceding works. The dispersion of the higher-energy Carlson-Goldman mode appears to be close to the phononic-like mode in a neutral Fermi superfluid with the same scattering length [19]. We recall that the other (low-velocity) gapless branch (red curve in Fig. 4 a) is computed here in the collisionless regime, since we assumed that the fermionic quasiparticle have an infinite lifetime. It is however reminiscent of the second gapless mode (gapless modes in charged condensates are also called Carlson-Goldman modes [11]) derived under hydrodynamic assumptions for the quasiparticle lifetime. This was already noticed in [12] for the case of the BCS limit (T c ≪ µ), where this sound branch is visible only in a tiny temperature range below T c . We recall however that the critical behavior [9,30] of the sound velocity is incorrectly predicted by Ref. [12]. Furthermore, we find no trace of the so-called "upper mode" introduced by the authors in the analytic continuation. This mode is an artifact of solving the truncated dispersion equation (see Eq. (2.25) in [12]) Re [det K (ω + i0 + )] = 0 instead of the equation (34). At stronger interactions, we note a broadening of the visibility range.

V. LIMIT OF LARGE PLASMA FREQUENCY
In three-dimensional superconductors (as opposed to layered materials), the concentration of carriers is typically rather high so that the plasma frequency appears several orders larger than the pair-breaking continuum edge. In the above subsections, we explored the regime when the plasma frequency is low enough to be in resonance with pair-breaking and gapless modes. Here, on the other hand, we focus on the regime of large ω p compared to both ∆ and T c . While this regime was previously considered only at q = 0 [6], or by using large but finite numerical value of ω p /T c , we show here how to analytically take the limit ω p → +∞ at fixed q, T and 1/k F a. This regime corresponds to most realistic superconductors and therefore is relevant for contemporary experiments.
In the limit ω p → +∞ and for complex frequencies low compared to the plasma frequency (|z| ≪ ω p ), one should redefine the inverse fluctuation propagator as where the rescaled matrix elements have a finite and nonzero limit when ω p → +∞ (we recall that α 0 = 3πω 2 p /8). The dispersion equation then transforms to Also the spectral weight functions (40), (41), (51) for (ω, ∆) ≪ ω p are determined by the expressions, which are independent of α 0 far from the plasma resonance, as expected. Since the double limit ω p → +∞ and 1/k F a → −∞ (i.e. the BCS limit) was already discussed in the literature [12], we show in Fig. 5 the collective excitations in the strongcoupling regime (more precisely at unitarity 1/k F a = 0). At the considered temperature, the excitation spectrum resembles the neutral spectrum, with both high-and low-velocity gapless branches (respectively black and red lines), a pair-breaking branch around 2∆. The main difference with the neutral case is ω IIb , which, as in the case of ω p finite, is responsible for a peak near 2∆ in the phase and density responses. This peak is absent in the neutral case and therefore characteristic of the charged system. We can see also that the gapless Carlson-Goldman branches can survive at strong coupling in a wider temperature range below T c with respect to the BCS case. Moreover, the lower-energy gapless mode is better resolved at stronger coupling. Fig. 6 shows the temperature evolution of density-density response at a fixed momentum q = 0.2 2m ∆| T =0 for the BCS regime with 1/k F a s = −1 and at unitarity. At low temperature, the peak linked to ω (2) IIb is well visible in the interval [2∆, ω (pp) 2 ]. When approaching T c , this peak disappears and instead a broad feature corresponding to the normal density response develops at energies comparable to ǫ F . FIG. 6: Contour plots of the spectral functions for the density-density response in the low-frequency region with respect to the plasma frequency when varying temperature at a fixed momentum q = 0.2 2m ∆| T =0 in the BCS regime with 1/k F a s = −1 (a) and at unitarity (b).

VI. CONCLUSIONS
In this paper, we investigated theoretically the dispersion and damping of collective excitations in charged and condensed Fermi gases at finite temperatures and in the BCS-BEC crossover. The treatment is performed within the Gaussian pair and density fluctuation (GPDF) method. The spectra of collective excitations are considered using two complementary methods: (i) exploration of the density-density, phase-phase and modulus-modulus response functions, (ii) determination of eigenfrequencies and damping factors of collective excitations from complex poles of the GPDF propagator analytically continued through its branch cuts. Comparison of results obtained by these two methods gives us a reliable and detailed picture of collective excitations.
At zero temperature, two collective excitation branches dominate the spectrum: the plasma mode and the pair-breaking mode. The gapped plasma branch continuously evolves to the gapless phononic branch of the neutral condensate as Coulomb interaction are turned off (that is, as the plasma frequency is sent to zero). When the plasma mode crosses the pairbreaking continuum edge, it exhibits a resonant anticrossing with the pair-breaking mode.
Correspondingly, the magnitude of the modulus-modulus spectral function, greatly increases on resonance, which can facilitate the experimental detection of pair-breaking modes.
At nonzero temperatures, a phononic branch may still exists in the charged system, coexisting with the gapped plasma and pair-breaking branch. This mode describes the motion of a minority of superconducting electrons in a majority of normal carrier. In preceding works devoted to BCS superconductors, it was shown to only existed in a close vicinity of T c . At stronger couplings, in the BCS-BEC crossover, the conditions are substantially more favorable for survival of the gapless mode.
The spectra of collective excitations considered in the present work can be a subject of experimental verification in spectral measurements of the density response of charged superfluids. The pair field response can also be experimentally detected, e. g. by tunneling experiments, similar to the Carlson-Goldman experiment [11]. Our method can easily be transposed to 2D or quasi-2D systems of condensed charged fermions, which makes it promising for the treatment of collective excitations in layered and high-T c superconductors.