Yvan Castin

Basic theory tools for degenerate Fermi gases

-The ideal Fermi gas

We consider in this section a non-interacting Fermi gas in a single spin component, at thermal equilibrium in the grand canonical ensemble. We review basic properties of such a gas. We concentrate on the degenerate regime ρλ 3 ≫ 1, with ρ the density and λ the thermal de Broglie wavelength defined as (1)

λ 2 = 2πh 2 mk B T 1 .
1. Basic facts. -In first quantization, the Hamiltonian of the non-interacting system is the sum of one-body terms,

(2)

H = N i=1 ĥ0 (i)
where N is the operator giving the number of particles, and ĥ0 is the one-body Hamiltonian given here by

(3) ĥ0 = p 2 2m + U ( r )
where p is the momentum operator, m is the mass of a particle and U ( r ) is a position dependent external potential. We shall consider two classes of external potential. Either where ω α is the oscillation frequency of a particle in the trap along one of the d dimensions of space.

The system is made of indistinguishable fermions, all in the same spin state. To implement the antisymmetry of the many-body wavefunction, we move to second quantization. Introducing the orthonormal basis of the eigenmodes φ λ of the one-body Hamiltonian, with eigenenergy ǫ λ , we expand the field operator as (5) ψ( r

) = λ φ λ ( r )a λ
where a λ annihilates one particle in the mode φ λ and obeys anticommutation relations with the creation and annihilation operators:

{a λ , a λ ′ } = 0 (6) {a λ , a † λ ′ } = δ λλ ′ . (7)
The Hamiltonian then reads in second quantized form:

(8) H = λ ǫ λ a † λ a λ ,
and the operator giving the number of particles:

(9) N = λ a † λ a λ .
The system is assumed to be at thermal equilibrium in the grand-canonical ensemble, with a many-body density operator given by (10) σ = Ξ -1 e -β(H-µ N )

where the factor Ξ ensures that σ has unit trace, β = 1/k B T and µ is the chemical potential. The choice of the grand-canonical ensemble allows to decouple one mode from the other and makes the density operator Gaussian in the field variables, so that Wick's theorem may be used to calculate expectation values, which are thus (possibly complicated) functions of the only non-zero quadratic moments

(11) a † λ a λ = n λ = n(ǫ λ ) ≡ 1 exp[β(ǫ λ -µ)] + 1 .
This is the famous Fermi-Dirac law for the occupation numbers n λ . An elementary derivation of Wick's theorem can be found in an appendix of [START_REF] Castin | Simple theoretical tools for low dimension Bose gases[END_REF]. We note that, because of the fermionic nature of the system, the grand-canonical ensemble cannot be subject to unphysically large fluctuations of the number of particles, contrarily to the ideal Bose gas case. Using Wick's theorem we find that the variance of the particle number is [START_REF] Blatt | Theoretical Nuclear Physics[END_REF] Var N =

λ n λ (1 -n λ ),
which is indeed always below the Poisson value N . In the zero temperature limit, for µ not coinciding with an eigenmode energy ǫ λ , one finds a vanishing variance of the particle number, so that the grand-canonical ensemble becomes equivalent to the canonical one.

In this lecture, we shall be concerned with the degenerate limit, where the temperature is much smaller than the chemical potential: [START_REF] Dalibard | Collisional dynamics of ultra-cold atomic gases[END_REF] k B T ≪ µ.

