Absorption and emission of a collective excitation by a fermionic quasiparticle in a Fermi superfluid

We study the process of absorption or emission of a bosonic collective excitation by a fermionic quasiparticle in a superfluid of paired fermions. From the RPA equation of motion of the bosonic excitation annihilation operator, we obtain an expression of the coupling amplitude of this process which is limited neither to resonant processes nor to the long wavelength limit. We confirm our result by independently deriving it in the functional integral approach using the gaussian fluctuation approximation, and by comparing it in the long wavelength limit to the quantum hydrodynamic result. Last, we give a straightforward application of the coupling amplitude we obtain by calculating the lifetime of the bosonic excitations of arbitrary wave number. We find a mode quality factor that decreases from its maximum at low wave numbers and vanishes when the bosonic branch hits the pair-breaking continuum.

We study the process of absorption or emission of a bosonic collective excitation by a fermionic quasiparticle in a superfluid of paired fermions. From the RPA equation of motion of the bosonic excitation annihilation operator, we obtain an expression of the coupling amplitude of this process which is limited neither to resonant processes nor to the long wavelength limit. We confirm our result by independently deriving it in the functional integral approach using the gaussian fluctuation approximation, and by comparing it in the long wavelength limit to the quantum hydrodynamic result. Last, we give a straightforward application of the coupling amplitude we obtain by calculating the lifetime of the bosonic excitations of arbitrary wave number. We find a mode quality factor that decreases from its maximum at low wave numbers and vanishes when the bosonic branch hits the continuum of fermionic biexcitations.

I. INTRODUCTION
At low temperature, when a quantum fluid is close to its ground state, it is generally described in terms of its elementary excitations or quasiparticles. These eigenmodes of the many-body system remain weakly populated as long as the gas is weakly excited. The quasiparticles are then weakly interacting, so that to first approximation one considers them as an ideal gas. Going beyond this description to include the coupling between quasiparticles is often a theoretical challenge but is essential to account for the ergodicity of the system and the dissipative phenomena such as the damping of the quasiparticles [? ], the loss of macroscopic coherence [? ? ] or the influence of temperature on the viscosity [? ] and the equation of state of the fluid [? ]. At zero temperature, taking the interactions among quasiparticles into account to compute the quantum fluctuations in the elementary modes can also improve the description of the ground state of the fluid [? ].
The quantum fluid we consider is a gas of ultracold spin-1/2 atomic fermions in which van der Waals interactions occur only in the s-wave channel between opposite spin fermions and favor the formation of ↑ / ↓ pairs. Below the critical temperature, those pairs condense into a macroscopically occupied wave function called a pair condensate, and the gas becomes superfluid. This phenomenon is nowadays commonly observed in the laboratory [? ? ? ? ]. The elementary excitations of this system are of two kinds : like in any superfluid, there exists a bosonic branch [? ? ] that in the long wavelength limit describes sound waves [? ] whose energy is proportional to the wave number. Since the elementary constituents of the superfluid are pairs of fermions, there also exists a fermionic branch of excitation of the relative motion of the pairs, with a gapped spectrum [? ]. Several important experimental results on the ground state and the elementary excitations of this system were obtained recently: the measurement of its zero temperature equation of state [? ] and of the dispersion relation of the fermionic [? ] and bosonic [? ? ] excitations, also at high energy [? ]. Experiments can now probe rather precisely the physics of the gas at non-zero temperatures [? ], but its description with analytical methods remains a theoretical challenge. Among the existing approaches to determine the non-zero temperature equation of state we mention a phenomenological model [? ] based on the thermal occupation of the bosonic and fermionic branch, and, at unitarity, effective field theories exploiting the extra symmetries of the system [? ? ]. The description of the interactions among quasiparticles relies on quantum hydrodynamics [? ], a low-energy effective theory limited to the leading order in temperature, except for the coupling between three bosonic excitations [? ]. In this framework, the resonant three-and four-body couplings between phonons [? ? ] and between phonons and fermionic quasiparticles [? ] were obtained.
In this article, we compute the amplitude of the inelastic process of absorption or emission of a collective bosonic excitation by a fermionic quasiparticle within two distinct microscopic approaches leading to the same result : the Random Phase Approximation (RPA) [? ] and the functional integral approach in the gaussian fluctuation approximation [? ? ]. Our result is in agreement with quantum hydrodynamics in its validity domain but is not limited to resonant processes nor to the long wavelength limit. We present a straightforward application: the computation of the damping rate of the collective excitations of arbitrary wave number due to their inelastic coupling to the fermionic quasiparticles in the collisionless regime. This rate is the main contribution to the lifetime of the collective excitations at the usual experimental temperatures [? ] when the concave collective branch [? ] forbids the Beliaev-Landau damping mechanism [? ]. It can be measured by Bragg spectroscopy, a technique recently applied to Fermi gases [? ]. With the microscopic expression of the coupling between three bosonic excitations obtained by Ref. [? ] and the present result, all three-body couplings between quasiparticles are known microscopically, which paves the way to a study of the dressing of quasiparticles by the interactions among them.

