H Kurkjian

Y Castin

A Sinatra

Three-phonon and four-phonon interaction processes in a pair-condensed Fermi gas

We study the interactions among phonons and the phonon lifetime in a pair-condensed Fermi gas in the BEC-BCS crossover in the collisionless regime. To compute the phonon-phonon coupling amplitudes we use a microscopic model based on a generalized BCS Ansatz including moving pairs, which allows for a systematic expansion around the mean field BCS approximation of the ground state. We show that the quantum hydrodynamic expression of the amplitudes obtained by Landau and Khalatnikov apply only on the energy shell, that is for resonant processes that conserve energy. The microscopic model yields the same excitation spectrum as the Random Phase Approximation, with a linear (phononic) start and a concavity at low wave number that changes from upwards to downwards in the BEC-BCS crossover. When the concavity of the dispersion relation is upwards at low wave number, the leading damping mechanism at low temperature is the Beliaev-Landau process 2 phonons ↔ 1 phonon while, when the concavity is downwards, it is the Landau-Khalatnikov process 2 phonons ↔ 2 phonons. In both cases, by rescaling the wave vectors to absorb the dependence on the interaction strength, we obtain a universal formula for the damping rate. This universal formula corrects and extends the original analytic results of Landau and Khalatnikov [ZhETF 19, 637 (1949)] for the 2 ↔ 2 processes in the downward concavity case. In the upward concavity case, for the Beliaev 1↔ 2 process for the unitary gas at zero temperature, we calculate the damping rate of an excitation with wave number q including the first correction proportional to q 7 to the q 5 hydrodynamic prediction, which was never done before in a systematic way.

Introduction

In several many-body systems the low-lying collective excitations are phonons. At low temperature, interactions among phonons determine their lifetime, correlation time and mean-free path. Therefore, they play a central role in transport phenomena, such as thermal conduction in dielectric solids, in hydrodynamic properties, such as temperature dependent viscosity and attenuation of sound in liquid helium [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF][START_REF] Khalatnikov | Relaxation Phenomena in Superfluid Helium[END_REF][START_REF] Maris | Phonon Dispersion and the Propagation of Sound in Liquid Helium-4 below 0.6K[END_REF][START_REF] Sherlock | The angular spreading of phonon beams in liquid helium 4: upward phonon dispersion[END_REF][START_REF] Damian | Pressure dependence of phonon interactions in liquid 4 He[END_REF][START_REF] Smith | Direct measurements of phononphonon scattering in liquid He 4[END_REF], and in macroscopic coherence properties of degenerate gases, as they determine the intrinsic coherence time of the condensate both in bosonic and fermionic gases [START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF][START_REF] Kurkjian | Brouillage thermique d'un gaz cohérent de fermions[END_REF]. The dominant phonon decay channel differs among physical systems, and depends in particular on the curvature of the phonon excitation branch [START_REF] Tucker | Four-phonon scattering in superfluid 4 He[END_REF][START_REF] Adamenko | Theory of scattering between two phonon beams in superfluid helium[END_REF][START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]. For a convex dispersion relation at low wave-number, Beliaev-Landau 2 ↔ 1 processes involving three quasiparticles dominate [START_REF] Beliaev | Application of the Methods of Quantum Field Theory to a System of Bosons[END_REF][START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF], while for a concave dispersion relation at low wave-number, Beliaev-Landau interactions are not resonant, and the Landau-Khalatnikov 2 ↔ 2 processes involving four quasiparticles dominate the dynamics.

Compared to other many-body systems, cold atomic gases offer the unique possibility to control and tune some microscopic parameters, in particular the interaction strength. For a cold atomic Fermi gas in two spin states ↑ and ↓, interactions occur only in the s-wave between fermions of opposite spins. We restrict ourselves to the universal regime of interactions of negligible range b, fully characterized by a single parameter, the s-wave scattering length a, |a| b. Using an external magnetic field, one can adjust this parameter of the interaction strength close to a socalled Feshbach resonance. This feature has allowed cold Fermi gases experiments to study the crossover between the Bardeen-Cooper-Schrieffer (BCS) regime 1/a → -∞, where the superfluid pairs (k ↑, -k ↓) are localized in Fourier space close to the Fermi surface, and the Bose-Einstein Condensate (BEC) regime 1/a → +∞, of dilute tightly bound dimers that behave as bosonic particles [START_REF] O'hara | Observation of a Strongly Interacting Degenerate Fermi Gas of Atoms[END_REF][START_REF] Bourdel | Measurement of the Interaction Energy near a Feshbach Resonance in a 6 Li Fermi Gas[END_REF][START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF][START_REF] Bartenstein | Collective Excitations of a Degenerate Gas at the BEC-BCS Crossover[END_REF][START_REF] Zwierlein | Condensation of Pairs of Fermionic Atoms near a Feshbach Resonance[END_REF][START_REF] Zwierlein | Vortices and superfluidity in a strongly interacting Fermi gas[END_REF][START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF][START_REF] Sidorenkov | Second sound and the superfluid fraction in a Fermi gas with resonant interactions[END_REF]. For the unpolarized gas of spin-1/2 fermions, the excitation spectrum consists of two branches: a gapped fermionic branch describing the excitation of the internal degrees of freedom of the ↑↓ pairs, and a bosonic branch describing the collective motion of the pair center of mass, which has a phononic behavior at small wave vectors [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF][START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF][START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF][START_REF] Klimin | Pair Excitations and Parameters of State of Imbalanced Fermi Gases at Finite Temperatures[END_REF][START_REF] Randeria | Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the Unitary Fermi Gas[END_REF][START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]. If the density of excitations is sufficiently small, the elementary excitations are long-lived weakly interacting quasiparticles.

In cold gases of spin-1/2 fermions, the dispersion relation of phonons changes from convex in the BEC regime [START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF] to concave in the BCS regime [START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF]; it is generally expected at unitarity to be convex at low wave number [START_REF] Salasnich | Extended Thomas-Fermi density functional for the unitary Fermi gas[END_REF][START_REF] Mañes | Effective theory for the Goldstone field in the BCSBEC crossover at T = 0[END_REF][START_REF] Rupak | Density functional theory for non-relativistic fermions in the unitarity limit[END_REF][START_REF] Salasnich | Viscosity-Entropy Ratio of the Unitary Fermi Gas from Zero-Temperature Elementary Excitations[END_REF][START_REF] Klimin | Pair Excitations and Parameters of State of Imbalanced Fermi Gases at Finite Temperatures[END_REF][START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF], although the predicted values of the curvature parameter significantly differ and there is no experimental verdict yet. A complete experimental investigation of the interactions between phonons, in both the concave and convex cases, seems in any case possible within this same physical system. In particular, experiments should observe a sharp increase in the phonon lifetime by switching from the convex to the concave case through a tuning of the scattering length. Experimental studies of collective excitations in cold gases have been performed both in the spectral domain with Bragg scattering [START_REF] Steinhauer | Excitation Spectrum of a Bose-Einstein Condensate[END_REF][START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF] and in the time domain [START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF][START_REF] Kinast | Damping of a Unitary Fermi Gas[END_REF]. Those studies combined with the possibility to realize homogeneous systems in flat-bottom traps [START_REF] Gaunt | Bose-Einstein Condensation of Atoms in a Uniform Potential[END_REF] make the measurement of the decay caused by either the 2 ↔ 1 or the 2 ↔ 2 processes a concrete prospect.

To date, an exhaustive theoretical study of phonon interactions and a general expression of the phonon damping rates at low temperature for any interaction strength is missing. The most natural example of a convex dispersion relation is the Bogoliubov excitation branch of the weakly interacting Bose gas. In this case, Beliaev obtained the coupling amplitudes in a microscopic framework and, from those, the damping rate associated to the 1 → 2 processes [START_REF] Beliaev | Application of the Methods of Quantum Field Theory to a System of Bosons[END_REF]. The contribution of the 2 → 1 Landau processes, which exist only at nonzero temperature, was obtained later in [START_REF] Hohenberg | Microscopic theory of superfluid helium[END_REF][START_REF] Szépfalusy | On the dynamics of continuous phase transitions[END_REF][START_REF] Pitaevskii | Landau damping in dilute Bose gases[END_REF]. In superfluid spin-1/2 Fermi gases, reference [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF] calculated the zero-temperature Beliaev 1 → 2 damping rate by a phenomenological low-energy effective theory going one step beyond quantum hydrodynamics in order to include a curvature in the phonon dispersion relation. However, the authors omitted to include corrections to the phonon coupling amplitude, which makes their treatment inconsistent. In the bosonic case of phase II of 4 He, whose dispersion relation was originally believed to be concave, Landau and Khalatnikov obtained the coupling amplitude for the 2 ↔ 2 processes, including the contribution of virtual off-resonant 3-phonon processes and introducing by hand corrections to hydrodynamics in the form of a cubic term in the excitation spectrum. However, the fact that they computed the damping rates only in the low and high wave-number limits makes this study incomplete and motivates us to revisit it. Studies of low-energy 2 ↔ 2 processes in 4 He were put on hold after 1970, when it was discovered that the dispersion relation of 4 He is in fact convex at low wave-number [START_REF] Maris | Phonon Dispersion and the Propagation of Sound in Liquid Helium-4 below 0.6K[END_REF][START_REF] Sherlock | The angular spreading of phonon beams in liquid helium 4: upward phonon dispersion[END_REF] at least at not too high pressure [START_REF] Sherlock | The angular spreading of phonon beams in liquid helium 4: upward phonon dispersion[END_REF][START_REF] Damian | Pressure dependence of phonon interactions in liquid 4 He[END_REF]. Studies on processes in the wave-number ranges where the dispersion relation becomes concave again, particularly challenging due to the existence of a small convex region, are still an active research topic [START_REF] Adamenko | Theory of scattering between two phonon beams in superfluid helium[END_REF]. To our knowledge, the study of the 2 ↔ 2 processes in cold Fermi gases is mainly at the unitary limit and in the collisional hydrodynamic regime where the relevant measure for dissipation is the shear viscosity of the gas [START_REF] Rupak | Shear viscosity of a superfluid Fermi gas in the unitarity limit[END_REF][START_REF] Enss | Viscosity and scale invariance in the unitary Fermi gas[END_REF][START_REF] Mannarelli | Phonon contribution to the shear viscosity of a superfluid Fermi gas in the unitarity limit[END_REF]. In this article, we restrict to temperatures low enough to be in the collisionless regime. There, we present a complete study of the interaction processes between phonons in cold Fermi gases, for any interaction strength in the BEC-BCS crossover, therefore in both the concave and convex cases.

We describe the unpolarized pair-condensed Fermi gas using a microscopic semiclassical model, based on a variational state including moving pairs. This model allows for a systematic expansion of the Hamiltonian in terms of two canonically conjugate fields β and β * that are weak if the density of excitations in the gas is low. The microscopic model and the principle of the method are exposed in section 2.

In section 3 we consider the expansion of the Hamiltonian to the quadratic order in the fields β and β * . By diagonalizing the quadratic problem, we find the spectrum of the excitations [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF][START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF][START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF][START_REF] Klimin | Pair Excitations and Parameters of State of Imbalanced Fermi Gases at Finite Temperatures[END_REF][START_REF] Randeria | Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the Unitary Fermi Gas[END_REF][START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]. Beyond the purely spectral results known in the literature, we derive here explicitly the quasiparticle modes and the annihilation/creation amplitudes, which we then quantize. This is rather elegantly done thanks to the Hamiltonian character of the equations in the semiclassical model.

Section 4 is devoted to the calculation of the coupling amplitudes for three-phonon and four-phonon interaction processes. We give the quantum hydrodynamics [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] predictions for the amplitudes of the 2 ↔ 1 and the 2 ↔ 2 processes, and compare them on the energy shell to a microscopic derivation, based for 2 ↔ 1 on the fermionic model of section 3, and for 2 ↔ 2 on a model of bosons with finite range interactions designed to have a concave dispersion relation at low wave numbers. Our microscopic test of the hydrodynamic predictions is particularly valuable for the 2 ↔ 2 processes, as (i) these processes involve nonresonant 2 ↔ 1 and 3 ↔ 0 processes, whose amplitudes are not correctly given by hydrodynamics, (ii) genuine hydrodynamics leads to a divergence in the 2 ↔ 2 amplitude that Landau and Khalatnikov had to regularize "by hand" by introducing a curvature in the phonon dispersion relation.

In section 5 we present a direct application of the previous results by calculating the damping rate of phonons in a Fermi gas in the whole BEC-BCS crossover in the collisionless regime. To this aim we use the effective on-shell coupling amplitudes for 2 ↔ 1 and 2 ↔ 2 processes calculated in section 4 within a Master equation approach. The three main results of section 5 are that (i) both for a convex and concave dispersion, by introducing properly normalized dimensionless quantities, we obtain a universal curve, giving the damping rate as a function of the wave vector, for any interaction strength in the BEC-BCS crossover, (ii) we give all the analytic asymptotic behaviors of the phonon damping rates, for cq k B T and cq k B T where c is the speed of sound and T is the temperature, correcting in particular the result of Landau and Khalatnikov [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] in the concave case, and (iii) for the unitary gas and at zero temperature only, we calculate the first correction to the hydrodynamic prediction for the phonon damping rate, using the effective field theory of Son and Wingate [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF][START_REF] Escobedo | Effective field theory and dispersion law of the phonons of a non-relativistic superfluid[END_REF], and improving the result of reference [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF].

The microscopic approach

Interaction Hamiltonian in the s-wave

We consider a gas of fermions in two internal states ↑ and ↓ on a cubic lattice of step l with periodic boundary conditions in the finite-size volume [0, L] 3 . Opposite spin fermions have an on-site interaction characterized by a coupling constant g 0