The Fermi-Dirac distribution has then the shape presented in Fig. 1. We shall repeatedly use the fact that the occupation number of holes 1n λ for ǫ λ < µ and the occupation number of particles n λ for ǫ λ > µ are narrow functions of ǫ λµ, of energy width ∼ k B T . More precisely, a standard low-T expansion technique proceeds as follows [START_REF] Diu | Physique Statistique[END_REF]. One writes the Fermi-Dirac formula as [START_REF] Castin | Coherent atomic matter waves[END_REF] n(ǫ) = n T =0 (ǫ) + sign (ǫµ) exp(β|ǫ -µ|) + 1 , singling out as a second term a narrow function of ǫµ. In the case of a continuous density of states ρ(ǫ), an integral over ǫ appears in thermodynamic quantities. One then extends the integration over ǫ to ] -∞, +∞[ for the second term of the right-hand side, neglecting a contribution O[exp(-µ/k B T )], after having expanded in powers of ǫµ the functions of ǫ multiplicating this second term, such as ρ(ǫ).

We illustrate this technique for the calculation of the density ρ in a spatially homogeneous system in the thermodynamic limit: [START_REF] Cohen-Tannoudji | Course at Collège de France[END_REF] ρ hom (µ, T ) =

d d k (2π) d n k = +∞ 0 dǫ ρ(ǫ)n(ǫ)
where we introduced the free space density of states ( 16)

ρ(ǫ) = d d k (2π) d δ(ǫ -h2 k 2 /2m) = d 2ǫ (2π) -d V d ( √ 2mǫ/h) with V d (u)
is the volume enclosed by the sphere of radius u in dimension d. We go to the limit T → 0 for a fixed chemical potential, and we take µ > 0 in order to have a non-zero density at zero temperature. From the rewriting Eq.( 14) of the Fermi-Dirac formula we obtain (17) ρ hom (µ, T ) = ρ hom (µ, 0) +

µ 0 d(δǫ) ρ(µ + δǫ) -ρ(µ -δǫ) exp(βδǫ) + 1 + O e -βµ ,
where we cut the integration over ǫ to 2µ, introducing an exponentially small error. The zero temperature contribution is readily evaluated by the integral in k space: (18) ρ hom (µ, 0) =

h2 k 2 /2m<µ d d k (2π) d = V d (1) (2π) d 2mµ h2 d/2
.

The thermal contribution is obtained as a low-T expansion by expanding ρ(µ ± δǫ) in powers of δǫ: only odd powers of δǫ contribute; in the resulting integrals over δǫ, we extend the upper bound from µ to +∞, paying the price of an exponentially small error. Finally a series expansion with even powers of T is obtained:

(19) ρ hom (µ, T ) = ρ hom (µ, 0) 1 + π 2 24 d(d -2) k B T µ 2 + O (k B T /µ) 4
The bidimensional case d = 2 deserves a particular discussion. One finds that all the powers of T in the expansion have vanishing coefficients, since the density of states is constant, ρ(ǫ) = m/(2πh 2 ). The integral over ǫ in Eq.( 15) can be performed exactly, (20) ρ d=2 hom (µ, T ) = mk B T 2πh 2 ln 1 + e βµ .

For a fixed and negative chemical potential, one sees that the density tends exponentially rapidly to zero when T → 0. For µ > 0 a more convenient rewriting in the low temperature limit is

(21) ρ d=2 hom (µ, T ) = mµ 2πh 2 1 + k B T µ ln 1 + e -βµ .
This explicitly shows that the deviation of the spatial density from the zero temperature value ρ hom (µ, 0) is exponentially small for T → 0 if d = 2.

1 . 2.
Coherence and correlation functions of the homogeneous gas. -We consider here a gas in a box with periodic boundary conditions. The eigenmodes are thus plane waves φ k ( r ) = exp(i k • r )/L d/2 , with a wavevector equal to 2π/L times a vector with integer components. We shall immediately go to the thermodynamic limit, setting L to infinity for a fixed µ and T .

The first-order coherence function of the field is defined as the following thermal average: (22) g 1 ( r ) = ψ † ( r ) ψ( 0 ) . In the thermodynamic limit, it is given by the Fourier transform of the momentum distribution of the gas, which is an isotropic function: (23) g 1 ( r ) = φ(r) ≡

d d k (2π) d n k e i k• r .
It is easy to see that no long-range order exists for this coherence function, even at zero temperature. At T = 0, n k = 1 for k < k 0 and n k = 0 for k > k 0 , where the wavevector k 0 is defined by (24) h2

k 2 0 2m = µ
and is simply the Fermi wavevector, in the present case of zero temperature. In the dimensions for d = 1 to d = 3 this leads to φ 1D (r) = sin k 0 r πr (25)

φ 2D (r) = k 0 2πr J 1 (k 0 r) (26) φ 3D (r) = - 1 2πr
∂ r φ 1D (r), [START_REF] Michael | Methods of modern mathematical physics 1: functional analysis[END_REF] where J 1 (x) is a Bessel function, behaving ascos(x + π/4)(2/πx) 1/2 for x → +∞. One then gets an algebraic decay of φ at large distances, in all dimensions, because of the discontinuity of the momentum distribution. Note that Eq.( 27) holds at all temperatures.

At low but non-zero temperature, the Fermi surface is no longer infinitely sharp but has a smooth variation with an energy width ∼ k B T . This corresponds to a momentum width δk 0 such that (28)

ǫ k0+δk0 -ǫ k0 ≃ h2 k 0 δk 0 m ≡ k B T.
We then expect a decreasing envelope for the oscillating function φ(r) with a decay faster than algebraic with a length scale 1/δk 0 . This can be checked by the following approximate calculation in 1D: we restrict to k > 0 using parity, we single out from n k the zero-temperature contribution θ(k 0k), where θ is the Heaviside function, and we linearize the kinetic energy dispersion relation around k 0 , writing ǫ k0+δk ≃ µ + hk 0 δk/m. Extending the integration over δk to -∞ we obtain [START_REF] Heiselberg | [END_REF] φ 1D (r)φ T =0 1D (r) ≃ -

2δk 0 π sin(k 0 r) +∞ 0 du sin(uδk 0 r) exp(u) + 1
where u originates from the change of variable u = |δk|/δk 0 . Integrating by parts over u, taking the integral of the sine function, gives a fully integrated term exactly compensating the T = 0 contribution, so that one is left with (30) φ 1D (r) ≃ sin(k 0 r) 2πr +∞ 0 du cos(uδk 0 r) cosh 2 (u/2) .

After extension to an integral over the whole real line (thanks to the parity of the integrand), one may close the integration contour at infinity and sum the residues over all the poles with Im u > 0 to get: (31) φ 1D (r) ≃ δk 0 sin(k 0 r) sinh(πδk 0 r) .

We have recovered a result derived in [START_REF] Efetov | [END_REF]. This approximate formula can also be used to calculate the 3D case, by virtue of Eq. [START_REF] Michael | Methods of modern mathematical physics 1: functional analysis[END_REF].

The correlation function of the gas, also called the pair distribution function, is defined as the following thermal average: (32) g 2 ( r ) = ψ † ( r ) ψ † ( 0 ) ψ( 0 ) ψ( r ) .

By virtue of Wick's theorem, g 2 is related to the function φ as (33)

g 2 ( r ) = φ 2 (0) -φ 2 (r)
where φ(0) is simply the mean particle density ρ. A key feature is that g 2 vanishes in the limit r → 0 (quadratically in r): this is a direct consequence of the Pauli exclusion principle, and expresses the fact that one cannot find (with a finite probability density) two fermions with the same spin state at the same location. This spatial antibunching is intuitively expected to lead to reduced density fluctuations for an ideal degenerate Fermi gas, as compared to an ideal Bose gas. This is quantified in the next subsection.

1

.

3. Fluctuations of the number of fermions in a given spatial zone. -We define the operator giving the number of particles within the sphere of radius R:

(34) NR = r<R d d r ψ † ( r ) ψ( r ),
the gas being spatially homogeneous and taken in the thermodynamic limit. We now show that the number NR has subpoissonian fluctuations for a degenerate ideal Fermi gas.

The mean value of the number of particles within the sphere is simply

(35) NR = ρV d (R)
where ρ = φ(0) is the mean density and V d (R) the volume enclosed by the sphere of radius R in dimension d. From Wick's theorem we find a variance (36) Var NR = NR -

r<R d d r r ′ <R d d r ′ φ 2 (| r -r ′ |).
This immediately reveals the subpoissonian nature of the fluctuations.

We quantify this subpoissonian nature in the limit of a large radius R. Replacing φ by its expression as a Fourier transform of n k , see Eq.( 23), we obtain (37) Var NR = NR -

d d k (2π) d d d k ′ (2π) d n k n k ′ F ( k -k ′ )
where F is the modulus squared of the Fourier transform of the characteristic function of the sphere:

(38) F ( q ) = r<R d d r e i q• r 2 .
One can show that F is an integrable function peaked in q = 0 and with a width ∼ 1/R. From the Parseval-Plancherel identity the integral over the whole momentum space of F is (2π) d V d (R). When Rδk 0 ≫ 1, where δk 0 is the momentum width of the finite temperature Fermi surface, see Eq.(28), we can replace F by a Dirac distribution:

(39) F ( q ) ≃ (2π) d V d (R)δ d ( q ).
This leads to

(40) Var NR ≃ V d (R) d d k (2π) d n k (1 -n k ).
To easily calculate the resulting integral to leading order in T , we note that (41)

∂ µ n k = βn k (1 -n k )
so that the T = 0 value of the density only is required, Eq.( 18). We finally obtain

(42) Var NR ≃ NR d k B T 2µ .
The above result does not apply when Rδk 0 ≪ 1. In the range of radius R ≪ (δk 0 ) -1 we may simply perform a T = 0 calculation. The technique to take advantage of the translational and rotational invariance of φ 2 (| rr ′ |) is to introduce the purely geometrical quantity ( 43)

K(X) = r<R d d r r ′ <R d d r ′ δ(| r -r ′ | -X), so that (44) Var NR = ρV d (R) - 2R 0 dX K(X)φ 2 (X).
An exact calculation of K(X) is possible: for X ≤ 2R, one obtains

K 1D (X) = 4(R -X/2) (45) K 2D (X) = πX 4R 2 arccos X 2R -X (4R 2 -X 2 ) 1/2 (46) K 3D (X) = π 2 3 X 2 (X + 4R)(X -2R) 2 . ( 47
)
When combined with Eqs. (25,[START_REF] Michael | Methods of modern mathematical physics 1: functional analysis[END_REF], this leads to exact expressions for the variance in 1D and 3D, in terms of the complex integral Ci and Si functions. The large k 0 R expressions are then readily obtained. The asymptotic expansion of the 2D has to be worked out by hand, singling out the contribution of the low X quadratic expansion of K(X). We obtain:

d = 1 : Var NR = γ + 1 + ln 4k 0 R π 2 + O(1/(k 0 R) 2 ) (48) d = 2 : Var NR = k 0 R π 2 [ln(4k 0 R) + C] + O(ln(k 0 R)/k 0 R) (49) d = 3 : Var NR = 1 2π 2 ln(A 3 k 0 R) (k 0 R) 2 - 1 24π 2 ln(B 3 k 0 R) + O(1/k 0 R) (50) with C = lim u→+∞ π 12 u 3 2 F 3 ( 3 2 , 3 2 ; 2, 5 2
, 3; -u 2 )ln u ≃ 0.656657 (51) . 4. Application to the 1D gas of impenetrable bosons. -Consider a 1D homogeneous gas of bosons interacting via a delta potential gδ(x 1x 2 ) in the limit g → +∞. In this so-called impenetrable boson limit, the interaction potential can be replaced by the contact condition that the many-body wavefunction vanishes when two bosons are in the same point. On the fundamental domain of ordered positions of the N bosons, x 1 < . . . < x N , one then realizes that the eigenwavefunction for bosons coincides with an eigenwavefunction of N non-interacting fermions, with the same eigenenergy. Out of the fundamental domain, the bosonic and fermionic wavefunctions may differ by a sign.

A 3 = exp(γ + 2 ln 2 -1/2) = 4.321100 . . . (52)
As a consequence, the spatial distribution of the impenetrable bosons, being sensitive to the modulus squared of the wavefunction, will be the same as for the ideal Fermi gas. The discussion on the variance of the number of particles in an interval [-R, R] in subsection 1

. 3 immediately applies. On the contrary, the first order coherence function for the impenetrable bosons g B 1 (x) is sensitive to the phase of the wavefunction and will differ from the g 1 fermionic function. More details, obtained with the Jordan-Wigner transformation, are given in [START_REF] Castin | Simple theoretical tools for low dimension Bose gases[END_REF] and references therein. At zero temperature, we simply recall the mathematical fact, showing the absence of condensate, that at large distances [4] (54)

g B 1 (x) ∼ Aρ |k F x| 1/2
where k F = πρ, ρ being the density of bosons, and A = 0.92418 . . .. We also use the fact (that one can check numerically [START_REF] Castin | Simple theoretical tools for low dimension Bose gases[END_REF]) that the power law behavior of g B 1 (x) is the same as the one of the simpler function ( 55)

G(x) = e iπ N[0,x]
where N[0,x] is the operator giving the number of fermions in the interval from 0 to x (here x > 0) and the expectation value is taken in the ground state of the ideal Fermi gas. In this way, one relates the long range behavior of g B 1 to the counting statistics of an ideal Fermi gas: G(x) is the difference between the probabilities of having an even number and an odd number of fermions in the interval [0, x]. One may assume that the probability distribution of the number n of fermions in the interval has a Gaussian envelope [5]. Using the Poisson formula

(56) +∞ n=-∞ f (n) = +∞ n=-∞ f (2πn)
where f (k) = dx f (x) exp(-ikx) is the Fourier transform of the arbitrary function f , and restricting to the leading terms n = 0 and n = 1, we obtain

(57) G(x) ≃ 2 cos(πρx)e -π 2 (Var N[0,x] )/2 .
From the asymptotic expression Eq.( 48) we obtain

(58) G(x) ∝ cos(πρx) (ρx) 1/2
which explains the 1/x 1/2 decay of g B 1 (x).

1 . 5.
In a harmonic trap. -In a harmonic trap, simple approximate results can be obtained in the limit where the chemical potential is much larger than the trapping frequencies ω α times h, in a temperature range to be specified, in the so-called semi-classical approximation. We illustrate this semi-classical approximation for two quantities, a thermodynamic one, the entropy, and a local observable, the density.

1

. 5.1. Semi-classical calculation of the entropy. The entropy S is a physically appealing way to evaluate to which extent a system is cold, as it is a constant in thermodynamically adiabatic transformations, contrarily to the temperature.

To calculate thermodynamic quantities in the grand-canonical ensemble, it is convenient to start from the grand-potential (59) Ω(µ, T ) = -k B T ln Tr e -β(H-µ N ) .

Then one uses the fact that

(60) dΩ = -S dT -N dµ
where N = N is the mean number of particles. Since the various field eigenmodes decouple in the grand-canonical ensemble, one finds

(61) Ω(µ, T ) = k B T dǫρ(ǫ) ln[1 -n(ǫ)]
where n(ǫ) is the Fermi-Dirac formula and ρ(ǫ) is the density of states for a particle in the trap:

(62) ρ(ǫ) = n δ ǫ - d α=1 (n α + 1/2)hω α
where the sum is over all vectors with non-negative integer coordinates. In the semiclassical limit one replaces the sum by an integral

(63) ρ(ǫ) ≃ ρ sc (ǫ) = [0,+∞[ d d d n δ(ǫ - α n α hω α )
where we dropped the ground mode energy ǫ 0 under the assumption ǫ ≫ ǫ 0 . By the change of variables n α = u α ǫ/hω α one obtains

(64) ρ sc (ǫ) = ǫ d-1 α hω α [0,+∞[ d d d u δ 1 - α u α .
The integral can be calculated by integrating over u d and by performing the change of variables

τ 1 = u 1 , τ 2 = u 1 + u 2 , . . . , τ d-1 = u 1 + . . . + u d-1 . We finally obtain (65) ρ sc (ǫ) = ǫ d-1 (d -1)!(hω) d
where ω = ( α ω α ) 1/d is the geometric mean of the trap frequencies.

We replace ρ by ρ sc in Eq.( 61), and we integrate by parts, taking the derivative of the ln[1n(ǫ)], to obtain the semi-classical approximation:

(66) Ω ≃ Ω sc = -(hω) -d +∞ 0 dǫ ǫ d d! n(ǫ).
In the degenerate regime, we use the technique sketched around Eq.( 14), based on the fact that n(ǫ)-n T =0 (ǫ) is a narrow function of the energy, to obtain the low-T expansion:

(67) Ω sc (µ, T ) = - µ d+1 (d + 1)!(hω) d 1 + π 2 6 d(d + 1) k B T µ 2 + . . .
Calculating S = -∂ T Ω and N = -∂ µ Ω, we obtain the semi-classical approximation for the number of particles and for the entropy per particle:

N ≃ µ d d!(hω) d (68) S/N ≃ dπ 2 3 k B T /T F (69) with the Fermi energy (70) k B T F ≡ hω(N d!) 1/d .
What are the validity conditions of this approximate formula for the entropy? The temperature should be T ≪ T F but, contrarily to the thermodynamic limit case, the temperature should remain high enough, to allow the approximation of the density of states (which is a sum of Dirac peaks) by a smooth function. The typical distance between the Dirac peaks at the Fermi energy should be smaller than the thermal energy width ∼ k B T of the Fermi surface. In the case of an isotropic harmonic trap, ω α = ω for all d dimensions, the Dirac peaks are regularly spaced by hω so that we impose (71) k B T > hω.

In the case of incommensurable trap frequencies, the peaks are irregularly spaced and we estimate a typical spacing from the inverse of the smoothed density of states at the Fermi surface: we require

(72) k B T > 1 ρ sc (µ)
,

with µ ≃ k B T F . This leads to the condition (73) k B T ≫ k B T F dN = hω (d!) 1/d dN 1-1/d .
One sees that this is a much weaker condition than in the isotropic case, for d > 1, since we are here in the large N limit. This fact is illustrated by a numerical example in 3D, calculating the entropy for an isotropic trap and an anisotropic trap, see Fig. 2.

1 . 5.2.
Semi-classical calculation of the density. The semi-classical approximation for the density profile in the external potential U ( r ) is strictly equivalent to the so-called local density approximation (LDA), as can be revealed by the writings:

ρ sc ( r ) = d d p (2πh) d n[U ( r ) + p 2 /2m] (74) ρ LDA ( r ) = ρ hom [µ -U ( r ), T ] (75)
where ρ hom is defined in Eq. [START_REF] Cohen-Tannoudji | Course at Collège de France[END_REF] and n(ǫ) is the Fermi-Dirac formula Eq. (11). The philosophies of the two approximations of course differ. In the semi-classical case, one relies on a phase space density of the particles, which is an approximation to the Wigner representation W ( r, p ) of the one-body density operator of the gas [6]:

(76) ρ = 1 exp[β( ĥ0 -µ)] + 1 ,
where the one-body Hamiltonian ĥ0 is defined in Eq.(3). One assumes 69). The chemical potential is fixed to µ = ǫ0 + 47.6h ω, where ǫ0 is the trap ground mode energy. In order to test the accuracy of Eq.( 69), in a situation where the number of particles N rather than the chemical potential µ is known, as it is the case in experiments with atomic gases, we use the exact, numerical value of N to evaluate Eq.( 69) rather than its approximation Eq.( 68).

(77) W ( r, p ) ≡ d d u (2πh) d r - u 2 |ρ| r + u 2 e i p• u/h ≃ 1 (2πh) d 1 exp{β[U ( r ) + p 2 2m -µ]} + 1 . 0 0.
Here N is about 2 × 10 4 .

It remains to integrate this approximation of W ( r, p ) over p to obtain ρ sc ( r ).

In the local density approximation, the gas is considered as a collection of quasimacroscopic pieces that have the same properties as the homogeneous gas of temperature T and chemical potential µ LDA = µ -U ( r ). Note that this last equation intuitively expresses the fact that the chemical potential is uniform in a gas at thermal equilibrium.

These approximations assume that the external potential varies weakly over the correlation length of the homogeneous gas, ∼ 1/k 0 where k 0 ( r ) is the local Fermi wavevector, such that h2 k 2 0 /2m = µ LDA . In a harmonic trap, the scale of variation of the potential may be taken as the spatial extension of the gas, the so-called Thomas-Fermi length R α along direction α such that

(78) µ = 1 2 mω 2 α R 2 α .
The condition k 0 ( r )R α ≫ 1, and correspondingly the semi-classical approximation and the LDA, will fail close to the borders of the cloud where k 0 diverges.

At T = 0, as is apparent from Eq.( 18), these approximations predict the density profiles:

(79) ρ( r ) = V d (1) (2π) d 2m(µ -U ( r )) h2 d/2
.

Integrating this density profile over the domain |r α | ≤ R α gives an expression of the total number of particles as a function of the chemical potential: this expression coincides with the righthand side of Eq.( 68), whatever the dimension of space d, which illustrates the consistency of the various approaches.

It is instructive to see if the LDA for the density profile can be used at T = 0 strictly, or if a lower bound on T is required, as was the case for the entropy. An answer to this question is obtained by calculating exactly the density of the gas and its second order derivatives in the trap center and comparing to the LDA. The value of ρ( 0 ) is simple to obtain from a numerical summation over the harmonic oscillator quantum numbers n α , since the value in x = 0 of the 1D harmonic oscillator eigenstate wavefunction of quantum number n ≥ 0 is known exactly: for n odd, φ n (0) = 0, and for n even, (80) φ

n (0) = (-1) n/2 (n!) 1/2 2 n/2 (n/2)!π 1/4 a -1/2 ho
where the harmonic oscillator length is a ho = (h/mω) 1/2 . This may be obtained from properties of the Hermite polynomials H n , using φ n (x) = N n H n (x/a ho ) exp(-x 2 /2a 2 ho ) with a positive normalisation factor such that N -2 n = a ho π 1/2 2 n n!, or from the recursion relation φ n (0) = -(1 -1/n) 1/2 φ n-2 (0), obtained from the a and a † formalism. The second order derivatives of the density in r = 0 are also easily obtained: they involve φ ′ n (0) = √ 2nφ n-1 (0)/a ho and the second order derivative of φ n , which is related to φ n by Schrödinger's equation, φ ′′ n (0) = -(2n + 1)φ n (0)/a 2 ho . We apply this check in 1D. At zero temperature, explicit expressions can be obtained. For a chemical potential n + 1/2 < µ/hω < n + 3/2, where n is an even non-negative integer, we obtain

ρ(0) = (n + 1)φ 2 n (0) = 1 π 1/2 a ho (n + 1)! 2 n [(n/2)!] 2 (81) ρ ′′ (0) = -2a -2 ho ρ(0). ( 82
)
For a chemical potential n + 3/2 < µ/hω < n + 5/2, where n is still even, the x = 0 density assumes the same value (since φ n+1 (0) = 0) whereas the second order derivative changes sign:

ρ(0) = 1 π 1/2 a ho (n + 1)! 2 n [(n/2)!] 2 (83) ρ ′′ (0) = 2a -2 ho ρ(0). ( 84
)
The T = 0 exact value of the density as a function of the chemical potential is compared to the LDA in Fig. 3a: the LDA does not reproduce the steps but nicely interpolates between them. By using Stirling's formula for the exact result, one finds that the relative error in the LDA is O(1/N ), where N is the particle number. If one cleverly approximates the value of a step in ρ(0) by taking the median value of the chemical potential, µ = (n + 3/2)hω, one even finds that the relative error in the LDA vanishes as O(1/N 2 ). The LDA, on the contrary, gives a totally incorrect prediction for ρ ′′ (0) (not shown in the figure): one has

(85) ρ ′′ LDA (0) = - 1 π 2 a 4 ho ρ LDA (0)
which tends to zero in the large N limit, whereas the exact result oscillates with a diverging amplitude. All this can be understood from the fact that the exact density profile at T = 0 is an oscillating function of x, oscillating at smaller and smaller scales in the large N limit. This fact can be revealed numerically [7], but more elegantly results from the exact summation formula [8]:

(86) m k=0 φ 2 k (x) = (m + 1)φ 2 m (x) -m(m + 1)φ m+1 (x)φ m-1 (x).
The LDA of course misses these oscillations. What happens at finite temperature ? We expect that, at k B T ∼ hω, these oscillations are sufficiently reduced to make the LDA more accurate, while the temperature is low enough to allow the use of the T = 0 limit of the LDA. This is confirmed by the k B T = hω data in Fig. 3: the plateaus in ρ(0) as a function of µ are hardly visible (see a), so are the oscillations of ρ(x) as a function of x [9], and the oscillations of ρ ′′ (0) are so strongly reduced that the LDA approximation becomes acceptable (see b).

-Two-body aspects of the interaction potential

In real gases, there are interactions among the particles. In fermionic atomic gases it is even possible to reach the maximally interacting regime allowed by quantum mechanics in a gas, the so-called unitary regime, without altering the stability and lifetime of the gas. This section reviews possible models that can be used to represent the interaction potential and recalls basic facts of two-body scattering theory. We restrict to a threedimensional two spin-state, single species Fermi gas; the p-wave interactions among atoms of a common spin state are neglected, usually an excellent approximation away from pwave Feshbach resonances. The s-wave interactions are on the contrary supposed to be strong, the scattering length being much larger than the potential range, in the vicinity of a Feshbach resonance. -Central density of a 1D harmonically trapped ideal Fermi gas, in the grand canonical ensemble, as a function of the chemical potential µ. (a) Density in the trap center x = 0: exact result for T = 0 (black solid line with steps), exact numerical sum for kBT = hω (red smooth solid line), vs the T = 0 LDA (dashed line). (b) Second order derivative of the density in the trap center, d 2 ρ/dx 2 (x = 0), for kBT = hω: exact numerical sum (solid line) vs the T = 0 LDA (dashed line).

2

. 1. Which model for the interaction potential ? . -The detailed description of the interaction of two atoms is involved. At large enough interatomic distances, in particular much larger than the size of the electronic orbitals of an atom, one can hope to represent this interaction by a position dependent potential V ( r ), which includes the van der Waals interaction term (87)

V (r) ≃ -C 6 r 6 , a simple formula that actually neglects retardation effects and long range dipole-dipole magnetic interactions. Forgetting about these complications, we can use Schrödinger's equation and the C 6 coefficient to construct a length b, called the van der Waals length, that we shall consider as the 'true' range of the interaction potential:

(88) h2 mb 2 = C 6 b 6 .
For alkali atoms, b is in the nanometer range. At short interatomic distances, however, this simple picture of a scalar interaction potential has to be abandoned, and one has to include the various Born-Oppenheimer potentials curves coming from the QED Hamiltonian for two atomic nuclei at fixed distance, including all the electronic and nuclear spin states, and the motional couplings among the Born-Oppenheimer curves due to the finite mass of an atom.

Fortunately we are dealing with gaseous systems: the mean interparticle distance is much larger than the potential range b:

(89) ρ 1/3 b ≪ 1
where ρ is the mean density. We shall also consider atomic dimers, but these dimers shall be very weakly bound, with a size on the order of the scattering length a, which is much larger than b in the vicinity of a Feshbach resonance. As a consequence, the central postulate for the theory of quantum gases is that the short-range details of the atomic interactions are unimportant, only the low-momentum behavior of the scattering amplitude between two atoms is relevant. For fermions of equal masses, this postulate has proved robust in its confrontation to experiments, even in the unitary regime.

A practical consequence of this postulate is that any short range model for the interaction leading to almost the same scattering amplitude f k as the true interaction, in the typical relative momentum range k of atoms in a quantum gas, is an acceptable model. We then put the constraint on any acceptable model for the interaction:

(90)

f model k ≃ f k
for the relevant values of the relative momenta k populated in the gas. We insist here that we impose similar scattering amplitudes over some momentum range, not just equal scattering lengths. Typical values of k can be the following ones:

(91)

k typ ∈ {a -1 , k F , λ -1 }
where a is the s-wave scattering length between opposite spin fermions, k F = (6π 2 ρ ↑ ) 1/3 is the Fermi momentum expressed in terms of the density of a single spin component ρ ↑ = ρ ↓ , and λ is the thermal de Broglie wavelength Eq.( 1). The choice of the correct value of k typ is left to the user and depends on the physical situation. The first choice k typ ∼ a -1 is well suited to the case of a condensate of dimers (a > 0) since it is the relative momentum of two atoms forming the dimer. The second choice k typ ∼ k F is well suited to a degenerate Fermi gas of atoms (not dimers). The third choice k typ ∼ λ -1 is relevant for a non-degenerate Fermi gas, a case not discussed in this lecture.

A condition for the gas not to be sensitive to the microscopic details of the potential is then that

(92) k typ b ≪ 1.
As we shall see, the scattering amplitude in this momentum range should be described by a very limited number of parameters, which allows to use very simple models in the many-body problem. The situation is more subtle for bosons, or for mixtures of fermions with widely different masses, where the Efimov phenomenon takes place (see the lecture by Gora Shlyapnikov in this volume): at the unitary limit, an effective 3-body attractive interaction occurs and accelerates the atoms; the wavevector k can then, for some type of Feshbach resonances (the broad Feshbach resonances to be defined below), reach the range ∼ b -1 , in which case the details of the true interaction may become important. We shall not consider this case here.

2

.

2. Reminder of scattering theory. -Suppose for simplicity that two particles of mass m interact via the short potential V ( r 1r 2 ) that tends to zero at infinity faster than 1/r 3 12 . Moving to the center of mass frame, with a center of mass with vanishing momentum, one gets the Schrödinger equation at energy E:

(93) Eφ( r ) = - h2 m ∆ r + V ( r ) φ( r ).
At negative energies E < 0, one looks for the discrete values of E such that φ is square integrable: these eigenstates correspond to two-body bound states, that we call dimers. At positive energies E > 0, one sets

(94) E = h2 k 2 0 m
with k 0 > 0, one looks for scattering states that obey the boundary conditions suited to describe a scattering experiment:

φ( r ) = φ 0 ( r ) + φ s ( r ) (95) φ 0 ( r ) = r | k 0 ≡ e i k0• r (96) φ s ( r ) ∼ f k0 ( n ) e ik0r r for r → +∞. ( 97
)
The wavefunction φ 0 represents the free wave coming from infinity, taken here to be a plane wave of wavevector k 0 of modulus k 0 . The function φ s ( r ) represents the scattered wave which, at infinity, depends on the distance as an outgoing spherical wave and possibly depends on the direction of observation n ≡ r/r through the scattering amplitude f k0 .

It is important to keep in mind that the Schrödinger equation with the above boundary conditions is formally solved by

(98) |φ = [1 + G(E + i0 + )V ]|φ 0 where G(z) = (z -H) -1 is the resolvent of the Hamiltonian, z a complex number. Using the identity G = G 0 + G 0 V G, where G 0 (z) = (z -H 0 ) -1
is the resolvent of the free, kinetic energy Hamiltonian, one gets

(99) |φ = [1 + G 0 (E + i0 + )T (E + i0 + )]|φ 0
where we have introduced the T operator (or T -matrix)

(100)

T (z) = V + V G(z)V.
Eq.( 99) allows to obtain the equivalence in momentum space of the asymptotic conditions Eqs.(95,96,97) in real space:

(101)

k |φ = (2π) 3 δ 3 ( k -k 0 ) + 1 E + i0 + -h2 k 2 /m k|T (E + i0 + )| k 0 ,
where the matrix element of T is an a priori unknown smooth, regular function of k. The link between the k-space and the r-space points of view is completed by the following relation:

(102) f k0 ( n ) = - m 4πh 2 k 0 n|T (E + i0 + )| k 0 ,
which follows from the expression of the kinetic energy Green's function for E > 0:

(103) r |G 0 (E + i0 + )| r ′ = - m 4πh 2 e ik0| r-r ′ | | r -r ′ | .
Clearly, the central object is the scattering amplitude. It obeys the optical theorem

(104) σ scatt (k 0 ) = d 2 n |f k0 ( n )| 2 = 4π k 0 Im f k0 ( k 0 /k 0 ),
where σ scatt is the total scattering cross section for distinguishable particles. This theorem is a direct consequence of the conservation of probability: the matter current  in the stationary state φ has a vanishing divergence, that is a zero flux through a sphere of radius r → +∞: using the asymptotic form of φ, one gets the announced theorem. Here we shall consider model potentials that scatter only in the s-wave: φ s is isotropic and so is the scattering amplitude, which reduces to an energy dependent complex number:

(105) f k0 ( n ) = f k0 .
The optical theorem simplifies to

(106) |f k0 | 2 = Im f k0 k 0 .
Using Im z -1 = -Im z/|z| 2 , for a complex number z, one realizes that this simply imposes:

(107)

f k0 = - 1 u(k 0 ) + ik 0
where u(k 0 ) is at this stage an arbitrary real function. The scattering amplitude takes its maximal modulus when u is negligible as compared to k 0 , this constitutes the unitary limit (108)

f k0 ≃ - 1 ik 0 .
A last point concerns the two-body bound states: their eigenergies correspond to poles of the resolvent G, which are also poles of the T matrix. A simple way to find the dimer eigenenergies is therefore to look for poles in the analytic continuation of the scattering amplitude to negative energies, setting k 0 = iq 0 , q 0 > 0, so that E = -h 2 q 2 0 /m < 0. Note that the determination +iq 0 is chosen for k 0 = √ mE/h so that the matrix elements of G 0 , see Eq.( 103), assume the correct value (that tends to zero at infinity). Furthermore one can easily have access to the asymptotic behavior of the dimer wavefunctions φ j , including the correct normalization factor. One uses the closure relation (109)

d 3 k 0 (2π) 3 |φ k0 φ k0 | + j |φ j φ j | = I
where I is the identity operator. One takes the matrix elements of this identity between r | and | r ′ , for large values of r and r ′ so that the asymptotic expression Eq.( 97) may be used. After angular average (for a s-wave scattering), one gets factors e ±ik0(r+r ′ ) and e ik0(r-r ′ ) . Using the optical theorem one shows that the coefficient of e ik0(r-r ′ ) vanishes.

Assuming that u(k) is an even analytical function of k one can extend the integral over

k 0 to ] -∞, +∞[, setting f * k0 = f -k0
, and use a contour integral technique by closing the contour with an infinite half-circle in the Im z > 0 part of the complex plane. From Eq.(107) the residue of f k in the pole k = iq j is obtained in terms of u ′ (iq j ). We finally obtain the large r behavior of any s-wave correctly normalized dimer wavefunction:

(110) φ j (r) ∼ q j 1 -iu ′ (iq j ) 1/2 e -qj r √ 2πr . 2 .
3. Effective range expansion and various physical regimes. -For the quantum gases, it is very convenient to introduce the usual low-k expansion of the scattering amplitude:

(111)

f k0 = - 1 a -1 + ik 0 -k 2 0 r e /2 + . . .
where the parameter a is called the scattering length, the parameter r e is called the effective range of the interaction, not to be confused with the true range b. As discussed below, this expansion is interesting if the omitted terms . . . in the denominator of Eq.( 111) are indeed negligible when Eq.( 92) is satisfied. We take the opportunity to point out that the general allowed values for the effective range can go from -∞ to +∞. This can be demonstrated on a specific example, the square well potential, with V (r) = -V 0 for r < b, V (r) = 0 otherwise, the true range b being fixed and V 0 being a variable positive quantity. Setting κ 0 = (mV 0 /h 2 ) 1/2 , we obtain the exact formulas:

a = b - tan κ 0 b κ 0 (112) r e = b - b 3 3a 2 - 1 κ 2 0 a . ( 113
)
This easily shows that for all non-zero values of κ 0 such that the scattering length vanishes, the effective range r e tends to -∞. When κ 0 tends to zero, the scattering length tends to zero (more precisely, 0 -), quadratically in κ 2 0 , as is evident also from the Born approximation; in the expression for r e , the last two terms diverge to infinity, with opposite sign; a systematic series expansion in powers of κ 0 b leads to ( 114)

r e = - 2b 2 5a + 176 175 b + O(κ 2 0 b 3 ) for κ 0 → 0 + .
One then sees that r e tends to +∞ for κ 0 → 0 + since a → 0 -in this limit. Note that the leading term in Eq.( 114) can also be obtained in the Born approximation.

As already mentioned, a very favorable case is when the . . . in the denominator of Eq.( 111) is negligible as compared to the sum of the first three terms, in the low momentum regime k 0 b ≪ 1: one can then characterize the true interaction by only two parameters, the scattering length and the effective range. In what follows, it is implicitly assumed that this favorable case is obtained, otherwise more realistic modeling of the interaction should be performed. One can then identify interesting limiting cases. The zero-range limit: This is a limit where the only relevant parameter of the true interaction is the scattering length a, i.e. we shall assume that the typical momentum

k typ satisfies (115) k 2 typ |r e | ≪ |a -1 + ik typ |,
where r e is the effective range of the true interaction. In practice, in present experiments and for typical atomic densities, this is true for the so-called broad Feshbach resonances, where r e ∼ b: the condition Eq.( 115) is then a direct consequence of Eq.( 92). For narrow Feshbach resonances, |r e | can be considerably larger than b and Eq.( 115) may be violated for typical densities. What type of model can be chosen in this zero-range limit ? A natural idea is to use a strictly zero-range model, having both a vanishing true range and effective range, and characterized by the scattering length as the unique parameter. This 'ideal world' idea corresponds to the Bethe-Peierls model of subsection 2 . 5, where the interaction is replaced by contact conditions.

In practice, the contact conditions are not always convenient to implement, in an approximate or numerical calculation. It is therefore also useful to introduce a finite range model, usually with a true range of the order of its effective range r model e . One then usually puts the model range to zero at the end of the calculation. Note that this is not always strictly speaking necessary: it may be sufficient that the model is also in the zero-range limit, (116)

k 2 typ |r model e | ≪ |a -1 + ik typ |.
In section 3 on BCS theory we shall restrict to the zero-range limit. We shall then approximate the true interaction by an on-site interaction in a Hubbard lattice model, see subsection 2 . 6: the effective range of the model is then of the order of the lattice spacing, and one may easily take the limit of a vanishing lattice spacing in the BCS theory. The unitary limit: What is the condition to reach the unitary limit in a gas, in a regime where the low-k expansion holds ? One should have

(117) |a -1 |, k 2 typ |r e | ≪ k typ .
The condition on a is in principle possible to fulfill, as a diverges close to a Feshbach resonance. The condition on r e shows that the zero-range limit Eq.( 115) is a necessary condition to reach the unitary limit. This is why the broad Feshbach resonances are favored in present experiments (over the narrow ones) to reach the strongly interacting regime. To make it simple, one may say that the unitary limit is simply the zero-range limit with infinite scattering length.

. 4.

A two-channel model . -We have the oversimplified view of a magnetically induced Feshbach resonance, see Fig. 4: two atoms interact via two potential curves, V open (r) and V closed (r). At short distances these two curves correspond to the electronic spin singlet and triplet, respectively, for atoms of electronic spin 1/2, so that they experience different shifts in presence of an external magnetic field. At large distances, V open (r) tends to zero whereas V closed (r) tends to a strictly positive value V ∞ .

The atoms are supposed to come from infinity in the internal state corresponding to the curve V open (r), the so-called open channel. The atoms are ultracold and have a low kinetic energy,

(118) E ≪ V ∞ .
In real life, the two curves are actually weakly coupled. Due to this coupling, the atoms can have access to the internal state corresponding to the curve V closed (r), but only at short distances; at long distances, the external atomic wavefunction in this so-called closed channel will be an evanescent wave that decays exponentially with r since E < V ∞ . Now assume that, in the absence of coupling between the channels, the closed channel supports a bound state of energy E mol , denoted in this text as the molecular state. Assume also that, by applying a judicious magnetic field, one sets the energy of this molecular state close to zero, that is to the dissociation limit of the open channel. In this case one may expect that the scattering amplitude of two atoms may be strongly affected, by a resonance effect, when a weak coupling exists between the two channels. A quite quantitative discussion is obtained thanks to a very simple two-channel model. This model is the most realistic of the models presented here, and we shall consider it as a 'reality' against which the other models should be confronted. Assuming that the particles are free in the open channel, and can populate only the molecular state in the closed channel, the only non-trivial ingredient is a coupling between the molecule with center of mass momentum zero, |mol , and a pair of atoms of opposite momenta and spin states || k , with

(119) || k ≡ a † k↑ a † -k↓ |0 ,
where the creation and annihilation operators obey the free space anticommutation relation

(120) {a kσ , a † k ′ σ ′ } = (2π) 3 δ( k -k ′ )δ σσ ′ .
This coupling in the momentum representation is defined as

V |mol = d 3 k (2π) 3 Λχ(k)|| k (121) V || k = Λχ(k)|mol . (122)
We have introduced a real coupling constant Λ, and a real and isotropic cut-off function χ(k), equal to 1 for k = 0 and to 0 for k = ∞, to avoid ultraviolet divergencies. It is convenient to take ( 123)

χ(k) = e -k 2 σ 2 /2
where the length σ is a priori of the order of the true potential range b. The two-channel model is thus defined by the three parameters E mol , Λ and σ. It scatters in the s-wave channel only, since χ is isotropic. An eigenstate of energy E of the two-channel Hamiltonian is of the form

(124) |φ = β|mol + d 3 k (2π) 3 α( k )|| k .
Including the molecular state energy and the kinetic energy of a pair of atoms, the stationary Schrödinger equation leads to

h2 k 2 m α( k ) + Λχ(k)β = Eα( k ) (125) E mol β + d 3 k (2π) 3 Λχ(k)α( k ) = Eβ. (126)
We restrict to a scattering problem, E > 0. Instructed by the previous discussion around Eq.(101), we correctly solve for α in terms of β in Eq.( 125):

(127) α( k ) = (2π) 3 δ( k -k 0 ) + Λχ(k) E + i0 + -h2 k 2 /m β,
where k 0 is the wavevector of the incident wave, and E = h2 k 2 0 /m. Insertion of this solution in Eq.( 126) leads to a closed equation for β that is readily solved. From Eq.(101) and Eq.( 102) we obtain the scattering amplitude (128)

f k0 = m 4πh 2 Λ 2 χ 2 (k 0 ) E mol -E + d 3 k (2π) 3 Λ 2 χ 2 (k) E+i0 + -h 2 k 2 /m .
Let us work a little bit on this scattering amplitude. First, using the identity

(129) 1 X + i0 + = P 1 X -iπδ(X)
where P denotes the principal value, one can check that f k0 is indeed of the form Eq.( 107), so that it obeys the optical theorem. Second, the scattering length is readily obtained by taking the k 0 → 0 limit:

(130) -a -1 = 4πh 2 m E mol Λ 2 - d 3 k (2π) 3 mχ 2 (k) h2 k 2 .
This shows that the location of the resonance (where a -1 = 0) is not E mol = 0 but a positive value E 0 mol , this shift being due to the coupling between the two channels. For E mol larger than this value, a < 0. For E mol below this value, a > 0.

Another general result, not specific to a Gaussian choice for χ(k) (but assuming that χ(k) varies quadratically around k = 0), concerns the zero range limit σ → 0 for fixed values of k 0 , Λ and the scattering length a (note that the molecular energy is then not fixed):

(131) f k0 → - 1 a -1 + ik 0 + 4πh 2 m E Λ 2 .
This can be seen formally as resulting from the identity +∞ 0 dX P1/(X 2 -1) = 0. This shows that the two-channel model has a well defined limit, in the zero true-range limit, characterized only by the scattering length and by a non-zero and negative effective range, (132)

r 0 e = - 8πh 4 m 2 Λ 2 .
This limit may be described by generalized Bethe-Peierls contact conditions [10]. This illustrates dramatically the difference between the true range and the effective range of an interaction. This also shows that, in this limit, the two-channel model can support a weakly bound state, that is a dimer with a size much larger than the true range σ.

Setting k 0 = iq 0 , one finds that the right hand side of Eq.( 131) has a pole with q 0 real and positive if and only if a > 0:

(133)

q 0 = 1 -1 -2r 0 e /a r 0 e .
Finally, an exact calculation of the scattering amplitude for the Gaussian choice Eq.( 123) is possible: the trick put forward by Mattia Jona-Lasinio is to calculate the analytic continuation to imaginary k 0 = iq 0 , so that one no longer needs to introduce the principal part distribution. One gets

(134) (f iq0 ) -1 = -a -1 - 1 2
r 0 e q 2 0 e -q 2 0 σ 2 + q 0 erfc(q 0 σ),

where erfc is the complementary error function, that tends to unity in zero, erfc(x) = 1 -2π -1/2 x + O(x 3 ). The scattering length is given by

(135) -a -1 = 4πh 2 E mol mΛ 2 - 1 σ √ π .
The calculation of the effective range for a finite σ is then straightforward:

(136)

r e = r 0 e + 4σ √ π 1 - √ πσ 2a .
We shall then classify the Feshbach resonances as follows: for a = ∞, the ones with r e ≃ r 0 e < 0, much larger than the true range in absolute value, constitute narrow Feshbach resonances. The ones with r e ≃ σ constitute broad Feshbach resonances. This classification is not only theoretical but also reflects the magnetic width ∆B of the resonance. If one assumes that E mol is an affine function of the magnetic field B, with a slope -µ, one finds from Eq.( 130) that the scattering length a of the two-channel model varies with B as

(137) a = a bg ∆B B 0 -B ,
where B 0 is the magnetic field value of the resonance center, and we have singled out the value a bg of the background scattering length (that is the value of the scattering length of the true interaction far from the resonance), motivated by the fact that a = a bg [1 + ∆B/(B 0 -B)] in a more complete theory including the direct interaction in the open channel [11]. This defines the width ∆B, such that

(138) µ∆B = 2h 2 ma bg r 0 e .
So, when a bg ∼ σ ∼ b, comparing µ∆B to the natural energy scale h2 /mb 2 amounts to comparing |r 0 e | to b.

. 5.

The Bethe-Peierls model . -In this model, the true interaction is replaced by contact conditions on the wavefunction. For two particles in free space, in the center of mass frame, the contact condition is that there exists a constant A such that

(139) φ( r ) = A r -1 -a -1 + O(r)
in the limit r → 0. For r > 0 the wavefunction is an eigenstate of the free Hamiltonian:

(140)

Eφ( r ) = - h2 m ∆ r φ( r ).
In mathematical terms, this defines the domain of the Hamiltonian. The Hilbert space remains however the same, with the usual scalar product, because a 1/r divergence in 3D is square integrable. An equation valid for all values of r can be written using the theory of distributions:

(141)

Eφ( r ) = - h2 m ∆ r φ( r ) + 4πh 2 a m φ reg δ( r )
where the so-called regular part of φ is (142)

φ reg = - A a = ∂ r [rφ( r )] r=0 = -a -1 lim r→0 rφ( r ).
This establishes the equivalence with the regularized pseudo-potential model [START_REF] Blatt | Theoretical Nuclear Physics[END_REF][START_REF] Dalibard | Collisional dynamics of ultra-cold atomic gases[END_REF][START_REF] Castin | Coherent atomic matter waves[END_REF][START_REF] Cohen-Tannoudji | Course at Collège de France[END_REF]. The last form of φ reg is convenient to use in the limit a = ∞, since a -1 simplifies with the factor a in front of the Dirac δ( r ).

The solution of Eq.( 140) with the boundary conditions Eq.( 97) and the contact condition is simply:

(143) φ( r ) = e i k0• r + f k0 e ik0r r
with the scattering amplitude (144)

f k0 = - 1 a -1 + ik .
This shows that the model has both a zero true range and zero effective range. It may be applied to mimic the effect of the true interaction potential when the condition Eq.( 115) is satisfied. Following the discussion above Eq.( 110), one finds that the model supports a dimer if and only if a > 0. The pole of f iq0 is then q 0 = 1/a, resulting in a dimer energy (145)

E 0 = - h2 ma 2 .
Since the model has a zero true range, it turns out that the dimer wavefunction coincides everywhere with its asymptotic form Eq.( 110),

(146) φ 0 (r) = 1 √ 2πa e -r/a r .
The advantage of the Bethe-Peierls model for the N -body problem is that it introduces the minimal number of parameters (the scattering length). In the unitary limit, it is a zero-parameter model, which is ideal to describe universal states. This advantage actually really exists if one restricts to analytical solutions of the model (numerical solutions tend to introduce extra parameters like a discretization step or a momentum cut-off). Analytical solutions in free space are available up to N = 3 in the unitary limit [18,19]. For particles trapped in harmonic potentials, analytical solutions are available for two particles at any a for isotropic traps [16] and cylindrically symmetric traps [17]. For three particles in an isotropic harmonic trap, the problem was solved analytically only in the unitary limit a = ∞ [20,21]. In the many-body problem, for the universal unitary gas in an isotropic harmonic trap, a scaling solution can be found for the many-body wavefunction in the time dependent case [22] and an exact mapping to free space zero energy eigenstates can be constructed [23,24,25].

It is worth mentioning that, for N = 3 bosons, the Bethe-Peirls model does not define a self-adjoint operator and has to be supplemented by extra contact conditions; this is related to the Efimov phenomenon [18,19]; the problem of course persists for N > 3, making the model non satisfactory. Fortunately this lecture is devoted to equal mass spin 1/2 fermions, for which no Efimov effect appears; the corresponding unitary gas is then called universal, to emphasize this aspect that no extra parameter has to be added to the Bethe-Peirls model. Note that this absence of Efimov effect is not simply due to the fact that the Pauli principle prevents one from having three fermions of spin 1/2 in the same point of space: if the spin up and spin down fermions have widely different masses, an Efimov effect may appear [26].

We take the opportunity to mention another trap to avoid, even for fermions. With the formulation of the model in terms of binary interactions with the pseudo-potential V PP ( r ) = gδ( r )∂ r (r •), where g = 4πh 2 a/m, it is easy to 'forget' that contact conditions are imposed on the many-body wavefunction. One may then be tempted to perform a variational calculation with a N -body trial wavefunction which does not satisfy the contact conditions Eq.(139), that is which is not in the domain of the Hamiltonian [START_REF] Michael | Methods of modern mathematical physics 1: functional analysis[END_REF]. Whereas in the low a-limit the corresponding prediction for the energy may make sense as a perturbative prediction, it may become totally wrong in the large a limit.

To illustrate this statement, let us consider two identical particles in a harmonic isotropic trap interacting via the pseudo-potential. After separation of the center of mass coordinates, one is left with the following Hamiltonian for the relative motion

(147) H = - h2 2µ ∆ r + 1 2 µω 2 r 2 + gδ( r )∂ r (r •),
where µ = m/2 is the reduced mass. Let us take as a trial wavefunction |φ t the Gaussian of the ground state of the harmonic oscillator, (148)

φ t (r) = e -r 2 /(4a 2 ho ) (2π)
E n | ψ n |φ t | 2 .
From [16] and for a finite a one finds that the overlap squared scales as 1/n 3/2 whereas the energy E n scales as n, for large n. The correct mean energy of φ t is thus infinite.

If one however falls in the trap, one gets the wrong mean energy (149)

E wrong = φ t |H|φ t = 3 2 hω + 2 π 1/2
hω a a ho .

In the low a limit this can be compared to the expansion of the exact minimal positive eigenenergy [16] (150)

E exact = 3 2 hω + hω 2 π 1/2 a a ho + 2 π (1 -ln 2) a a ho 2 + . . . .
One sees that the wrong variational calculation gives the correct linear correction in a, that is the correct perturbative result. One may then be tempted to use the wrong variational calculation to predict the sign of the coefficient of a 2 ; one would naively assume E exact ≤ E wrong and one would wrongly predict a negative sign for the coefficient of a 2 [28]. In the opposite limit of an infinite scattering length, one directly realizes the absurdity of the assumption E exact ≤ E wrong : for a → -∞, one finds E wrong → -∞, whereas E exact → hω/2, with the exact wavefunction ∝ e -r 2 /4a 2 ho /r. Note that the same pathology occurs for the many-body problem. Consider for example a two spin-component spatially homogeneous Fermi gas, with interactions modeled by the pseudo-potential. If one uses the Hartree-Fock variational ansatz, which is not in the domain of the Hamiltonian, one correctly obtains, in the weakly attractive regime k F a → 0 -, the mean field shift of the ground state energy, but one wrongly predicts a negative sign for the coefficient of a 2 in the low k F a expansion of the ground state energy, by wrongly arguing that the exact ground state energy has to be below the 'variational' Hartree-Fock energy. As can be checked on the systematic expansion [START_REF] Heiselberg | [END_REF] the coefficient of a 2 is actually positive, because 11 > 2 ln 2.

One should not leave this section with the impression that any variational calculation is impossible within the Bethe-Peierls model. To be safe, one simply has to take a variational ansatz which is in the domain of the Hamiltonian. An example of a successful variational calculation is the derivation of the virial theorem for a unitary gas in a non necessarily isotropic harmonic trap; this derivation, due to Frédéric Chevy, is reported in [25], and uses the fact, for the universal unitary gas 1/a = 0, that the function of the rescaled coordinates Ψ( r 1 /λ, . . . , r N /λ), where λ > 0, is in the domain of the Hamiltonian if the wavefunction Ψ( r 1 , . . . , r N ) is in the domain of the Hamiltonian.

2

.

6. The lattice model . -The last model we consider is at the intersection of two very popular classes of models in condensed matter physics, the separable potentials and the lattice models. It is slightly simpler than the two-channel model, but it does not apply to narrow Feshbach resonances in a situation where the effective range term cannot be neglected in the scattering amplitude. It is more tractable than the Bethe-Peierls model for variational treatments (like the BCS method, see section 3) or for exact numerical treatments since the domain of the Hamiltonian is the one of usual quantum mechanics. In particular, since it is a sort of Hubbard model, the machinery of quantum Monte Carlo methods applicable to the Hubbard model [30,31,32,33] could be applied to it.

In this model, the spatial coordinates r of the particles are discretized on a cubic grid of step l. As a consequence, the components of the wavevector of a particle have a meaning modulo 2π/l only, since the function r → exp(i k • r ) defined on the grid is not changed if a component of k is shifted by an integer multiple of 2π/l. We shall therefore restrict the wavevectors to the first Brillouin zone:

(151) k ∈ D ≡ - π l , π l 3 .
This also shows that the lattice structure in real space provides a cut-off in momentum space.

The interaction between opposite spin particles takes place when two particles are at the same lattice site, as in the Hubbard model. In first quantized form, it is represented by a discrete delta potential:

(152) V = g 0 l 3 δ r1, r2 .
The coupling constant g 0 is a function of the grid spacing l. It is adjusted to reproduce the scattering length of the true interaction. The scattering amplitude of two atoms on the lattice with vanishing total momentum is given in [34], and a detailed discussion is presented in [START_REF] Castin | Simple theoretical tools for low dimension Bose gases[END_REF]. We give here only the result, generalized to an arbitrary even dispersion relation k → E k for the kinetic energy on the lattice:

(153)

f k0 = - m 4πh 2 1 g -1 0 -D d 3 k (2π) 3 1 E+i0 + -2E k .
Adjusting g 0 to recover the correct scattering length then gives (154) (157)

1 g 0 = 1 g - D d 3 k (2π) 3 1 2E k , with g = 4πh 2 a/m.
r model e = l 2π 3 k / ∈[-1/2,1/2[ 3 d 3 k k 4 = 0.336 868 47 . . . l.
In the l → 0 limit, for a fixed a (not fixed g 0 ) the lattice model then reproduces the scattering amplitude of the Bethe-Peierls model, (158) lim

l→0 f k0 = - 1 a -1 + ik 0 ,
and admits a dimer for a > 0. Note that g 0 indeed assumes negative values that tend to 0 -linearly with l in this limit: the interaction is attractive, whatever the sign of the scattering length ( 1 ).

Studying the weakly interacting regime: In the opposite case of l ≫ |a|, the bare coupling constant is (159) g 0 ≃ 4πh 2 a m , so that g 0 has the sign of the scattering length. For a > 0, the model interaction is repulsive, for a < 0 it is attractive. In both cases, it does not support any two-body bound state. Its scattering length can be calculated in the Born approximation, since the term in a in the denominator of Eq.( 155) is small as compared to unity. This makes it an appropriate model to study the weakly interacting regime, as was done for Bose gases in [34]: the Hartree-Fock theory, for example, which relies on the Born approximation, may be applied, and its accuracy may be supplemented by keeping higher order terms in a perturbative expansion of the interaction. In the context of fermions, it may be used in the weakly attractive or weakly repulsive regime; in the latter case, it may be an alternative to the celebrated hard sphere model (with hard spheres of diameter a > 0). It is instructive to apply the philosophy around Eq.( 90) to understand why the lattice model with l ≫ |a| is restricted to the weakly interacting regime, for a quantum degenerate gas of atoms with k typ = k F .

( 1 ) More generally, the lattice model admits a dimer if g0 < g c 0 , where g c 0 is the a → ∞ limit of Eq.( 155). This condition is equivalent to g0 < 0 and a > 0.

• In an effort to have f model k ≃ f k , one first adjusts the scattering length of the model to the correct value a by tuning g 0 .

• If one further adjusts the lattice spacing l to have r model e = r e : this assumes a broad Feshbach resonance with r e ∼ b, since the narrow Feshbach resonance has r e < 0; one then finds l ∼ b; since k typ b ≪ 1 (otherwise microscopic details of the interaction will come in, and the . . . in Eq.( 111) is not negligible), having |a| ≪ l leads to k F |a| ≪ 1, which is the weakly interacting regime.

• Alternatively, one may consider the more favorable case of a true interaction being in the zero-range limit Eq.( 115). The model should also be in the zero-range limit Eq.( 116), which can be written here

(160) k F l ≪ 1 k F a + i ,
and this is sufficient (r model e ∝ l may differ a lot from r e ). Since we wish to have |a| ≪ l, the above equation leads to k F |a| ≪ |(k F a) -1 + i|, a condition that may be satisfied only in the weakly interacting limit k F |a| ≪ 1.

• The same result may be obtained in the following much faster way: since we are dealing here with fermions, there will be at most one fermion per spin component per lattice site, so that k F l < (6π 2 ) 1/3 ; the condition |a| ≪ l then immediately leads to k F |a| < 1. This reasoning is however specific to the lattice model, whereas the previous reasoning is easily generalized to the hard sphere model (where r model e ∝ a).

The 'true' Hubbard model: To be complete, we consider a second case for the dispersion relation, the one of the true Hubbard model: this makes the link with condensed matter physics, and this also shows that a universal quantum gas may be described by an attractive Hubbard model in the limit of a vanishing filling factor. The Hubbard model is usually presented in terms of the tunneling amplitude between neighboring sites, here t = -h 2 /2ml 2 , and of the on-site interaction U = g 0 /l 3 . The dispersion relation is then

(161) E k = h2 ml 2 α [1 -cos(k α l)]
where the summation is over the three dimensions of space. Close to k = 0, the Hubbard dispersion relation by construction reproduces the free space one h2 k 2 /2m. The explicit version of Eq.( 154) is obtained from Eq.(155) by replacing the numerical constant K by K Hub = 3.175911 . . .. At the unitary limit this leads to U/t = 7.913552 . . ., which corresponds to an attractive Hubbard model since t < 0, lending itself to a Quantum Monte Carlo analysis for equal spin populations with no sign problem [30,31,32]. We have also calculated the effective range of the Hubbard model, which remarkably is negative:

(162) r Hub e ≃ -0.3056l.

.

7. Application of Bethe-Peierls to a toy model: two macroscopic branches. -In this subsection, we move to the problem of N interacting fermions, N/2 with spin component ↑ and N/2 with spin component ↓. In order to have a global view of what the BCS theory will be useful for, it is instructive to start with a purely qualitative model [START_REF] Pricoupenko | [END_REF].

Consider a matter wave of isotropic wavefunction φ(r) in a hard wall spherical cavity of radius R, in presence of a point-like scatterer of fixed position in the center of the cavity. The effect of the point-like scatterer is characterized by the scattering length a and is treated by imposing the Bethe-Peierls contact condition Eq.( 139).

What is the link between this model and the many-body problem of N interacting fermions ? The wavefunction φ(r) describes the relative motion of a, let us say, spin ↑ atom, with respect to the nearest spin ↓ atom modelized by the point-like scatterer. The cavity represents (i) the interaction effect of the other N/2 -1 spin ↓ atoms and (ii) the Pauli blocking effect of the other N/2 -1 spin ↑ atoms. Interaction effect (i): the radius of the cavity should then be of the order of the mean interparticle separation in the gas,

(163) R ∝ 1 k F .
Pauli blocking effect (ii): in the case g = 0, the zero point energy of φ should be of the order of the Fermi energy so that one has also the choice (163). This explains also the choice of hard walls for the cavity; in the case of bosons, imposing that the derivative ∂ r φ vanishes in r = R [47], or taking a cubic cavity with periodic boundary conditions would be appropriate. Finally, the total energy of the gas is related to ǫ by (164) E ∝ N ǫ.

Specific but arbitrary proportionality coefficients are used in Eqs.(163,164) to produce the Fig. 5 to come, as detailed in [START_REF] Pricoupenko | [END_REF]. Let us proceed with the calculation of the eigenenergies of the matter wave in the cavity. An eigenmode of the cavity with energy ǫ solves Schrödinger's equation for 0

< r < R (165) - h2 m ∆φ(r) = ǫφ(r)
with the Bethe-Peierls contact condition Eq.(139). If ǫ > 0, we set ǫ = h2 k 2 /m and k solves (166) tan kR = ka.

If ǫ < 0, we set ǫ = -h 2 κ 2 /m and κ solves (167) tanh κR = κa.

-10 -5 0 5 10 In the cavity, there is an infinite number of discrete modes. To each of this mode we associate a distinct macroscopic state of the gas. In figure 5 we have plotted the energy per particle and the pressure P = -∂ V E, where V = N/ρ is the total volume of the gas, for the first two branches, as functions of -1/k F a. We have taken -1/k F a as the abscissa because it allows an almost direct mapping with the B field axis in a real experiment around the location B 0 of the Feshbach resonance, see Eq.(137).

-1/k F a -2 -1 0 1 2 E/N [E F ] (a)
The first excited branch is metastable. It starts with a weakly repulsive Fermi gas on the extreme left of the figure and has a larger energy than the ideal Fermi gas, indicating effective repulsion. When 0 < a ≪ k -1 F , three-body collisions (not included in the toy model) lead in a real gas to the formation of dimers φ 0 : in the language of the toy model, the system starts populating the ground branch. This suggests a first experimental way to produce a gas of dimers.

The ground branch continuously connects the weakly attractive Fermi gas (on the right) to a gas of dimers (on the left). This provides a second experimental way to produce a condensate of dimers, by adiabatic following of the ground branch. The sharp decrease of the total energy on the left part of the figure reflects the 1/a 2 dependence of the dimer binding energy ǫ 0 . The pressure is always less than the Fermi pressure of the ideal Fermi gas, indicating effective attraction with respect to the ideal Fermi gas. Note that P drops very rapidly on the left side, due to the absence of interaction between the molecules in the toy model. In real life, this interaction occurs with a scattering length that was calculated by solution of the 4-body problem [48]. As discussed in other lectures of this volume, the molecular condensates have been obtained experimentally [39,40,41,42].

The regime of infinite scattering length |a| = +∞ is universal in the toy model: one finds that the energy of the corresponding unitary gas has to be proportional to the ground state energy of the ideal Fermi gas:

(168)

E unitary = ηE ideal 0 where E ideal 0 = 3N E F /5, η is a numerical constant depending on the branch, and E F is the Fermi energy of the ideal gas. In the present state of the art of the field, it is accepted that this universality property also holds for a real ground state Fermi gas. Approximate fixed node Monte Carlo calculations [35,36] give for the ground branch η ≃ 0.4, in disagreement with early experiments [37] but in good agreement with recent experimental measures [38,39,43] and with recent exact Quantum Monte Carlo results [33].

The toy model allows to easily solve the following paradox:

• A common saying is that a gas with a positive scattering length a > 0 experiences effective repulsion, and that it experiences effective attraction for a < 0.

• Another saying is that for |a| = +∞ the energy of the gas is universal and does not depend on the sign of a.

• However, if a → +∞ (scenario 1), one expects that the universal state has effective repulsion, whereas for a → -∞ (scenario 2), one expects that the universal state experiences attraction. How can it be that there is no dependence on the sign of a in the unitary limit?

The answer provided by the toy model is simple: scenario 1 and scenario 2 are predicted to lead in the unitary limit to two different universal states, belonging respectively to the first excited branch and to the ground branch, one experiencing effective repulsion, the other effective attraction.

-Zero temperature BCS theory: study of the ground branch

In this section, we consider the many-body problem of a two-component Fermi gas with an interaction characterized only by the s-wave scattering length a.

To be able to use the BCS theory for all values of k F a, while getting sensible results, we restrict to the zero temperature case: the usual static BCS theory is indeed unable to get a fair approximation to the critical temperature on the bosonic side of the Feshbach resonance (a > 0, k F a ≪ 1), for reasons that will become clear at the end of this section. An improved finite temperature theory can be found in [44,[START_REF] Randeria | Bose-Einstein Condensation[END_REF].

For simplicity we also assume that the two spin states ↑ and ↓ have the same number of particles, so that they have a common chemical potential µ. The case of imbalanced chemical potentials and particle numbers is presently the object of experimental and theoretical studies, and may lead to a variety of observed and/or predicted phenomena, as the not yet observed non-demixed BCS phases with spatially modulated order parameter (the so-called LOFF phases) [49,50] or the already observed spatial demixing of the two spin components [51,43,52,53].

Going beyond the toy model of previous section, the zero temperature BCS theory predicts that the particles arrange in pairs and gives several properties of the gas of pairs:

• The existence of long range order in the pair coherence function:

g pair 1 ( r ) = Ψ † ( r ) Ψ( 0 )
where Ψ( r ) = ψ↓ ( r ) ψ↑ ( r ) annihilates a pair of particles with opposite spin in r. In the thermodynamic limit, the zero temperature BCS theory predicts that g pair 1 has a non-zero limit for large r: the pairs form a condensate. Except in the regime k F a → 0 + , the pairs are not bosons, so that this condensate is not simply a Bose-Einstein condensate.

• The BCS theory predicts an energy required to break a pair.

• The time-dependent BCS theory also predicts collective modes for the motion of the pairs, like sound waves, associated to the famous Anderson-Bogoliubov phonon branch [START_REF] Bogoliubov | New Method in the Theory of Superconductivity[END_REF][START_REF] Anderson | [END_REF]. In a trapped system, the time dependent BCS theory predicts the equivalent of superfluid hydrodynamic modes.

In what follows, we shall use the lattice model of subsection 2 . 6. In second quantized form we therefore take the grand-canonical Hamiltonian (169)

H = k∈D σ=↑,↓ h2 k 2 2m -µ a † kσ a kσ +l 3 r,σ U ( r ) ψ † σ ( r ) ψσ ( r )+g 0 l 3 r ψ † ↑ ( r ) ψ † ↓ ( r ) ψ↓ ( r ) ψ↑ ( r ).
The external potential U ( r ) may be accompanied by a rotational term if one wishes to study vortices. Here the a kσ 's obey the usual discrete anticommutation relations Eqs. (6,7). The field operators ψσ ( r ) obey anticommutation relations mimicking the ones in continuous space in the limit l → 0:

(170) { ψσ ( r ), ψ † σ ′ ( r ′ )} = l -3 δ r r ′ δ σσ ′ .
As discussed in subsection 2 . 6, this lattice model is very close to the usual Hubbard model of condensed matter physics. What is unusual here is that we take a quadratic dispersion relation, and more important the model makes sense (to describe a real atomic gas with continuous positions) in the low filling factor limit, an unusual limit in solid state physics. We note that the original Hubbard model (with the usual cosine dispersion relation Eq.( 161) and with a filling factor of the order of unity) may be realized in a real experiment with atoms in an optical lattice [56,57,58,59]. where α is a complex number and a † creates a boson in some normalized one-body wavefunction φ( r ), ( 172)

a † = r l 3 φ( r ) ψ † ( r ).
By analogy, the BCS theory [60], which is a variational theory, takes as a trial state a coherent state of pairs:

(173) |Ψ BCS = N e γC † |0
where N is a normalization factor, γ is a complex number and C † creates two particles in the normalized pair wave function φ( r 1 , r 2 ):

(174)

C † = r1, r2 l 6 φ( r 1 , r 2 ) ψ † ↑ ( r 1 ) ψ † ↓ ( r 2 ).
Note that in general C and C † do not obey bosonic commutation relations. For simplicity, we omit terms creating two particles in the same spin state. This shall be sufficient for our purposes as we restrict here to the case of balanced spin state populations.

3

. 1.2. A more convenient form from the Schmidt decomposition. To briefly review the main properties of this ansatz, it is convenient to introduce the Schmidt decomposition of γ|φ , which is, in the physics of entanglement, routinely applied to the state vector of two arbitrary quantum particles:

(175) γ|φ = α Γ α |A α ⊗ |B α
where the coefficients Γ α are real numbers, the set of |A α is an orthonormal basis here for a spinless particle, and the set of |B α is also an orthonormal basis for a spinless particle. Note that in general these two basis are distinct. As |φ is normalized to unity, one has α Γ 2 α = |γ| 2 . Then the pair creation operator takes the simple expression

(176) γ C † = α Γ α ĉ † α↑ ĉ † α↓
where ĉ † α↑ is the creation operator of a particle in the state |A α | ↑ and ĉ † α↓ is the creation operator of a particle in the state |B α | ↓ . Note that these creation operators and the associated annihilation operators obey the usual fermionic anticommutation relations.

The key advantage of the Schmidt decomposition is to allow a rewriting of the BCS state, more familiar to the reader:

(177) |Ψ BCS = N α 1 + Γ α ĉ † α↑ ĉ † α↓ |0 .
To obtain this expression, we have used the fact that any binary product ĉ † α↑ ĉ † α↓ commutes with any ĉ † β↑ ĉ † β↓ , and the fact that the series expansion of exp(γĉ † α↑ ĉ † α↓ ) terminates to first order in γ since (ĉ † ) 2 = 0 for fermions. The form Eq.( 177) shows that one can consider each pair of modes {α ↑, α ↓} as independent, since each factor in the product commutes with any other factor. The normalization factor is then readily calculated, N = α 1/ 1 + Γ 2 α . It can be absorbed in the following rewriting, that may be the one directly familiar to the reader:

(178) |Ψ BCS = α U α -V α ĉ † α↑ ĉ † α↓ |0 ,
where U α and V α are the real numbers (179)

U α = 1 1 + Γ 2 α and V α = -Γ α 1 + Γ 2 α .
Note that they satisfy U 2 α + V 2 α = 1, and the minus sign in V α was introduced to ensure consistency with coming notations.

3

. 1.3. As a squeezed vacuum: Wick's theorem applies. A third interesting rewriting of the BCS state can be obtained from the identity

(180) e θ(b † c † -cb) |0 = cos θ|0 + sin θ b † c † |0
where θ is a real number, b and c are two fermionic annihilation operators with standard anticommutation relations ( 2 ). Then |Ψ BCS = S|0 where the unitary operator is

(181) S = α e θα(ĉ † α↑ ĉ † α↓ -h.c.)
( 2 ) This identity can be proved by a direct expansion of the exponential in powers of θ, or by obtaining a differential equation satisfied by the left hand side considered as a function of θ.

and the angles θ α are such that (182)

U α = cos θ α V α = -sin θ α .
For bosons, S would be interpreted as a squeezing operator [START_REF] Gardiner | Quantum Noise[END_REF]. The BCS state is therefore the equivalent for fermions of the squeezed vacuum for bosons.

Calculating the expectation value in the BCS state of a product of operators ĉ and ĉ † is therefore equivalent to calculating the expectation value in the vacuum state of the product of the transformed operators S † ĉS and S † ĉ † S. These transformed operators have a linear expression in terms of the original ĉ and ĉ † , since one has

S † ĉα↑ S = U α ĉα↑ -V α ĉ † α↓ (183) S † ĉ † α↓ S = V α ĉα↑ + U α ĉ † α↓ . ( 184
)
The linearity of these transformations is evident since S can be formally considered as the time evolution operator for a quadratic Hamiltonian, which generates indeed linear equations of motion of the creation and annihilation operators. As a consequence, Wick's theorem can be applied to calculate expectation values in the BCS state, since it applies for the vacuum.

3

. 1.4. Some basic properties of the BCS state. Before moving to the energy minimization within the BCS ansatz, we calculate some interesting quantities. Since Wick's theorem applies, the expectation value of any quantity is a function of the only non-zero two-point averages:

ĉ † α↑ ĉα↑ = ĉ † α↓ ĉα↓ = Γ 2 α 1 + Γ 2 α ( 185
) ĉα↓ ĉα↑ = ĉ † α↑ ĉ † α↓ = Γ α 1 + Γ 2 α , ( 186
)
where the expectation values are taken in the BCS state.

For the total number of particles in the BCS state, we obtain the mean value and the variance

N = α 2Γ 2 α 1 + Γ 2 α (187) Var N = α 4Γ 2 α (1 + Γ 2 α ) 2 . ( 188
)
One then immediately obtains the inequality

(189) Var N ≤ 2 N ,
an inequality that becomes an equality in the limit where all Γ α ≪ 1. This shows that, in the large N limit, the fact that the BCS state has not a well defined number of particles is not a problem for most practical purposes, since the relative fluctuations are O(1/ √ N ). For the expectation value of the commutator of C and C † in the BCS state we obtain (190) [C,

C † ] = 1 -α 2Γ 4 α /(1 + Γ 2 α ) α Γ 2 α .
We see that, in the limit where all Γ 2 α ≪ 1, this expectation value is close to unity. In this limit, one may consider that C † creates a bosonic pair: the BCS state becomes a Glauber coherent state for this bosonic pair, one obtains Poissonian fluctuations in the number N /2 of pairs, which explains the upper limit of Eq.( 189). This also shows that the BCS ansatz will have no difficulty to predict the formation of a Bose-Einstein condensate of dimers in the a → 0 + limit [START_REF] Leggett | [END_REF].

A last interesting property is that the vacuum of an arbitrary quadratic Hamiltonian (quadratic in the quantum field) is actually a BCS state. We defer the proof of this statement to §3 . 3. A similar property is that the vacuum of a set of arbitrary operators b α obeying anticommutation relations is a BCS state: in other words, the BCS state is a quasi-particle vacuum. We refer to [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] for a detailed description of this aspect.

3

.

2. Energy minimization within the BCS family.

-We now proceed with the minimization of the energy of the lattice model Hamiltonian within the family of not necessarily normalized BCS states. We define

(191) E[Φ] = Ψ BCS |H|Ψ BCS Ψ BCS |Ψ BCS
where the BCS state |Ψ BCS is of the form Eq.(173) parametrized by the unnormalized pair wavefunction Φ = γφ and by the factor N now considered as an independent variable. If one performs a variation of Φ around a minimizer of E[Φ], Φ = Φ 0 + δΦ, this induces a variation δ|Ψ BCS of the BCS state around |Ψ 0 BCS , (192)

δ|Ψ BCS = r1, r2 l 6 δΦ( r 1 , r 2 ) ψ † ↑ ( r 1 ) ψ † ↓ ( r 2 )|Ψ 0 BCS ≡ δ X|Ψ 0 BCS .
The corresponding first order variation of the energy function E[Φ] has to vanish for all possible values of δΦ,

(193) δE = (δ Ψ BCS |)(H -E[Φ 0 ])|Ψ 0 BCS + c.c. = 0 ∀ δΦ,
where we assumed that the minimizer |Ψ 0 BCS is normalized to unity. The second step is to introduce the quadratic Hamiltonian deduced from the full Hamiltonian H by performing incomplete Wick's contractions. This recipe has to be applied to the interaction term only, since the other terms of the Hamiltonian are quadratic.

The quartic on-site interaction is replaced by the quadratic expression:

g 0 ψ † ↑ ψ † ↓ ψ↓ ψ↑ → ψ † ↑ ψ † ↓ g 0 ψ↓ ψ↑ + h.c. + ψ † ↑ ψ↑ g 0 ψ † ↓ ψ↓ + ↑↔↓, (194)
where the expectation values are taken in the state |Ψ 0 BCS . The first line of the righthand side consists in a pairing term, that creates/annihilates a pair of particles; it involves the following pairing field, also called 'gap' for historical reasons that will become clear:

(195) ∆( r ) ≡ g 0 ψ↓ ( r ) ψ↑ ( r ) .
The second line of the righthand side is simply the Hartree mean field term: each spin-up particle sees a mean field potential which is g 0 times the density of spin-down particle. Note that, in the lattice model, this mean field term does not diverge at the unitary limit, since g 0 remains bounded in this limit, see Eq.( 155). The quadratic Hamiltonian associated by Wick's contractions to the full Hamiltonian is the so-called BCS Hamiltonian, and it has the structure

(196) H ≡ r, r ′ ,σ l 3 ψ † σ ( r )h r, r ′ ψσ ( r ′ ) + l 3 r ψ † ↑ ( r ) ψ † ↓ ( r ) ∆( r ) + h.c. + E adj .
Here the matrix h is the sum of the one-body part of the Hamiltonian H and of the Hartree mean field terms,

(197) h r, r ′ = [kin] r, r ′ + U ( r ) -µ + g 0 ψ † ↑ ( r ) ψ↑ ( r ) δ r, r ′
where [kin] is the matrix representing the kinetic energy operator. We have assumed for simplicity that the one-body Hamiltonian is spin independent and that the two spin density profiles are identical, so that the matrix h does not depend on the spin state.

Finally, the additive constant E adj is adjusted to the value

(198) E adj = -g 0 r l 3 ψ † ↑ ψ↑ 2 + | ψ↓ ψ↑ | 2
to ensure that the mean energy of H and of H coincide in the minimizer

|Ψ 0 BCS , (199)
H = H .
Using Wick's theorem and the fact that the operator δ X is a quadratic function of the quantum field, one can then show the marvelous property

(200) (δ Ψ BCS |)H|Ψ 0 BCS = (δ Ψ BCS |)H|Ψ 0 BCS ,
whatever the variation δ|Ψ BCS of the BCS state within the BCS family.

There is the following key consequence. Considering the quadratic Hamiltonian H as a given Hamiltonian, it turns out that the BCS state |Ψ 0 BCS , obeying Eq.( 193), also obeys the condition to be a stationary state of the energy functional 193) is thus equivalent to the condition

(201) E[Φ] = Ψ BCS |H|Ψ BCS Ψ BCS |Ψ BCS , noting that one has E[Φ 0 ] = E[Φ 0 ]. Condition Eq.(
(202) δE = (δ Ψ BCS |)(H -E[Φ 0 ])|Ψ 0 BCS + c.c. = 0 ∀ δΦ,
That is, to minimize E[Φ] within the BCS family, we shall take the BCS state that minimizes the energy of H.

We are left with the study of a quadratic Hamiltonian. As reminded to the reader in §3 .

3, the quadratic Hamiltonian H may be diagonalized by a canonical transformation, the so-called Bogoliubov transformation. This allows to show that the ground state |ground of H is in fact a BCS state (see §3

.

3), and to calculate the expectation values ψ † σ ψσ ground and ψ↓ ψ↑ ground in |ground . Since |Ψ 0 BCS and |ground actually coincide, one is left with the self-consistency conditions:

ψ † σ ψσ ground = ψ † σ ψσ (203) ψ↓ ψ↑ ground = ψ↓ ψ↑ . (204)
Since the expectation values in the ground state of H are non-linear functionals of the coefficients ψ † σ ψσ and ψ↓ ψ↑ of H, this constitutes a non-linear self-consistent problem, the so-called BCS equations for the density and the gap function.

Another crucial consequence is to realize that one may write the self-consistency conditions taking an excited eigenstate of the quadratic BCS Hamiltonian H rather than the ground state [START_REF]One may wonder if condition Eq.(202) really implies that |Ψ 0 BCS is an eigenstate of H. This is not evident a priori since the variation δ|ΨBCS does not explore the whole Hilbert space but only the manifold tangent to the BCS family in |Ψ 0 BCS . One can show however that the answer is yes[END_REF]: in this case, the variational BCS calculation is used to predict excited states of the gas! This shows that BCS theory immediately predicts elementary excitations of the gas. As we shall see in the thermodynamic limit, these excitations correspond to pair breaking and are characterized by an energy spectrum with a gap.

3

.

3. Reminder on diagonalization of quadratic Hamiltonians. -We consider here a quadratic Hamiltonian of the form Eq. (196), where, at this stage, h is an arbitrary (spin independent) hermitian matrix and ∆( r ) an arbitrary complex field. We wish to put this Hamiltonian in canonical form by performing a Bogoliubov transformation.

The key property of a quadratic Hamiltonian is that it leads in Heisenberg picture to linear equations of motion. Collecting the field operators and their hermitian conjugate at all points of the lattice in big vectors, we get

(205) ih∂ t ψσ ψ † -σ = L σ ψσ ψ † -σ
where -σ stands for the spin component opposite to σ and where we have introduced two matrices

(206) L ↑ = h diag(∆) diag(∆ * ) -h * and L ↓ = h diag(-∆) diag(-∆ * ) -h * .
These two matrices are hermitian: they can be diagonalized in an orthonormal basis. Furthermore they obey the following symmetry property: if the vector (u, v) is an eigenstate of L ↑ with the energy ǫ, then

• the vector (-v * , u * ) is an eigenvector of L ↑ , with the eigenvalue -ǫ,

• the vector (u, -v) is an eigenstate of L ↓ , with the eigenvalue ǫ,

• the vector (v * , u * ) is an eigenvector of L ↓ , with the eigenvalue -ǫ.

Assuming for simplicity that all the eigenenergies are non-zero, we can collect the eigenvectors in pairs of opposite eigenenergies. Expanding ( ψσ , ψ † -σ ) in the eigenbasis of L σ , and expressing the negative energy modes of L ↑ and all the eigenmodes of L ↓ in terms of the positive energy modes of L ↑ , noted as (u α , v α ) , we obtain

ψ↑ ( r ) = α b α↑ u α ( r ) -b † α↓ v * α ( r ) (207) ψ↓ ( r ) = α b α↓ u α ( r ) + b † α↑ v * α ( r ) (208)
where the positive energy eigenvectors of L ↑ are normalized in a way mimicking the continuous space limit, with Dirac's notation r |u α = u α ( r ):

(209) u α |u β + v α |v β = l 3 r u * α ( r )u β ( r ) + v * α ( r )v β ( r ) = δ αβ ,
The coefficients in the expansion of the field operators are themselves operators that are easy to express using the orthonormal nature of the eigenbasis of L σ :

b α↑ = l 3 r u * α ( r ) ψ↑ ( r ) + v * α ( r ) ψ † ↓ ( r ) (210) b α↓ = l 3 r u * α ( r ) ψ↓ ( r ) -v * α ( r ) ψ † ↑ ( r ). ( 211
)
Using Eq.(209) again, one can check that the b's and their hermitian conjugates obey the standard fermionic anticommutation relations, The last step is to express the Hamiltonian in normal form, using the b α 's. One first realizes that the quadratic Hamiltonian can be expressed as a quadratic form defined by the matrices L σ , (212)

H = Tr h + E adj + 1 2 l 3 σ ψ † σ , ψ-σ L σ ψσ ψ † -σ
where the extra constant term, given by the trace of the matrix h, originates from the anticommutator of ψ † σ with ψσ . Inserting the modal decomposition Eqs.(207,208) and using the fact that (u α , v α ) is an eigenvector of L ↑ , etc, we obtain the canonical form:

(213) H = E 0 + α,σ ǫ α b † ασ b ασ ,
where we recall that all the ǫ α > 0 by construction. The immediate expression for E 0 is E 0 = Tr h + E adjα ǫ α ; it involves two quantities that are not well behaved in the zero lattice spacing limit. Expressing the matrix L ↑ as a sum of ±ǫ α times dyadics involving the eigenvectors as ket and bra, and projecting over the upper left block leads to

(214) h r, r ′ = l 3 α ǫ α [u α ( r )u * α ( r ′ ) -v * α ( r )v α ( r ′ )] ,
where we used the orthonormal nature of the eigenbasis of L ↑ , as defined in Eq.( 209).

Taking the trace of this expression in the spatial basis of the lattice leads to the more convenient form for E 0 :

(215)

E 0 = E adj -2 α ǫ α v α |v α ,
where E adj is given by Eq.(198). The canonical form Eq.( 213), with all the ǫ α > 0 by construction, immediately shows that the ground state |Ψ 0 of H is the vacuum of all the b α 's, with an eigenenergy E 0 . Please remember that we are dealing here with the grand canonical Hamiltonian; this is why E 0 < 0 for the ideal gas g 0 = 0. To show that this ground state is a BCS state, we take for |Ψ 0 an ansatz of the BCS type Eq.( 173

[b ασ , e γC † ] = γe γC † [b ασ , C † ]. (217)
Since exp(γC † ) is an invertible operator, Eq.(216) reduces to

(218) b ασ + γ[b ασ , C † ] |0 = 0.
From the expressions Eqs.(210,211) we calculate the commutator with C † and obtain the linear system

v * α ( r ) + γ r ′ l 3 φ( r ′ , r )u * α ( r ′ ) = 0 (219) v * α ( r ) + γ r ′ l 3 φ( r, r ′ )u * α ( r ′ ) = 0. ( 220
)
The solution may be written in matrix form as follows. Since the lattice has discrete positions, we consider φ( r, r ′ ) as the element of a square matrix [φ] on row r and column r ′ . Similarly we form the square matrix [U ] (respectively [V ]) such that the element on row r and column α is u α ( r ) (respectively v α ( r )). Then one has the explicit writing

(221) γ [φ] = γ t [φ] = -l -3 ([U ] † ) -1 [V ] † ,
where t M is the transposed matrix of a matrix M . Since the set of all the eigenvectors of L ↑ forms an orthonormal basis, one has l 3 (U † U + V † V ) = Id and t U V = t V U , where Id is the identity matrix. This ensures that the solution [φ] is indeed a symmetric matrix

(222) φ( r 1 , r 2 ) = φ( r 2 , r 1 ).
This symmetry property could be expected from the fact that the quadratic Hamiltonian H is invariant by an arbitrary rotation R of the spin degree of freedom, and that its ground state is not degenerate since we have assumed ǫ α strictly positive. The condition that C † = RC † R † for all spin rotations indeed imposes that φ is a symmetric function of r 1 and r 2 . Using Wick's theorem one can easily calculate expectation values of observables in |Ψ 0 using the modal decompositions Eqs. (207,208). Two important examples are the mean density and the anomalous average

ρ ↑ ( r ) = ρ ↓ ( r ) = α |v α ( r )| 2 , (223) ψ↓ ( r ) ψ↑ ( r ) = - α u α ( r )v * α ( r ). (224)

3

. 4. Summary of BCS results for the homogeneous system. -

3 .
4.1. Gap equation in the thermodynamical limit. In the case of a spatially homogeneous system (cubic box with periodic boundary conditions) explicit predictions are easily extracted from the BCS theory. For a spatially uniform gap parameter ∆, assumed to be real non negative, and for a uniform density profile ρ ↑ = ρ ↓ the spectral decomposition of L ↑ is readily obtained: from the translational invariance one seeks eigenvectors of the form

(225) (u k ( r ), v k ( r )) = L -3/2 e i k• r (U k , V k ).
Restricting to the positive part of the spectrum, one finds the eigenenergies

(226) ǫ k = h2 k 2 2m -μ 2 + ∆ 2 1/2 = h2 k 2 2m -μ + i∆
where μ = µg 0 ρ ↑ is the chemical potential minus the Hartree mean field term. The amplitudes U k , V k , chosen to be real, are normalized as U 2 k + V 2 = 1 and are given by

(227) (U k + iV k ) 2 = h2 k 2 2m -μ + i∆ ǫ k .
To make the link with the quantities Γ α of section 3 . 1.2, we take for |A α = |A k the plane wave of wavevector k and for |B α = |B k the plane wave of wavevectork. This gives

(228) Γ k = - V k U k = - ∆ h2 k 2 2m -μ + ǫ k
, and a pair wavefunction [START_REF]This expression coincides with Eq.(221[END_REF] (229)

γφ( r 1 , r 2 ) = 1 L 3 k∈D Γ k e i k•( r1-r2) .
Remember that the regime where all the Γ 2 k ≪ 1 corresponds to the case where the BCS state is a condensate of bosonic pairs.

Taking the real and imaginary parts of Eq.( 227), and going directly to the thermodynamic limit, one obtains from Eqs.(223,224)

ρ = D d 3 k (2π) 3 1 - h2 k 2 2m - μ ǫ k (230) ∆ = -g 0 D d 3 k (2π) 3 ∆ 2ǫ k , ( 231
)
where ρ is the total density. This makes explicit the self-consistent character of the conditions Eqs. (203,204). The last equation Eq.( 231) is called the gap equation.

3

. 4.2. In the limit of a vanishing lattice spacing. An important point is to discuss the dependence of these predictions with the lattice spacing l, since the relevant regime of a continuous gas is described by the lattice model only in the l → 0 limit. The integral giving ρ converges at large k, since the integrand is O(1/k 4 ), so that one may directly set l → 0 in Eq.( 230). On the contrary the integral in Eq.( 231) has an ultraviolet divergence when integrated over the whole momentum space. However g 0 tends to zero in the l → 0 limit, compensating exactly this divergence. Dividing Eq.( 231) by g 0 and expressing 1/g 0 from Eq.( 154) leads to a more convenient expression for the gap equation, in which one may take directly l → 0:

(232) ∆ g = ∆ D d 3 k (2π) 3 m h2 k 2 - 1 2ǫ k .
Furthermore, we note that the BCS theory becomes simpler when l → 0: since the coupling constant g 0 vanishes in this limit, the Hartree mean field term also vanishes, so that μ → µ and the matrix h in the quadratic Hamiltonian H reduces to the matrix of the kinetic energy minus µ. The a priori unknown coefficients of H thus reduce to the anomalous average ψ↓ ψ↑ , and Eq.( 204) remains as the only self-consistent equation, the other one Eq.( 203) giving explicitly the density as a function of µ and of the anomalous average. In general, one has in practice to use the BCS theory in this l → 0 limit, to be in the continuous gas regime. As we shall see in section 3 . 4.4, a notable exception is the weakly attractive regime k F a → 0 -.

3

. 4.3. BCS prediction for an energy gap. As sketched in the last paragraph of section 3 . 2, the excited states of the BCS Hamiltonian H correspond to a variational BCS approximation for excited states of the full Hamiltonian of the gas. Here we come back to this statement, and show how to derive one of the most celebrated predictions of BCS theory, the minimum energy required to break a pair.

First we revisit briefly section 3 . 2: in the process of energy minimization within the BCS family of states, we ended up looking for a BCS state that is the ground state of a quadratic Hamiltonian of the form Eq. (196). Taking here for simplicity the continuous limit l → 0, for a spatially homogeneous solution, we can omit the Hartree mean field term and perform a variational calculation that immediately restricts to BCS states being the ground state of the following Hamiltonian, (233)

H λ = k∈D σ (E k -µ) a † kσ a kσ + l 3 r λ ψ † ↑ ( r ) ψ † ↓ ( r ) + h.c.
where

E k = h2 k 2 /
H g λ = k∈D (E k -µ) 1 - E k -µ |E k -µ + iλ| + g 0 L 3 ψ↓ ψ↑ g λ 2 ,
with the anomalous average ψ↓ ψ↑

g λ = -L -3 k∈D λ/(2|E k -µ + iλ|).
In these expressions, the appex g stands for "ground state". It remains to require that the first order derivative of H g λ with respect to λ vanishes, which reduces, after explicit calculation of the derivative, to the following equation satisfied by the optimal value λ g of λ:

(235) g 0 ψ↓ ψ↑ g λg = λ g .

One has recovered, in a slightly different way, the self-consistency condition Eq.( 204), and one finds that λ g satisfies the gap equation (with μ = µ). This optimal λ g was called earlier ∆.

We can now perform a similar variational calculation, taking as a variational ansatz an excited state of H λ . To preserve the symmetry between the two spin components ↑ and ↓, we consider the excited state ( 236)

|Ψ e λ = b † q ↑ b † -q↓ |Ψ g λ ,
where q is a fixed wavevector. Two Bogoliubov modes of the Hamiltonian H λ are now in the excited state, with occupation number one, rather than being in the vacuum state of the corresponding b's. The mean energy of H in this excited state is simply, in the thermodynamic limit, (237)

H e λ = H g λ + 2(E q -µ) 2 + 2λg 0 ψ↓ ψ↑ g λ |E q -µ + iλ| + O(L -3 ).
It differs from H g λ by a term O(1), whereas the full energy is O(L 3 ), so that the optimal value λ e of λ for the minimization of H e λ differs from λ g = ∆ by a small term δλ. We can then simply expand Eq.(237) in powers of δλ. The key point now is that the ground state energy varies only to second order in δλ since it is minimal in ∆; since the second order derivative of

H g λ is O(L 3 ), δλ is O(1/L 3
) and the contribution to the total energy of the corresponding second order variation of H g λ is O(1/L 3 ) and negligible. In the remaining terms of Eq.(237) we can simply set λ = ∆. Using Eq.( 235), we are left with ( 238)

H e λe = H g ∆ + 2ǫ q + O(1/L 3 )
where ǫ q is the BCS spectrum Eq.(226).

To summarize, one considers excited states of H as variational states, with a small number of Bogoliubov excitations. One should then in principle solve again the selfconsistency conditions Eqs.(203,204), but in the thermodynamic limit the excitation of a few Bogoliubov modes has a small effect on the density and gap parameter, an effect that we have shown to be negligible on the energy. One can then directly consider that the elementary excitations of the quadratic Hamiltonian H associated to the ground energy BCS state actually give the energy of the elementary excitations of BCS theory [START_REF]) + µeN . Here the chemical potentials µ and µe differ by O(1/L 3 ) for a few Bogoliubov excitations, so that one may expand to first order in δµ = µe -µ[END_REF].

When ∆ > 0, we see that the minimal value of the BCS spectrum ǫ k with respect to k is non-zero: it has a gap E gap = min k ǫ k . When μ > 0, which contains the usual regime of the BCS theory, the regime of condensation of Cooper pairs (see below), the gap is E gap = ∆, hence the name of ∆. When μ < 0, which contains the regime of a Bose-Einstein condensate of dimers, the gap is E gap = μ2 + ∆ 2 1/2 . Finally, we get from BCS theory a prediction for the miminal energy 2E gap required to break a pair, that is to get from a condensate of N/2 pairs, a condensate of N/2 -1 pairs and two unpaired atoms [START_REF] De Gennes | Superconductivity of Metals and Alloys, section 4.3[END_REF].

3

. 4.4. BCS predictions in limiting cases. We shall not discuss the full solution of the BCS equations but we briefly review simple limiting cases. We introduce the Fermi momentum k F of the ideal Fermi gas, such that k 3 F /(6π 2 ) = ρ ↑ = ρ ↓ . Limit k F a → 0 -: the gap parameter ∆ tends exponentially to zero in this limit, as we shall see. For such a small value of ∆, one may set ∆ = 0 in Eq.( 230), to obtain that μ is the Fermi energy of the ideal gas h2 k 2 F /2m, so that

(239) µ ≃ h2 2m (6π 2 ρ ↑ ) 2/3 + g 0 ρ ↑ .
If one takes the mathematical limit l → 0 the Hartree mean field term disappears since g 0 → 0 in this limit, as already mentioned. However in the present weakly attractive limit one can choose |a| ≪ l ≪ k -1 F so that the continuous space limit and the zero-range condition Eq.( 116) for the lattice model are obtained while g 0 ≃ g, see Eq.( 155). The on-site interaction potential is then treatable in the Born approximation, which makes the BCS approach more accurate. One then indeed recovers the first term gρ ↑ of a systematic expansion of µ in powers of k F a [START_REF] Heiselberg | [END_REF]. Now turning to the gap equation, one finds a gap parameter [START_REF] Randeria | Bose-Einstein Condensation[END_REF] (240) ∆ ≃ 8e -2 μ e -π/(2kF |a|) .

Since μ = h2 k 2 F /2m > 0 in this weakly attractive limit, |Γ k | can assume extremely large values, for h2 k 2 /2m < μ: the pairs that are condensed are not bosons. Note the nonanalytic dependence of the gap on the small parameter k F |a|, which indicates that the BCS state in the thermodynamic limit cannot be obtained by a perturbative treatment of the interaction potential. This non-analytic dependence can be readily seen from Eq.(232), whose integrand diverges in k = √ 2mμ/h, in the limit ∆ → 0. Replacing the k integration variable by the ideal gas mode energy E = h2 k 2 /2m, and approximating the ideal gas density of modes by a constant in the energy interval of half-width δ ≪ μ around μ, one gets a contribution (241

) μ+δ μ-δ dE ∆ 2 + (E -μ) 2 ∆→0 ∼ 2 ln 2δ ∆ .
The same technique applied over e.g. the interval E ∈ [0, 4μ] leads to Eq.( 240), when one also takes into account (to lowest order in ∆) the difference between the exact integral in the gap equation and the approximate one.

Limit k F a → 0 + : the coefficients Γ k tend uniformly to zero in this limit, because μ < 0 and ∆ ≪ |μ|, so that the pair creation operator C † obeys approximately bosonic commutation relations. This was expected physically, since the ground state of the gas is a condensate of almost bosonic dimers. Since two atoms in a dimer are at a relative distance ∼ a, one should ensure l ≪ a in the lattice model, so that one takes the mathematical limit l → 0, which implies μ → µ. Let us check that the BCS theory correctly predicts a condensate of such dimers. We first simplify the gap equation Eq.( 232) by using ǫ k ≃ -µ + h2 k 2 /2m. After division by ∆, we get an equation for the chemical potential

(242) 1 g ≃ - d 3 k (2π) 3 1 h2 k 2 m -2µ - m h2 k 2 which leads to (243) µ ≃ - h2 2ma 2 .
This is minus half the binding energy of a dimer, exactly what was expected (keep in mind that µN = µ mol N mol where N mol is the total number of dimers and is equal to N/2, so that the molecular chemical potential µ mol is twice the atomic one). The next step is to expand the integrand of Eq.(230) to leading order in ∆ to calculate the gap parameter [START_REF] Randeria | Bose-Einstein Condensation[END_REF]:

(244) ∆ ≃ 2 √ 3π (k F a) 3/2 h2 ma 2 ≪ h2 ma 2 .
Note that in this molecular BEC regime, ∆ is not proportional to the energy required to break a pair. One sees from Eq.( 226) that the gap in ǫ k is ≃ |µ|, since µ is negative and much larger in absolute value than ∆. The energy to break a pair is then 2|µ|, which is indeed the binding energy of a dimer. Finally, by performing the Fourier transform in Eq.( 229), approximating V k /U k to leading order in ∆, one obtains the pair wavefunction

(245) φ( r 1 -r 2 ) ≃ 1 L 3/2 φ 0 (| r 1 -r 2 |)
where φ 0 is the wavefunction of the bound state of two atoms given by Eq.( 146), and where we used

γ 2 = k (V k /U k ) 2 .
Limit k F |a| = +∞: the numerical solution of the gap equation Eq.( 232) and of the density equation Eq.(230) in the limit l → 0 gives the BCS estimate of the numerical coefficient η of Eq.(168) for the ground branch. This estimate is an upper bound [START_REF] Engelbrecht | [END_REF]:

(246) η ≤ η BCS = 0.5906 . . .

A much better estimate was obtained by approximate (fixed node) Monte Carlo calculations [35,36], η ≃ 0.4, a value confirmed by recent exact quantum Monte Carlo calculations [33]. Early measurements of η were in contradiction with these theoretical values [37], but more precise measurement performed in Innsbruck [38], in Paris [39] and at Rice University [43] are consistent with them.

3

. 5. Derivation of superfluid hydrodynamic equations from BCS theory. -A key point of the BCS theory is to predict a spectrum of elementary excitations having a gap, see Eq.( 226). One could then be tempted to infer strong predictions on the thermodynamic properties, e.g. that the entropy obeys an activation law O(e -Egap/kB T ) at low temperature, where E gap is the energy gap. However, in the BEC limit, where the gas is simply a Bose-Einstein condensate of bosonic dimers at low temperature, this prediction of an activation law is obviously wrong, since one knows that the relevant excitation spectrum is the Bogoliubov spectrum, which has no gap.

In reality, the spatially homogeneous gas, whatever the considered regime (BEC or BCS), is expected to have a branch of collective excitations with no gap, behaving as sound waves at low momenta. These collective excitations correspond to coherent oscillations of the density of the pairs, which are not gapped, to be distinguished from the pair breaking excitations which are gapped. The key point briefly addressed in this section is that the time-dependent BCS theory contains such a branch of collective, nonpair-breaking oscillations [START_REF]We consider here the case of neutral particles. For charged particles the situation is different, because[END_REF]. This is similar to what happens for the weakly interacting Bose gas: the Bogoliubov spectrum is obtained from a linearization of the time dependent Gross-Pitaevskii equation.

Here, rather than performing an exact linearization of the time-dependent BCS equations, which leads to the so-called RPA approach [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF][START_REF] Côté | [END_REF]71,72,73,74], we go through a sequence of simple approximations allowing one to derive superfluid hydrodynamic equations from the time-dependent BCS theory. This has several advantages: it is more physical, it easily applies to harmonic trapping [START_REF]The RPA approach can be applied in an isotropic harmonic trap, with some non negligible numerical effort[END_REF], and it is applicable to the nonlinear time dependent regime, to study analytically low-energy collective excitations of a trapped superfluid Fermi gas [START_REF] Baranov | [END_REF]77], but also the expansion of the gas after the trap was switched off [78] and the response of the gas to a rotation of the harmonic trap [79,80] even in the non-linear regime [80].

Coming back to thermodynamical aspects, one may fear that the straightforward finite temperature extension of the BCS theory, with a variational density operator which is Gaussian in the field operators [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF], is not able to calculate the critical temperature T c with a good accuracy, because it does not correctly include the collective excitations. It turns out that this simple finite temperature BCS theory correctly gives T c in the BCS limit k F a → 0 -within a numerical factor [81], but is indeed totally wrong in the BEC limit k F a → 0 + where it does not reproduce at all Einstein's prediction for the critical temperature of the ideal Bose gas. Refinements of the BCS theory have been developed to recover Einstein's prediction and to obtain a calculation of T c approximately valid over the whole range of values of k F a [44,[START_REF] Randeria | Bose-Einstein Condensation[END_REF].

3

. In practice, the equation of motion for Φ is not useful. One rather moves to the Heisenberg picture with respect to the time dependent Hamiltonian H(t), as we did in section 3 . 3. Since this Hamiltonian is quadratic, the equations of motion for the fields are linear, of the form Eq.( 205), where L ↑ is now time dependent. Since these equations are linear, we can solve them by evolving in Eqs.(207,208) the mode functions (u α , v α ) while keeping constant the quantum coefficients bασ where σ =↑ or ↓:

(252) ih∂ t u α v α = L ↑ (t) u α v α .
These equations on the mode functions are effectively non-linear since coefficients of L ↑ , involving the mean density ψ † ↑ ψ↑ and the anomalous average ψ↓ ψ↑ , depend on the mode functions; since the bασ are constants of motion, this dependence is still given by Eqs.(223,224), now taking the time dependent u α 's and v α 's.

3

. 5.2. Semi-classical approximation. The physical situation that we have in mind here is a gas in the BCS regime (a < 0, µ > 0) in a time dependent harmonic trap. The trap may be rotating, with an angular velocity Ω(t), in which case one moves to the rotating frame to eliminate the time dependence of the trapping potential due to rotation; this introduces an additional term -Ω(t) • L in the one-body Hamiltonian, where L is the angular momentum operator of a particle. In this case one may expect that quantum vortices form [82] for a high enough rotation frequency; the semi-classical approximation that we present is however restricted to a vortex-free gas.

We wish to treat the equation of motion of (u α , v α ) in the semi-classical approximation. The general validity condition of a semi-classical approximation is that the applied potentials vary slowly over the coherence lengths of the gas, a coherence length being the typical width of a correlation function of the gas.

A first correlation function here is ψ † ↑ ( r )ψ ↑ ( r ′ ) . Using the homogeneous gas expression of this function, one sees that it is the Fourier transform of the function V 2 k of width the Fermi momentum, ∆k = k F , where h2 k 2 F /2m = µ. The associated coherence length is 1/∆k = k -1 F . A second correlation function is ψ↓ ( r ) ψ↑ ( r ′ ) . For the homogeneous gas this is the Fourier transform of the function U k V k , which according to Eq.( 227) has a momentum width ∆k = m|∆|/h 2 k F . The associated coherence length 1/∆k corresponds in BCS theory to the length of a Cooper pair, (253) l BCS = h2 k F m|∆| .

Since ∆ < µ in the BCS regime, we keep l BCS as the relevant coherence length. A first typical length scale of variation of the matrix elements in Eq.(252) originates from the position dependence of |∆|: in the absence of rotation and for an isotropic trap, this is expected to be the Thomas-Fermi radius R TF of the gas, defined as µ = mω 2 R 2 TF /2, where ω is the atomic oscillation frequency. This assumes that the scale of variation of the modulus of the gap is the same as the one of the density, a point confirmed in the adiabatic approximation to come. The condition that the mean field Hartree term and the harmonic potential have a weak relative variation over l BCS for typical values of the position also leads to the condition (254) l BCS ≪ R TF .

For an isotropic harmonic trap this is equivalent to the condition (255) |∆| ≫ hω.

In the general time dependent case, however, this is not the whole story, since the phase of ∆ can also become position dependent. In the rotating case for example, with a rotation frequency of the order of ω, the phase of ∆ may vary as ≃ mωxy/h; this introduces a local wavevector mωR TF /h ≃ k F , making a semi-classical approximation hopeless. We eliminate this problem with a gauge transform of the u's and v's [80]: ũα ( r, t) ≡ u α ( r, t)e -iS( r,t)/h (256) ṽα ( r, t) ≡ v α ( r, t)e +iS( r,t)/h , (257) where S is defined as An efficient frame to perform a semi-classical approximation in a systematic way for the evolution of a wave in a potential is to use the Wigner representation [6] of the density operator of the wave. Here the wave is represented by the two-component spinor (ũ α , ṽα ) so that we introduce the corresponding density operator We introduce the Wigner representation of σ [6], assuming for simplicity a continuous position space:

(263) W ( r, p, t) = Wigner{σ} ≡ d 3 x (2πh) 3 rx/2|σ| r + x/2 e i p• x/h . Since σ is the density operator of a two-component spinor, the Wigner distribution W is a two by two matrix. Within this representation, the key quantities of BCS theory, the total density, the modulus of the gap function and the total matter current have the following expressions in the general case of a rotating frame: ρ( r, t) = 2 where Ω is the angular velocity of the rotating frame, and the so-called velocity field is defined as (267) v( r, t) ≡ grad S( r, t) m .

We have introduced here the total matter current ( r, t), that obeys by construction (268) ∂ t ρ + div  = 0,

In the rotating frame, in a many-body state invariant by exchange of the spin states ↑ and ↓, it is very generally given by (269) ( r, t) = h im ψ † ↑ ( r, t) grad ψ↑ ( r, t)c.c.ρ( r, t) Ω(t) × r, Within BCS theory, these last two relations still hold [START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF] and lead to Eq.( 266). The semi-classical expansion then consists e.g. in

(270) Wigner{V ( ˆ r )σ} = [V ( r ) + ih 2 ∂ r V • ∂ p + . . .]W ( r, p, t),
where V is a position-dependent potential. The successive terms in this expansion we call zeroth order, first order, etc, in the semi-classical approximation. Assuming that the momentum derivative of W is ≃ W/∆p, where ∆p is the momentum width of W , we recover the fact that the first order term is small as compared to the zeroth order one if the variation of V over the coherence length h/∆p, which is of the order of the spatial derivative of V times the coherence length, is small as compared to V .

Here we shall need only the equations of motion of the Wigner distribution W up to zeroth order in the semi-classical approximation:

(271) ih∂ t W ( r, p, t) ≃ L ↑ 0 ( r, p, t), W ( r, p, t)

where the two by two matrix L ↑ 0 is equal to that may be called effective chemical potential for reasons that will become clear later. At time t = 0, the gas is at zero temperature. By introducing the spectral decomposition of L ↑ (t = 0) one can then check that where (U 0 , V 0 ) is the eigenvector of L ↑ 0 ( r, p, t = 0) of positive energy and normalized to unity. At time t, according to the zeroth order evolution Eq.( 271), each two by two matrix W remains a pure state, with components U and V solving (277) ih∂ t U ( r, p, t) V ( r, p, t) = L ↑ 0 ( r, p, t) U ( r, p, t) V ( r, p, t) .

(272) L ↑ 0 (

3

. 5.3. Adiabatic approximation. We then introduce the so-called adiabatic approximation: the vector (U, V ), being initially an eigenstate of L ↑ 0 ( r, p, t = 0), will be an instantaneous eigenvector of L ↑ 0 ( r, p, t) at all later times t. This approximation holds under the validity condition of the adiabaticity theorem [START_REF] Messiah | Quantum Mechanics[END_REF], discussed for our specific case in [80], generically requiring that the energy difference between the two eigenvalues of L ↑ 0 ( r, p, t) be large enough. As this energy difference can be as small as the gap parameter, this imposes a minimal value to the gap, not necessarily coinciding with the one of Eq.( 255), as shown in [80]. In this adiabatic approximation, one can take where | + ( r, p, t) , having real components (U inst , V inst ), is the instantaneous eigenvector with positive eigenvalue of the matrix L ↑ 0 defined in Eq.( 272). (U inst , V inst ) are simply the amplitudes on the plane wave exp(i p • r/h) of the BCS mode functions of a spatially uniform BCS gas of chemical potential µ eff and of gap parameter |∆( r, t)|. Using Eq.( 264) and Eq.( 265) with the approximate Wigner distribution Eq.( 278), one further finds that this fictitious spatially uniform BCS gas is at equilibrium at zero temperature so that the zero temperature equations of state may be used, relating the chemical potential to the density, µ 0 [ρ], and the gap to the density, ∆ 0 [ρ], where the functions µ 0 and ∆ 0 may be calculated from the spatially homogeneous BCS theory. In particular, the equation of state relating the chemical potential to the density leads to (279) µ eff ( r, t) = µ 0 [ρ( r, t)] which leads, together with Eq.( 273), to one of the time dependent hydrodynamic equations, the Euler-type one. Also, U inst and V inst are even functions of p, so that the integral in the right hand side of Eq.(266) vanishes and Eq.(268) reduces to the hydrodynamic continuity equation in a rotating frame, (280) ∂ t ρ( r, t) + div ρ( r, t) v ( r, t) -Ω(t) × r = 0.

Under the adiabatic approximation, the superfluid hydrodynamic equations are thus derived from BCS theory, remarkably without having to postulate the expression of the superfluid current in terms of the gradient of the phase of the order parameter. We refer to [80] for a discussion of the validity condition of the adiabatic approximation. The presence of superfluid hydrodynamic equations allows to conclude in the existence of collective modes in the BCS theory, in a transparent way, by a simple linearization around steady state. Although we have considered here a trapped system, we note that the formalism may also be developed in the spatially homogeneous case. To create a sound wave of wavevector q, one may apply an external potential varying with spatial harmonics e ±i q• r , e.g. with a Bragg pulse [START_REF] Meacher | [END_REF]85]. The semi-classical approximation holds if this external potential varies slowly at the scale of the Cooper pair size l BCS , that is ql BCS ≪ 1. Linearization of hydrodynamic equations in the linear response regime predicts a linear dispersion relation ω q = c s q, with a sound velocity c s such that (281)

mc 2 s = ρ dµ 0 dρ [ρ],
of the order of hk F /m in the regime a < 0. The semi-classicality condition ql BCS ≪ 1 then results in hω q ≪ ∆. This is satisfactory, as for hω q > 2∆ one observes in the full RPA theory a coupling of the sound wave to the elementary, pair breaking excitations, which may distort the dispersion relation and even damp the sound wave, see e.g. [72,73,74] and references therein, a phenomenon not included in superfluid hydrodynamics. We acknowledge useful suggestions from Félix Werner, Mattia Jona-Lasinio, Christoph Weiss, Alice Sinatra.

c Società Italiana di Fisica 1 a

 1 zero external potential, supplemented by periodic boundary conditions in a box of size L, defined by r ∈ [0, L[ d , or a harmonic potential mimicking the traps used in experiments:

Fig. 1 .

 1 Fig. 1. -Fermi-Dirac distribution in the degenerate regime. Solid line: occupation numbers n(ǫ) for the particles. Dashed line: occupation numbers 1 -n(ǫ) for the holes. The temperature is kBT = µ/10.

B 3 =

 3 exp(γ + 2 ln 2 -5/12) = 4.696621 . . . (53) and 2 F 3 is a hypergeometric function, γ = 0.5772156649 . . . is Euler's constant.

1

 1

  Fig. 3. -Central density of a 1D harmonically trapped ideal Fermi gas, in the grand canonical ensemble, as a function of the chemical potential µ. (a) Density in the trap center x = 0: exact result for T = 0 (black solid line with steps), exact numerical sum for kBT = hω (red smooth solid line), vs the T = 0 LDA (dashed line). (b) Second order derivative of the density in the trap center, d 2 ρ/dx 2 (x = 0), for kBT = hω: exact numerical sum (solid line) vs the T = 0 LDA (dashed line).

Fig. 4 .

 4 Fig. 4. -Oversimplified view of a magnetically induced Feshbach resonance. The interaction between two atoms is described by two Born-Oppenheimer curves. Solid line: open channel potential curve. Dashed line: closed channel potential curve. When one neglects the coupling Λ between the two curves, the closed channel has a molecular state of energy E mol with respect to the dissociation limit of the open channel. Note that the energy spacing between the solid curve and the dashed curve was greatly exaggerated, for clarity.

Fig. 5 .

 5 Fig.5. -In the toy model with a scatterer in a hard wall spherical cavity: (a) energy per particle and (b) pressure of the gas in the first two branches, as functions of -1/kF a. kF , EF and PF are respectively the Fermi wavevector, the Fermi energy and the pressure of the T = 0 ideal Fermi gas with the same density as the interacting gas.

3 . 1 . 3 . 1 . 1 .

 31311 The BCS ansatz . -A coherent state of pairs. Let us recall the Glauber coherent state of a bosonic field,(171) |α = e -|α| 2 /2 e αa † |0 ,

  ) to solve the equations (216) b ασ |Ψ 0 = 0 for all mode index α and spin component σ =↑, ↓. The commutator of a b ασ with the pair creation operator C † is a linear combination of the ψ † 's so it commutes with C † . As a consequence,

  , t) = |∆( r, t)| e 2iS( r,t)/h . The time dependent BCS equations Eq.(252) are modified, ∆ being replaced by |∆| and the single particle Hamiltonian h being replaced by the gauge transformed Hamiltonian (259)h = e -iS/h he +iS/h + ∂ t S.

α

  |ũ α ũα | |ũ α ṽα | |ṽ α ũα | |ṽ α ṽα | .Note that the matrix elements of σ in position space are related to the previously mentioned correlation functions of the gas, up to the gauge transform, which makes the discussion consistent:| r |σ 22 | r ′ | = | ψ † ↑ ( r ) ψ↑ ( r ′ ) | (261) | r |σ 12 | r ′ | = | ψ↓ ( r ) ψ↑ ( r ′ ) |. (262)

hD d 3

 3 p W 22 ( r, p, t) (264) |∆|( r, t) = -g 0 hD d 3 p W 12 ( r, p, t) (265)  ( r, t) = ρ v ( r, t) -Ω(t) × r -2 m hD d 3 p p W 22 ( r, p, t),(266)

  = 0) = θ[L ↑ (t = 0)]where θ(x) is the Heaviside function. Since L ↑ 0 (t = 0) is the classical limit of the operator L ↑ (t = 0), the leading order semi-classical approximation for the corresponding Wigner function is, in a standard way, given by (275)W ( r, p, t = 0) ≃ 1 (2πh) 3 θ[L ↑ 0 ( r, p, t = 0)].This implies that the 2 × 2 matrix W is proportional to a pure spin-1/2 state |ψ ψ| with (276) |ψ( r, p, t = 0) = U 0 ( r, p ) V 0 ( r, p )

  3 θ[L ↑ 0 ( r, p, t)] = 1 (2πh) 3 | + ( r, p, t) +( r, p, t)|

  3/4 a where a ho = h/mω. This trial wavefunction does not obey the contact conditions, so it is not in the domain of the Hamiltonian. It is therefore mathematically incorrect to calculate the mean energy of φ t by representing the action of H on φ t by the differential operator Eq.(147). One rather has to calculate the overlap ψ n |φ t of |φ t with each exact eigenstate ψ n of H of eigenenergy E n , and calculate the sum n

	3/2	,
	ho

  This formula is reminiscent of the technique of renormalization of the coupling constant.A natural case to consider is the one of the usual parabolic dispersion relation, E k = h2 k 2 /2m. A more explicit form of Eq.(154) is then and that may be expressed analytically in terms of the dilog special function. The effective range of the model is easily calculated with the complexification technique k 0 = iq 0 ; it is positive and proportional to l:

	with a numerical constant given by
	(156)	K =	12 π	0	π/4	dθ ln(1 + 1/ cos 2 θ) = 2.442 749 607 806 335 . . . ,
	(155)					g 0 =	4πh 2 a/m 1 -Ka/l

  2m and the real non-negative number λ is now the only variational parameter. The expectation value of the full Hamiltonian H in the ground state |Ψ g λ of

	H λ is readily evaluated from section 3	. 3 and from Eqs.(226,227),
	(234)

  5.1. Time dependent BCS theory. This is a direct generalization of the static case of section 3 .2. The exact N -body Schrödinger equation can be obtained from extremalisation over the time-dependent state vector |Ψ(t) of the following action, We have integrated by parts over t to get rid of (d/dt)|δΨ and we have used Eq.(248) to show that the fully integrated term vanishes. The condition that δS vanishes for all ket variations indeed leads to Schrödinger's equation.For time dependent BCS theory, one forces the ket to be of the BCS form Eq.(173), |Ψ = |Ψ BCS . The variation of the ket is performed within the BCS family, by a variation of Φ = γφ and of N , so that |δΨ = δ(|Ψ BCS ). We can then introduce the now time dependent quadratic Hamiltonian H(t) constructed in section 3It remains to check that the ket evolving this way remains of the BCS form. To this end, one shows, with the reasoning having led to Eq.(217), that ψσ |Ψ BCS (t) is equal to a linear combination of the ψ † σ ′ ( r ) acting on |Ψ BCS (t) . One then sees that H|Ψ BCS is of the form "constant plus polynomial of degree exactly two in ψ † " acting on |Ψ BCS (t) , which can be reproduced by an appropriate time dependence of Φ.

	(247)	S =	t f ti	dt	ih 2		Ψ(t)|	d dt	|Ψ(t) -c.c. -Ψ(t)|H(t)|Ψ(t)
	for fixed initial |Ψ(t i ) and final |Ψ(t f ) values of the state vector. For an arbitrary variation |δΨ(t) of the ket |Ψ(t) , subject to the condition
	(248)					|δΨ(t i ) = |δΨ(t f ) = 0,
	we calculate the first order variation of the action,
				t f
	(249)		δS =
									. 2 such that, for all
	variations within BCS family,
	(250) (δ Ψ One then sees that δS in Eq.(249) identically vanishes if the ket evolves with the quadratic
	Hamiltonian,
	(251)				ih	d dt	|Ψ BCS (t) = H(t)|Ψ BCS (t) .

ti dt δΨ(t)| ih d dt -H(t) |Ψ(t) + c.c. . BCS (t)|)H(t)|Ψ BCS (t) = (δ Ψ BCS (t)|)H(t)|Ψ BCS (t) .

  r, p, t) =

		p 2 2m -µ eff ( r, t) |∆|( r, t)	|∆|( r, t) -p 2 2m + µ eff ( r, t)	.
	We have introduced the position and time dependent function,
	(273)	µ

eff ( r, t) ≡ µ -U ( r, t) -1 2 mv 2 ( r, t) + m v ( r, t) • ( Ω(t) × r ) -∂ t S( r, t),