II. COLD FERMI GAS
We consider a gas of neutral atomic fermions of mass m interacting via the van der Waals force. The atoms are equally distributed in two internal states labelled by ↑ and ↓ and evolve in a homogeneous cubic space of volume L 3 . If the gas is cold and dilute enough, the atoms interact significantly only two by two, in the s-wave channel, which means that the interactions among same spin fermions are negligible. In this regime, the exact form of the potential does not influence the macroscopic physics, allowing us to select an effective contact potential V (r, r ′ ) = g 0 l 3 δ r,r ′ where δ is the Kronecker symbol, and where, rather than introducing a cutoff in Fourier space, we choose to discretise space into a cubic lattice of step l. The bare coupling constant g 0 is related to the s-wave scattering length a, the experimentally accessible parameter, through the renormalisation relation [? ] In second quantisation, the grand canonical Hamiltonian of the system readŝ H = l 3 r,σ=↑/↓ψ † σ (r) where P is the momentum operator, whose eigenfunctions e ik·r have eigenvalues k.

III. BCS THEORY AND FERMIONIC QUASIPARTICLES
In this section, we recall how the fermionic excitation branch is described in the Bardeen-Cooper-Schrieffer (BCS) theory [? ]. We introduce the BCS Hamiltonian (4) where in the interaction term we have replaced the quadratic quantum field g 0ψ↓ψ↑ by a self-consistent order parameter ∆. The resulting quadratic Hamiltonian is readily diagonalised in Fourier space by the Bogoliubov rotation: where D = 2π L Z 3 ∩ [−π/l, π/l[ 3 is the set of wave vectors of the first Brillouin zone compatible with the periodic boundary conditions. This change of basis introduces fermionic operatorsγ † kσ (σ =↑ / ↓) that create excitations whose energy ǫ k depends on ∆ and on the dispersion relation of the normal gas ξ k = 2 k 2 /2m − µ: The U k and V k coefficients of the Bogoliubov rotation are given by and the self-consistent relation giving the value of ∆ is The ground state of the BCS Hamiltonian (4) is a state in which all fermions are paired with a probability |V k | 2 of finding the k ↑ /−k ↓ pair. We have introduced the Fourier transform of the fermionic field operatorâ kσ = 1 L 3/2 k ψ σ (r)e ik·r . The action of the operatorγ † kσ on the BCS ground state (9) is to destroy the kσ/ − kσ ′ pair and replace it by an unpaired kσ fermion. This is why the fermionic excitations are referred to as pair-breaking excitations. In terms of thosê γ excitations, the BCS Hamiltonian is diagonal

IV. DESCRIPTION OF THE BOSONIC EXCITATIONS IN THE RPA
The idea of the RPA is to linearise, in the Heisenberg picture, the equations of motion of the quadratic fermionic field operatorsψ † σψ σ ,ψ † ↑ψ † ↓ andψ ↓ψ↑ by performing incomplete Wick contractions on the quartic terms: abĉd →âb ĉd + âb ĉd −âĉ bd − âĉ bd +âd bĉ + âd bĉ − âb ĉd + âĉ bd − âd bĉ (11) where the average value . . . is taken in the BCS ground state (9). Contrarily to Anderson [? ], we choose to express the resulting linear system in the basis of the fermionic quasiparticle operatorsγ, by introducing the variableŝ In these notations, the q coordinate is the centre-of-mass wave vector of the pair of quasiparticles and k is the wave vector of its relative motion. In this basis, the RPA equations take the remarkably simple form [1] Before giving the expressions of the source termsŜ q k and Y q k , let us briefly comment the linear system (14-16). The choice of writing the equations in Fourier space has decoupled operators with a different center-of-mass wave vector q. Then the choice of singling out the quasiparticle basis has decoupled the evolution (16) of the quasiparticle-hole operatorsγ †γ , and reduced it to the trivial one given by the BCS Hamiltonian (10). These operators still enter into the equations of motion ofŝ q k andŷ q k as source termŝ where we complete the notations introduced in (17) We temporarily put these source terms aside to focus on the homogeneous system inŝ q k andŷ q k . This system coincides exactly with the one obtained by Ref.
[? ] in a semi-classical approach [? ]. We briefly recall the results obtained by this reference on the collective excitation branch. The eigenenergy of the collective mode is determined by the implicit equation [? ] where the I σσ ′ are sums over the internal degrees of freedom of the pairs At low q, this energy is phononic with a speed of sound c given by the hydrodynamic formula mc 2 = ρdµ/dρ. The bosonic operators associated to the collective eigenmodes are expressed in terms of the fermionic quasiparticle-pair operatorŝ with the coefficients equivalent to those (44,45) of Ref. [? ] given that that sets the value of N q . Note that in the RPA we directly obtain quantumb q operators, without need of the quantisation procedure described in Refs.

V. COUPLING BETWEEN BOSONIC AND FERMIONIC QUASIPARTICLES IN THE RPA
With the work we have done in the previous section on the RPA equations, it is rather straighforward to derive the coupling between bosonic and fermionic quasiparticles. It is in fact contained in the source terms (18) and (19). To see this, we contract the RPA system of equations (14) and (15) to form the equation of motion of the bosonic operators (25) The first term in the right hand side of this equation represents the free evolution of the bosonic quasiparticle annihilation operator, at an angular frequency ω q . The second describes its coupling to the fermionic operatorŝ γ, with a coupling amplitude [2] that we choose independent of the system size L 3 . We interpret the equation of motion (29) as an Heisenberg equation derived from the fictitious Hamiltonian where the operatorsb q are bosonic and commute with the fermionic operatorsγ kσ . The amplitude A kq now clearly appears as the coupling amplitude of the non-linear process represented in Fig. 1, where a bosonic excitation is absorbed or emitted by a fermionic quasiparticle. In Sec. VII, we use our knowledge of A kq to compute the limitation of the bosonic quasiparticle lifetime caused by the process of Fig. 1. Before that, we explain in the following section how our result (30) can be obtained in a totally different approach based on functional integration.

VI. COMPARISON TO THE FUNCTIONAL INTEGRAL APPROACH
In this section we explain how the coupling between bosonic and fermionic quasiparticles can be obtained in the functional integral formalism [? ? ? ] applied to Fermi gases [? ]. The starting point is to introduce the action where the imaginary time τ = it varies from 0 to 1/k B T , the Grassmann fields ψ σ andψ σ are the analogues of the quantum field operatorsψ σ andψ † σ respectively and H[ψ ↑ , ψ ↑ ,ψ ↓ , ψ ↓ ] is the Hamiltonian (3) where the quantum fields have been replaced by their Grassmann equivalent. The partition function is expressed as a functional integral of this action where the integral Dψ spans all the possible configurations of the Grassmann fields. This functional integral of a quartic action cannot be performed analytically, forcing us to look for an approximation. Two disctinct yet equivalent approaches have been proposed: (i) In the gaussian fluctuation approach, one introduces an auxiliary complex field ∆(r, τ ), integrates out the where z q is the eigenenergy and M , sometimes called the gaussian fluctuation matrix, is a 2 × 2 matrix whose coefficients, after correction [? ] of the sign mistake in Ref. [? ], are given by M 22 (z, q) = M 11 (−z, −q), M 12 = M 21 and where we permit ourselves the short-hand notation a ± = a k± q 2 and we introduce the Fermi-Dirac occupation numbers of the fermionic quasiparticles Notice that in (35)-(36) the variable z can take in principle only the pure imaginary values corresponding to the Matsubara frequencies z = i 2πnk B T , n ∈ N. To relate the functional integral approach to the RPA, we do a small change of variables These new sums coincide at zero temperature with those we introduced in the RPA [3] Moreover, they fulfill a similar eigenvalue equation We chose the variables (38) and (39) also for a practical reason, as they substantially simplify the long wavelength calculations (see Sec. VII A). Since Eq.(41) has no solution on the imaginary axis, it is natural to extend the definition of the M σσ ′ functions to the complex plane by setting where the factor 2 is chosen so that Γ corresponds to the damping rate of the mode. This is not a problem at zero temperature (that is when n F k = 0, ∀k ∈ D) and this leads to the eigenvalue equation (21) already encountered in the framework of the RPA, which for fixed q possesses a unique real solution ω q < min k (ǫ k+ q 2 + ǫ k− q 2 ). At non zero temperature, when n F k = 0, ∀k ∈ D, the situation is more complicated: the function z → 2 has a branch cut along the whole real axis, corresponding to the continuum k → ǫ k+ q 2 − ǫ k− q 2 (which in the RPA we view as the eigenenergy of the quasiparticle-hole operatorγ † k− q 2γ k+ q 2 and which clearly reaches zero 0 when q = 0) appearing in the denominators of the last lines of (38) and (39). This function has then a discontinuity when crossing the real axis and no root outside of it : Eq.(41), as such, has no solution.
However, as in the case of a discret state coupled to a continuum [? ], the inverse function z → (M ++ (z, q)M −− (z, q) − z 2 [M +− (z, q)] 2 ) −1 retains the memory of its pole at zero temperature in the shape of its variations in the vicinity of the branch cut. In the perturbative regime, when the coupling to the continuum is weak, these variations are fast around z = ω q and denote the existence of a pole z q in the analytic continuation of the inverse function across the branch cut, for Imz < 0. Here we restrict to this perturbative regime and look for the new solution z q of the analytically continuated Eq.(41) in the vicinity of the zero temperature solution ω q : Clearly, this regime is reached for sufficiently low temperatures. More precisely, we assume that the deviations to the zero temperature spectrum behave as We justify this hypothesis (which the final result will also confirm) by the fact that M σσ ′ depends on temperature through the fermionic occupation numbers n F k ≃ e −ǫ k /kBT at low temperature. To expand M σσ ′ in powers of |z q − ω q |, we split the contributions of the second and third lines of (38) and (39). The terms of the second lines have no branch cut that reaches 0, and are therefore straightforwardly continued to the lower half-plane : we simply expand in powers of |z q − ω q | inside the summation. The terms of the third lines are those whose branch cut reaches 0 but they are already of order O(e −β∆ ) (because of the n F k+ q 2 − n F k− q 2 factor): it is enough to obtain their value for |z q − ω q | → 0, which we do by approaching the branch cut from above setting z = ω q + iη with η → 0 + . In this article, we focus on the imaginary part Γ q = O(e −β∆ ) of the correction to the spectrum. To obtain it, it is enough to expand the real and imaginary parts of M σσ ′ up to first non-zero order in e −∆/kBT : where in Eq. (46) the K σσ ′ and J σσ ′ contributions stem respectively from the second and third lines of (38) and (39). Explicitly, we have Inserting the expansions (45) and (46) in the eigenvalue equation (41) yields the expression of the damping rate where all the J σσ ′ and I σσ ′ functions are evaluated in (ω q , q). To identify the normalisation constant N q of the RPA, we have used the relation which is shown in the framework of the RPA by replacing expressions (26) and (27) in the normalisation condition (28) and then using the zero temperature eigenvalue equation (21). We conclude this calculation by replacing in (51) the J σσ ′ functions by their expressions (49) and (50) This expression of Γ q makes an elegant connection with the RPA expression of the coupling amplitude A kq : it is the damping rate one obtains by applying the Fermi golden rule to the population of the bosonic excitations of wave vector q using the RPA Hamiltonian (31) and adding up the damping rates due to the ↑ and ↓ fermionic excitations, here identical. Expression (53) is valid only to leading order in e −∆/kBT , so we should keep only the leading term in the fermionic occupation numbers Our perturbative approximation therefore assumes that the gas of fermionic quasiparticles is non-degenerate.

VII. PHONON DAMPING
In this section, we compute an explicit expression of the phonon damping rate (53) due to the coupling to the fermionic quasiparticles. We first concentrate on the long wavelength limit, in which our result can be compared to other approaches.

A. Long wavelength limit
In the long wavelength limit, q → 0 (that is both 2 q 2 /m∆ → 0 and ω q /k B T → 0), we first simplify in Eq. (53) the argument of the Dirac delta function expressing energy conservation (55) where we denote the cosine of the angle between k and q by There exists a value of u satisfying the energy conservation constraint (55) if and only if Physically, this constraint means that the absorption of a phonon is allowed only if the slope dǫ k dk = ξ k ǫ k k m of the fermionic branch is greater in absolute value than the speed of sound c. This condition may be interpreted as a Landau criterion if we view the fermionic quasiparticle as an impurity travelling in the superfluid at the group velocity dǫ k dk . Eliminating the norm of k in favor of ξ = ξ k /∆ with the change of variable ξ k = 2 k 2 /2m − µ, this constraint takes a polynomial form Since mc 2 is a function of µ/∆ through the BCS equation of state [? ], the coefficients of this polynomial depend The integration domain X over the variable ξ = ξ k /∆ is shown in the long wavelength limit as a function of µ/∆. The BEC limit is found at µ/∆ → −∞ and the BCS limit at µ/∆ → +∞. Without any constraint, the integration domain is physically restricted to ξ > −µ/∆, that is outside the hatched region. For all values of µ/∆, there exists an integration domain X1 ⊂ R + (upper part of the graph) in the increasing part of the BCS branch; this interval is bounded from below by ξ1 (full orange line) and not bounded from above. On the BCS side, for µ/∆ > mc ≃ 2.45 (or 1/kFa < −0.594), there exists a second interval X2 ⊂ R − in the decreasing part of the BCS branch, bounded from below by ξ2 (full blue line) and from above by ξ3 (full red line). Note that the vicinity of the minimum of the BCS branch (at ξ = 0) is always outside the integration domain.
in fact on only one parameter. In Fig. 2 we show as a function of µ/∆ the set X of values of ξ for which the energy conservation constraint can be fulfilled. On the Bose-Einstein Condensate (BEC) side µ/∆ → −∞ where the BCS branch k → ǫ k is a strictly increasing function of k, X is of the form [ξ 1 , +∞[ with a lower bound ξ 1 > 0. When µ > 0, the BCS branch is a decreasing function for k < (2mµ/ 2 ) 1/2 ; when the slope of the branch in this region becomes larger than the speed of sound c, that is for µ/∆ = 2.45, there appears a second contribution to X , of the form [ξ 2 , ξ 3 ] with ξ 2 < ξ 3 < 0. We then expand in powers of q expression (30) of the coupling amplitude [4] as well as the factor containing the occupation numbers Last, we compute the rate Γ q in the thermodynamic limit d 3 k, in a spherical frame with z-axis along q, and using the Dirac delta function to integrate out the colatitude u (61) We have introduced the notation ǫ = ξ 2 + 1 and the Fermi velocity v F related to the density by 3π 2 ρ = (mv F / ) 3 . The ratio Γ q /ω q tends to a non-zero constant (t 2 −1) 2 . In the BEC limit the dashed lines indicate the oblique asymptote for 1/kFa To the left of the vertical dotted line (for 1/kFa < −0.594) the domain X contains a contribution of the decreasing part of the BCS branch.
at low q. This allows us to represent on Fig. 3 the long wavelength limit of Γ q /ω q as a function of the interaction strength for different values of temperature. As expected, this value is exponentially small in temperature because of the e −∆ǫ/kBT factor in the integrand. At very low temperature, this exponential factor is extremely peaked around its maximum for ξ 1 = min{|ξ|, ξ ∈ X } (full orange line on Fig. 2), and we can replace the non exponential factor in the integrand by its value in ξ = ξ 1 to obtain In this expression there appears an effective gap ∆ ′ = ∆ ξ 2 1 + 1 that reflects the fact that the absorptionemission process is not resonant for excitations located at the minimum of the fermionic branch.

B. Comparison to other approaches
Several other methods exist to obtain the coupling amplitude A kq in the long wavelength limit. In their description of superfluid helium-4 with quantum hydrodynamics, Landau and Khalatnikov [? ] derive an expression for the coupling between phonons and rotons by treating the phonons as a semiclassical hydrodynamic perturbation acting on a generic gapped roton Hamiltonian. Their expression can be generalized [? ] to ultracold atomic Fermi gases where the fermionic branch k → ǫ k plays the role of the rotons: This expression for A coincides with our result (59) for resonant processes only (i.e. processes obeying the energy conservation relation (55) which is used to eliminate the angular variable u), as expected for an effective low energy theory such as quantum hydrodynamics. This is suitable for a calculation involving on-shell energies, such as that of the damping rate Γ q . However, for quantities depending on off-shell coupling amplitudes, such as energy shifts of the dispersion relation, the hydrodynamic result (63) cannot be used. Note that Ref.
[? ] uses the quantum hydrodynamic theory to calculate the damping rate of phonons due to the scattering of phonons on BCS excitations. This is the dominant process at the lowest temperatures because it can resonantly couple to the excitations at the bottom of the fermionic branch: the associated damping rate behaves as e −∆/kBT , without an effective gap as in (62). However, in practice it can be much smaller than the rate Γ q at the experimentally accessible temperatures. Obtaining four-body couplings (between two BCS excitations and two collective excitations) in a microscopic theory requires us to extend our description beyond the RPA, or beyond the gaussian fluctuation approximation in the functional integral framework, and this is beyond the scope of the present paper.
Also the functional integral approach has been used previously to calculate the damping rate Γ q in the long wavelength limit [? ], and a different result was reported than what we obtain in the present treatment, Eq. (61). One of the reasons for this difference is that in Ref. [? ] the coefficients K σσ ′ arising from the non-singular part of M σσ ′ are neglected, based on the fact that the denominator ω q ± (ǫ k+ q 2 + ǫ k− q 2 ) is not resonant. Here we show that this contribution cannot be neglected, as it is essential to retrieve the RPA result and modifies significantly the value of Γ q , also at low q.
C. Beyond the long wavelength limit Next, we turn to the general result (53) beyond the long wavelength limit. As mentioned above, this limit requires two conditions to be met: cq/k B T ≪ 1 and 2 q 2 /m∆ ≪ 1. Relaxing the first condition is straightforward as only the expansion of the fermionic occupations numbers (60) uses the inequality cq/k B T ≪ 1. Using energy conservation and the relations (1 − n F k )/n F k = e βǫ k and (1 + n B q )/n B q = e β ωq for fermionic and bosonic occupation numbers respectively, one finds [? ]: with the bosonic occupation numbers This allows us to write (61) for arbitrary values of cq/k B T as (66) Relaxing the remaining condition 2 q 2 /m∆ ≪ 1 is more difficult, as no analytic expression is available for the collective branch q → ω q outside the long wavelength limit. Similarly as in Refs. [? ? ], we numerically solve ω q from the implicit equation (21) and use this result to evaluate the integral in Eq. (53). Details on how to satisfy the energy conservation constraint can be found in Appendix A.
In figures 4 and 5 we plot the inverse quality factor Γ q /ω q as a function of wave number q, for several interaction strengths, at a temperature ∆/k B T = 5. For all interaction strengths the inverse quality factor decreases quadratically at low q, starting from its limiting value (61) for q → 0. Figure 4 shows the damping rate in the BCS regime (1/k F a = −0.47) and at unitarity (1/k F a = 0). For these two cases, the collective branch q → ω q disappears when hitting the continuum of fermionic biexcitations k → ǫ k+ q 2 + ǫ k− q 2 at q/(2mµ) 1/2 ≃ 2.0 and 2.3, respectively. In these points, the inverse quality factor becomes vanishingly small along a vertical tangent. The suppression of Γ q is due to a divergence in the normalization constant N q (52) which, in turn, comes from the denominators ( ω q − ǫ k+ q 2 − ǫ k− q 2 ) 2 in the integral quantities K σσ ′ (47,48), which becomes a second order pole as ω q enters the pair breaking continuum. Physically, this suppression relates to the fact that when the branch approaches the continuum, the emission and absorption processes that are considered here are dominated by the level repulsion exerted by the continuum. Figure 5 shows the interesting case, on the BEC side, that occurs when the the collective branch q → ω q is not singly connected [? ], but consists of two pieces [0, q sup ]∪[q inf , +∞[ with q sup < q inf , as shown in the inset. Now Γ q goes to zero in the two points q sup , q inf where the collective branch reaches the continuum. For large q the collective excitations acquire an energy 2 q 2 /4m equal to the kinetic energy of a molecule of mass 2m and momentum q. In Ref. [? ] this is interpreted as a tightly bound molecule that is ejected out of the pair condensate. These high energy excitations can only be absorbed by fermionic quasiparticles whose wave number k is large enough so that the maximal energy difference for aligned wavevector ǫ k+ q 2 − ǫ k− q 2 exceeds ω q . This explains the strong decrease of the inverse quality factor for large q.

D. Towards an experimental observation
Last, we discuss the experimental observability of the damping rate Γ q we predict. Reaching the low temperature regime where our theory is applicable is not anymore a serious limitation: experiments can reach temperatures of a few hundredths [? ] of the Fermi temperature T F , which, choosing for the gap the value ∆/k B T F ≃ 0.69 pre- ] when it reaches the lower edge of the pair-breaking continuum (dotted curve in the inset) at qsup/(2mµ) 1/2 ≃ 3.0 and reappears as it re-emerges from the continuum at q inf /(2mµ) 1/2 ≃ 4.8, after which it stays close to the continuum edge. The horizontal dotted line shows the analytic q → 0 limit. The damping rate Γq goes to zero at qsup and q inf , and remains small in the second part of the branch (q > q inf ).
dicted by BCS theory at unitarity (the value measured by Ref. [? ] ∆/k B T F ≃ 0.44 is of the same magnitude), corresponds to values of ∆/k B T of the order of 10. Similarly, typical experiments have a large enough precision to detect a lifetime of order 1/Γ q ; for example, choosing q/(2mµ) 1/2 = 0.5, 1/a = 0 and ∆/k B T = 5 (parameters of the full line in Fig.4), we have ω q ≃ 0.5∆ and Γ q /ω q ≃ 0.015, which, taking for the Fermi temperature the typical value T F ≈ 1 µK, corresponds to a lifetime 1/Γ q ≈ 1.5 ms. Much longer lifetimes (and quality factors much larger than ω q /Γ q ≃ 67) have been observed in experiments on the low-energy modes of a cold Bose gas [? ]. Similar measurements of the frequency broadening of the collective mode have been performed recently in paired Fermi gas using Bragg excitations [? ], but their precision is limited by the spatial inhomogeneity of the gas. We are confident that they will soon reach the precision required to measure Γ q thanks to the many ameliorations offered by the flat-bottom potentials [? ].
To conclude on the observability of the damping of the collective excitations by absorption-emission by the fermionic quasiparticles, we still need to compare it to the others three-or four-body damping mechanisms acting on the collective excitations, that are: (i) the Landau-Beliaev processes between three collective excitations [? ? ], (ii) the Landau-Khalatnikov process between four collective excitations [? ? ] and (iii) the process of scattering of a collective excitation on a fermionic quasiparticle [? ]. When allowed Beliaev-Landau damping is the dominant phenomenon. For example, keeping q/(2mµ) 1/2 = 0.5, 1/a = 0 and ∆/k B T = 5, we obtain Γ Beliaev−Landau q /ω q ≃ 0.14; the ratio Γ q /ω q of the process we are interested in is thus one order of magnitude smaller, yet we believe still measurable with a good precision. In fact the situation is even more favorable: the collective branch in superfluid Fermi gases has the particularity of being concave at low wave number on the BCS side [? ? ] (for 1/k F a < −0.14 according to RPA) which energetically forbids Beliaev-Landau damping. In this case, only processes (ii) and (iii) remain, which, as far as the damping of collective excitations is concerned, are dominated by the process of absorption-emission by aγ quasiparticle over a large temperature domain [? ]. Taking 1/k F a = −0.47 (parameters of the dashed line in Fig.4, in the near BCS regime where Beliaev-Landau damping is still forbidden) and still ∆/k B T = 5 and q/(2mµ) 1/2 = 0.5, we obtain 8,1×10 −3 for the inverse quality factor of the absorptionemission process, 2.9×10 −4 for Landau-Khalatnikov process (ii) and 6.6 × 10 −3 for the scattering process (iii). The process studied in this article is therefore the main limitation to the lifetime of the collective excitations on the BCS side of BEC-BCS crossover at the usual experimental temperatures. We are thus optimistic about an experimental observation in the near future.

VIII. CONCLUSION
The coupling amplitude for the absorption or emission of a collective excitation by a fermionic quasiparticle in a Fermi superfluid is obtained using two complementary methods, namely the Random Phase Approximation and the functional integral formalism. In the limit of long wavelengths our result agrees with the prediction of quantum hydrodynamics for resonant processes. This result is a necessary ingredient in order to take interactions between quasiparticles into account in any description of a Fermi superfluid at low (but non-zero) temperatures. We also calculate the damping rate of collective excitations, resulting from the coupling to fermionic excitations beyond the long wavelength limit in the collisionless regime. We show that the inverse quality factor decreases (at first quadratically) as the wave number of the collective mode increases, and vanishes when the collective mode rejoins the continuum of fermionic biexcitations. The damping rate we compute can be measured directly using Bragg spectroscopy on the BCS side of the crossover where the collective branch is concave, at the temperatures currently reached by the experiments. The main difficulty in evaluating the integral in expression (53) of the damping rate consists in finding the integration domain K over the norm of k. The energy conservation constraint ω q + ǫ k− q 2 − ǫ k+ q 2 = 0 can be rewritten in dimensionless form A2) and the dimensionless quantities with x > 0, m > 0, x q > 0 and u ∈ [−1, 1]. The function E is taken to depend explicitly on x and u which represent the integration variables linked to the vector k, but not on m and x q which are constant parameters fixed by the interaction strength and the wave number q respectively. The dimensionless eigenenergy ω q /∆ > 0 is connected to x q via the zero-temperature dispersion relation for the collective mode q → ω q , obtained by numerically solving (21). Our approach consists in finding the value of u that satisfies (A1) for fixed x,m and x q . This value, if it exists, is unique since E is a monotonic function of u: For x = 0 the constraint (A1) can never be satisfied, so this case can be eliminated. The case x q = 0 was discussed in section VII A. The remaining condition m = x 2 + x 2 q /4 is independent of u so that dE/du is either identically zero or has a constant sign for u ∈ [−1, 1]. More precisely, Moreover, since we have E(x, 0) = ω q /∆, ∀x > 0, the function u → E(x, u) has a root in [−1, 1] if and only if These conditions determine the set K. Close to x = m − x 2 q /4 neither condition is met, because in this point E(x, u) is constant and strictly positive. This implies that the integration domain splits up as in the low q limit into q /4, the upper do-main K 1 is of the form [x 1 , ∞[ with a lower bound x 1 that can be determined by numerically solving E(x, −1) = 0 m − x 2 q /4. Contrarily to K 1 , it does not exist for all values of m and x q (it clearly disappears for m − x 2 q /4 < 0). To know whether it exists, we look for the minimum x min of E(x, 1) in ]. If at that point the function is of negative sign E(x min , 1) < 0, we look for the boundaries of K 2 such that x 2 < x min < x 3 . Otherwise K 2 is the empty set.
[1] Here we neglected the Wick contractions in the Hartree- ↑ whose contributions to the spectrum of the bosonic excitations [? ] and to the coupling between quasiparticles [? ] vanish when the continuus space limit l → 0 is taken at fixed scattering length a.
[2] To establish Eq.(30) we have used the relations that are demonstrated from the definitions (26) and (27) of M kq and N kq using the self-consistent relation on the order parameter (8) to recognize the sums (22) and (23).
[3] At zero temperature we analytically continue the M σσ ′ functions to the real axis by setting z → ω, which we can do without trouble as long as ω stays below the gapped k → ǫ k+ q 2 + ǫ k− q 2 continuum.