V(r, r ) = g 0 δ r,r l 3 (1)
The discretisation of space automatically provides a momentum cut-off, restricting the single-particle wave vectors to the first Brillouin zone FBZ = [-π/l, π/l[ 3 of the lattice, while allowing one to keep the simplicity of contact interactions. The Hamiltonian of the system in the grand canonical ensemble of chemical potential µ is given by:

Ĥ = l 3 r,σ=↑/↓ ψ † σ (r) - 2 2m ∆ r -µ ψσ (r) + g 0 l 3 r ψ † ↑ (r) ψ † ↓ (r) ψ↓ (r) ψ↑ (r) (2)
where the discrete Laplacian operator ∆ r has the eigenfunctions e ik•r with eigenvalues -k 2 and the field operator of the fermions has the discrete anticommutation relations { ψσ (r), ψ † σ (r )} = δ σσ δ rr /l 3 , with σ, σ =↑ or ↓. The desired regime of zero-range interactions is obtained by taking the continuous space limit l → 0 for a scattering length a fixed to the experimental value. In practice, one eliminates the bare coupling constant g 0 using the following relation, derived from the open space (L = ∞) scattering theory applied to the potential (1) [START_REF] Castin | Simple theoretical tools for low dimension Bose gases[END_REF][START_REF] Castin | Basic Theory Tools for Degenerate Fermi Gases[END_REF]:

1 g 0 = m 4π 2 a - FBZ d 3 k (2π) 3 m 2 k 2 (3)

Ground state in the BCS approximation

The idea behind the theory of Bardeen-Cooper-Schrieffer (BCS) [START_REF] Bardeen | Theory of Superconductivity[END_REF] is to search for an approximation of the ground state of the Hamiltonian (2) within the family of states

|ψ BCS = k∈D U k -V k â † k↑ â † -k↓ |0 (4)
where

D = 2π L Z 3 ∩ [-π/l, π/l[ 3 is
the ensemble of wavevectors of the first Brillouin zone satisfying the periodic boundary conditions, the operator âkσ is a Fourier component of the field operator ψσ (r) and annihilates a fermion of wavevector k and spin σ, V k is the probability amplitude of finding a kσ fermion in |ψ BCS and

U k = 1 -|V k | 2 .
The locus of the minimizers of the classical energy functional ψ BCS | Ĥ|ψ BCS is the circle

U k = U 0 k V k = V 0 k e iφ (5)
where φ ∈ [0, 2π[ and the values of reference U 0 k and V 0 k are chosen real.

To study the low-energy behavior of the system, we select on this circle the state of phase φ = 0

|ψ BCS 0 = k U 0 k -V 0 k â † k↑ â † -k↓ |0 (6)
as the origin of our expansion. In this symmetry-breaking state, we define the order parameter of the BCS theory

∆ ≡ g 0 ψ BCS 0 | ψ↓ ψ↑ |ψ BCS 0 = - g 0 L 3 k∈D U 0 k V 0 k , (7)
This quantity replaces g 0 or a as the most natural parameter characterizing the interaction strength in the BCS theory.

It can be used to express the coefficients U 0 k , V 0 k in the form

V 0 k = 1 2 1 - ξ k k and U 0 k = 1 2 1 + ξ k k ( 8
)
with the energies

ξ k = 2 k 2 2m -µ + g 0 ρ 2 (9) k = ∆ 2 + ξ 2 k ( 10
)
With the useful relation U 0 k V 0 k = ∆/2 k , the gap equation ( 7) takes its more usual form:

1 g 0 = - 1 L 3 k∈D 1 2 k (11)
In the state (6), the average total density ρ, or average total number of particles N per unit volume, is given by

ρ ≡ N L 3 ≡ k 3 F 3π 2 = 2 L 3 k∈D (V 0 k ) 2 (12)
where k F is the Fermi wave number of the ideal gas of density ρ. Combined with equation [START_REF] Kurkjian | Brouillage thermique d'un gaz cohérent de fermions[END_REF] this leads to the BCS equation of state which relates ρ to µ and ∆.

Performing a Bogoliubov rotation on the particle creation and annihilation operators, one finally defines the creation and annihilation operators of fermionic excitations:

γk↑ = U 0 k âk↑ + V 0 k â † -k↓ (13) γ
-k↓ = -V 0 k â † k↑ + U 0 k â-k↓ (14)
The BCS ground state [START_REF] Smith | Direct measurements of phononphonon scattering in liquid He 4[END_REF] is the vacuum of these operators, that annihilate k ↑ and -k ↓ quasiparticles of energy k . For completeness we give the corresponding modal expansion of the fermionic fields:

ψ ↑ (r) ψ † ↓ (r) = 1 L 3/2 k∈D γk↑ U 0 k V 0 k e ik•r + γ † k↓ -V 0 k U 0 k e -ik•r (15)

Ansatz of moving pairs

To study the fluctuations of the system, we view the symmetry-breaking BCS ground state as the quasiparticle vacuum and we construct a coherent state of pairs of quasiparticles, or squeezed vacuum, from it: [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF] where the independent variational parameters z k+q,k are the coherent state complex amplitudes and N(t) is the normalisation factor. This Ansatz is inspired from section 9.9b of reference [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF]. Contrarily to the BCS ground state [START_REF] Smith | Direct measurements of phononphonon scattering in liquid He 4[END_REF] in which all the pairs are at rest it includes pairs of quasiparticles with nonzero center-of-mass wavevector q. This will allow us to study the gapless, phononlike, collective modes of the system and their interactions. We do not include in the Ansatz elementary fermionic excitations obtained by the direct action of γ † kσ on the BCS ground state |ψ 0 BCS since they are gapped, and we are interested here in excitation energies much smaller than the gap. Note that the quasiparticle annihilation and creation operators γkσ and γ † kσ are linear combinations of the particle annihilation and creation operators âkσ and â † kσ , see equations [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF][START_REF] O'hara | Observation of a Strongly Interacting Degenerate Fermi Gas of Atoms[END_REF], so that considering a coherent state of quasiparticle pairs ( 16) is equivalent to considering a coherent state of pairs of particles, that is in fine a general time-dependent BCS Ansatz. In the end, all these Ansatz are Gaussian in the fermionic field operators, but the form ( 16) leads to much simpler calculations in the weakly excited regime as it immediately identifies the z k+q,k as small parameters of the expansion.

|ψ(t) = N(t) exp         k,q z k+q,k (t)γ † k+q↑ γ † -k↓         |ψ 0 BCS
References [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] and [START_REF] Kurkjian | Cohérence, brouillage et dynamique de phase dans un condensat de paires de fermions[END_REF] explain how to apply the variational principle to the state [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF]. One cleverly introduces the variables

β k,k = -z(1 + z † z) -1/2 k,k (17)
with z the matrix (z k,k ) k,k ∈D , so that the equations of motion take a Hamiltonian form

i dβ k ,k dt = ∂E ∂β * k ,k (18)
-i dβ * k ,k dt = ∂E ∂β k ,k (19)
One evaluates the associated classical Hamiltonian

E ≡ ψ| Ĥ|ψ (20)
through Wick's theorem, that relates E to averages of operators that are bilinear in γkσ . Finally one obtains an elegant expression for these averages by performing a Schmidt decomposition of the matrix z and expressing the result in terms of the field β, in matrix notation:

γ-k↓ γk ↑ = -β(1 -β † β) 1/2 k ,k (21) γ
† k↑ γ † -k ↓ = -β(1 -β † β) 1/2 * k,k (22) γ
† k ↑ γk↑ = ββ † k,k (23) γ
† -k ↓ γ-k↓ = β † β k ,k . (24)
3. Quadratic order of the classical Hamiltonian: normal variables of the bosonic branch and phonon operators

Linearized equations of motion

To linearize the equations of motion ( 18) and ( 19), we assume that the state ( 16) differs only slightly from the BCS ground state, so that

∀k, k , |z k,k | 1 and |β k,k | 1, (25)
We expand the energy functional in powers of the field β

E = E 0 + E 2 + O(β 3 ) ( 26
)
where

E 0 = ψ BCS 0 | Ĥ|ψ BCS 0 is a constant and E 2 is bilinear in β.
There are no linear terms in this expansion since the BCS ground state (corresponding to β = 0) is a minimizer of E. As in the Random Phase Approximation (RPA) [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF], the linearized equations of motion, which follow from [START_REF] Zwierlein | Condensation of Pairs of Fermionic Atoms near a Feshbach Resonance[END_REF] and [START_REF] Zwierlein | Vortices and superfluidity in a strongly interacting Fermi gas[END_REF] after replacement of E by E 2 , are decoupled according to the total wavevector q. For this reason, we rewrite the coordinates of the field β k 1 ,k 2 in the Anderson fashion, with the relative wavevector (k 1 + k 2 )/2 in the index and the center-of-mass one k 1k 2 in the exponent, keeping in mind that β k 1 ,k 2 actually involves the physical wave vectors k 1 and -k 2 , see equation [START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF], so that:

β q k ≡ β k+q/2,k-q/2 (27)
In terms of the vectors β q = (β q k ) k∈D and βq = ((β -q k ) * ) k∈D , the equations of motion take the matrix form i d dt

β q βq = L q β q βq , (28)
with an evolution operator L q that is symplectic

σ z L q σ z = (L q ) † ( 29
)
and particle-hole symmetric

σ x L q σ x = -(L q ) * (30)
We have introduced σ x = 0 1 1 0 and σ z = 1 0 0 -1 in block notations, and we termed (30) a particle-hole symmetry because σ x exchanges the field β with the complex-conjugate field β, which is reminiscent of the fact that, in the quantum world, the hole field operator is the hermitian conjugate of the particle field operator.

We give an explicit expression of the evolution operator L q in the sum-and-difference basis

y q k = β q k -βq k (31)
s q k = β q k + βq k (32)
We have i dy

q k dt = kq s q k + g 0 L 3 k ∈D W - kq W - k q + w + kq w + k q s q k (33) i ds q k dt = kq y q k + g 0 L 3 k ∈D W + kq W + k q -w - kq w - k q y q k ( 34
)
The coefficients W ± k,q and w ± k,q are cleverly chosen combinations of the coefficients

U 0 k and V 0 k W ± kq = U 0 k+q/2 U 0 k-q/2 ± V 0 k+q/2 V 0 k-q/2 (35)
w ± kq = U 0 k+q/2 V 0 k-q/2 ± U 0 k-q/2 V 0 k+q/2 (36)
and the energies kq are those of the continuum of two fermionic quasiparticles,

kq = k+q/2 + k-q/2 (37)
The equations of motion [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] contain two terms: first, an individual part that couples the amplitudes of same relative and center-of-mass wavevectors k and q and contains the trivial evolution of the operators γkσ under the BCS Hamiltonian ĤBCS = E 0 + kσ k γ † kσ γkσ , and second, a collective part that couples the normal amplitudes y q k and s q k to collective amplitudes of same total wavevector q. Our semi-classical equations of motion coincide with the quantum average in state [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF] of the RPA equations [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF]. This can be seen by using the definitions [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF] and [START_REF] O'hara | Observation of a Strongly Interacting Degenerate Fermi Gas of Atoms[END_REF] to express equations (78a-d) of reference [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF] in terms of the quasiparticle operators and by neglecting the expectation values of the γ † kσ γk σ operators, which, by virtue of equations ( 23) and [START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF], is justified at the linear order of the variational theory.

Branch of collective excitations

We now look for the eigenmodes of the linear system [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] with a positive energy ω q below the continuum k → k+q/2 + k-q/2 of two fermionic excitations, which is the spectrum we would obtain by restricting ourselves to the individual part of the system [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF]

0 < ω q < inf k ( k+q/2 + k-q/2 ) ( 38
)
The eigenvalue problem associated to [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] is solved in full generality in reference [START_REF] Kurkjian | Cohérence, brouillage et dynamique de phase dans un condensat de paires de fermions[END_REF]. We give the main steps of the solution in the simplifying continuous space limit l → 0 in Appendix A. The implicit equation on the mode eigenfrequency reads

I ++ (ω q , q)I --(ω q , q) = 2 ω 2 q I +-(ω q , q) 2 ( 39
)
where we introduced the integrals

I ++ (ω, q) = R 3 d 3 k        kq (W + kq ) 2 ( ω) 2 -( kq ) 2 + 1 2 k        (40)
I --(ω, q) = R 3 d 3 k        kq (W - kq ) 2 ( ω) 2 -( kq ) 2 + 1 2 k        (41)
I +-(ω, q) = R 3 d 3 k W + kq W - kq ( ω) 2 -( kq ) 2 (42)
We note that the same equation ( 39) can be obtained by the RPA; this was done by Anderson [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF] in the weak coupling limit and by reference [START_REF] Kurkjian | Cohérence, brouillage et dynamique de phase dans un condensat de paires de fermions[END_REF] in the general case. Although the linear system of the RPA contains extra terms proportional to the operators γ † k+q↑ γk↑ and γ † -k↓ γ-k-q↓ , which are absent in [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] because their expectation values in the Ansatz [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF] are quadratic in the field β, these extra terms can be treated as source terms and do not affect the spectrum. Equation [START_REF] Pitaevskii | Landau damping in dilute Bose gases[END_REF] can also be obtained by a Gaussian approximation of the action in a path integral framework [START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF][START_REF] Randeria | Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the Unitary Fermi Gas[END_REF] and a Greens functions approach associated with a diagrammatic approximation [START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF]. The conditions on q for the existence of a solution ω q are discussed in reference [START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF], while the concavity of the spectrum is studied in reference [START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]. Beyond those previous works, we construct here in the next section the quantum operators associated to the collective modes.

Construction of the normal variables of the collective branch

General case

Using the symmetries ( 29) and (30) of the evolution operator, we obtain the normal amplitudes exactly as in the bosonic case with the Bogoliubov theory [START_REF] Castin | Low-temperature Bose-Einstein Condensates in Time-dependent traps[END_REF]. The first step is to derive the eigenvector e + of energy

ω q > 0 e + (q) =        M q k N q k        k∈D (43)
To obtain the analytic expression of the coefficients M and N in the continuous limit l → 0, we solve the system [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] setting d/dt → -iω q and we use equations (A.6) and (A.7) to eliminate the collective amplitudes. We have

M q k -N q k = 2∆ kq W + kq -W - kq I ++ (ω q ,q) I +-(ω q ,q) N q ( 2 kq -2 ω 2 q ) (44)
M q k + N q k = 2∆ 2 ω 2 q W + kq -kq W - kq I ++ (ω q ,q) I +-(ω q ,q) ω q N q ( 2 kq -2 ω 2 q ) ( 45
)
where the dimensionless normalization constant N q will be determined by equation [START_REF] Kurkjian | Phase operators and blurring time of a pair-condensed Fermi gas[END_REF]. Important properties of M and N are their invariance by internal (k → -k) and external (q → -q) parity, where internal and external refers to the structure of the pairs:

M q k = M q -k = M -q k ( 46
)
N q k = N q -k = N -q k (47)
This is a consequence of the parity invariance of the problem and of the s-wave pairing. Due to the particle-hole symmetry [START_REF] Rupak | Density functional theory for non-relativistic fermions in the unitarity limit[END_REF], one can associate to the eigenvector e + another eigenvector e -of energy -ω q by a multiplication by σ x :

e -(q) = σ x        k∈D [(M q k ) 2 -(N q k ) 2 ] = 1 (51)
To obtain the amplitudes b q of the collective modes, we project the field β:

b q = d + (q) • β q βq = k∈D [M q k β q k -N q k (β -q k ) * ] (52) b * -q = d -(q) • β q βq = k∈D [ -N q k β q k + M q k (β -q k ) * ] (53)
This is the first central result of this paper. The equality b * q = (b q ) * suggested by our notation is a consequence of the invariance by external parity [START_REF] Castin | Basic Theory Tools for Degenerate Fermi Gases[END_REF][START_REF] Bardeen | Theory of Superconductivity[END_REF].

Conversely, to express the classical field in terms of the phonon amplitudes, we expand it on the eigenvectors:

β q βq = b q e + (q) + b * -q e -(q) + . . . (54)
where in the ellipsis . . . we omitted the component of the field on the other eigenmodes of total wavevector q. One can show [START_REF] Kurkjian | Cohérence, brouillage et dynamique de phase dans un condensat de paires de fermions[END_REF] that in the continuous limit those omitted modes are nothing else than the continuum k → k+q/2 + k-q/2 of fermionic-quasiparticle biexcitations with center-of-mass wavevector q. Projecting the vectorial equation ( 54), we obtain:

β k+q/2,k-q/2 = M q k b q + N q k b * -q + . . . (55)
β * k-q/2,k+q/2 = N q k b q + M q k b * -q + . . . (56)
After quantization (see Sec.3.4 below), M and N appear in these expressions as the coefficients of a new Bogoliubov rotation rearranging the fermionic bilinear operators appearing in equations [START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF][START_REF] Sidorenkov | Second sound and the superfluid fraction in a Fermi gas with resonant interactions[END_REF] into bosonic quasiparticle operators. It comes on top of the rotation induced by the U 0 k and V 0 k coefficients rearranging the particle operators into fermionic quasiparticle operators. This new rotation acts on the pairs of fermions, hence the two indices of M and N.

Long-wave limit

In the long-wave limit q → 0, we recall the expansion of the energy ω q obtained by reference [START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]:

ω q = q→0 cq        1 + γ 8 q mc 2 + O q mc 4        . (57)
In this expression, c is the speed of sound, derived, as in any superfluid, from the equation of state via the hydrodynamic formula

mc 2 = ρ dµ dρ , (58)
where the derivative with respect to ρ is taken for a fixed scattering length a. The coefficient γ of the cubic order of the expansion ( 57) is a rational fraction [START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF] of the variables

x = ∆ µ and y = d∆ dµ = R 3 d 3 k ξ k 3 k R 3 d 3 k ∆ 3 k (59)
Note that y can be related to x by the use of the BCS equation of state [START_REF] Beliaev | Application of the Methods of Quantum Field Theory to a System of Bosons[END_REF] in the last equality of (59). The q → 0 expansion of the sum [START_REF] Escobedo | Effective field theory and dispersion law of the phonons of a non-relativistic superfluid[END_REF] and difference (45) of M q k and N q k then read:

N q (M q k -N q k ) = ∆ k + O(q 2 ) (60) N q (M q k + N q k ) = - cq k 2∆ dW - k0 dµ + O(q 3 ) (61)
N q /L 3 = cq 2 dρ dµ + O(q 3 ) (62)
Note that to lowest order these expressions coincide with the coefficients of the zero-energy mode e n and of the anomalous mode e a of the zero-momentum subspace evolution operator L 0 obtained in reference [START_REF] Kurkjian | Phase operators and blurring time of a pair-condensed Fermi gas[END_REF].

As an application of the modal expansion [START_REF] Galitskii | The Energy Spectrum of a Non-ideal Fermi Gas[END_REF][START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF] and of the above low-q limit (60,61) of the corresponding amplitudes, we derive a modal expansion for the density ρ(r, t) = σ ψ † σ (r) ψσ (r) [START_REF] Castin | Bose-Einstein Condensates in Atomic Gases: Simple Theoretical Results[END_REF] and for the phase θ(r, t) of the order parameter, such that

g 0 ψ↓ (r) ψ↑ (r) = |g 0 ψ↓ (r) ψ↑ (r) |e iθ(r,t) , (64)
where the expectation value is taken in the time-dependent variational Ansatz [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF]. Expanding the Ansatz to first order in the amplitudes z k+q,k (t) -β k+q,k (t), using Wick's theorem to calculate the expectation value in the BCS ground state and replacing the fermionic fields ψσ (r) by their modal expansion from equation [START_REF] Bourdel | Measurement of the Interaction Energy near a Feshbach Resonance in a 6 Li Fermi Gas[END_REF], we obtain to leading order the deviations from the BCS ground state values:

δρ(r, t) = 1 L 3 k,q w + kq [β k+ q 2 ,k-q 2 (t)e iq•r + c.c.] + O(β 2 ) (65) δθ(r, t) = -g 0 2i∆L 3 k,q W + kq [β k+ q 2 ,k-q 2 (t)e iq•r -c.c.] + O(β 2 ) ( 66
)
where the coefficients w + kq and W + kq are defined by equations [START_REF] Kinast | Damping of a Unitary Fermi Gas[END_REF][START_REF] Gaunt | Bose-Einstein Condensation of Atoms in a Uniform Potential[END_REF]. Next, we replace the field β by its expansion [START_REF] Galitskii | The Energy Spectrum of a Non-ideal Fermi Gas[END_REF][START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF] on the collective modes to obtain, using the symmetry property [START_REF] Bardeen | Theory of Superconductivity[END_REF],

δρ(r, t) = 1 L 3/2 q B q (b q + b * -q )e iq•r + O(β 2 ) (67) δθ(r, t) = -i L 3/2 q C q (b q -b * -q )e iq•r + O(β 2 ) ( 68
)
with the coefficients

B q = 1 L 3/2 k w + kq (M q k + N q k ) (69) C q = g 0 2∆L 3/2 k W + kq (N q k -M q k ) (70)
In the low-q limit, also using the density [START_REF] Beliaev | Application of the Methods of Quantum Field Theory to a System of Bosons[END_REF], the gap equation [START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF] and the thermodynamic expression (58) of the sound velocity, we obtain the simple scaling laws

B q = qρ 2mc 1/2 [1 + O(q 2 )] (71)
C q = 2mc qρ 1/2 [1 + O(q 2 )] (72)
This result will be compared to the quantum hydrodynamic theory in section 4.1.2.

Quantization of the normal variables

To quantize the amplitudes of the bosonic modes obtained in our semi-classical approach, we remember that, at the linear order of the small amplitude approximation [START_REF] Combescot | Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover[END_REF], the field β q k reduces toγ-k+q/2↓ γk+q/2↑ , that is to the expectation value of a quasiparticle pair operator, see equation [START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF]. We therefore perform the substitution

β q k = β k+q/2,k-q/2 → -γ -k+q/2↓ γk+q/2↑ (73)
which transforms the amplitude b q in a quantum operator: bq

= - k∈D M q k γ-k+q/2↓ γk+q/2↑ -N q k γ † k-q/2↑ γ † -k-q/2↓ (74)
This intuitive action can in fact be justified by several arguments. First, the quantum operator (74) can also be obtained in the framework of the RPA by diagonalizing the homogeneous system for the quasiparticle pair creation and annihilation operators γk↑ γk ↓ and γ † k ↓ γ † k↑ , and treating the operators γ † kσ γk σ , whose dynamics is trivial in the RPA and whose average value is sub-leading in the variational theory, as source terms. Second, one can apply the quantization procedure based on bosonic images as described in chapter 11 of reference [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] where the quantum version of the field β is taken to be a bosonic field operator B (not to be confused with the bosonic excitation operators bq ). The expression of the bosonic image of a two-body fermionic operator such as γ-k+q/2↓ γk+q/2↑ in terms of the field B, which would provide an exact version of the substitution (73), is not simple in the general case since it involves an infinite series expansion in powers of B. Fortunately, in the limit of a weakly excited gas, one can neglect the operator population of the bosonic images BB † , and thus justify the substitution (73), once we assimilate the twobody fermionic operators and their bosonic images. Last, we highlight the bosonic nature of the operator bq (74) when the gas is weakly excited. The commutator of bq and b †

q reads: bq , b † q -1 = k∈D N q k 2 γ † k-q/2↑ γk-q/2↑ + γ † -k-q/2↓ γ-k-q/2↓ -M q k 2 γ † k+q/2↑ γk+q/2↑ + γ † -k+q/2↓ γ-k+q/2↓ (75)
It differs from unity by fermionic quasiparticle population operators, that are exactly zero in the BCS ground state and of second order in the field z = O(β) in a quasiparticle coherent state such as [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF].

Beyond the quadratic order: Interactions between phonons and comparison to hydrodynamics

The central actor in this work is the phonon-phonon interaction, that eventually causes phonon damping. It is derived in this section, in the form of three-phonon or four-phonon coupling amplitudes in phonon-interaction Hamiltonians that are cubic or quartic in the phonon field. Although the results involve sometimes lengthy formulas, the calculations are conceptually very simple: starting on one side from the microscopic energy functional [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF] or the easier-to-handle bosonic model (E.1) and on the other side from the quantum hydrodynamic Hamiltonian (87), one performs an expansion up to order three or four in the fields, and one inserts in the cubic or quartic terms the expansion of the fields on the phonon normal modes involving as linear coefficients the phonon creation and annihilation operators. That is all for the three-phonon processes in subsection 4.1. In the four-phonon processes of subsection 4.2, there is the additional subtlety that the three-phonon processes treated to second order have to be added to the direct four-phonon ones to get the correct effective four-phonon coupling as already understood in reference [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF]. An important conclusion of this section will be that the three-phonon and four-phonon coupling amplitudes of the microscopic model and of quantum hydrodynamics coincide in the low wave number limit on the energy shell (that is for processes conserving the unperturbed phonon energy), although they strongly differ outside the energy shell.

Three-phonon processes

In this subsection, we study the processes involving three bosonic quasiparticles: both the 2 ↔ 1 Beliaev-Landau processes and the off-resonant 3 ↔ 0 processes. Our goal is to identify these processes in the expansion of the Hamiltonian in powers of the excitation field, and to extract the associated matrix elements. We shall use two independent theories and compare them. First, we will use our microscopic variational theory and expand the energy functional E (20) up to order three in β,

E = E 0 + E 2 + E 3 + O(β 4 ) (76)
and we will insert the expansion (54) of the field on the collective eigenmodes in the trilinear term E 3 . We shall perform this microscopic calculation in the case 2 ↔ 1 only. Second, we will use the quantum hydrodynamics of Landau and Khalatnikov [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF]. This mesoscopic theory treats the pairs of fermions at large spatial scale as bosonic fields. It has the advantage of relying on the exact equation of state. It will be applied to the 2 ↔ 1 and 3 ↔ 0 processes. The comparison of the two theories will allow us to discuss the validity of quantum hydrodynamics.

Microscopic approach

The idea behind the microscopic calculation of the three-phonon coupling amplitudes is simple: we express the classical Hamiltonian [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF] in terms of the amplitudes b q of the collective modes (see Sec.3.3), isolate the terms containing the creation and annihilation amplitudes b * q and b q corresponding to the processes we study, and extract their coefficient. For the three-phonon processes, we focus on the cubic terms gathered in E 3 .

General case. The cubic part of E can be written as

E 3 = g 0 L 3 k,k ,q∈D T q k,k β k +q/2,k -q/2 γ † k+q/2↑ γk-q/2↑ +β k -q/2,k +q/2 γ † -k-q/2↓ γ-k+q/2↓ + c.c. ( 77
)
where we introduced the tensor

T q k,k = W - kq w + k q -w + kq W - k q -w - kq W + k q -W + kq w - k q 2 . ( 78
)
We use the relations [START_REF] Anderson | Random-Phase Approximation in the Theory of Superconductivity[END_REF][START_REF] Marini | Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions[END_REF] and we insert the expansion (54) in the expression (77) to obtain the following result1 :

E 3 = mc 2 (ρL 3 ) 1/2 q 1 ,q 2 ,q 3 ∈D δ q 1 +q 2 ,q 3 A 2↔1 micro (q 1 , q 2 ; q 3 )b * q 1 b * q 2 b q 3 + c.c. + . . . (79)
The vectors q 1 , q 2 and q 3 are the three wavevectors involved in the Beliaev-Landau process, with q 3 the wavevector of the phonon which decays into phonons of wavevectors q 1 and q 2 or which originates from the merging of two such phonons. In the ellipsis . . ., we omit the three-body processes involving non-bosonic excitations and the terms proportional to b q 1 b q 2 b q 3 or b * q 1 b * q 2 b * q 3 describing the off-resonant 3 ↔ 0 processes which we do not study with the microscopic theory. We have factorized the quantity mc2 /(ρL 3 ) 1/2 and defined the dimensionless coupling amplitude A 2↔1 micro (q 1 , q 2 ; q 3 ) of the process b * q 1 b * q 2 b q 3 , such that it is finite and nonzero in the thermodynamic limit:

mc 2 (ρL 3 ) 1/2 A 2↔1 micro (q 1 , q 2 ; q 3 ) = g 0 L 3 k,k ∈D T q 1 k,k M q 1 k M q 3 k-q 2 /2 M q 2 k-q 3 /2 + N q 2 k+q 3 /2 N q 3 k+q 2 /2 +T q 1 k,k N q 1 k M q 2 k+q 3 /2 M q 3 k+q 2 /2 + N q 3 k-q 2 /2 N q 2 k-q 3 /2 + T q 3 k,k N q 3 k N q 2 k+q 1 /2 M q 1 k-q 2 /2 + M q 3 k M q 2 k+q 1 /2 N q 1 k-q 2 /2 + q 1 ↔ q 2 (80)
where the notation +q 1 ↔ q 2 means that one must add to the terms already present in (80) those obtained by exchanging q 1 and q 2 while leaving q 3 unchanged. Although the principle of our microscopic approach to calculate the 2 ↔ 1 coupling amplitude is fairly simple, we obtain a relatively cumbersome expression (80). This is due to the resummation over the internal wave vectors k and k of the colliding pairs. A much more tractable expression is obtained in the long-wave limit in the next subsection.

Long-wave limit. At low temperature, the low-wavevector phonons dominate the kinetics of the gas. This motivates a specific study of A 2↔1 in the limit q 1 , q 2 , q 3 → 0, where a comparison to hydrodynamics is meaningful. We perform the expansion of the expression (80) in this limit and for a continuous space, l → 0, which allows us to integrate over the internal degrees of freedom k and k of the pairs and to obtain an expression of the coupling amplitude that depends only of the external wavevectors q 1 , q 2 and q 3 . The microscopic calculation is done in detail in reference [START_REF] Kurkjian | Cohérence, brouillage et dynamique de phase dans un condensat de paires de fermions[END_REF]; the underlying fermionic nature of the problem makes it rather tedious. We give here only the final result, expressed in terms of the energies ω q instead of the wavevectors, with the shorthand notation ω i ≡ ω q i , i = 1, 2, 3. The two sets of variables are connected by the dispersion relation ( 57) and momentum conservation relations. In the limit ω i → 0, we obtain:

A 2↔1 micro (q 1 , q 2 ; q 3 ) = mc 2 -1/2 1 2 3/2 √ ω 1 ω 2 ω 3 2J(x, y) (ω 1 + ω 2 -ω 3 ) + mc 2 2 A(x, y)(ω 1 + ω 2 -ω 3 )(ω 2 1 + ω 2 2 + ω 2 3 ) + B(x, y) ω 3 3 -ω 3 1 -ω 3 2 + C(x, y)ω 1 ω 2 ω 3 + O(ω 5 ) (81)
where the rational fractions A, B, C and J of the variables x and y defined in [START_REF] Kurkjian | Landau-Khalatnikov phonon damping in strongly interacting Fermi gases[END_REF] are given in Appendix B. The denominator (ω 1 ω 2 ω 3 ) 1/2 leads to a divergence of the coupling amplitude in the long-wave limit whenever the process is not on the energy shell, that is when it does not obey the energy conservation relation

ω 3 = ω 1 + ω 2 . ( 82
)
One can check that the Bose-Einstein Condensate (BEC) limit (a → 0 + ) of our result (81) coincides with the prediction of Bogoliubov theory [START_REF] Giorgini | Damping in dilute Bose gases: A mean-field approach[END_REF][START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF] for a weakly interacting gas of bosonic dimers 2 .

Resonant processes. To conclude this microscopic study, we evaluate the coupling amplitude on the energy shell, that is for processes satisfying the energy conservation (82). Such processes are allowed in the limit ω 3 → 0 for a positive γ parameter only, that is when the dispersion relation q → ω q is convex at low q. We collect the rational fractions B and C to form the thermodynamic quantity 83)), is plotted as a function of the interaction strength 1/k F a. The black solid line is obtained from the equation of state measured in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. In the BEC (k F a → 0 + ) and BCS (k F a → 0 -) limits the dashed black lines are the asymptotic behaviors given by the Lee-Huang-Yang corrections [START_REF] Lee | Many-Body Problem in Quantum Mechanics and Quantum Statistical Mechanics[END_REF][START_REF] Abrikosov | Concerning a Model for a Non-Ideal Fermi Gas[END_REF][START_REF] Galitskii | The Energy Spectrum of a Non-ideal Fermi Gas[END_REF] (that include the first correction to the mean field energy, scaling as (k F a) 2 on the BCS side and as ρa ρa 3 on the BEC side). For comparison the prediction of the BCS equation of state ( 12) is plotted as a dash-dotted blue line. Note that 1 + Λ F tends to 1 in the BEC limit, as for a weakly interacting Bose gas. It is equal to 8/9 in the unitary limit, and tends to the same value in the BCS limit, since in both cases µ ∝ ρ 2/3 .

2B(x, y) + 2 3 C(x, y) = 1 + ρ 3 d 2 µ dρ 2 dµ dρ -1 ≡ 1 + Λ F (83)
Its expression as a function of x and y is given in Appendix B and can also be obtained by differentiating twice the BCS equation of state [START_REF] Beliaev | Application of the Methods of Quantum Field Theory to a System of Bosons[END_REF] with respect to µ. This leads to an elegant expression of the on-shell coupling amplitude:

A 2↔1 OnS (q 1 , q 2 ; q 3 ) = 3 (1 + Λ F ) 3 q 1 q 2 q 3 32m 3 c 3 + O q 7/2 3 ( 84
)
The thermodynamic quantity 1 + Λ F is related to the so-called Grüneisen parameter u = ∂ ln c ∂ ln ρ by 1 + Λ F = 2 3 (1 + u). Its occurrence as a factor in equation ( 84) is the only difference between this fermionic formula and its equivalent for a gas of weakly interacting bosons obtained by the Bogoliubov theory (see equations (D8) and (D9) of reference [START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF]). It is plotted as a function of the interaction strength in the BEC-BCS crossover in figure 1. It tends to 1 in the BEC limit (x → 0 + , y ∼ -4/x) like in a gas of weakly interacting bosons where µ ∝ ρ, and it is equal to 8/9 at unitarity (x = y) and in the BCS limit (x → 0, y → 0) since in both cases µ ∝ ρ 2/3 .

Hydrodynamic approach

We now compare our microscopic result (81) to the irrotational quantum hydrodynamics of Landau and Khalatnikov. This theory performs a large-scale description of the gas in terms of two hermitian quantum fields ρ(r, t) and v(r, t), directly neglecting the exponentially weak density of fermionic quasiparticles at low T . The velocity field v(r, t) is assumed to be irrotational and is written as the (discrete) gradient of the phase field operator φ(r, t),

v(r, t) = m ∇ φ(r, t), (85)
canonically conjugated to the density field operator ρ(r, t):

[ ρ(r, t), φ(r , t)] = i δ r,r l 3 . ( 86
)
The dynamics of these fields is governed by the hydrodynamic Hamiltonian

Ĥhydro = l 3 r 2 2m (∇ φ) • ρ (∇ φ) + e 0,0 ( ρ) , (87)
where e 0,0 is the bare energy density, that shall be renormalized by the eigenmodes zero-point energy as described in [START_REF] Kurkjian | Brouillage thermique d'un gaz cohérent de fermions[END_REF] to give rise to the true energy density e 0 in the ground state, related to the zero-temperature chemical potential by µ = de 0 /dρ. The procedure to follow is standard and similar to the microscopic approach. We linearize the equations of motion for weak spatial fluctuations of the density and phase fields:

ρ(r, t) = ρ0 + δ ρ(r, t) (88) φ(r, t) = φ0 (t) + δ φ(r, t) (89)
where ρ0 will be replaced by the mean density ρ. We then expand the fields on the eigenmodes of the linearized dynamics:

δ ρ(r, t) = ρ 1/2 L 3/2 q∈D * q 2mc 1/2 ( bq + b † -q ) e iq•r (90) δ φ(r, t) = -i ρ 1/2 L 3/2 q∈D * mc 2 q 1/2 ( bq -b † -q ) e iq•r (91)
where the bosonic operators bq are the hydrodynamic counterpart of those in equation ( 74). The corresponding bosonic excitations have a purely linear spectrum ω hydro q

= cq where the sound velocity c at density ρ is still given by the hydrodynamic expression [START_REF] Landau | The theory of the viscosity of Helium II: II. Calculation of the viscosity coefficient[END_REF]. Next, we insert in the cubic part of the Hamiltonian (87)

Ĥ(3) hydro = l 3 r 2 2m ∇δ φ • δ ρ∇δ φ + 1 6 d 2 µ dρ 2 (δ ρ) 3 , (92)
the modal expansions ( 90) and ( 91) to obtain

Ĥ(3) hydro = mc 2 ρL 3 1/2
q 1 ,q 2 ,q 3 ∈D * δ q 1 +q 2 ,q 3 A 2↔1 hydro (q 1 , q 2 ; q 3 ) b †

q 1 b † q 2 bq 3 + h.c. + δ q 1 +q 2 +q 3 ,0 A 3↔0 hydro (q 1 , q 2 , q 3 ) b † q 1 b † q 2 b † q 3 + h.c. ( 93
)
where the coupling amplitudes of the 2 ↔ 1 and 3 ↔ 0 processes are given by A 2↔1 hydro (q 1 , q 2 ; q 3 ) = 3 q 1 q 2 q 3 32m 3 c 3 (3Λ F +u 12 +u 13 +u 23 ) (94)

A 3↔0 hydro (q 1 , q 2 , q 3 ) = 1 3 3 q 1 q 2 q 3 32m 3 c 3 (3Λ F +u 12 +u 13 +u 23 ) (95)
in terms of the parameter Λ F defined in equation ( 83) and of the cosine of the angle between the wave vectors q i and q j , u i j = q i • q j q i q j (96)

The 2 ↔ 1 amplitude clearly differs from the amplitude (81) obtained by the microscopic approach (in particular it does not diverge in the large wavelength limit). Physically, this may be surprising at first sight. The collective-mode component of the field in the microscopic theory is indeed expected to correspond to the field of the quantum hydrodynamic theory: for example, the modal expansions of the density and phase fields (67,68) in the microscopic theory seem to match, after quantization of the normal modes b q , the modal expansions (90,91) of quantum hydrodynamics (apart from an overall factor 2 in equation (68) due to the fact that the phase θ(r, t) of the pairing field is conjugated to the pair density, rather than to the particle density), so that the operators bq in equation ( 74) indeed seem to correspond to those in equations (90,91).

This discrepancy between the three-phonon coupling amplitudes can be attributed to two reasons: (i) the most obvious one is the difference of the Hamiltonians of the microscopic and hydrodynamic models, (ii) the somehow more hidden one is the fact that the modal expansions (67,68) in the microscopic model are valid only to first order in β, whereas the modal expansions (90,91) in quantum hydrodynamics are exact, with the consequence that the operators bq coincide in the two points of view only to first order in β. The first effect is due to the fact that quantum hydrodynamics is an effective field theory only valid at small wave numbers. The second effect is due to the fact that quantum hydrodynamics directly use the density field ρ(r, t) and the phase field φ(r, t) as the canonically conjugated variables of the theory (this corresponds to a modulus-phase point of view), whereas the microscopy theory uses the pair fields β k,k (t) and β * k,k (t) as conjugate variables, in which case the field density ρ(r, t) and the phase field θ(r, t) are only specific subfields that are nonlinear functions of β and β * , see equations [START_REF] Castin | Bose-Einstein Condensates in Atomic Gases: Simple Theoretical Results[END_REF]64). To qualitatively assess the relative importance of the two effects in the coupling amplitude discrepancies, the simplest way is to study the weakly interacting Bose gas case, which leads to less tedious calculations and reproduces the fermionic case in the BEC limit k F a → 0 + . This study is performed in Appendix C. We then conclude that the leading effect is effect (ii), that is the slight difference between the operators bq in the microscopic and hydrodynamic theories, as it explains the out-of-the-energy-shell divergence of the microscopic model coupling amplitude at vanishing wave numbers q i . On the contrary, effect (i) only gives a negligible correction O(q 7/2 i ) to the ≈ q 3/2 i hydrodynamic coupling amplitude. One can check however that the microscopic-theory and the hydrodynamic results for the 2 ↔ 1 coupling amplitude agree when the energy is conserved, see equation (84) and equation (94) written on the energy shell,

A 2↔1 hydroOnS (q 1 , q 2 ; q 3 ) = 3(1 + Λ F ) 3 q 1 q 2 q 3 32m 3 c 3 (97)
In the hydrodynamic theory on the energy shell the three wave vectors q 1 , q 2 and q 3 must indeed be colinear in the same direction in order to satisfy the equality |q 1 + q 2 | = q 1 + q 2 imposed by energy conservation for a linear spectrum (this is the well-known equality case in the triangular inequality). This is a satisfactory and important conclusion as, we shall see it, the Beliaev-Landau phonon damping rates Γ q of Sec.5.2 only involve the on-shell coupling amplitudes. In this respect, the microscopic and hydrodynamic theories lead as expected to the same low-energy physics.

Four-phonon processes

We now consider the 2 ↔ 2 four-phonon process. When the excitation branch q → ω q is concave at low wave numbers, this is the resonant process involving the minimal number of phonons because the 1 ↔ 2 and 1 ↔ 3 processes are now forbidden by energy conservation. Also, this process is more intriguing on a theoretical point of view since it involves virtual nonresonant 1 ↔ 2 or 3 ↔ 0 intermediate processes. In this case, the equivalence of hydrodynamics with the microscopic approach is not obvious, as hydrodynamics does not correctly describe the 1 ↔ 2 processes off-shell. In this section, we give the hydrodynamic prediction for the 2 ↔ 2 effective coupling amplitude that includes the virtual processes, then we validate the result with a microscopic model. Since the fermionic microscopic model would be quite cumbersome, we use a boson model with finite range interactions, such that the excitation branch is concave at low q.

Transition amplitude

We need to calculate the transition amplitude between an initial state of energy E i , which is an arbitrary Fock state of bosonic quasiparticles,

|i = |(n q ) q∈D (98)
and a final state of energy E f where two phonons of wave vectors q 1 and q 2 were annihilated and replaced by phonons of wavevectors q 3 and q 4 :

| f = b † q 3 b † q 4 bq 1 bq 2 n q 1 n q 2 (1 + n q 3 )(1 + n q 4 ) |i (99)
Whatever the specific model, the Hamiltonian can be expanded as

Ĥ = E 0 + Ĥ2 + Ĥ3 + Ĥ4 + . . . ( 100
)
where E 0 is a constant, Ĥ2 is the free quasiparticle Hamiltonian and Ĥ3 , Ĥ4 are the third order and fourth order terms. Ĥ3 cannot directly couple |i to | f , so we calculate the coupling to second order in perturbation theory, which amounts to treating Ĥ4 to first order and Ĥ3 to second order to construct an effective Hamiltonian [START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF]:

f | Ĥ2↔2,eff |i = f | Ĥ4 |i + λ f | Ĥ3 |λ λ| Ĥ3 |i E i -E λ ≡ A i→ f (101)
There exist 6 intermediate states |λ , labeled from I to VI, that can be accessed at zero temperature that is when all the modes q q 1 , q 2 are initially empty. These states correspond to the creation and reabsorption of a virtual phonon by a three-phonon nonresonant process. They are represented by the diagrams on the left part of figure Ĥ2↔2,eff = mc 2 ρL 3 q 1 ,q 2 ,q 3 ,q 4 ∈D δ q 1 +q 2 ,q 3 +q 4 A 2↔2,eff (q 1 , q 2 ; q 3 , q 4 ) b †

q 3 b † q 4 bq 1 bq 2 (102)
By construction, the matrix element of Ĥ2↔2,eff between |i and | f is the transition amplitude A i→ f . We thus have the relation

A i→ f = n q 1 n q 2 (1 + n q 3 )(1 + n q 4 ) 4mc 2 ρL 3 A 2↔2,eff (103)
where the factor 4 is a counting factor.

From now one we restrict to on-shell processes, such that

ω 1 + ω 2 = ω 3 + ω 4 , (104)
with the shorthand notation ω i = ω q i , i = 1, 2, 3, 4. In this case, a simplification occurs between each of the λ=I,II,III,IV,V,VI diagrams and its λ =I',II',III',IV',V',VI' counterpart on the right column of figure 2, which formally reduces the problem to zero temperature. Each diagram and its counterpart have indeed opposite energy denominators, with numerators that differ only through the factor involving the occupation number of the intermediate phonon q, (1 + n q ) in the left column (where the intermediate phonon is first created) and n q in the right column (where the intermediate phonon is first annihilated). Taking I and I' as an example, one has E i -

E I = ω 1 + ω 2 -ω q 1 +q 2 and E i -E I = ω q 3 +q 4 -ω 3 -ω 4 = -(E i -E I ); in the matrix element f | Ĥ3 |I I| Ĥ3 |i one has the factor (1 + n q 1 +q 2 )
and in f | Ĥ3 |I I | Ĥ3 |i one has the factor n q 1 +q 2 . Collecting the diagrams by pairs, we obtain an effective coupling amplitude identical to the zero-temperature one, that is with n q = 0 for all q q 1 , q 2 . In terms of the direct on-shell 2 ↔ 2 coupling amplitude A 2↔2,dir OnS (that is related to Ĥ4 in the same way as A 2↔2,eff is related to Ĥ2↔2,eff ) and of the amplitudes A 2↔1 and A 3↔0 introduced in section 4.1, the effective on-shell amplitude can finally be written as follows:

A 2↔2,eff OnS (q 1 , q 2 , q 3 , q 4 ) = A 2↔2,dir OnS (q 1 , q 2 , q 3 , q 4 ) + A 2↔1 (q 1 , q 2 ; q 1 + q 2 )A 2↔1 (q 3 , q 4 ; q 1 + q 2 ) ω1 + ω2 -ω1+2 + 9A 3↔0 (q 3 , q 4 , -q 1q 2 )A 3↔0 (q 1 , q 2 , -q 1q 2 ) -( ω1 + ω2 + ω1+2 ) + A 2↔1 (q 3 , q 1q 3 ; q 1 )A 2↔1 (q 1q 3 , q 2 ; q 4 ) ω1 -ω3 -ω1-3 + A 2↔1 (q 4 , q 3q 1 ; q 2 )A 2↔1 (q 3q 1 , q 1 ; q 3 ) ω3 -ω1 -ω3-1 + A 2↔1 (q 4 , q 1q 4 ; q 1 )A 2↔1 (q 1q 4 , q 2 ; q 3 ) ω1 -ω4 -ω1-4 + A 2↔1 (q 3 , q 4q 1 ; q 2 )A 2↔1 (q 4q 1 , q 1 ; q 4 ) ω4 -ω1 -ω4-1

where the energies were rescaled by mc 2 , ω ≡ ω mc 2 with the notation ω i± j ≡ ω q i ±q j .

(106)

Note that the coupling due to diagram VI (whose contribution is written here immediately after the one of diagram I) has a counting factor equal to 9.

Effective amplitude in hydrodynamics

At low temperature, we need the large wavelength limit of the on-shell expression (105). If one wants to get it from quantum hydrodynamics, one must introduce in quantum hydrodynamics a correction grasping some element of microscopic physics. Indeed, the hydrodynamic excitation spectrum is purely linear, which causes the energy denominators in (105) to vanish, and the effective coupling amplitude to diverge, when the wave vectors are colinear. This is an artefact of hydrodynamics. In reality the spectrum has a nonzero γ curvature parameter, here γ < 0, so that the energy denominators do not vanish. For almost colinear wave vectors, it is then natural to regularize the coupling amplitude by replacing the hydrodynamic dispersion relation ω hydro q = cq with the expansion (57) as done by Landau and Khalatnikov [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF]. We shall bring a microscopic justification to this: for colinear wave vectors, our bosonic microscopic model (see subsection 4.2.4) gives a large wavelength equivalent of formula (105) that indeed agrees with the Landau-Khalatnikov modified hydrodynamics.

First we determine the direct 2 ↔ 2 process amplitude from the quartic terms of the Hamiltonian:

Ĥ(4) hydro = 1 24 d 3 µ dρ 3 l 3 r δ ρ4 (107)
As we did previously, we insert the expansions (90) and (91) into Ĥ(4) hydro and we keep only the 2 ↔ 2 terms b †

q 3 b † q 4 bq 1 bq 2 : Ĥ2↔2,dir hydro = mc 2 ρL 3
q 1 ,q 2 ,q 3 ,q 4 ∈D δ q 1 +q 2 ,q 3 +q 4 A 2↔2,dir hydro (q 1 , q 2 ; q 3 , q 4 ) b †

q 3 b † q 4 bq 1 bq 2 (108)
to obtain the rescaled direct 2 ↔ 2 coupling amplitude A 2↔2,dir hydro (q 1 , q 2 ; q 3 , q 4 ) = Σ F 16 4 q 1 q 2 q 3 q 4 m 4 c 4 (109

)
where we introduced

Σ F ≡ ρ 3 mc 2 d 3 µ dρ 3 (110)
Second we combine the amplitude (109) with our previous expressions for the 2 ↔ 1 (94) and 3 ↔ 0 (95) amplitudes to obtain A 2↔2,eff hydro corr,OnS (q 1 , q 2 ; q 3 , q 4 ) =

1 16 4 ω 1 ω 2 ω 3 ω 4 m 4 c 8 Σ F + (ω 1 + ω 2 ) 2 A 1234 + ω 2 1+2 B 1234 (ω 1 + ω 2 ) 2 -ω 2 1+2 + (ω 1 -ω 3 ) 2 A 1324 + ω 2 1-3 B 1324 (ω 1 -ω 3 ) 2 -ω 2 1-3 + (ω 1 -ω 4 ) 2 A 1423 + ω 2 1-4 B 1423 (ω 1 -ω 4 ) 2 -ω 2 1-4       (111)
where the index "corr" means that one goes beyond the hydrodynamic approximation for the dispersion relation in the denominators of expression (111) by using the cubic approximation (57), and we introduced the coefficients

A i jkl = (3Λ F + u i j )(1 + u kl ) + (3Λ F + u kl )(1 + u i j ) + (1 + u i j )(1 + u kl ) ( 112
)
B i jkl = (3Λ F + u i j )(3Λ F + u kl ) (113) • • q S q 1 q 2 q 3 q 4 I • • q 1 q 2
q S q S q 3 q 4 I'

• • q t1 q 1 q 3 q 2 q 4 II • • q t1 q 2 q 4 q 1 q 3 q t1 II' • • -q t1 q 2 q 4 q 1 q 3 III • • -q t1 q 1 q 3 q 2 q 4 -q t1 III' • • q t2 q 1 q 4 q 2 q 3 IV • • q t2 q 2 q 3 q 1 q 4 q t2 IV' • • -q t2 q 2 q 3 q 1 q 4 V • • -q t2 q 1 q 4 q 2 q 3 -q t2 V' • • -q S q 3 q 4 q 1 q 2 VI • • q 1 q 2 -q S q 4 q 3 -q S VI'
Figure 2: Second order diagrams for the 4-phonon process (q 1 , q 2 ) → (q 3 , q 4 ) with q 1 + q 2 = q 3 + q 4 , considered as two successive three-phonon processes. The incoming wave vectors q 1 and q 2 and the emergent ones q 3 and q 4 are plotted as a solid line with an arrow. On the left column the diagrams include a virtual intermediate phonon, plotted as a dashed line with an arrow. On the right column they include a real intermediate phonon, plotted as a double line with an arrow. On any given row, the two diagrams have the same intermediate phonon: in I and I' q S = q 1 + q 2 = q 3 + q 4 , in II and II' q t1 = q 1q 3 = q 4q 2 , in III and III' -q t1 , in IV and IV' q t2 = q 1q 4 = q 3q 2 , in V and V' -q t2 , in VI and VI' -q S .

Effective amplitude in a weakly interacting Bose gas with nonzero range interactions

To study the 2 ↔ 2 process with a microscopic approach and understand how the divergence in hydrodynamics must be regularised, we shall not use the fermionic variational theory of section 2, which would be very heavy to manipulate due to the internal degrees of freedom of the pairs. We rather use a bosonic model with a large enough interaction range b so that the Bogoliubov excitation branch is concave at low q. In Appendix E, we describe this model in detail and we calculate the 2 ↔ 2 effective coupling amplitude in the weakly interacting limit from Bogoliubov theory. On shell, in the limit of small wave vectors, we then recover the prediction (111) of corrected quantum hydrodynamics specialised to the bosonic equation of state (E.11), that is with Λ F = Σ F = 0 and (m, c, ρ) → (m B , c B , ρ B ), where m B is the boson mass, c B the sound velocity in the Bose gas and ρ B the density of the bosons. The curvature of the spectrum, which eliminates the divergence of the 2 ↔ 2 effective coupling amplitude, is present from the start in the microscopic model and does not need to be added by hand. This provides a microscopic justification to the prescription of Landau and Khalatnikov [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] described in the beginning of subsection 4.2.3. The result was not obvious from the start, and it is due to a subtle cancellation between the large wavelength divergences of the direct coupling term (E.20) and of the second-order-perturbation-theory virtual coupling term. The details of the microscopic calculations are thus very different from the hydrodynamic ones, where the direct coupling term is zero for the considered equation of state (E.11).

Application: phonon damping in the BEC-BCS crossover

A general master-equation expression of the damping rates

To calculate the damping rate, we view the phonon mode of wave vector q as a harmonic oscillator coupled to the reservoir formed by the other quasiparticle modes [START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF], assumed to be at thermal equilibrium at temperature T . We thus split the low-energy effective Hamiltonian as

Ĥ = ω q b † q bq +         q q ω q b † q bq         + ( R † bq + b † q R) + . . . (114)
The first and the second terms, originating from Ĥ2 , describe the free evolution of the mode q and of the reservoir, respectively. The third term, originating from the part of Ĥ3 or Ĥ2↔2,eff involving bq or b † q , gives the coupling between the reservoir and the mode q. The ellipsis . . . includes (i) higher order nonlinear processes, (ii) processes that do not involve the mode q, and (iii) terms of the form R b † q bq , where R is an operator of the reservoir, that shift the energy of the mode q. In the Born-Markov approximation, valid in the weak coupling and collisionless limit of low temperature T and low wave number q at fixed nonzero cq/k B T ,3 one gets for the equation of motion of the mean number of excitations nq in the mode q [START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF]:

d dt nq = -Γ q ( nq -nq ) ( 115
)
where nq is the thermal equilibrium population of the mode:

nq = 1 exp ω q k B T -1 (116)
The damping rate Γ q is given by

Γ q = +∞ -∞ dt 2 e -iω q t Tr R [ R, R † (t)]σ eq R (t) (117)
where Tr R is the trace over the states of the reservoir, σ eq R is the thermal equilibrium density operator of the reservoir, [ Â, B] is the commutator of operators  and B, and the time evolution of the operator R † (t) of the reservoir is calculated in the interaction picture with the Hamiltonian q q ω q b † q bq .

Convex case: Beliaev-Landau damping

For a convex dispersion relation, the bilinear terms in R lead to resonant processes, so they give the leading contribution to Γ q at low temperature. We thus insert in (117) the following expression for R † :

R † = mc 2 (ρL 3 ) 1/2         q 2 ,q 3 A 2↔1 (q 2 , q 3 ; q) b † q 2 b † q 3 δ q,q 2 +q 3 + q 1 ,q 2
2A 2↔1 (q 2 , q; q 1 ) b † q 1 bq 2 δ q 2 +q,q

1         (118)
We split the contribution of the Beliaev process q ↔ (q , (qq )) (first contribution in ( 118)) from the one of the Landau process (q, q ) ↔ q + q (second contribution in (118)):

Γ Bel q = (mc 2 ) 2 2π 2 2 ρ R 3 d 3 q |A 2↔1
OnS (q , qq ; q)| 2 δ(ω q + ω q-q -ω q )(1 + nq-q + nq ) (119)

Γ Lan q = (mc 2 ) 2 π 2 2 ρ R 3 d 3 q |A 2↔1
OnS (q , q; q + q)| 2 δ(ω q+q -ω q -ω q )(n qnq +q ) (120) and in the low temperature limit we use the leading order expression (84) of the on-shell coupling amplitude 2 ↔ 1, or equivalently the hydrodynamic on-shell expression (97). We integrate over the wave vector q in spherical coordinates of polar axis q to obtain the low temperature equivalents 4

Γ Bel q ∼ T →0 9(1 + Λ F ) 2 32π mc 2 ρ mc 3 k B T mc 2 5
ΓBel ( q) (121)

Γ Lan q ∼ T →0 9(1 + Λ F ) 2 32π mc 2 ρ mc 3 k B T mc 2 5 ΓLan ( q) (122)
The wave numbers are here rescaled by the typical thermal wave number as follows:

q = cq k B T (123)
and the T → 0 limit is taken at fixed q. The functions ΓBel and ΓLan are universal functions of q, that can be expressed in terms of the Bose functions g α (z) = +∞ n=1 z n /n α , also called polylogarithms Li α (z), and of the Riemann ζ function, ζ(α) = g α (1), 5 ΓBel ( q) = q5 30 -4π 4 15 q + 48 ζ(5)g 5 (e -q) -24 qg 4 (e -q) + 4 q2 [ζ(3)g 3 (e -q)] (124) ΓLan ( q) = ΓBel ( q) -q5 30 + 8π 4 15 q (125)

This leads to the following limiting behaviors:

ΓBel = q→0 q4 6 + q6 360 + O( q8 ) (126) ΓBel = q→+∞ q5 30 +4ζ(3) q2 - 4π 4 15 q + 48ζ(5)+O( q2 e -q) (127)
The corresponding ones for ΓLan can be trivially deduced from equation (125). The two damping rates and their sum are plotted as functions of q in figure 3.

4. If one directly uses the quantum hydrodynamic linear dispersion relation ω q = cq in the argument of the Dirac distribution in equations (119,120), we face the polar integral 1 -1 duδ(u -1) (u is the cosine of the angle between q and q and the energy is conserved for u = 1 only), which seems to be equal to 1/2. If one correctly takes into account the strict convexity of the dispersion relation, through the curvature parameter γ > 0, we obtain the correct value 1 -1 duδ(uu 0 ) = 1 with the energy conserving root u 0 ∈] -1, 1[ arbitrarily close to 1 in the low temperature limit. Obviously this value does not depend on γ because the Dirac distribution always has a unit mass, which explains why the results (121,122) are γ-independent, but in a subtle way it utilises the sign of γ: if γ was negative, there would be no energy conserving root u 0 in the interval of integration, and the resulting integral would vanish, indicating the absence of Beliaev-Landau damping for a concave dispersion relation.

5. One just has to series expand 1/(e q -1) and 1/(e q + q -1) in powers of the variable e -q and to exchange summation and integration. We recall that ζ(4) = π 4 /90. 125) rates (dashed lines). The rates are rescaled as in equations ( 121) and (122).

Concave case: Landau-Khalatnikov damping

For a concave dispersion relation, the mode q is damped at sufficiently low temperature by the 2 ↔ 2 process only; in equation ( 117), we thus keep in R † the terms b † q 3 b † q 4 bq 2 , that originate from the quartic Hamiltonian (102). More explicitly, we replace R † in equation ( 117) by the following expression: R † = mc 2 ρL 3 q 2 ,q 3 ,q 4 2A 2↔2 eff (q, q 2 ; q 3 , q 4 ) b † q 3 b † q 4 bq 2 δ q+q 2 ,q 3 +q 4 (128)

This leads to

Γ 2↔2 q = (mc 2 ) 2 4π 5 2 ρ 2 R 6 d 3 q 2 d 3 q 3 |A 2↔2,eff
OnS (q, q 2 ; q 3 , q 4 )| 2 δ(ω q 3 + ω q 4 -ω q 2 -ω q ) nq 2 (1 + nq 3 )(1 + nq 4 ) -(1 + nq 2 )n q 3 nq 4

(129) where the vector q 4 is expressed in terms of the other wave vectors through momentum conservation:

q 4 = q + q 2 -q 3 (130)
In what follows, we explain how to obtain a low-temperature equivalent for the rate Γ 2↔2 q at fixed q = cq/k B T . We take as the polar z axis of the spherical coordinates the direction of q; the vectors q i then have coordinates (q i , θ i , φ i ). We consider a temperature T , controlled by the small parameter

≡ k B T mc 2 1, (131)
so low that the typical wave numbers are much smaller than mc/ and the bosonic branch is populated in its quasilinear part only. In this case, the coupling amplitude A 2↔2 OnS , that would be divergent at vanishing angles for a linear dispersion relation as already pointed out in section 4.2.3, is extremely peaked around θ 2 = θ 3 = 0 with a width of order in θ 2 and θ 3 , 6 and a height 1/ 2 times larger than the typical amplitude at nonzero angles, as we shall see.

6. The peak width is of order because the intermediate-virtual-phonon energy change when the θ i vary from 0 to is of the same order as the cubic correction to the hydrodynamic dispersion relation, knowing that the typical wave numbers are ≈ k B T/ c. Taking as an example the intermediate phonon q S = q + q 2 , one finds that c[|q + q 2 | -(q + q 2 )] ∼ -cqq 2 q+q 2 θ 2 2 ≈ k B T θ 2 2 ; this deviation is of the same order as the cubic term ≈ k B T |γ| 2 in the equation ( 57) when θ 2 2 ≈ |γ| 2 , hence the change of variable (132).

Using conservation of energy and momentum (130), one also finds that θ 4 = O( ) over the width of the peak. We thus rescale the wave numbers as in equation ( 123) and the polar angles as

θi = θ i |γ| 1/2 (132)
with γ < 0 the curvature parameter ( 57), then we perform a Taylor expansion of the coupling amplitude (111) for → 0 at fixed rescaled quantities:

A 2↔2 OnS (q, q 2 ; q 3 , q 4 ) = →0 3(1 + Λ F ) 4 2 ( q q2 q3 q4 ) 1/2 |γ| A 2↔2 red ( q, q2 , q3 , θ2 , θ3 ) + O( 2 ) ( 133
)
The Σ F -dependent term of the A 2↔2,dir direct amplitude does not contribute at this order. This property, combined with the clever rescaling (132) of the polar angles, allowed us in equation ( 133) to pull out the factors 1 + Λ F and γ that depend on the interaction strength. This leads to a universal A 2↔2 red reduced amplitude:

A 2↔2 red ( q, q2 , q3 , θ2 , θ3 ) = 1 q q2 θ2 2 ( q+ q2 ) 2 + 3 4 - 1 q q3 θ2 3 ( q-q3 ) 2 + 3 4 - 1 q( q + q2 -q3 ) θ2 4 ( q3 -q2 ) 2 + 3 4 ( 134
)
The first, second and third terms in (134) originate from the second, third and fourth terms in (111). In the last two terms, we carefully distinguished the cases q > q 3 and q < q 3 , q > q 4 and q < q 4 before taking the → 0 limit. Next, the implicit relation issued from energy conservation,

q4 = q + q2 -q3 - 2 |γ| 8 q3 + q3 2 -q3 3 -q3 4 + O( 4 ) (135)
is iterated once and combined with a spherical geometry calculus projecting relation (130) over q. This gives the following expression for θ4 :

θ2 4 = q2 θ2 2 -q3 θ2 3 -1 4 q3 + q3 2 -q3 3 -( q + q2 -q3 ) 3 q + q2 -q3 + O( 2 ) ( 136
)
Note that q3 < q + q2 according to (135). Also, due to the rotational invariance around q, the integrand of (129) depends on the azimuthal angles only through their difference φ ≡ φ 2 -φ 3 . The last step is to integrate the Dirac distribution ensuring energy conservation. This is conveniently done in a polar representation of the rescaled angles:

θ2 = R cos α θ3 = R sin α (137)
We also write the energy difference between the initial state and the final state as

ω q 3 + ω q 4 -ω q 2 -ω q = mc 2 3 |γ| 2 uR 2 + v + O( 5 ) (138) with u = q( q3 sin 2 α -q2 cos 2 α) + q2 q3 (1 -sin 2α cos φ) q + q2 -q3 (139) v = 1 4 q3 + q3 2 -q3 3 -( q + q2 -q3 ) 3 (140)
In the form (138), the Dirac distribution is readily integrated over R. We can finally express the 2 ↔ 2 damping rate in terms of a universal function Γ2↔2 of the rescaled wave number q:

Γ 2↔2 q mc 2 ∼ →0 81(1 + Λ F ) 4 256π 4 |γ| k B T mc 2 7 mc ρ 1/3 6 Γ2↔2 ( q) (141)
This is one of the central results of this paper. The function Γ2↔2 ( q) is given by a quadruple integral 7

Γ2↔2 ( q) = ∞ 0 d q2 q+ q2 0 d q3 q q3 2 q3 3 ( q + q2 -q3 ) |v| [1 + f ( q2 )] f ( q3 ) f ( q + q2 -q3 ) f ( q) × π/2 0 dα π 0 dφ sin α cos α Θ - v u v u A 2↔2 red q, q2 , q3 , v u 1/2 cos α, v u 1/2 sin α 2 ( 142
)
where we introduced the Heaviside function Θ(x ≥ 0) = 1, Θ(x < 0) = 0 and the reduced Bose function f (x) = 1/(e x -1). The occupation numbers in (129) were rewritten using the property (1 + nq i )/n q i = e ω q i /k B T and the conservation of energy. The function Γ2↔2 is plotted in figure 4. Its low-and high-q behaviors can be obtained analytically:

Γ2↔2 ( q) = q→0 16π 5 135 q3 + O( q4 ) (143) Γ2↔2 ( q) = q→∞ 16πζ(5) 3 q2 + O( q) ( 144
)
These limiting behaviors disagree with the results of Landau and Khalatnikov in reference [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] (see their equations (7.6) and (7.12) in the version [START_REF] Landau | The theory of the viscosity of Helium II: II. Calculation of the viscosity coefficient[END_REF]), even for the order in q of the leading terms. This discrepancy is due to the fact that these authors have neglected in the coupling amplitude at low q and at high q the contribution of the diagrams II to V of figure 2 (the diagram VI is nonresonant). This is not justified as already pointed out in reference [START_REF] Tucker | Four-phonon scattering in superfluid 4 He[END_REF], where it was observed numerically that the neglected diagrams are comparable in magnitude with the kept one. We find that at leading order, the neglected diagrams actually interfere destructively with the diagram I, which makes the exact results (143,144) sub-leading with respect to the Landau and Khalatnikov predictions by two orders in q.

Phonon damping beyond hydrodynamics

In this section we concentrate on the zero temperature unitary Fermi gas. Our goal is to calculate the phonon damping rates beyond the hydrodynamic result (127), which we rewrite in the form

Γ Bel q unitary gas hydro = 2 9π mc 2 mc ρ 1/3 3 q5 30 ( 145
)
where we introduced the notations q ≡ q mc and ωq ≡ ω q mc 2 (146)

To this aim, we first calculate the amplitude A 2↔1 OnS (q 1 , q 2 ; q 3 ) to first order in the curvature parameter γ of the excitation spectrum, that is we go in expression (84) one step beyond the hydrodynamic order in the low wave number limit. Then we calculate the Beliaev rate from (119) including first order corrections in γ also in the excitation spectrum. Last, we question the accuracy of the Fermi golden rule delta distribution in (119) and we investigate processes of higher order than Beliaev.

7. Let us give the explicit intervals of integration [α min , α max ] and [φ min , φ max ] imposed by the Heaviside function. We set χ = q 2 q 3 +q(q 3 sin 2 α-q 2 cos 2 α) q 2 q 3 sin 2α

,

A = q 2q 3 + q 2q 2 , C = 1 -q 2q 3 + q 2q 2 , A = A (1+A 2 ) 1/2 and C = C (1+A 2 ) 1/2 . (i)
If v > 0 and q < q 3 , 2α min = acos Aacos C and 2α max = acos A + acos C . (ii) If v > 0 and q > q 3 , α min = 0 and 2α max = acos A + acos C . (iii) If v < 0 and q < q 3 , α min = 0 and α max = π/2. (iv) If v < 0 and q > q 3 , 2α min = acos Cacos A and α max = π/2. In the cases (i) and (ii) (v > 0), φ min = 0, φ max = π if χ < -1 and φ max = acos χ otherwise. In the cases (iii) and (iv) (v < 0), φ max = π, φ min = 0 if χ > 1 and φ min = acos χ otherwise. The integral over φ can be calculated analytically. One faces Φ 0 dφ i b i a i +cos φ 2 that can be expressed in terms of the primitive F(Φ) = Φ 0 dφ 1 a+cos φ and of its derivative with respect to a. In order to make the limiting behaviors more apparent, the rate is divided by q2 and plotted as a function of argsh( q)= ln( q + 1 + q2 ). The oblique tangent at q = 0 corresponds to the low-q limiting behavior (see equation ( 143)). The horizontal straight line is the q → ∞ limit of Γ2↔2 / q2 (see equation ( 144)) and the asymptotic curve at large q is a fit of Γ2↔2 / q2 by an affine function of 1/ q. The black symbols are numerical calculations of the integral (129). In the numerics, we used a dispersion relation of the Bogoliubov form (149) restricted to an interval [0, q max ] over which its first order derivative is positive. The used dispersion relation coincides with [START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF] up to the third order in q. The discs are obtained by a linear extrapolation of the value of the integral for = 0, and the error bars halfwidth is given by the difference with a quadratic extrapolation.

If a ∈] -1, 1[, F(Φ) = 2 (1-a 2 ) 1/2 argth 1-a 1+a 1/2 tan Φ 2 . If |a| > 1, F(Φ) = 2 (a 2 -1) 1/2 atan a-1 (a 2 -1) 1/2 tan Φ 2 . If a = 1, F(Φ) = tan Φ 2 .

of a Beliaev process 2 ↔ 1 for the unitary gas

In this subsection we calculate the amplitude A 2↔1 OnS (q 1 , q 2 ; q 3 ) including the first correction in γ from the results of Son and Wingate using conformal invariance in an effective field theory [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF]. Son and Wingate constructed a Lagrangian that includes the first correction to hydrodynamics, expressed in terms of a phase field φ. For details on the formalism, we refer to their work [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF], and to reference [START_REF] Escobedo | Effective field theory and dispersion law of the phonons of a non-relativistic superfluid[END_REF] where it was used to calculate nonzero-temperature corrections to the sound velocity.

In short, we first quadratize the Lagrangian to determine the phonon excitation spectrum ω q and the corresponding eigenmode amplitudes. The spectrum is of the cubic form [START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF] with the curvature parameter γ = -64 45c 0 c 1 + 3 2 c 2 , where the dimensionless constant c 0 appears in the leading (hydrodynamic) order of the Lagrangian, and the dimensionless constants c 1 and c 2 appear in the next-to-leading order and parameterize the first corrections beyond hydrodynamics. The constant c 0 is linked to the already measured [START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF] Bertsch parameter ξ B relating the chemical potential µ to the Fermi wave number k F or the total density ρ:

µ = ξ B 2 k 2 F 2m and ρ = k 3 F 3π 2 (147)
The constants c 1 and c 2 , on the contrary, have not been measured yet; they have not been calculated either with good precision from a microscopic theory, only estimates are available, see [START_REF] Mannarelli | Phonon contribution to the shear viscosity of a superfluid Fermi gas in the unitarity limit[END_REF] and references therein. Second, we use the cubic terms in the Lagrangian to obtain the coupling between the phonon modes. This leads to the on-shell Beliaev coupling amplitude to first order beyond hydrodynamics:

A 2↔1 OnS (k, k ; q) = - √ 2 3 ( ωq ωk ωk ) 1/2 1 - 7γ 32 ( ω2 q + ω2 k + ω2 k ) + o( ω2 q ) (148)
Note that the same linear combination γ = -64 45c 0 c 1 + 3 2 c 2 as in the spectrum appears in this result. Details of the calculation are given in Appendix D.

Phonon damping rate for the T = 0 unitary gas

Up to cubic order in q, we can rewrite the excitation spectrum in a Bogoliubov form

ωq = q 1 + γ 4 q2 1/2 (149)
This allows us to recycle the result (A14) in Appendix A of reference [START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF] to perform the angular integration in (119),

Γ Bel q unitary gas = mc 2 π mc ρ 1/3 3 q 0 d ǩ|A 2↔1 OnS (k, k ; q)| 2 ǩ q ωq -ωk 1 + γ( ωq -ωk ) 2 1/2 (150)
where we acknowledge in the notation the fact that the coupling amplitude (148), where ωk = ωqωk , depends on the moduli of the wavevectors only. By performing the change of variable κ ≡ ǩ/ q, we express [ Γ Bel q unitary gas γ 2 / q] as an integral over κ between 0 and 1 of a function depending only on κ and on the small parameter γ q2 . By expanding this function to the sub-leading order, that is to third order in γ q2 , and performing the integral over κ, we obtain the provisional result

Γ Bel q unitary gas prov = 2 9π mc 2 mc ρ 1/3 3 q5 30 1 - 25 112 γ q2 + o( q2 ) (151)
The sub-leading term in our provisional result (151) is qualitatively different from the one in [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF], in particular it predicts a reduction of the damping rate with respect to the hydrodynamic prediction at low wavevectors rather than an increase. The disagreement is unexpected as the calculation of reference [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF] is performed with the same methodology as ours and in the same spirit. We think it is due to fact that the dependence of A 2↔1 on γ was finally neglected in [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF] while it gives a contribution of the same order as the dependence on γ in the spectrum. Before accepting the result (151), we should ask ourselves which kind of correction originates from the fact that the considered one-phonon state is unstable, of width Γ q /2, and as a consequence, that the free phonon energy is not exactly conserved in the one-to-two phonon decay process, contrarily to the constraint imposed by the Dirac delta function in equation ( 119). In order to estimate the order of magnitude of this effect, one can replace the Dirac delta function by a Lorentzian of half-width Γ/2,

δ(ω q + ω q-q -ω q ) → Γ q /2 (ω q + ω q-q -ω q ) 2 + Γ 2 q /4 (152)
where Γ q can be identified to its leading term in q 5 . One finds then that this effect introduces a correction to Γ q that is of order q 7 in the limit q → 0. 8 To get the exact contribution of this broadening effect to Γ q , we of course do not use the heuristic prescription (152) but we do a rigorous calculation in Appendix F. In short, using the resolvent of the Hamiltonian, we write at the Beliaev order a self consistent equation for the complex energy z q = ω qi Γ q /2 of the q phonon. If one replaces z q by its zeroth order approximation ω q + iη (η → 0 + ) in the implicit part of the equation, one recovers exactly (151). If one performs one self consistence iteration, that is if one replaces -2 Im z q / by its usual hydrodynamic approximation (145), one obtains the final result at the Beliaev order:

Γ Bel q unitary gas = 2 9π mc 2 mc ρ 1/3 3 q5 30        1 - 25 112 γ q2 + 4 √ 3ξ 3/2 B 243γ q2 + o( q2 )        (153)
where ξ B is the Bertsch parameter (147) of the unitary gas.

In order to conclude, and to make the result (153) rigorous, we have to verify that no other process of higher order than the Beliaev 1 → 2 process, gives a contribution of order q7 . A natural candidate is the cascade process represented in figure 5, which combines two 1 → 2 processes induced by the cubic Hamiltonian Ĥ3 , giving rise to a an effective coupling amplitude 1 → 3, which is of second order with a virtual phonon 9 . We estimate the corresponding change in the complex energy of the phonon q by treating the effective 1 → 3 coupling to second order in perturbation theory.

8. Using the note 4, we consider the polar integral J = 1 -1 du η /π η 2 +(u-u 0 ) 2 where η ≈ Γ q / q ≈ q4 and 1u 0 = 1cos θ 0 ≈ q2 as in (F.26). Then J -1 ≈ q2 gives a correction O( q7 ) to Γ Bel q coming from the finite energy width of the initial state. 9. The attentive reader has probably noticed that the quartic Hamiltonian Ĥ4 provides a direct 1 → 3 coupling to first order. The corresponding amplitude is in q 2 . This direct coupling contributes to Γ q at the order q 9 , that is here negligible even by taking into account the effect of a small denominator in q 3 . For this reason we do not mention it in the main text. q H 3 q 3 q-q 3 H 3 q 2 q 1 =q-q 2 -q 3

Figure 5: Effective coupling process, of second order in Ĥ3 involving the emission of a virtual phonon. As one can deduce from simple power counting (see the text), in the case of a convex dispersion relation and when the phonons q 1 , q 2 and q 3 are emitted within a small angle O( q/mc) with respect to q, this process contributes to the order ( q/mc) 7 to the T = 0 damping rate of a phonon of small wave vector q. Its contribution has then to be added to the result (153) as we verify in Appendix F by an explicit calculation.

At first sight, the result is O( q9 ) hence negligible. One has to integrate over two independent emitted phonon wave vectors, for example q 2 and q 3 , the third one q 1 being imposed by momentum conservation. As the wave numbers q i are of order q, this gives a factor q 6 . It comes then the product of four matrix elements of Ĥ3 , as the effective 1 → 3 coupling contains two matrix elements and it is treated to the second order, which adds a global factor (q 3/2 ) 4 = q 6 . Finally, there are three energy denominators, one from second order perturbation theory and the other two from the effective 1 → 3 coupling. As the phonon energies are of order cq, this provides a factor q 3 in the denominator. The whole thing is O( q9 ) as announced.

The previous reasoning however neglects the enhancement effect of the small denominators of order q 3 , that occurs when the wave vectors q i are emitted forward, with small angles θ i with respect to q, which already played a crucial role in section 5.3. In the limiting case where q and the q i are all aligned in the same direction, momentum conservation imposes q = q 1 + q 2 + q 3 (154) so that the energy difference ω q -( ω q 1 + ω q 2 + ω q 3 ) is not of order q but rather of order q 3 , taking into account the cubic term in the dispersion relation [START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF]. This conclusion extends to all the energy denominators as far as the emission angles are O( q). One can check indeed that q 1 = |qq 2q 3 | depends in relative value to second order in the emission angles θ 2 and θ 3 with coefficients of order q0 , in the same way as the true dispersion relation deviates in relative value from that of hydrodynamics to second order in q.

Let us then refine the naive estimation O( q9 ) of the previous paragraph, by considering the integration over q 2 and q 3 within cones of angular aperture O( q) around q. Each cone occupies a solid angle O( q2 ) so that we lose a factor q 4 in the integration over polar angles. On the other hand we gain a factor q -2 for each energy denominator, that is a global factor q -6 . We then predict a change in the complex energy of order q7 , that is the same order as the correction to hydrodynamics appearing in equation ( 153).

In Appendix F, we then explicitly calculate the contribution to the damping rate Γ q that comes from the effective coupling 1 → 3 treated to second order. We find that its expression, a rather tedious quintuple integral, leads to

Γ 1→3 q unitary gas = 2 9π mc 2 mc ρ 1/3 3 q5 30        - 2 √ 3ξ 3/2 B 567γ q2 + o( q2 )        (155)
In the same appendix, we verify that no other process, of arbitrarily high order in Ĥ3 , Ĥ4 , etc., contributes to Γ q to the order q7 , even accounting for the enhancement due to small denominators in q 3 . We then can sum the contributions (153) and ( 155) to obtain

Γ q unitary gas = 2 9π mc 2 mc ρ 1/3 3 q5 30        1 - 25 112 γ q2 + 22 √ 3ξ 3/2 B 1701γ q2 + o( q2 )        (156)
an exact expansion, plotted in figure 6, to be counted among the successes of this paper.

Conclusion

We presented a complete study of interaction processes between phonons in cold Fermi gases at low temperature, for any zero-range interaction strength between fermions, therefore in both the concave and convex cases for the phonon dispersion relation at low wave number q. We clarified the conditions of validity of a low-energy effective theory such as hydrodynamics by comparing it to a microscopic approach which takes into account the internal degrees of freedom of the pairs. Those effective theories correctly predict the phonon coupling amplitudes only on the energy shell.

One of the main contributions of this study is the microscopic derivation of the 2 ↔ 2 coupling amplitude. Indeed, its expression in second order perturbation theory includes nonresonant 2 ↔ 1 and 3 ↔ 0 processes. Since these processes can be resonant in the quantum hydrodynamic treatment, where the excitation spectrum is linear, Landau and Khalatnikov had to introduce "by hand" a curvature to the dispersion relation to avoid the divergence of the coupling amplitude. At first sight this procedure is risky since, as we just said, quantum hydrodynamics does not predict correctly the coupling amplitudes for nonresonant (off-shell) processes. To provide a microscopic test of this procedure would be particularly cumbersome within our fermionic microscopic approach, so we rather considered a model of weakly interacting bosons with finite range interactions designed to have a concave dispersion relation. We found that although the microscopic expressions of the coupling amplitudes differ from that of hydrodynamics, when one sums up all the nonresonant three-phonon processes to second order and the direct 2 ↔ 2 processes to first order, and one restricts on-shell with respect to the effective 2 ↔ 2 interaction, the Landau-Khalatnikov prescription and the microscopic result for the effective 2 ↔ 2 coupling amplitude agree.

As a second result of this paper, we gave universal formulas for the damping rates of both the 2 ↔ 1 and 2 ↔ 2 processes at low temperature, universal in the sense that the introduced reduced-rate functions q → Γ( q) do not depend on the atomic species, on the interaction strength, or even on the temperature. The most interesting result is the analytic derivation, as a function of q, of the phonon damping rate in the concave case, given to leading order in temperature by the 2 ↔ 2 processes, see (141) and figure 4, that is the subject of a Letter [START_REF] Kurkjian | Landau-Khalatnikov phonon damping in strongly interacting Fermi gases[END_REF]. In the limiting cases cq k B T and cq k B T , our result disagrees with the one of reference [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] and is sub-leading by two orders in cq/k B T . This is due to the failed assumption in [START_REF] Landau | Teoriya vyazkosti Geliya-II[END_REF] that some interaction diagrams are negligible while in reality they destructively interfere with the supposedly leading diagram.

Finally we also calculated, at zero temperature and in the unitary limit, the first correction ∝ q 7 to the hydrodynamic prediction ∝ q 5 for the single phonon decay rate. Our calculation allows to refine the prediction of [START_REF] Bighin | Beliaev damping of the Goldstone mode in atomic Fermi superfluids[END_REF] by (i) the actual inclusion of the beyond-hydrodynamics expression of the coupling amplitude [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF], (ii) the inclusion of a finite width Γ q in energy conservation in the Fermi golden rule, that is of a purely imaginary term of order q 5 in the energy denominator of perturbation theory and (iii) the inclusion of other processes, of higher order than the Beliaev process, in particular of the effective coupling 1 → 3 treated to second order.

All our predictions can be tested in state-of-art experiments on cold fermionic gases. In particular, a discussion of the observability of the Landau-Khalatnikov 2 → 2 damping, using cold atoms trapped in a flat-bottom potential, is done in [START_REF] Kurkjian | Landau-Khalatnikov phonon damping in strongly interacting Fermi gases[END_REF], with a proposed experimental protocol for the selective excitation of a phonon standing wave in the box in a thermodynamic-limit regime, that is at a wave number q much larger than the inverse of the system size but still in the linear part of the excitation spectrum.

Appendix A. Implicit equation on the collective excitation spectrum

We explain here how to obtain from the linear system [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] the implicit equation ( 39) on the collective mode angular frequency ω q at wave vector q, in the zero lattice spacing limit l → 0. We introduce the collective amplitudes

Y q = g 0 L 3 k∈D W + kq y q k S q = g 0 L 3 k∈D W - kq s q k (A.1)
y q = g 0 L 3 k∈D w - kq y q k s q = g 0 L 3 k∈D w + kq s q k (A.2)
and we solve the 2 × 2 linear system [START_REF] Katz | Beliaev Damping of Quasiparticles in a Bose-Einstein Condensate[END_REF][START_REF] Chevy | Transverse Breathing Mode of an Elongated Bose-Einstein Condensate[END_REF] to express the unknowns y q k and s q k in terms of the collective amplitudes. Next, we replace the result in (A.1,A.2) and come up with the homogeneous system

              Σ W + W + -1 Σ ω W + W - -Σ W + w - Σ ω W + w + Σ ω W -W + Σ W -W --1 -Σ ω W -w - Σ W -w + Σ w -W + Σ ω w -W - -Σ w -w --1 Σ ω w -w + Σ ω w + W + Σ w + W - -Σ ω w + w - Σ w + w + -1                             Y q S q y q s q               = 0 (A.3)
where we introduced the notations

Σ ab = g 0 L 3 k∈D kq a kq b kq ( ω q ) 2 -( kq ) 2 (A.4) Σ ω ab = g 0 L 3 k∈D ω q a kq b kq ( ω q ) 2 -( kq ) 2 (A.5)
with a and b that can be any of W + , W -, w + , w -. The system (A.3) simplifies in the continuous limit l → 0. Since g 0 → 0 and w ± kq = k→+∞ O(1/k 2 ), all the Σ tend to 0 in the third and fourth lines, and therefore we must have

y q = s q = 0 (A.6)
Next we divide the first two lines of (A.3) by g 0 . The gap equation [START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF] ensures that all the divided matrix elements have a finite nonzero limit. The system therefore reduces to its 2 × 2 upper left block, which, in the thermodynamic limit, we write as: I ++ (ω q , q) ω q I +-(ω q , q) ω q I +-(ω q , q) I --(ω q , q) Y q S q = 0 (A.7)

where the integrals I depend on the eigenfrequency ω q and on the wavevector q, see equations [START_REF] Rupak | Shear viscosity of a superfluid Fermi gas in the unitarity limit[END_REF][START_REF] Enss | Viscosity and scale invariance in the unitary Fermi gas[END_REF][START_REF] Mannarelli | Phonon contribution to the shear viscosity of a superfluid Fermi gas in the unitarity limit[END_REF]. The sought implicit equation on the eigenfrequency [START_REF] Pitaevskii | Landau damping in dilute Bose gases[END_REF] is the condition that the determinant of the 2 × 2 matrix in equation (A.7) is zero.

and hermitian conjugate, which brings an extra term A 2↔1 via H 2 to the Beliaev-Landau coupling amplitude. At the end, the Beliaev-Landau coupling amplitude in the ( ρ, θ) point of view reads:

A 2↔1 ρθ = A 2↔1 δψ + A 2↔1 via H 2 (C.8)
with

A 2↔1 via H 2 (q 1 , q 2 ; q 3 ) = ω3 d 3 8 (d 1 d 2 -s 1 s 2 ) + ω3 s 3 8 (s 1 d 2 + s 2 d 1 ) - ω2 d 2 8 (d 1 d 3 + s 1 s 3 ) - ω1 d 1 8 (d 2 d 3 + s 2 s 3 ) + d 3 8 (s 1 s 2 ω2 + s 2 s 1 ω1 ) - s 3 8 (d 2 s 1 ω1 + d 1 s 2 ω2 ) (C.9)
with ωi = ω q i /mc 2 , knowing that mc 2 = µ at this order. From the Bogoliubov expressions ω i = [ Ěi ( Ěi + 2)] 1/2 and s i = 1/d i = [ Ěi /( Ěi + 2)] 1/4 , one has ωi d i = ( Ěi + 2)s i and ˇ i s i = Ěi d i , where Ěi = E q i /mc 2 and E q = 2 q 2 /2m. Using as well the conservation of momentum q 3 = q 1 + q 2 and its squares q 2 3 = (q 1 + q 2 ) 2 , q 2 2 = (q 3q 1 ) 2 and q 2 1 = (q 3q 2 ) 2 , we obtain

A 2↔1 ρθ (q 1 , q 2 ; q 3 ) = 1 8 (s 3 d 1 d 2 q1 • q2 + s 2 d 1 d 3 q1 • q3 + s 1 d 2 d 3 q2 • q3 ) - 1 16 s 1 s 2 s 3 ( q2 1 + q2 2 + q2 3 ) (C.10)
with qi = q i /mc. We can now compare the Beliaev-Landau coupling amplitudes A 2↔1 δψ and A 2↔1 ρθ among themselves and to the quantum hydrodynamic coupling amplitude. With the expressions of s i and d i given above, it is found from equation (C.4) that A 2↔1 δψ diverges when all the q i tend to zero (out of the energy shell). On the contrary, even out of the energy shell, A 2↔1 ρθ becomes equivalent to the quantum hydrodynamic result (94) (here Λ F = 0) that tends to zero, but it does not coincide with it (it differs by relative corrections O(q 2 i )) due to the difference between the microscopic Hamiltonian (C.1) and the quantum hydrodynamic Hamiltonian.

As a final check, one can also directly insert the expansion (C.7) (turned into an expansion of δ ρ and θ) in the terms of the Hamiltonian (40) of reference [START_REF] Mora | Extension of Bogoliubov theory to quasicondensates[END_REF] that are cubic in δ ρ and θ. One then directly recovers the Beliaev-Landau coupling amplitude (C.10) among the bρθ q . Furthermore, if one neglects everywhere in this Hamiltonian (40) the kinetic energy stored in the density-variation δ ρ of the field, at the origin of the so-called quantum pressure term, one recovers exactly the quantum hydrodynamic Beliaev-Landau coupling (94) (with here Λ F = 0).

Appendix D. Beliaev coupling beyond hydrodynamics

We give here details of the calculation allowing one to obtain, for the unitary Fermi gas, the first beyondhydrodynamic correction to the Beliaev phonon coupling amplitude from the Son and Wingate Lagrangian [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF], see equation (148).

For convenience we move to a Hamiltonian formalism, introducing the field Π that is canonically conjugate to the phase field φ and that represents (up to a sign) density fluctuations, to obtain (here = m = 1)

H (0) 2 = 2µ -1/2 15c 0 Π 2 + 5 2 c 0 µ 3/2 1 2 (grad φ) 2 (D.1)
H (0) 3 = 4µ -2 3(15c 0 ) 2 Π 3 - 1 2 Π(grad φ) 2 (D.2) H (2) 2 = - 16c 1 µ -3/2 (15c 0 ) 2 (grad Π) 2 -c 2 µ 1/2 (∆φ) 2 (D.3) H (2) 3 = - 96c 1 µ -3 (15c 0 ) 3 Π(grad Π) 2 + 2c 2 µ -1 15c 0 Π(∆φ) 2 (D.4)
where the index 2 or 3 refers to the expansion order of the Hamiltonian in Π and φ while the exponent (0) or (2) refers to the expansion order in powers of spatial gradients, the zeroth order being the standard hydrodynamics. Note that Π and grad φ are of the same order. At the hydrodynamic level, the Hamiltonian depends on a single constant c 0 that must be determined from a microscopic theory. It is linked as follows to the Bertsch parameter (147):

c 0 = 2 5/2 15π 2 ξ 3/2 B (D.5)
The first correction to hydrodynamics involves two other dimensionless constants c 1 and c 2 . By following the procedure already used in this paper, we first use the quadratized Hamiltonian H (0) 2 + H (2) 2 to determine the excitation spectrum as in [START_REF] Son | General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas[END_REF]:

ω q = 2µ 3 1/2 q        1 -π 2 (2ξ B ) 1/2 c 1 + 3 2 c 2 q k F 2 + o(q 2 )        ≡ cq        1 + γ 8 q mc 2 + o(q 2 )        (D.6)
and the modal expansion of the fields in the quantization volume L 3 :

Π(r) = 1 L 3/2 q 0 Π q (b q + b * -q )e iq•r (D.7) φ(r) = 1 L 3/2 q 0 φ q (b q -b * -q )e iq•r (D.8)
with the amplitudes

Π q = 1 √ 2 A q B q 1/4
(D.9)

φ q = i √ 2 A q B q -1/4
(D.10)

A q B q = 75 8 (c 0 µq) 2 1 + 8µ -1 15c 0 c 1 - 3 2 c 2 q 2 + o(q 2 ) (D.11)
Indeed Π q φ q is an eigenvector of eigenvalue -iω q of the matrix 0 -A q B q 0 that appears in the Hamiltonian equations of motion. Note that different linear combinations of the constants c 1 and c 2 appear in the spectrum (D.6) and in the modal amplitudes (D.11). By inserting the modal decomposition in the cubic Hamiltonian H (0) 3 + H (2) 3 , and isolating the Beliaev terms 2 ↔ 1 as we did in (93), we obtain the on-shell Beliaev coupling amplitude to first order beyond hydrodynamics as given in (148).

The Bose gas is in the weakly interacting regime (ρ B a 3 B ) 1/2 1, where ρ B = N B /L 3 is the density of the bosons and the s-wave scattering length a B is given by 4π 2 a B /m B = Ṽ0 in the Born approximation. Following the Bogoliubov theory [START_REF] Bogoliubov | A new method in the theory of superconductivity I[END_REF] in its U(1) symmetry preserving version [START_REF] Gardiner | Particle-number-conserving Bogoliubov method which demonstrates the validity of the time-dependent Gross-Pitaevskii equation for a highly condensed Bose gas[END_REF][START_REF] Castin | Low-temperature Bose-Einstein Condensates in Time-dependent traps[END_REF][START_REF] Castin | Bose-Einstein Condensates in Atomic Gases: Simple Theoretical Results[END_REF][START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF], we split the bosonic field operator as ψ(r) = e i θ0 n1/2 0 φ 0 (r) + Λ(r) (E.4)

where θ0 is the condensate phase operator, n0 is the number of bosons in the condensate mode φ 0 (r) = 1/L 3/2 and the non-condensed field operator Λ(r), orthogonal to the condensate mode, conserves the number of particles. Within the subspace with fixed total number of bosons N B , we eliminate n0 through the relation

n0 = N B -l 3 r Λ † (r) Λ(r) (E.5)
To describe the 2 ↔ 2 processes, we expand the Hamiltonian in powers of Λ up to order 4:

ĤB = ĤB0 + ĤB2 + ĤB3 + ĤB4 + . . . (E.6) We obtain 10 ĤB0 = Ṽ0 N 2 B 2L 3 (E.7) ĤB2 = l 3 r Λ † (r) - 2 2m B ∆ r Λ(r) + ρ B l 6 r,r V(r -r ) Λ † (r) Λ(r ) + 1 2 Λ † (r) Λ † (r ) + Λ(r) Λ(r ) (E.8) ĤB3 = ρ 1/2 B 2 l 6 r,r V(r -r ) Λ † (r) + Λ † (r ) Λ(r ) Λ(r) + Λ † (r ) Λ † (r) Λ(r) + Λ(r ) (E.9) ĤB4 = l 6 2 r,r V(r -r ) Λ † (r) Λ † (r ) Λ(r ) Λ(r) - Ṽ0 2L 3        l 3 r Λ † (r) Λ(r)        2 - 1 L 3        l 3 r Λ † (r) Λ(r)                l 6 r,r V(r -r ) Λ † (r) Λ(r ) + 1 2 Λ † (r) Λ † (r ) + Λ(r) Λ(r )         (E.10)
0th order. From ĤB0 we get the equation of state to leading order:

µ B = ρ B Ṽ0 ≡ 2 2mξ 2 (E.11)
where we introduced the healing length ξ and the chemical potential µ B of the Bose gas. We then apply relation [START_REF] Landau | The theory of the viscosity of Helium II: II. Calculation of the viscosity coefficient[END_REF] to obtain the sound velocity m B c 2 B = µ B (E.12) 2nd order. ĤB2 can be diagonalised by a Bogoliubov transformation:

Λ(r) = 1 L 3/2 q 0 U B q bq e iq•r + V B q b † q e -iq•r (E.13)
where the bq are annihilation operators of bosonic quasiparticles and the amplitudes U B q and V B q are given by

U B q + V B q ≡ S q =          2 q 2 2m B 2 q 2 2m B + 2ρ B Ṽ(q)          1/4 (E.14) U B q -V B q ≡ D q = 1 S q (E.15)
10. We are here in the large N B limit and we neglect 1 as compared to n0 and to N B . In the equations (E.8) and (E.10), there was a cancellation of the Hartree contribution with a chemical potential-type contribution originating from the expansion of the condensate interaction energy Ṽ0 n2 0 /2L 3 in powers of the number of non-condensed particles. that is quadratic in the phonon creation b † and annihilation b operators. It is nevertheless coupled to the two-phonon, three-phonon, etc. continua by the rest of the Hamiltonian V = Ĥ3 + Ĥ4 + . . ., that contains cubic, quartic, etc. terms, when written in the normal order for the b.

As a consequence, the discrete state will in general get diluted in the continua, giving rise to a complex pole z q in the analytic continuation of the resolvent Ĝ(z) = (z -Ĥ) -1 of the full Hamiltonian [START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF]. This pole can be written as

z q = ω q -i Γ q 2 (F.1)
where ω q is the angular eigenfrequency of the phonon and Γ q its damping rate 11 .

Appendix F.2. Perturbative calculation and power counting The q → 0 limit corresponds to the weak coupling limit. This is apparent in the quantum hydrodynamic theory as the modal expansion of the velocity and density fluctuations v and δ ρ involves coefficients that tend to zero as q 1/2 , so that the matrix elements of Ĥp between phonon Fock states behave as q p/2 : Ĥp ≈ q→0 q p/2 (F.2)

One can then attempt a perturbative calculation of Γ q starting from the exact expression obtained by the projector method [START_REF] Cohen-Tannoudji | Processus d'interaction entre photons et atomes[END_REF] q|

Ĝ(z)|q = 1 z -q| Ĥeff (z)|q (F.3)
then expanding in powers of V the matrix element of the effective Hamiltonian in the state |q :

q|H eff (z)|q = q| Ĥ2 |q + q| V Q Q z Q -Q Ĥ2 Q Q V|q + q| V Q Q z Q -Q Ĥ2 Q Q V Q Q z Q -Q Ĥ2 Q Q V|q + . . . (F.4)
where Q = 1 -|q q| projects orthogonally to |q . The result (153) corresponds to the second term on the right-hand side of the equation, the Beliaev term, in which the contribution of V reduces to that of Ĥ3 . In all the following terms of order 3 in V or higher, z can be approximated in the denominator by z = ω q + iη, η → 0 + , (F.5) the displacement of z by an imaginary part that is O(q 5 ) in these terms does not contribute to Γ q at the order q 7 . The same conclusion does not hold for the Beliaev term, as we shall see in subsection Appendix F.3. Let us consider in (F.4) the term T n of order n in V and let us try to give an upper bound to its order in q as done in the simple reasoning below equation (153).

Taking (F.2) into account, we must retain in the V factors as many as possible contributions of Ĥ3 and, if not, those of Ĥ4 . As Ĥ3 changes the parity of the phonon number, while Ĥ4 preserves it, we will retain only Ĥ3 if n is even, but it will be necessary to keep at least one factor Ĥ4 if n is odd. The minimal number of independent wave vectors of the virtual phonons is n/2 for n even, (n + 1)/2 for n odd, and the integration over each independent wave vector provides a factor q 3 . Finally, one has to count a factor 1/q for each of the n -1 energy denominators. All together, this leads to

T n = O(q 2n+1 ) if n even (F.6) T n = O(q 2n+3 ) if n odd (F.7)
11. We cannot calculate here the corrections to ω q due to the coupling with the continua, even at the order 2 in V, because we rely on an effective Hamiltonian. If we try, we would encounter principal part integrals with ultraviolet divergences that, once inserted in the unperturbed value in Ĥ2 would give, following the ideas of renormalisation, the true value, which would however remain uncalculated and unknown in the absence of a microscopic model [START_REF] Kurkjian | Brouillage thermique d'un gaz cohérent de fermions[END_REF]. Within this context, we do not find it useful to distinguish between the non perturbed phonon eigenenergy for Ĥ2 , which we should rigorously note (2) q , and the true energy q . By the way, in order to obtain the scaling laws with q of the perturbative terms of order n in V, we assume that the wave numbers of the virtual phonons are O(q), hence that the ultraviolet cutoff is set at a wave number Aq, A

1. This has the big advantage that the real correction to the angular eigenfrequency (2) q due to Ĥ3 to the leading order is O(q 5 ) and it affects neither the γ parameter in equation ( 57) nor the damping rate Γ q to the order q 7 . Figure F.7: Diagram of order 2 in V that contributes to the damping rate Γ q of the phonon q to the order q 7 . It represents the Beliaev process whose leading order is q 5 . On the figure, each vertex represents an action of Ĥ(±)

3 .

One would then conclude with no surprise, that the processes of higher order than Beliaev, that is of order n ≥ 3 in V, give a contribution to Γ q that is O(q 9 ) hence negligible. The simple reasoning above forgets the possibility that the energy denominators are of order q 3 . For this case to occur, the contribution to the energy denominator of the linear part of the spectrum must have a chance to cancel out so the considered processes must have the possibility to be resonant. In order to have all the energy denominators of order q 3 , all the independent wave vectors must be emitted forward, within a narrow cone of angular aperture O(q) with respect to the direction of q. The angular integration then brings a factor q 2 for each independent wave vector, that is q n for n even and q n+1 for n odd, but we gain a factor q -2(n-1) thanks to n -1 small denominators.

This leads to the refined upper bounds

T n = O(q n+3 ) if n even (F.8) T n = O(q n+6 ) if n odd (F.9)

Figure F.8: The two diagrams of order 4 in V that may contribute to the damping rate Γ q of the phonon q to the order q 7 thanks to the "small denominator" effect. The one on the left (I) is a Beliaev process with a single loop correction (itself of the Beliaev nature) to the virtual phonons angular eigenfrequency. The one on the right (II) is a Beliaev process enriched by an interaction between the virtual phonons. Each vertex corresponds to an action of Ĥ(±) 3 .

The wave numbers q i must be rescaled by q, meaning that qi ≡ q i q are kept fixed when taking the q → 0 limit, with the constraint q2 + q3 ≤ 1 (F.28)

following from (154) and from the positivity of q 1 . We give as an intermediate result the expression of an energy denominator in the Ŵ(-) matrix elements, q -( q-q 3 + q 3 ) mc 2 = q→0 γ q3 2

3 4 q3 (1 -q3 ) - q3 R 2 sin 2 α 1 -q3 + O( q5 ) (F.29)
and the energy denominator in the 3 phonon subspace, q -( q 1 + q 2 + q 3 ) mc 2 = + iη][ 3 4 q3 (1 -q3 ) -q3 X sin 2 α 1-q3

+ iη] (F. [START_REF] Kinast | Damping of a Unitary Fermi Gas[END_REF] In the integration over the radius R we have performed the change of variables X = R 2 . This allows to see that the integral over R, that is over X, in (F.35) has a zero imaginary part in the limit η → 0 + , and hence Proof: We first note that the theorem would be trivially proved if the integral over X was taken over R -. Indeed, for all η > 0 and for all X < 0, X [a(X)

+ iη][b(X) + iη][c(X) + iη] < X a(X)b(X)c(X) (F.38)
where the upper bound has a finite integral over R -, so that the dominated convergence theorem can be used to exchange the order of the integral and limit η → 0 + . We can then extend the integration domain in (F.37) to the whole R without changing the result. On the other hand we can use the following identity valid for any η > 0 not necessarily infinitesimal, concerning both the real and the imaginary parts: Θ(u 0 A 2 q-2v) hence the correction (155) to be added to the result (153).

To conclude, one might ask whether there is a physical interpretation to the contributing diagram of order 4 in V, which is the left one (I) in figure F.8. In the calculation of the complex energy shift ∆ Bel q of a phonon q induced by a Beliaev process, this diagram takes into account the effect of a modification of the angular eigenfrequency of the virtual phonons, modification itself induced by a Beliaev process at this order. More quantitatively, let us introduce the Beliaev self-energy (F.10), where we note now more rigorously (2) k the unperturbed eigenenergy of a phonon with wave vector k in Ĥ2 , to distinguish it from the exact energy q (see the note 11). We verified that the change in the complex energy of the phonon q originating from the left diagram (I) is exactly δz (I) q = 1 2 q 3 ( q 3 , qq 3 | Ĥ3 |q ) 2

[ q + iη -( (2) q 3 + (2) q-q 3 )] 2

∆ Bel q 3 ( q + iη -(2) q-q 3 ) + ∆ Bel q-q 3 ( q + iη - (2) q 3 ) (F.45)

The reader will notice that the argument z of ∆ Bel q-q 3 is here the exact energy q minus the non perturbed energy of the virtual phonon that is spectator in the process of diagram (I), that is the one that does not participate to the loop. The sum of the correction (I) and of the usual perturbative Beliaev shift of the q phonon energy can then be written, to the order 4 in Ĥ3 , as a Beliaev shift for a renormalized complex dispersion relation:

∆ Bel
q ( q + iη) + δz (I) q 1 2 q 3 ,q 4 ( q 3 , q 4 | Ĥ3 |q ) 2 δ q 3 +q 4 ,q q + iη -[ (2) q 3 + ∆ Bel q 3 ( q + iη -(2) q 4 ) + (2) q 4 + ∆ Bel q 4 ( q + iη -(2) q 3 )] (F.46)

By the way, the formulation (F.45) provides a second calculation method, alternative to the one of equation (F.34).

With the rescalings of subsection 5.3 suitable to the "small denominators", and by using (F.18) with (q, k) successively equal to (q 3 , k) and (q, q 3 ), where k is the integration variable appearing in the expression taken from (F.10) of ∆ Bel q 3 ( q + iη -(2) q-q 3 ), we find ∆ Bel q 3 ( q + iη -(2) q-q 3 ) mc Here q3 = q 3 /q, k = k/q 3 , θ3 = θ 3 /(γ 1/2 q), θ 3 is the non oriented angle between the vectors q and q 3 . The integral over θ = θ/(γ 1/2 q3 ), where θ is the angle between k and q, has been explicitly performed over the interval [0, A/ q3 ] with the same cut-off parameter A as in equation (F.34). After the insertion of (F.47) in equation (F.45), we obtain Im δz (I) q mc 2 = -

Figure 1 :

 1 Figure1:The factor 1 + Λ F , which is the Fermi gas equation-of-state dependent part in the 2 ↔ 1 phonon coupling amplitude (see equation (83)), is plotted as a function of the interaction strength 1/k F a. The black solid line is obtained from the equation of state measured in[START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. In the BEC (k F a → 0 + ) and BCS (k F a → 0 -) limits the dashed black lines are the asymptotic behaviors given by the Lee-Huang-Yang corrections[START_REF] Lee | Many-Body Problem in Quantum Mechanics and Quantum Statistical Mechanics[END_REF][START_REF] Abrikosov | Concerning a Model for a Non-Ideal Fermi Gas[END_REF][START_REF] Galitskii | The Energy Spectrum of a Non-ideal Fermi Gas[END_REF] (that include the first correction to the mean field energy, scaling as (k F a) 2 on the BCS side and as ρa ρa 3 on the BEC side). For comparison the prediction of the BCS equation of state (12) is plotted as a dash-dotted blue line. Note that 1 + Λ F tends to 1 in the BEC limit, as for a weakly interacting Bose gas. It is equal to 8/9 in the unitary limit, and tends to the same value in the BCS limit, since in both cases µ ∝ ρ 2/3 .

Figure 3 :

 3 Figure3: Low-temperature leading behavior of the total phonon damping rate (solid line) as a function of the rescaled wave number q = cq/k B T for a convex dispersion relation. This is the sum of the Beliaev (124) and Landau (125) rates (dashed lines). The rates are rescaled as in equations (121) and (122).

Figure 4 : 2 ↔

 42 Figure4: 2 ↔ 2 damping rate as a function of the rescaled wave number q= cq/k B T . In order to make the limiting behaviors more apparent, the rate is divided by q2 and plotted as a function of argsh( q)= ln( q + 1 + q2 ). The oblique tangent at q = 0 corresponds to the low-q limiting behavior (see equation (143)). The horizontal straight line is the q → ∞ limit of Γ2↔2 / q2 (see equation (144)) and the asymptotic curve at large q is a fit of Γ2↔2 / q2 by an affine function of 1/ q. The black symbols are numerical calculations of the integral (129). In the numerics, we used a dispersion relation of the Bogoliubov form (149) restricted to an interval [0, q max ] over which its first order derivative is positive. The used dispersion relation coincides with[START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF] up to the third order in q. The discs are obtained by a linear extrapolation of the value of the integral for = 0, and the error bars halfwidth is given by the difference with a quadratic extrapolation.

Figure 6 :

 6 Figure6: We show in relative value the first correction (156) to the hydrodynamic prediction for the phonon damping rate of the unitary gas at zero temperature. We took for the curvature parameter the RPA value γ 0.083769 and the Bertsch parameter ξ B = 0.376 experimentally measured[START_REF] Mark | Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas[END_REF].

  One must then keep in Ĥ3 only the terms in b † b † b, that constitute the raising part Ĥ(+) 3 of the cubic Hamiltonian, or the terms in b † bb , that constitute its lowering part Ĥ(-) 3 . The other terms in b † b † b † and bbb induce nonresonant processes. In Ĥ4 , we shall only keep terms b † b † bb , that are the most favorable potentially resonant terms, as they add a single independent wave vector.

6 ( 2 3X sin α cos α( 1 - 3 (Figure F. 9 : 2 q3 3 (

 6213923 Figure F.9: Thick line: contour integration in the complex plane that should be used to apply the residue theorem, in order to obtain the identity (F.39). The real numbers a 0 , a 1 , b 0 , b 1 , c 0 , c 1 and η are positive. The radius ρ of the half-circle tends to infinity. The crosses correspond to the poles of the integrand in the case they are all simple, but the result holds also in the general case.

  Im I (II) = 0 (F.36) Indeed all the factors in the integrand denominator in (F.[START_REF] Kinast | Damping of a Unitary Fermi Gas[END_REF]) have the form α -βX + iη, with α > 0 and β > 0. As a consequence, the integral over X ∈ [0, +∞[ can be evaluated using the following elementary theorem:Theorem: We define a(X) = a 0a 1 X, b(X) = bb 1 X and c(X) = c 0c 1 X for all X ∈ R. If all the coefficients a 0 , a 1 , b 0 , b 1 , c 0 , c 1 are > 0 and η ∈ R, then lim η→0 + R + dX Im X [a(X) + iη][b(X) + iη][c(X) + iη]= 0 (F.37)

A 2 / q2 0 ( 1 -

 201 X) + iη][b(X) + iη][c(X) + iη] = 0 (F.39)To show it, we use the residue theorem by closing downwards the integration contour following a half-circle of radius ρ → +∞ in the lower half-plane; the behavior of the integrand as 1/X 2 for large |X| allows for it, which leads to a vanishing O(1/ρ) contribution of the half-circle. As the integrand poles are all in the upper half-plane, they are not part in the limit η → 0 + : Im q3 )dX 2 q3 (vu 0 X + iη)[3 4 q3 (1 -q3 ) + iη] u 0 A 2 q-2 + iη)ln(v + iη)

1 -

 1 Θ is the Heaviside function and ln is the usual branch of the complex logarithm, with a branch cut along the real negative axis. We are then left withIm I (I) q3 )(1 -q2 -q3 ) 2 q2 2 q3 Θ(u 0 A 2 q-2v) →

3 iπ 3 1-q3 2 (

 332 ln iη + 3 q3 (1 -q3 ) + 3 q3 3 k(1 -k) -4 q3 X 3 1-q3 iη + 3 q3 (1 -q3 ) -4 q3 X F.[START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF]

  2, with the virtual intermediate phonon plotted as a dashed line. To these six intermediate states λ=I-VI correspond six other diagrams λ=I'-VI', shown on the right part of figure2, where the intermediate phonon has the same wave vector but is annihilated and recreated rather than being created then annihilated. These intermediate states exist only at nonzero temperature since the intermediate phonon must preexist in state |i .4.2.2. Effective 2 ↔ 2 coupling amplitudeThe effective coupling amplitude A 2↔2,eff of the 2 ↔ 2 process is defined by the following writing of the effective Hamiltonian:

Everywhere in this article the coupling amplitudes A m↔n are symmetric functions of their m first arguments, the wavevectors q 1 , q 2 , . . . , q m of the incoming phonons, as well as of their n last arguments, the wavevectors q m+1 , q m+2 , . . . , q m+n of the outgoing phonons. With this property the amplitudes are uniquely determined.

The BEC limit corresponds to x → 0 + and y ∼ -4/x[START_REF] Kurkjian | Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover[END_REF]; in Appendix B we give the expressions (B.12) of the rational fractions in this limit. It is then straightforward to check that our expression of A 2↔1 micro in the BEC limit coincides with equation (16) of reference[START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF] expanded to third order in q and expressed in terms of the Bogoliubov excitation energies of a gas of NB = N /2 bosonic dimers of mass m B = 2m, and speed of sound c.

When T and q tend to zero, the wave numbers q i of the intermediate phonons contributing to Γ q tend to zero, and so do in general the phonon coupling amplitudes appearing in the operator R. Still a validity condition of the result (117) is that ω q Γ q th , with cq th = k B T , that is the angular frequency ω q of the considered mode must greatly exceed the typical relaxation rate of the thermal phonons, which defines the so-called collisionless regime (as opposed to the hydrodynamic regime).

The Taylor contribution of order n is O(q 2n+5 ) as the nth order derivative of the self-energy, that involves an energy denominator to the power n + 1, is of order q 8 /q 3(n+1) accounting for the "small denominators" effect.

We omitted here the θ3 -independent logarithmic real term, in between square brackets in equation (F.47), because its contribution to the final result is clearly real in the limit η → 0 + . The dependence on the cut-off parameter A then disappears. The integral over X 3 = θ2 3 of the bit with the logarithm is calculated by parts, by taking the derivative of the logarithm; one finds that its contribution is zero in the limit η → 0 + by using the identity (F.42). Finally, only the term iπ matters, and it leads exactly to the result (155).

Appendix B. x and y rational fractions in the three-phonon process amplitudes

The rational fractions of the variables x and y that appear in the expression (81) of the microscopic phonon coupling amplitude are given by J(x, y) = 3xy 2xy + 2 (B.1)

A(x, y) = Appendix C. Low wave number out-of-the-energy-shell discrepancy between the microscopic and quantum hydrodynamic Beliaev-Landau coupling amplitudes

In this appendix, as announced in the penultimate paragraph of section 4.1, we clarify the physical origin of the discrepancy of the Beliaev-Landau 2 phonons ↔ 1 phonon coupling amplitudes of the microscopic model and of quantum hydrodynamics in the low wave number limit out of the energy shell: the former indeed diverges whereas the latter tends to zero. To avoid tedious calculations in the fermionic microscopic model, we rather use a bosonic microscopic model, the weakly interacting spatially homogeneous Bose gas, of grand canonical Hamiltonian Ĥ at chemical potential µ, with

We then easily make the link between this microscopic theory and quantum hydrodynamics through the modulusphase representation of the bosonic field operator [START_REF] Mora | Extension of Bogoliubov theory to quasicondensates[END_REF],

where the phase field θ(r) and the density field ρ(r) are canonically conjugated.

In the Bogoliubov U(1)-symmetry breaking formalism, one splits the bosonic field as ψ(r) = ψ 0 + δ ψ(r), where ψ 0 > 0 is the classical field energy minimiser, and one expands the Hamiltonian in powers of δ ψ, Ĥ = H 0 + Ĥ2 + Ĥ3 +. . .. The field δ ψ(r) is then expanded over the normal modes of Ĥ2 (Bogoliubov expansion):

where we omitted for simplicity the contribution of the anomalous modes associated to symmetry breaking and we explicitly indicated that the quasi-particle creation and annihilation operators are relative to a δ ψ-expansion of Ĥ. Insertion of the expansion (C.3) in Ĥ3 leads to the reduced Beliaev-Landau coupling amplitude in the Cartesian δ ψ point of view:

with s i = s q i , s q = U q + V q , d i = d q i and d q = U q -V q .

In the modulus-phase point of view, adapting reference [START_REF] Mora | Extension of Bogoliubov theory to quasicondensates[END_REF] to a symmetry breaking approach, we expand ψ(r) in powers of the density fluctuation δ ρ(r) = ρ(r) -ψ 2 0 and of the phase θ(r) so that δ ψ(r) = δ ψ1 (r) + δ ψ2 (r) + . . ., with

To quadratic order in δ ρ and θ, one replaces δ ψ in Ĥ2 with δ ψ1 ; the normal mode expansion of the corresponding quadratic Hamiltonian in (δ ρ, θ) gives (again neglecting for simplicity the symmetry-breaking related anomalous modes):

U q e iq•r bρθ q + V q e -iq•r ( bρθ q ) † (C.7)

with the same amplitudes (U q , V q ) as in equation (C.3) but with operator-valued coefficients bρθ q , ( bρθ q ) † that coincide with bδψ q , ( bδψ q ) † , as δ ψ1 and δ ψ do, only to first order in the fluctuations. To cubic order in δ ρ and θ, two different contributions are obtained: (i) the obvious one coming from the replacement of δ ψ with δ ψ1 in Ĥ3 , which leads to the Beliaev-Landau coupling amplitude A 2↔1 δψ , this time among the bρθ q , and (ii) the more indirect one, coming from the replacement of δ ψ with δ ψ1 + δ ψ2 in Ĥ2 and extraction of the corresponding cubic contributions of the type δ ψ † 2 δ ψ1

Appendix E. Three-phonon and four-phonon processes in a bosonic model

Here we describe and we implement the microscopic bosonic model used in section 4.2.4 to test the 2 ↔ 2 effective phonon coupling predicted by quantum hydrodynamics. We consider spinless bosons on a cubic lattice of lattice constant l, with a large enough interaction range b so that their excitation spectrum is concave at low wave number q. The Hamiltonian of the lattice model reads

with the interaction potential

of Fourier transform Ṽ(q) = Ṽ0 e -q 2 b 2 /2 with Ṽ0 = (2π) 3/2 b 3 V 0 (E.3)

We have introduced the angular eigenfrequencies of the Bogoliubov quasiparticles

The dispersion relation q → ω q is then concave in the vicinity of q = 0 under the condition, assumed to be satisfied in what follows, b > ξ (E.17) 3rd order. We insert the modal expansion (E.13) in the cubic Hamiltonian (E.9) that we write in the form (93) with the constants m, c and ρ replaced with m B , c B and ρ B respectively, and with coupling amplitudes now given by

with the notations S i ≡ S q i and D i ≡ D q i . 4th order. We insert the modal expansion (E.13) in the quartic Hamiltonian (E.10) that we write in the form ( 108) with (m, c, ρ) → (m B , c B , ρ B ) and a direct 2 ↔ 2 coupling amplitude given by

We considered here the general case where the q i and their opposite are two-by-two distinct, in which case only the first term of ĤB4 contributes. The other terms of ĤB4 , that originate from the expansion (E.5) of n0 in powers of the number of non-condensed particles, contribute to the equation of state beyond Bogoliubov theory [START_REF] Mora | Extension of Bogoliubov theory to quasicondensates[END_REF].

Effective coupling amplitude. We obtain the on-shell effective coupling amplitude of the bosonic model A 2↔2,eff B,OnS from the amplitudes (E.18,E.19,E.20) as prescribed by the equation (105). We then expand it in the limit of small wave vectors and express the result in terms of the angular frequencies ω B q i , i = 1, 2, 3, 4 and ω B q 1 +q 2 , ω B q 1 -q 3 , ω B q 1 -q 4 . We then recover the hydrodynamic prediction corrected by the Landau and Khalatnikov prescription, as explained in section 4.2.4.

Appendix F. Contribution of the Beliaev and higher order processes to the T = 0 decay of the single phonon state Appendix F.1. Presentation of the problem and link with the resolvent of the Hamiltonian We consider the damping rate Γ q of a phonon with wave vector q prepared at zero temperature in a spatially homogeneous gas of fermions. We assume that the phononic excitation branch q → ω q is convex in vicinity of q = 0. The general problem is to determine the behavior of Γ q in the limit q → 0, where one can use an effective low energy theory to describe the coupling among phonons. The dominant process is of course the 1 → 2 Beliaev one. By introducing a dimensionless parameter γ > 0 describing both the correction in q 3 to the hydrodynamic linear excitation spectrum, and the first correction to the hydrodynamic coupling amplitude 1 → 2, we have obtained the provisional result (151) up to the sub-leading order q 7 . In this appendix, we give some details on the obtention of the final result at the Beliaev order (153), and we examine the possible contributions of order q 7 of all the higher order processes, such as the cascade process 1 → 2 → 3 in figure 5, leading to the correction (155) and to the final result (156).

In the subspace with a fixed total momentum q, the one-phonon state |q is the only discrete state, as it is the only one completely characterized by a single wave vector value. It is an eigenstate of Ĥ2 , the part of the Hamiltonian The conclusion of this discussion is that only the order 4 in V has a chance to produce a correction in q 7 to the result Γ Bel q in equation ( 153). We analyze it in more detail in subsection Appendix F.4. We can also confirm that the order 2 in V gives the leading contribution to Γ q in q 5 . This however does not exonerate us from a rigorous study of its sub-leading contributions up to the order q 7 , which is done in subsection Appendix F.3.

Appendix F.3. Study to second order in V

In this subsection we restrict the effective Hamiltonian in equation (F.4) to second order in V and calculate the resulting value of Γ q up to the order q 7 .

To this order in V several different diagrams actually contribute, as V = Ĥ3 + Ĥ4 + Ĥ5 + . . .. The leading order in q can be estimated by power counting. If the two factors V are both equal to Ĥ3 , one has to integrate over at least one independent wave vector and the leading order is q 5 as in equation (F.8). If they are both equal to Ĥ4 , one has to integrate over at least two independent wave vectors and the leading order is q 9 without the "small denominators" effect, but not any better with the "small denominators". If one of the factors V involves Ĥp , with p ≥ 5, the contribution is even smaller and negligible. We can then restrict to V = Ĥ3 , and furthermore to V = Ĥ(±)

3 as the non resonant terms in b † b † b † and bbb of Ĥ3 involve two independent wave vectors and contribute in O(q 8 ).

There remains then a single diagram, the Beliaev one, represented in figure F.7. Let us introduce the corresponding Beliaev self-energy, for any complex number z with Im z > 0:

where we considered the thermodynamic limit and we have used the form (93) of the hydrodynamic Hamiltonian. To this order of approximation, the pole z q associated to the phonon q is the solution of the implicit equation

where the downward arrow denotes the analytic continuation of the self-energy from the upper half-plane Im z > 0 to the lower one Im z < 0 through the branch cut on the real positive axis. Let us separate z Bel q into the real and imaginary parts as in (F.1) and let us expand the right-hand side in powers of Γ q :

It is enough here to truncate the Taylor expansion to the order one; on the other hand the zeroth order is not enough 12 .

The values of the analytic continuation and of its derivative in q are obtained as the limits when η → 0 + of the noncontinuated functions in z = q + iη:

The second term is the ordinary perturbative result. It leads to the damping rate Γ Bel q,pert given by the zero temperature limit of equation ( 119),

and its value for the unitary gas has been calculated up to the order q 7 in equation ( 151). Once the imaginary part of equation (F.13) has been taken, we find a nonnegligible correction to this order, resulting from a first self consistent iteration of z q in the implicit equation (F.11):

[A 2↔1 hydro (k, qk; q)] 2 ( qkq-k + iη) 2 (F. [START_REF] Bartenstein | Crossover from a Molecular Bose-Einstein Condensate to a Degenerate Fermi Gas[END_REF] In order to obtain an equivalent of this derivative when q≡ q/mc → 0, we apply the same technique as in section 5.3. We use spherical coordinates of polar axis q and perform the rescalings k = kq and θ = γ 1/2 qθ on the modulus and polar angle of the vector k, with k < 1 to profit from the "small denominators" effect. In the integrand and hence in the coupling amplitude (94), we take the limit q → 0 with k and θ fixed, with the intermediate results

taking into account the fact that 1 + Λ F = 8/9 at unitarity. The integration with respect to θ over [0, +∞[ is simple after the change of variable X = θ2 , and we are left with

The dominated convergence theorem allows us to take the η → 0 limit in the integrand and to obtain the result (153), knowing that ρ( /mc) 3 = √ 3/(π 2 ξ 3/2 B ). As a side remark, the fact that the derivative with respect to z (F.20) is real in the limit η → 0 + , means that a change in the real dispersion relation q → ω q of order q 5 does not affect Γ q to the order q 7 , which the calculation technique of subsection 5.4.2 also shows.

Appendix F.4. Study to the order 4 in

V

We have seen in subsection Appendix F.2 that the order 4 in V may contribute to Γ q to the order q 7 . Let us write the corresponding correction δΓ q restricting to the leading order, that is keeping Ĥ(-)

3 in the first two factors V and Ĥ(+)

3 in the last two, this sequence being imposed by the presence of the projectors Q:

where we introduced, for ε = ±, the effective coupling operator 1 ↔ 3 to second order

In these expressions, one can use for Ĥ3 the quantum hydrodynamic approximation (93), provided one includes in q = ω q the cubic correction (57) to the linear hydrodynamic spectrum. As the matrix elements of Ĥ3 are real in the Fock basis, the matrices representing Ŵ(±) are obtained one from the other by transposition. After the insertion of a closure relation in the three-phonon subspace, one obtains

where, by virtue of equation ( 93),

A 2↔1 hydro (q 1 , q 2 ; q 1 + q 2 )A 2↔1 hydro (q 1 + q 2 , q 3 ; q) q + iη -( q 1 +q 2 + q 3 ) + A 2↔1 hydro (q 2 , q 3 ; q 2 + q 3 )A 2↔1 hydro (q 2 + q 3 , q 1 ; q) q + iη -( q 2 +q 3 + q 1 ) + A 2↔1 hydro (q 1 , q 3 ; q 1 + q 3 )A 2↔1 hydro (q 1 + q 3 , q 2 ; q) q + iη -(

The contribution to the matrix element of the first term between square brackets, that we note P(q 1 , q 2 |q 3 ), corresponds to the process in figure 5. The other two terms, P(q 2 , q 3 |q 1 ) and P(q 1 , q 3 |q 2 ), are deduced from the first one by circular permutation. As the rest of the summand in (F.23) is invariant by permutation of the three wave vectors, we can replace its numerator [P(q 1 , q 2 |q 3 ) + P(q 2 , q 3 |q 1 ) + P(q 1 , q 3 |q 2 )] 2 by 3[P(q 1 , q 2 |q 3 )] 2 + 6P(q 1 , q 2 |q 3 )P(q 1 , q 3 |q 2 ), which amounts to considering the two diagrams in figure F.8 and the integral expression

hydro (q 1 , q 2 ; q 1 + q 2 )A 2↔1 hydro (q 1 + q 2 , q 3 ; q)] 2 [ q + iη -( q 1 + q 2 + q 3 )][ q + iη -( q 1 +q 2 + q 3 )] 2

+

2A 2↔1 hydro (q 1 , q 2 ; q 1 + q 2 )A 2↔1 hydro (q 1 + q 2 , q 3 ; q)A 2↔1 hydro (q 1 , q 3 ; q 1 + q 3 )A 2↔1 hydro (q 1 + q 3 , q 2 ; q)

where we eliminated the sum over q 1 thanks to momentum conservation, q 1 = q -(q 2 + q 3 ), and took the thermodynamic limit.

In practice, the integration over q 2 and q 3 is performed in spherical coordinates of polar axis q. To evaluate the contribution to δΓ q of the emission cones of q 2 and q 3 of angular aperture O(q) around q, we rescale the polar angles θ i as follows: θi = θ i γ 1/2 q (F.26) then we let q = q/mc tend to zero with θi fixed. The calculation is similar to the one in section 5.3. One introduces the polar representation as in (137

enclosed by the contour, see the figure F.9, hence (F.39) and (F.37) 13 . One might think that the same reasoning applies to the contribution (F.34) and that Im I (I) = 0, in which case δΓ q = o( q7 ) and there would be no correction to add to the result (153). This is not the case because, in order to obtain a finite value of I (I) , we have this time to keep in the integral over X a finite value A 2 / q2 of the upper bound, which prevents the application of the theorem. Here, A > 0 is a cut-off constant whose precise value is not relevant, and the power law in q-2 comes from the fact that θ i ≤ π and then θi ≤ π/(γ 1/2 q) in (F.26), which implies X = R 2 = θ2 2 + θ2 3 = O( q-2 ). If simply replaces A 2 / q2 by +∞ in equation (F.34), one finds indeed that the integral over X diverges as α -2 when α tends to zero, which leads to an integral over α that diverges logarithmically for α = 0. 14 To progress, we cut the integration interval over α into two, a sub-interval [0, ν] for which the cut-off A 2 / q2 in the integration over X is necessary, and a sub-interval [ν, π/2] for which the integration over X can be extended to +∞:

The contribution of the second bit to Im I (I) is zero, because of the theorem. The contribution of the first bit does not depend on the value of ν. We can then calculate it for a value of ν sufficiently small to approximate each contribution depending on α in the integrand by its leading order in α, that is sin α α, cos α 1,

The dependence in φ then disappears and the integration over φ simply gives a factor 2π. After a simple calculation of the integral over α:

we shall concentrate on the integration over X. The fully integrated term α = ν of equation (F.41), after a multiplication by the factor X in the numerator of (F.34) and a division by the denominator (vu 0 X + iη), gives an integrand that is O(1/X 2 ); one can in this case replace the upper bound A 2 / q2 in the integration over X by +∞ and use the following variant of the theorem, in order to show that its contribution to Im I (I) is exactly zero:

The fully integrated term α = 0 of equation (F.41) leads on the contrary to an integral over X of nonzero imaginary 13. We supposed here that the leading contribution to δΓ q comes from the "bicone" configuration in which both q 2 and q 3 are in the forward emission cone of angular aperture O( q) with respect to the direction of q. One can imagine a more subtle scenario, called "unicone", in which only the vector q 3 would be in this cone, while q 2 would be at an angle ≈ q0 with q. In this case, only the denominator of P(q 1 , q 2 |q 3 ) is ≈ q3 , while that of P(q 1 , q 3 |q 2 ) is ≈ q. The crossed term P(q 1 , q 3 |q 2 )P(q 1 , q 2 |q 3 ) is then negligible with respect to [P(q 1 , q 2 |q 3 )] 2 . The global energy denominator q + iη -( q 1 + q 2 + q 3 ) of (F.23) is also ≈ q, which makes us lose a factor q2 , but this is exactly compensated by the loss of a factor q2 in the numerator in the polar integral dθ 2 sin θ 2 . The ensemble seems then to contribute to the same order q7 as the bicone configuration. The integration over the polar angle θ 3 rescaled as in (F.26), with the inclusion of the three-dimensional Jacobian q2 3 from the integration over q3 , however leads to which has a real limit -2 3 when η → 0 + . Hence the sought imaginary part in δΓ q may only come from the global energy denominator, which would provide a Dirac δ[ q -( q 1 + q 2 + q 3 )]; however, to the leading order in q that we consider here, q 3 = q3 q, with 0 < q3 < 1, and the dispersion relation is linear, so that the argument of the Dirac distribution cannot be zero unless q 2 and q are collinear and of the same direction, which is in contradiction with the hypothesis of having a q 2 outside the forward emission cone. The same arguments hold for the damping rate Γ 2↔2 q of section 5.3, and justifies that we only have considered there the "bicone" configuration.

14. This phenomenon does not occur in the contribution (F.35) because neither the factor containing the term X sin 2 α, dangerous when α → 0, nor the one containing the term X cos 2 α, dangerous when α → π/2, are taken to the square. As a consequence, the integral over X from 0 to +∞ diverges only as ln α or ln( π 2 -α) when α → 0 or α → π/2. Mathematically, by replacing in equation (F.35) the integral +∞ 0 dX f (X) by -0 -∞ dX f (X) from the identity (F.39), then by using the dominated convergence theorem, we can justify the exchange of the limit η → 0 + and the integration over q2 , q3 , φ and α, and prove the result (F.36).