Félix Werner

Yvan Castin

General relations for quantum gases in two and three dimensions: Two-component fermions

Keywords: numbers: 67.85.Lm, 67.85.-d, 34.50.-s, 31.15.ac I. GENERAL INTRODUCTION

We derive exact general relations between various observables for N spin-1/2 fermions with zerorange or short-range interactions, in continuous space or on a lattice, in two or three dimensions, in an arbitrary external potential. Some of our results generalize known relations between the large-momentum behavior of the momentum distribution, the short-distance behaviors of the pair distribution function and of the one-body density matrix, the norm of the regular part of the wavefunction, the derivative of the energy with respect to the scattering length or to time, and the interaction energy (in the case of finite-range interactions). The expression relating the energy to a functional of the momentum distribution is also generalized. Moreover, expressions are found (in terms of the regular part of the wavefunction) for the derivative of the energy with respect to the effective range re in 3D, and to the effective range squared in 2D. They express the fact that the leading corrections to the eigenenergies due to a finite interaction-range are linear in the effective range in 3D (and in its square in 2D) with model-independent coefficients. There are subtleties in the validity condition of this conclusion, for the 2D continuous space (where it is saved by factors that are only logarithmically large in the zero-range limit) and for the 3D lattice models (where it applies only for some magic dispersion relations on the lattice, that sufficiently weakly break Galilean invariance and that do not have cusps at the border of the first Brillouin zone; an example of such relations is constructed). Furthermore, the subleading short distance behavior of the pair distribution function and the subleading 1/k 6 tail of the momentum distribution are related to ∂E/∂re (or to ∂E/∂(r 2 e ) in 2D). The second order derivative of the energy with respect to the inverse (or the logarithm in the two-dimensional case) of the scattering length is found to be expressible, for any eigenstate, in terms of the eigenwavefunctions' regular parts; this implies that, at thermal equilibrium, this second order derivative, taken at fixed entropy, is negative. Applications of the general relations are presented: We compute corrections to exactly solvable two-body and three-body problems and find agreement with available numerics; for the unitary gas in an isotropic harmonic trap, we determine how the finite-1/a and finite range energy corrections vary within each energy ladder (associated to the SO(2, 1) dynamical symmetry) and we deduce the frequency shift and the collapse time of the breathing mode; for the bulk unitary gas, we compare to fixed-node Monte Carlo data, and we estimate the deviation from the Bertsch parameter due to the finite interaction range in typical experiments.

The experimental breakthroughs of 1995 having led to the first realization of a Bose-Einstein condensate in an atomic vapor [1][2][3] have opened the era of experimental studies of ultracold gases with non-negligible or even strong interactions, in dimension lower than or equal to three [4][5][START_REF]Coherent atomic matter waves[END_REF][START_REF]Lecture notes of Les Houches Spring School[END_REF][START_REF]Ultra-Cold Fermi gases[END_REF]. In these systems, the thermal de Broglie wavelength and the typical distance between atoms are much larger than the range of the interaction potential. This so-called zero-range regime has interesting universal properties: Several quantities such as the thermodynamic functions of the gas depend on the interaction potential only through the scattering length a, a length that can be defined in any dimension and that characterizes the low-energy scattering amplitude of two atoms. This universality property holds for the weakly repulsive Bose gas in three dimensions [START_REF] Lee | [END_REF] up to the order of expansion in (na 3 ) 1/2 corresponding to Bogoliubov the-ory [10,11], n being the gas density. It also holds for the weakly repulsive Bose gas in two dimensions [12][13][14][15], even at the next order beyond Bogoliubov theory [16]. For a much larger than the range of the interaction potential, the ground state of N bosons in two dimensions is a universal N -body bound state [17][18][19][20][21]. In one dimension, the universality holds for any scattering length, as exemplified by the fact that the Bose gas with zerorange interaction is exactly solvable by the Bethe ansatz both in the repulsive case [22] and in the attractive case [23][24][25].

For spin 1/2 fermions, the universality properties are expected to be even stronger. The weakly interacting regimes in 3D [26][27][28][29][30][31] and in 2D [32] are universal, as well as for any scattering length in the Bethe-ansatz-solvable 1D case [33,[START_REF] Gaudin | La fonction d'onde de Bethe[END_REF]. Universality is also expected to hold for an arbitrary scattering length even in 3D, as was recently tested by experimental studies on the BEC-BCS crossover using a Feshbach resonance, see [START_REF]Ultra-Cold Fermi gases[END_REF] and Refs.

therein and e. g. [START_REF] Partridge | [END_REF][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52], and in agreement with unbiased Quantum Monte Carlo calculations [53][54][55][56][57][58][59]; and in 2D, a similar universal crossover from BEC to BCS is expected when the parameter ln(k F a) [where k F is the Fermi momentum] varies from -∞ to +∞ [60][61][62][63][64][65][66][67]. Mathematically, results on universality were obtained for the Nbody problem in 2D [68]. In 3D, mathematical results were obtained for the 3-body problem (see, e.g., [69][70][71][72][73]). The universality for the fermionic equal-mass N -body problem in 3D remains mathematically unproven [215].

Universality is also expected for mixtures in 2D [64,74,75], and in 3D for Fermi-Fermi mixtures below a critical mass ratio [74,[76][77][78]. Above a critical mass ratio, the Efimov effect takes place, as it also takes place for bosons [79,80]. In this case, the three-body problem depends on a single additional parameter, the three-body parameter. The Efimov physics is presently under active experimental investigation [81][82][83][84][85][86][87]. It is not the subject of this paper (see [88]).

In the zero-range regime, it is intuitive that the shortrange or high-momenta properties of the gas are dominated by two-body physics. For example the pair distribution function g (2) (r 12 ) of particles at distances r 12 much smaller than the de Broglie wavelength is expected to be proportional to the modulus squared of the zeroenergy two-body scattering-state wavefunction φ(r 12 ), with a proportionality factor Λ g depending on the manybody state of the gas. Similarly the tail of the momentum distribution n(k), at wavevectors much larger than the inverse de Broglie wavelength, is expected to be proportional to the modulus squared of the Fourier component φ(k) of the zero-energy scattering-state wavefunction, with a proportionality factor Λ n depending on the many-body state of the gas: Whereas two colliding atoms in the gas have a center of mass wavevector of the order of the inverse de Broglie wavelength, their relative wavevector can access much larger values, up to the inverse of the interaction range, simply because the interaction potential has a width in the space of relative momenta of the order of the inverse of its range in real space.

For these intuitive reasons, and with the notable exception of one-dimensional systems, one expects that the mean interaction energy E int of the gas, being sensitive to the shape of g (2) at distances of the order of the interaction range, is not universal, but diverges in the zerorange limit; one also expects that, apart from the 1D case, the mean kinetic energy, being dominated by the tail of the momentum distribution, is not universal and diverges in the zero-range limit, a well known fact in the context of Bogoliubov theory for Bose gases and of BCS theory for Fermi gases. Since the total energy of the gas is universal, and E int is proportional to Λ g while E kin is proportional to Λ n , one expects that there exists a simple relation between Λ g and Λ n .

The precise link between the pair distribution function, the tail of the momentum distribution and the energy of the gas was first established for one-dimensional systems. In [22] the value of the pair distribution function for r 12 = 0 was expressed in terms of the derivative of the gas energy with respect to the one-dimensional scattering length, thanks to the Hellmann-Feynman theorem. In [89] the tail of n(k) was also related to this derivative of the energy, by using a simple and general property of the Fourier transform of a function having discontinuous derivatives in isolated points.

In three dimensions, results in these directions were first obtained for weakly interacting gases. For the weakly interacting Bose gas, Bogoliubov theory contains the expected properties, in particular on the short distance behavior of the pair distribution function [90][91][92] and the fact that the momentum distribution has a slowly decreasing tail. For the weakly interacting spin-1/2 Fermi gas, it was shown that the BCS anomalous average (or pairing field) ψ↑ (r 1 ) ψ↓ (r 2 ) behaves at short distances as the zero-energy two-body scattering wavefunction φ(r 12 ) [93], resulting in a g (2) function indeed proportional to |φ(r 12 )| 2 at short distances. It was however understood later that the corresponding proportionality factor Λ g predicted by BCS theory is incorrect [94], e.g. at zero temperature the BCS prediction drops exponentially with 1/a in the non-interacting limit a → 0 -, whereas the correct result drops as a power law in a.

More recently, in a series of two articles [95,96], explicit expressions for the proportionality factors Λ g and Λ n were obtained in terms of the derivative of the gas energy with respect to the inverse scattering length, for a spin-1/2 interacting Fermi gas in three dimensions, for an arbitrary value of the scattering length, that is, not restricting to the weakly interacting limit. Later on, these results were rederived in [97][98][99], and also in [100] with very elementary methods building on the aforementioned intuition that g (2) (r 12 ) ∝ |φ(r 12 )| 2 at short distances and n(k) ∝ | φ(k)| 2 at large momenta. These relations were tested by numerical four-body calculations [101]. An explicit relation between Λ g and the interaction energy was derived in [99]. Another fundamental relation discovered in [95] and recently generalized in [102,103] to fermions in 2D, expresses the total energy as a functional of the momentum distribution and the spatial density.

II. CONTENTS

Here we derive generalizations of the relations of [22,89,95,96,99,102,103] to two dimensional gases, and to the case of a small but non-zero interaction range (both on a lattice and in continuous space). We also find entirely new results for the first order derivative of the energy with respect to the effective range, as well as for the second order derivative with respect to the scattering length. We shall also include rederivations of known relations using our elementary methods. We treat in detail the case of spin-1/2 fermions, with equal masses in the two spin states, both in three dimensions and in two dimensions. The discussion of spinless bosons and arbitrary mixtures is deferred to another article, as it may involve the Efimov effect in three dimensions [START_REF] Werner | Sections III and IV of[END_REF].

This article is organized as follows. Models, notations and some basic properties are introduced in Section III. Relations for zero-range interactions are summarized in Table II and derived for pure states in Section IV. We then consider lattice models (Tab. III and Sec. V) and finite-range models in continuous space (Tab. IV and Sec. VI). In Section VII we derive a model-independent expression for the correction to the energy due to a finite range or a finite effective range of the interaction, and we relate this energy correction to the subleading short distance behavior of the pair distribution function and to the coefficient of the 1/k 6 subleading tail of the momentum distribution (see Tab. V). The case of general statistical mixtures of pure states or of stationary states is discussed in Sec. VIII, and the case of thermodynamic equilibrium states in Sec. IX. Finally we present applications of the general relations: For two particles and three particles in harmonic traps we compute corrections to exactly solvable cases (Sec. X A and Sec. X B). For the unitary gas trapped in an isotropic harmonic potential, we determine how the equidistance between levels within a given energy ladder (resulting from the SO(2, 1) dynamical symmetry) is affected by finite 1/a and finite range corrections, which leads to a frequency shift and a collapse of the breathing mode of the zero-temperature gas (Sec. X C). For the bulk unitary gas, we check that general relations are satisfied by existing fixed-node Monte Carlo data [START_REF] Astrakharchik | [END_REF][106][107] for correlation functions of the unitary gas (Sec. X D). We quantify the finite range corrections to the unitary gas energy in typical experiments, which is required for precise measurements of its equation of state (Sec. X E). We conclude in Section XI.

III. MODELS, NOTATIONS, AND BASIC PROPERTIES

We now introduce the three models used in this work to account for interparticle interactions and associated notations, together with some basic properties to be used in some of the derivations.

For a fixed number N σ of fermions in each spin state σ =↑, ↓, one can consider that particles 1, . . . , N ↑ have a spin ↑ and particles N ↑ + 1, . . . , N ↑ + N ↓ = N have a spin ↓, so that the wavefunction ψ(r 1 , . . . , r N ) (normalized to unity) changes sign when one exchanges the positions of two particles having the same spin [216].

A. Zero-range model

In this well-known model (see e.g. [79,80,[START_REF] Albeverio | Solvable models in quantum mechanics[END_REF][START_REF] Castin | Coherent atomic matter waves[END_REF][START_REF] Petrov | [END_REF][111][START_REF] Castin | Ultracold Fermi gases[END_REF][START_REF] Werner | [END_REF] and refs. therein) the interaction potential is replaced by boundary conditions on the N -body wavefunction: For any pair of particles i = j, there exists a function A ij , hereafter called regular part of ψ, such that [Tab. I, Eq. (1a)] holds in the 3D case and [Tab. I, Eq. (1b)] holds in the 2D case, where the limit of vanishing distance r ij between particles i and j is taken for a fixed position of their center of mass R ij = (r i + r j )/2 and fixed positions of the remaining particles (r k ) k =i,j different from R ij . Fermionic symmetry of course imposes A ij = 0 if particles i and j have the same spin. When none of the r i 's coincide, there is no interaction potential and Schrödinger's equation reads H ψ(r 1 , . . . , r N ) = E ψ(r 1 , . . . , r N ) with

H = - 2 2m N i=1
∆ ri + H trap , where m is the atomic mass and the trapping potential energy is

H trap ≡ N i=1 U (r i ), (1)
U being an external trapping potential. The crucial difference between the Hamiltonian H and the noninteracting Hamiltonian is the boundary condition [Tab. I, Eqs. (1a,1b)].

B. Lattice models

These models are used for quantum Monte Carlo calculations [53-56, 58, 114]. They can also be convenient for analytics, as used in [15,16,100,115] and in this work. Particles live on a lattice, i. e. the coordinates are integer multiples of the lattice spacing b. The Hamiltonian is

H = H kin + H int + H trap (2)
with, in first quantization, the kinetic energy

H kin = - 2 2m N i=1 ∆ ri , (3)
the interaction energy

H int = g 0 i<j δ ri,rj b -d , (4)
and the trapping potential energy defined by (1); i.e. in second quantization

H kin = σ D d d k (2π) d ǫ k c † σ (k)c σ (k) (5)
H int = g 0 r b d (ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ )(r) (6)
H trap = r,σ b d U (r)(ψ † σ ψ σ )(r). ( 7
)
Here d is the space dimension, ǫ k is the dispersion relation, ψ obeys discrete anticommutation relations { ψσ (r), ψ † σ ′ (r ′ )} = b -d δ rr ′ δ σσ ′ . The operator c † σ (k) creates a particle in the plane wave state |k defined by r|k = e ik•r for any k belonging to the first Brillouin

Three dimensions

Two dimensions ψ(r1, . . . , rN ) =

r ij →0 1 rij - 1 a Aij (Rij , (r k ) k =i,j ) + O(rij ) (1a) ψ(r1, . . . , rN ) = r ij →0 ln(rij /a) Aij (Rij , (r k ) k =i,j ) + O(rij) (1b) (A (1) , A (2) ) ≡ i<j k =i,j d d r k d d Rij A (1) * ij (Rij , (r k ) k =i,j )A (2) ij (Rij , (r k ) k =i,j ) (2)
(A (1) , HA (2) )

≡ i<j k =i,j d d r k d d Rij A (1) * ij (Rij , (r k ) k =i,j )Hij A (2) ij (Rij , (r k ) k =i,j ) (3)
TABLE I: Notation for the regular part A of the N -body wavefunction appearing in the contact conditions (first line, with Rij = (ri + rj )/2 fixed), for the scalar product between such regular parts (second line) and for corresponding matrix elements of operators Hij acting on Rij and on the r k , k = i, j (third line). 15,16,55,58,115].

zone D = -π b , π b d . The corresponding anticommuta- tion relations are {c σ (k), c † σ ′ (k ′ )} = (2π) d δ σσ ′ δ(k -k ′ ) if k and k ′ are
ǫ k = 2 k 2 2m [
Another choice, used in [54,[START_REF] Burovski | [END_REF], is the dispersion relation of the Hubbard model:

ǫ k = 2 mb 2 d i=1 [1 -cos(k i b)].
More generally, what follows applies to any ǫ k such that

ǫ k → b→0 2 k 2 2m sufficiently rapidly and ǫ -k = ǫ k .
A key quantity is the zero-energy scattering state φ(r), defined by the two-body Schrödinger equation (with the center of mass at rest)

- 2 m ∆ r + g 0 δ r,0 b d φ(r) = 0 (8)
and by the normalization conditions

φ(r) ≃ r≫b 1 r - 1 a in 3D (9) φ(r) ≃ r≫b ln(r/a) in 2D. (10)
A two-body analysis, detailed in Appendix A, yields the relation between the scattering length and the bare coupling constant g 0 , in three and two dimensions:

1 g 0 3D = m 4π 2 a - D d 3 k (2π) 3 1 2ǫ k (11) 1 g 0 2D = lim q→0 - m 2π 2 ln( aqe γ 2 )+ D d 2 k (2π) 2 P 1 2(ǫ q -ǫ k ) ( 12
)
where γ = 0.577216 . . . is Euler's constant and P is the principal value. This implies that (for constant b):

d(1/g 0 ) d(1/a) = m 4π 2 in 3D (13) d(1/g 0 ) d(ln a) = - m 2π 2 in 2D. ( 14
)
Another useful property derived in Appendix A is

φ(0) = - 4π 2 mg 0 in 3D (15)
φ(0) = 2π 2 mg 0 in 2D, (16)
which, together with (13,14), gives

|φ(0)| 2 = 4π 2 m d(-1/a) dg 0 in 3D (17) |φ(0)| 2 = 2π 2 m d(ln a) dg 0 in 2D. ( 18
)
In the zero-range limit (b → 0 with g 0 adjusted in such a way that a remains constant), it is expected that the spectrum of the lattice model converges to the one of the zero-range model, as explicitly checked for three particles in [115], and that any eigenfunction ψ(r 1 , . . . , r N ) of the lattice model tends to the corresponding eigenfunction of the zero-range model provided all interparticle distances remain much larger than b. For any stationary state, let us denote by 1/k typ the typical length-scale on which the zero-range model's wavefunction varies: e.g. for the lowest eigenstates, this is on the order of the mean interparticle distance, or on the order of a in the regime where a is small and positive and dimers are formed. The zerorange limit is then reached if k typ b ≪ 1. This notion of typical wavevector k typ can also be applied to the case of a thermal equilibrium state, since most significantly populated eigenstates then have a k typ on the same order; it is then expected that the thermodynamic potentials converge to the ones of the zero-range model when b → 0, and that this limit is reached provided k typ b ≪ 1. For the homogeneous gas, defining a thermal wavevector k T by

2 k 2 T /(2m) = k B T , we have k typ ∼ max(k F , k T ) for a < 0 and k typ ∼ max(k F , k T , 1/a) for a > 0.
For lattice models, it will prove convenient to define the regular part A by

ψ(r 1 , . . . , r i = R ij , . . . , r j = R ij , . . . , r N ) = φ(0) × A ij (R ij , (r k ) k =i,j ). ( 19
)
In the zero-range regime k typ b ≪ 1, when the distance r ij between two particles of opposite spin is ≪ 1/k typ while all the other interparticle distances are much larger than b and than r ij , the many-body wavefunction is proportional to φ(r jr i ), with a proportionality constant given by (19):

ψ(r 1 , . . . , r N ) ≃ φ(r j -r i ) A ij (R ij , (r k ) k =i,j ) (20)
where R ij = (r i + r j )/2. If moreover r ij ≫ b, φ can be replaced by its asymptotic form [START_REF] Lee | [END_REF]10); since the contact conditions [Tab. I, Eqs. (1a,1b)] of the zero-range model must be recovered, we see that the lattice model's regular part tends to the zero-range model's regular part in the zero-range limit.

C. Finite-range continuous-space models

Such models are used in numerical few-body correlated Gaussian and many-body fixed-node Monte Carlo calculations (see e. g. [5,65,101,[START_REF] Astrakharchik | [END_REF][116][117][118] and refs. therein). They are also relevant to neutron matter [119]. The Hamiltonian reads

H = H 0 + N ↑ i=1 N j=N ↑ +1 V (r ij ), (21)
H 0 being defined by (3) where ∆ ri now stands for the usual Laplacian, and V (r) is an interaction potential between particles of opposite spin, which vanishes for r > b or at least decays quickly enough for r ≫ b. The twobody zero-energy scattering state φ(r) is again defined by the Schrödinger equation -( 2 /m)∆ r φ + V (r)φ = 0 and the boundary condition [START_REF] Lee | [END_REF] or (10). The zero-range regime is again reached for k typ b ≪ 1 with k typ the typical relative wavevector [218]. Equation (20) again holds in the zero-range regime, where A now simply stands for the zero-range model's regular part.

IV. RELATIONS IN THE ZERO-RANGE LIMIT

We now derive relations for the zero-range model. For some of the derivations we will use a lattice model and then take the zero-range limit. We recall that we derive all relations for pure states in this section, the generalization to statistical mixtures and the discussion of thermal equilibrium being deferred to Sections VIII and IX.

A. Tail of the momentum distribution

In this subsection as well as in the following subsections IV B, IV D, IV E, IV G, we consider a many-body pure state whose wavefunction ψ satisfies the contact condition [Tab. I, Eqs. (1a,1b)]. We now show that the momentum distribution n σ (k) has a σ-independent tail proportional to 1/k 4 , with a coefficient denoted by C [Tab. II, Eq. ( 1)]. C is usually referred to as the "contact". We shall also show that C is related by [Tab. II, Eqs. (2a,2b)] to the norm of the regular part A of the wavefunction (defined in Tab. I). In 3D these results were obtained in [96] [219]. Here the momentum distribution is defined in second quantization by

n σ (k) = nσ (k) = c † σ (k)c σ (k)
where c σ (k) annihilates a particle of spin σ in the plane-wave state |k defined by r|k = e ik•r ; this corresponds to the normalization

d d k (2π) d n σ (k) = N σ . (22)
In first quantization,

n σ (k) = i:σ l =i d d r l d d r i e -ik•ri ψ(r 1 , . . . , r N ) 2 (23)
where the sum is taken over all particles of spin σ: i runs from 1 to N ↑ for σ =↑, and from N ↑ + 1 to N for σ =↓.

Three dimensions:

The key point is that in the large-k limit, the Fourier transform with respect to r i is dominated by the contribution of the short-distance divergence coming from the contact condition [Tab. I, Eq. (1a)]:

d 3 r i e -ik•ri ψ(r 1 , . . . , r N ) ≃ k→∞ d 3 r i e -ik•ri × j,j =i 1 r ij A ij (r j , (r l ) l =i,j ). (24)
A similar link between the short-distance singularity of the wavefunction and the tail of its Fourier transform was used to derive exact relations in 1D in [89]. From ∆(1/r) = -4πδ(r), we have d 3 r e -ik•r 1 r = 4π k 2 , so that

d 3 r i e -ik•ri ψ(r 1 , . . . , r N ) ≃ k→∞ 4π k 2 j,j =i e -ik•rj × A ij (r j , (r l ) l =i,j ). ( 25
)
One inserts this into (23) and expands the modulus squared. After spatial integration over all the r l , l = i, the crossed terms rapidly vanish in the large-k limit, as they are the product of e ik•(rj-r j ′ ) and of regular functions of r j and

r j ′ [220]. This yields n σ (k) ∼ k→∞ C/k 4 ,
with the expression [Tab. II, Eq. (2a)] of C in terms of the norm (A, A) defined in [Tab. I, Eq. ( 2)].

Two dimensions:

The 2D contact condition [Tab. I, Eq. (1b)] now gives

d 2 r i e -ik•ri ψ(r 1 , . . . , r N ) ≃ k→∞ d 2 r i e -ik•ri × j,j =i ln(r ij )A ij (r j , (r l ) l =i,j ). ( 26
)
Three dimensions Two dimensions (with real wavefunctions for line 7). Line 2 holds for any pure state; here A is the regular part of the wavefunction appearing in the contact condition, and (A, A) is its squared norm (defined in Tab. I). Lines 4 and 8 hold for any stationary state. Lines 9-11 hold at thermal equilibrium in the canonical ensemble. Line 12 holds for any time-dependence of scattering length and trapping potential, and any corresponding time-dependent statistical mixture. Many of the 3D relations were originally obtained in [95,96] (see text), while the 2D relation (5b) was obtained in [103] for the homogeneous system and in [102] (in a different form) for the general case.

C ≡ lim k→+∞ k 4 nσ(k) (1)
C = (4π) 2 (A, A) (2a) C = (2π) 2 (A, A) (2b) d 3 R g (2) ↑↓ R + r 2 , R - r 2 ∼ r→0 C (4π) 2 1 r 2 (3a) d 2 R g (2) ↑↓ R + r 2 , R - r 2 ∼ r→0 C (2π) 2 ln 2 r (3b) dE d(-1/a) = 2 C 4πm (4a) dE d(ln a) = 2 C 2πm (4b) E -Etrap = 2 C 4πma E -Etrap = lim Λ→∞ - 2 C 2πm ln aΛe γ 2 + σ d 3 k (2π) 3 2 k 2 2m nσ(k) - C k 4 (5a) + σ k<Λ d 2 k (2π) 2 2 k 2 2m nσ(k) (5b) d 3 R g (1) σσ R + r 2 , R - r 2 = r→0 Nσ - C 8π r + O(r 2 ) (6a) d 2 R g (1) σσ R + r 2 , R - r 2 = r→0 Nσ + C 4π r 2 ln r + O(r 2 ) (6b) 1 3 3 i=1 σ d 3 R g (1) σσ R + ru i 2 , R - ru i 2 = r→0 N 1 2 2 i=1 σ d 2 R g (1) σσ R + ru i 2 , R - ru i 2 = r→0 N - C 4π r - m 3 2 E -Etrap - 2 C 4πma r 2 + o(r 2 ) (7a) + C 4π r 2 ln r a -1 - m 2 2 (E -Etrap) r 2 + o(r 2 ) (7b) 1 2 d 2 En d(-1/a) 2 = 4π 2 m 2 n ′ ,E n ′ =En |(A (n ′ ) , A (n) )| 2 En -E n ′ (8a) 1 2 d 2 En d(ln a) 2 = 2π 2 m 2 n ′ ,E n ′ =En |(A (n ′ ) , A (n) )| 2 En -E n ′ (8b) d Ē d(-1/a) S = dF d(-1/a) T = 2 C 4πm (9a) d Ē d(ln a) S = dF d(ln a) T = 2 C 2πm (9b) d 2 F d(-1/a) 2 T < 0 (10a) d 2 F d(ln a) 2 T < 0 (10b) d 2 Ē d(-1/a) 2 S < 0 (11a) d 2 Ē d(ln a) 2 S < 0 (11b) dE dt = 2 C 4πm d(-1/a) dt + dHtrap dt (12a) dE dt = 2 C 2πm d(ln a) dt + dHtrap dt (12b)
From ∆(ln r) = 2πδ(r), one has d 2 r e -ik•r ln r = -2π k 2 and

d 2 r i e -ik•ri ψ(r 1 , . . . , r N ) ≃ k→∞ - 2π k 2 j,j =i e -ik•rj × A ij (r j , (r l ) l =i,j ). ( 27
)
As in 3D this leads to [Tab. II, Eq. (2b)].

B. Pair distribution function at short distances

The pair distribution function gives the probability density of finding a spin-↑ particle at r ↑ and a spin-↓ particle at r ↓ : g

(2) ↑↓ (r ↑ , r ↓ ) = ( ψ † ↑ ψ↑ )(r ↑ )( ψ † ↓ ψ↓ )(r ↓ ) = ( N k=1 d d r k )|ψ(r 1 , . . . , r N )| 2 N ↑ i=1 N j=N ↑ +1 δ (r ↑ -r i ) δ (r ↓ -r j ).
We set r ↑,↓ = R ± r/2 and we integrate over r i and r j :

g (2) ↑↓ R + r 2 , R - r 2 = N ↑ i=1 N j=N ↑ +1 k =i,j d d r k ψ r 1 , . . . , r i = R + r 2 , . . . , r j = R - r 2 , . . . , r N 2 (28)
Let us define the spatially integrated pair distribution function

[221] G (2) ↑↓ (r) ≡ d d R g (2) ↑↓ R + r 2 , R - r 2 , (29)
whose small-r singular behavior we will show to be related to C via [Tab. II, Eqs. (3a,3b)].

Three dimensions:

Replacing the wavefunction in ( 28) by its asymptotic behavior given by the contact condition [Tab. I, Eq. (1a)] immediately yields

G (2) ↑↓ (r) ∼ r→0 (A, A) r 2 . ( 30
)
Expressing (A, A) in terms of C through [Tab. II, Eq. (2a)] finally gives [Tab. II, Eq. (3a)].

In a measurement of all particle positions, the mean total number of pairs of particles of opposite spin which are separated by a distance smaller than s is N pair (s) =

r<s d d r G (2)
↑↓ (r), so that from [Tab. II, Eq. (3a)]

N pair (s) ∼ s→0 C 4π s, (31)
as obtained in [95,96].

Two dimensions:

The contact condition [Tab. I, Eq. (1b)] similarly leads to [Tab. II, Eq. (3b)]. After integration over the region r < s this gives

N pair (s) ∼ s→0 C 4π s 2 ln 2 s. (32)
C. First order derivative of the energy with respect to the scattering length

The relations [Tab. II, Eqs. (4a,4b)] can be derived straightforwardly using the lattice model, see Sec.V E. Here we derive them by directly using the zero-range model, which is more involved but also instructive.

Three dimensions:

Let us consider a wavefunction ψ 1 satisfying the contact condition [Tab. I, Eq. (1a)] for a scattering length a 1 . We denote by A (1) ij the regular part of ψ 1 appearing in the contact condition [Tab. I, Eq. (1a)]. Similarly, ψ 2 satisfies the contact condition for a scattering length a 2 and a regular part A

(2) ij . Then, as shown in Appendix B using the divergence theorem, the following lemma holds: 1) , A (2) ) (33) where the scalar product between regular parts is defined by [Tab. I, Eq. ( 2)]. We then apply (33) to the case where ψ 1 and ψ 2 are N -body stationary states of energy E 1 and E 2 . The left hand side of (33) then reduces to (E 2 -E 1 ) ψ 1 |ψ 2 . Taking the limit a 2 → a 1 gives

ψ 1 , Hψ 2 -Hψ 1 , ψ 2 = 4π 2 m 1 a 1 - 1 a 2 (A (
dE d(-1/a) = 4π 2 m (A, A) (34)
for any stationary state. Expressing (A, A) in terms of C thanks to [Tab. II, Eq. (2a)] finally yields [Tab. II, Eq. ( 4a)]. This result as well as [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] is contained in Ref. [95, 96][222]. We recall that here and in what follows, the wavefunction is normalized: ψ|ψ = 1.

Two dimensions:

The 2D version of the lemma (33) is

ψ 1 , Hψ 2 -Hψ 1 , ψ 2 = 2π 2 m
ln (a 2 /a 1 ) (A (1) , A (2) ), Finally we mention the equivalent form of relation [Tab. II, Eq. (5a)]:

(
E -E trap = lim Λ→∞ 2 C 4πm 1 a - 2Λ π + σ k<Λ d 3 k (2π) 3 2 k 2 2m n σ (k) . ( 37
)
Two dimensions:

The 2D version of ( 37) is [Tab. II, Eq. (5b)]. This was shown for a homogeneous system in [103] and in the general case in [102] [223]. This can easily be rewritten in the following forms, which resemble [Tab. II, Eq. (5a)]:

E -E trap = - 2 C 2πm ln aqe γ 2 + σ d 2 k (2π) 2 2 k 2 2m × n σ (k) - C k 4 θ(k -q)
for any q > 0, (38) where the Heaviside function θ ensures that the integral converges at small k, or equivalently

E -E trap = - 2 C 2πm ln aqe γ 2 + σ d 2 k (2π) 2 2 k 2 2m × n σ (k) - C k 2 (k 2 + q 2 )
for any q > 0. (39) To derive this we again use the lattice model. We note that, if the limit q → 0 is replaced by the limit b → 0 taken for fixed a, Eq. ( 12) remains true (see Appendix A); repeating the reasoning of Section V B then shows that [Tab. III, Eq. (3b)] remains true; as in 3D we finally get in the limit b → 0

E -E trap = - 2 C 2πm ln aqe γ 2 + σ d 2 k (2π) 2 2 k 2 2m × n σ (k) - C k 2 P 1 k 2 -q 2 (40)
for any q > 0; this is easily rewritten as [Tab. II, Eq. (5b)].

E. One-body density matrix at short-distances

The one-body density matrix is defined as g (1) σσ (r, r ′ ) = ψ † σ (r) ψσ (r ′ ) where ψσ (r) annihilates a particle of spin σ at point r. Its spatially integrated version G (1) σσ (r) ≡ d d R g (1) σσ R -

r 2 , R + r 2 (41)
is a Fourier transform of the momentum distribution:

G (1) σσ (r) = d d k (2π) d e ik•r n σ (k). ( 42
)
The expansion of G

σσ (r) up to first order in r is given by [Tab. II, Eq. (6a)] in 3D, as first obtained in [95], and by [Tab. II, Eq. (6b)] in 2D. The expansion can be pushed to second order if one sums over spin and averages over d orthogonal directions of r, see [Tab. II, Eqs. (7a,7b)] where the u i 's are an orthonormal basis [224]. Such a second order expansion was first obtained in 1D in [89]; the following derivations however differ from the 1D case [225].

Three dimensions:

To derive [Tab. II, Eqs. (6a,7a)] we rewrite (42) as

G (1) σσ (r) = N σ + d 3 k (2π) 3 e ik•r -1 C k 4 + d 3 k (2π) 3 e ik•r -1 n σ (k) - C k 4 . ( 43
)
The first integral equals -(C/8π)r. In the second integral, we use

e ik•r -1 = r→0 ik • r - (k • r) 2 2 + o(r 2 ). ( 44
)
The first term of this expansion gives a contribution to the integral proportional to the total momentum of the gas, which vanishes since the eigenfunctions are real.

σσ (r) = N σ + I(r) + J(r) with

I(r) = d 2 k (2π) 2 e ik•r -1 C k 4 θ(k -q) (45) J(r) = d 2 k (2π) 2 e ik•r -1 n σ (k) - C k 4 θ(k -q) (46)
where q > 0 is arbitrary and the Heaviside function θ ensures that the integrals converge.

To evaluate I(r) we use standard manipulations to rewrite it as I(r) = Cr 2 /(2π) +∞ qr dx[J 0 (x) -1]/x 3 , J 0 being a Bessel function. Expressing this integral with Mathematica in terms of an hypergeometric function and a logarithm leads for r → 0 to I(r) = Cr 2 /(8π)[γ -1ln 2 + ln(qr)] + O(r 4 ). To evaluate J(r) we use the same procedure as in 3D: expanding the exponential [see (44)] yields an integral which can be related to the total energy thanks to (38) [226].

F. Second order derivative of the energy with respect to the scattering length

We denote by |ψ n an orthonormal basis of N -body stationary states that vary smoothly with 1/a, and by E n the corresponding eigenenergies. We will derive [Tab. II, Eqs. (8a,8b)], where the sum is taken on all values of n ′ such that E n ′ = E n . This implies that for the ground state energy E 0 ,

d 2 E 0 d(-1/a) 2 < 0 in 3D ( 47
)
d 2 E 0 d(ln a) 2 < 0 in 2D. ( 48
)
Eq. ( 47) was intuitively expected [120]: Eq. (31) shows that dE 0 /d(-1/a) is proportional to the probability of finding two particles very close to each other, and it is natural that this probability decreases when one goes from the BEC limit (-1/a → -∞) to the BCS limit (-1/a → +∞), i.e. when the interactions become less attractive [227]. Eq. ( 48) also agrees with intuition [228].

For the derivation, it is convenient to use the lattice model (defined in Sec. III B): As shown in Sec.V F one easily obtains (60) and [Tab. III, Eq. ( 6)], from which the result is deduced as follows. |φ(0)| 2 is eliminated using (17,18). Then, in 3D, one uses

d 2 E n d(-1/a) 2 = d 2 E n dg 2 0 dg 0 d(-1/a) 2 + dE n dg 0 d 2 g 0 d(-1/a) 2 (49)
where the second term equals 2g 0 dE n /d(-1/a) m/(4π 2 ) and thus vanishes in the zero-range limit. In 2D, similarly, one uses the fact that d 2 E n /d(ln a) 2 is the zero-range limit of (d 2 E n /dg 2 0 ) • (dg 0 /d(ln a)) 2 .

G. Time derivative of the energy

We now consider the case where the scattering length a(t) and the trapping potential U (r, t) are varied with time. The time-dependent version of the zero-range model (see e.g. [121]) is given by Schrödinger's equation

i ∂ ∂t ψ(r 1 , . . . , r N ; t) = H(t) ψ(r 1 , . . . , r N ; t) (50)
when all particle positions are distinct, with

H(t) = N i=1 - 2 2m ∆ ri + U (r i , t) , (51)
and by the contact condition [Tab. I, Eq. (1a)] in 3D or [Tab. I, Eq. (1b)] in 2D for the scattering length a = a(t). One then has the relations [Tab. II, Eqs. (12a,12b)], where E(t) = ψ(t)|H(t)|ψ(t) is the total energy and H trap (t) = N i=1 U (r i , t) is the trapping potential part of the Hamiltonian. In 3D, this relation was first obtained in [96]. A very simple derivation of these relations using the lattice model is given in Section V G. Here we give a derivation within the zero-range model.

Three dimensions:

We first note that the generalization of the lemma (33) to the case of two Hamiltonians H 1 and H 2 with corresponding trapping potentials U 1 (r) and U 2 (r) reads: 1) , A (2) )

ψ 1 , H 2 ψ 2 -H 1 ψ 1 , ψ 2 = 4π 2 m 1 a 1 - 1 a 2 (A (
+ ψ 1 | N i=1 [U 2 (r i , t) -U 1 (r i , t)] |ψ 2 . ( 52
)
+ ψ(t)| N i=1 [U (r i , t + δt) -U (r i , t)] |ψ(t + δt) . ( 53
)
Dividing by δt, taking the limit δt → 0, and using the expression [Tab. II, Eq. (1a)] of (A, A) in terms of C, the right-hand-side of (53)

ψ 1 , H 2 ψ 2 -H 1 ψ 1 , ψ 2 =
2π 2 m ln(a 2 /a 1 )(A (1) , A (2) )

+ ψ 1 | N i=1 [U 2 (r i , t) -U 1 (r i , t)] |ψ 2 . ( 54
)

V. RELATIONS FOR LATTICE MODELS

In this Section, it will prove convenient to introduce an operator Ĉ by [Tab. III, Eqs. (1a,1b)] and to define C by its expectation value in the state of the system,

C = Ĉ ( 55
)
In the zero-range limit, this new definition of C coincides with the definition [Tab. II, Eq. ( 1)] of Section IV, as shown in Appendix C.

A. Interaction energy and Ĉ

The interaction part H int of the lattice model's Hamiltonian is obviously equal to g 0 dH dg 0 [see Eqs. (2,3,4)].

Rewriting this as 1 g 0 dH d(-1/g 0 )

, and using the simple expressions (13,14) for d(1/g 0 ), we get the relation [Tab. III, Eq. ( 2)] between H int and Ĉ, both in 3D and in 2D.

B. Total energy minus trapping potential energy in terms of momentum distribution and Ĉ

Here we derive [Tab. III, Eqs. (3a,3b)]. We start from the expression [Tab. III, Eq. ( 2)] of the interaction energy and eliminate 1/g 0 thanks to (11,12). The desired expression of H -H trap = H int +H kin then simply follows from the expression (5) of the kinetic energy.

Three dimensions

Two dimensions

Ĉ ≡ 4πm 2 dH d(-1/a) (1a) Ĉ ≡ 2πm 2 dH d(ln a) (1b) Hint = 4 m 2 Ĉ g0 (2)
H -Htrap = 2 Ĉ 4πma H -Htrap = lim q→0 - 2 Ĉ 2πm ln aqe γ 2 + σ D d 3 k (2π) 3 ǫ k nσ(k) -Ĉ 2 2mǫ k 2 (3a) + σ D d 2 k (2π) 2 ǫ k nσ(k) -Ĉ 2 2mǫ k P 2 2m(ǫ k -ǫq) (3b) C = (4π) 2 (A, A) (4a) C = (2π) 2 (A, A) (4b) dE d(-1/a) = 2 C 4πm (5a) dE d(ln a) = 2 C 2πm (5b) 1 2 d 2 En dg 2 0 = |φ(0)| 4 n ′ ,E n ′ =En |(A (n ′ ) , A (n) )| 2 En -E n ′ (6)
d 2 F dg 2 0 T < 0, d 2 E dg 2 0 S < 0 (7) R b 3 (ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ )(R) = Ĉ (4π) 2 |φ(0)| 2 (8a) R b 2 (ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ )(R) = Ĉ (2π) 2 |φ(0)| 2 (8b)
In the zero-range regime ktypb ≪ 1

R b 3 g (2) ↑↓ R + r 2 , R - r 2 ≃ C (4π) 2 |φ(r)| 2 , for r ≪ k -1 typ (9a) R b 2 g (2) ↑↓ R + r 2 , R - r 2 ≃ C (2π) 2 |φ(r)| 2 , for r ≪ k -1 typ (9b) nσ(k) ≃ C 2 2mǫ k 2 , for k ≫ ktyp (10)
TABLE III: Relations for spin-1/2 fermions for lattice models. Ĉ is defined in line 1 and C = Ĉ . Lines 2, 3 and 8 are relations between operators. Line 4 holds for any pure state [the regular part A being defined in Eq. ( 19) in the text]. Lines 5-6 hold for any stationary state. Line 7 holds at thermal equilibrium in the canonical ensemble. Lines 9-10 are expected to hold in the zero-range regime ktypb ≪ 1, where ktyp is the typical wavevector, for any stationary state or at thermal equilibrium.

C. Interaction energy and regular part

In the forthcoming subsections V D, V E and V F, we will use the following lemma: For any wavefunctions ψ and ψ ′ ,

ψ ′ |H int |ψ = g 0 |φ(0)| 2 (A ′ , A) (56)
where A and A ′ are the regular parts related to ψ and ψ ′ through (19), and the scalar product between regular parts is naturally defined as the discrete version of [Tab. I, Eq. ( 2)]:

(A ′ , A) ≡ i<j (r k ) k =i,j Rij b (N -1)d A ′ * ij (R ij , (r k ) k =i,j ) × A ij (R ij , (r k ) k =i,j ). ( 57
)
The lemma simply follows from

ψ ′ |H int |ψ = g 0 i<j (r k ) k =i,j b (N -2)d rj b d × (ψ ′ * ψ)(r 1 , . . . , r i = r j , . . . , r j , . . . , r N ). ( 58
)

D. Relation between Ĉ and (A, A)

Lemma ( 56) with ψ ′ = ψ writes

ψ|H int |ψ = g 0 |φ(0)| 2 (A, A). ( 59
)
Expressing ψ|H int |ψ in terms of C = ψ| Ĉ|ψ thanks to [Tab. III, Eq. ( 2)], and using the expressions (17,18) of |φ(0)| 2 , we get [Tab. III, Eqs. (4a,4b)].

E. First order derivative of an eigenenergy with respect to the coupling constant

For any stationary state, the Hellmann-Feynman theorem, together with the definition [Tab. III, Eqs. (1a,1b)] of Ĉ and the relation [Tab. III, Eqs. (4a,4b)] between C and (A, A), immediately yields [Tab. III, Eqs. (5a,5b)].

F. Second order derivative of an eigenenergy with respect to the coupling constant

We denote by |ψ n an orthonormal basis of N -body stationary states which vary smoothly with g 0 , and by E n the corresponding eigenenergies. We apply second order perturbation theory to determine how an eigenenergy varies for an infinitesimal change of g 0 . This gives:

1 2 d 2 E n dg 2 0 = n ′ ,E n ′ =En | ψ n ′ |H int /g 0 |ψ n | 2 E n -E n ′ , (60)
where the sum is taken over all values of n ′ such that E n ′ = E n . Lemma (56) then yields [Tab. III, Eq. ( 6)].

G. Time derivative of the energy

The relations [Tab. II, Eqs. (12a,12b)] remain exact for the lattice model. Indeed, dE/dt equals dH/dt from the Hellmann-Feynman theorem. In 3D, we can rewrite this quantity as dH trap /dt +d(-1/a)/dt dH/d(-1/a) , and the desired result follows from the definition [Tab. III, Eq. ( 1a)] of Ĉ. The derivation of the 2D relation [Tab. II, Eq. (12b)] is analogous.

H. On-site pair distribution operator

Let us define a spatially integrated pair distribution operator Ĝ( 2)

↑↓ (r) ≡ R b d (ψ † ↑ ψ ↑ ) R + r 2 (ψ † ↓ ψ ↓ ) R - r 2 .
(61) Using the relation [Tab. III, Eq. ( 2)] between Ĉ and H int , expressing H int in terms of Ĝ( 2) ↑↓ (0) thanks to the secondquantized form [START_REF]Coherent atomic matter waves[END_REF], and expressing g 0 in terms of φ(0) thanks to (15,16), we immediately get: Ĝ( 2)

↑↓ (0) = Ĉ (4π) 2 |φ(0)| 2 in 3D (62)
Ĝ( 2)

↑↓ (0) = Ĉ (2π) 2 |φ(0)| 2 in 2D. (63)
[Here, |φ(0)| 2 may of course be eliminated using (15,16).] These relations are analogous to the one obtained previously within a different field-theoretical model, see Eq. ( 12) in [97].

I. Pair distribution function at short distances

The last result can be generalized to finite but small r, see [Tab. III, Eqs. (9a,9b)] where the zero-range regime k typ b ≪ 1 was introduced at the end of Sec. III B. Here we justify this for the case where the expectation values g

(2) ↑↓ R + r 2 , R -r 2 = (ψ † ↑ ψ ↑ ) R + r 2 (ψ † ↓ ψ ↓ ) R -r 2
and C = Ĉ are taken in an arbitrary stationary state ψ in the zero-range regime; this implies that the same result holds for a thermal equilibrium state in the zerorange regime, see Section IX. We first note that the expression (28) of g

↑↓ in terms of the wavefunction is valid for the lattice model with the obvious replacement of the integrals by sums, so that

G (2) ↑↓ (r) ≡ Ĝ(2) ↑↓ (r) = R b d N ↑ i=1 N j=N ↑ +1 (r k ) k =i,j b (N -2)d × ψ r 1 , . . . , r i = R + r 2 , . . . , r j = R - r 2 , . . . , r N 2 . ( 64
)
For r ≪ 1/k typ , we can replace ψ by the short-distance expression (20), assuming that the multiple sum is dominated by the configurations where all the distances |r k -R| and r kk ′ are much larger than b and r:

G (2) ↑↓ (r) ≃ (A, A) |φ(r)| 2 . ( 65
)
Expressing (A, A) in terms of C thanks to [Tab. III, Eqs. (4a,4b)] gives the desired [Tab. III, Eqs. (9a,9b)].

J. Momentum distribution at large momenta

Assuming again that we are in the zero-range regime k typ b ≪ 1, we will justify [Tab. III, Eq. ( 10)] both in 3D and in 2D. We start from

n σ (k) = i:σ (r l ) l =i b d(N -1) ri b d e -ik•ri ψ(r 1 , . . . , r N ) 2 .
(66) We are interested in the limit k ≫ k typ . Since ψ(r 1 , . . . , r N ) is a function of r i which varies on the scale of 1/k typ , except when r i is close to another particle r j where it varies on the scale of b, we can replace ψ by its short-distance form (20): (67) where φ

ri b d e -ik•ri ψ(r 1 , . . . , r N ) ≃ φ(k) × j,j =i e -ik•rj A ij (r j , (r l ) l =i,j ),
(k) = k|φ = r b d e -ik•r φ(r).
Here we excluded the configurations where more than two particles are at distances b, which are expected to have a negligible contribution to (66). Inserting ( 67) into (66), expanding the modulus squared, and neglecting the crossproduct terms in the limit k ≫ k typ , we obtain

n σ (k) ≃ | φ(k)| 2 (A, A). (68)
Finally, φ(k) is easily computed for the lattice model: for k = 0, the two-body Schrödinger equation (A1) directly gives φ(k) = -g 0 φ(0)/(2ǫ k ), and φ(0) is given by (15,16), which yields [Tab. III, Eq. ( 10)].

K. Minorization of C by the order parameter (This subsection is supplementary to the published paper) None of the previous relations involve the macroscopic quantum properties of the spin-1/2 Fermi gas, such as superfluidity and off-diagonal long range order. For an arbitrary state |ψ in which the gas is pair-condensed, with a nonzero order parameter ∆(r) of arbitrary position dependence, one obtains an additional relation, in the form of the following minorization:

C ≥ m 2 4 r b d |∆(r)| 2
where d = 2 or 3 is the dimension of space. This inequality is straightforwardly obtained in a U (1) symmetry breaking point of view, where the order parameter is related to the pairing field in the lattice model in 3D [START_REF] Castin | Ultracold Fermi gases[END_REF] and in 2D [G. Tonini, F. Werner, Y. Castin, Eur. Phys. J. D 39, 283 (2006)] by

∆(r) ≡ g 0 ψ ↓ (r)ψ ↑ (r) .
We then split the operator Ô = (ψ ↓ ψ ↑ )(r) as the sum of its expectation value Ô and of fluctuations δ Ô. From the identity Ô † Ô = | Ô | 2 + (δ Ô) † (δ Ô) and the nonnegativeness of the last term in that identity, we obtain

(ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ )(r) ≥ | (ψ ↓ ψ ↑ )(r) | 2 = |∆(r)| 2 g 2 0 .
It remains to sum this inequality over r and to use [Tab. III, Eq. ( 2)] and the expression ( 6) of H int to obtain the announced minorization.

The generalisation to the U (1)-symmetry preserving case is straightforward. When the gas is pair-condensed, the two-body density operator ρ2 , defined by

r 1 , r 2 |ρ 2 |r ′ 1 , r ′ 2 ≡ ψ † ↑ (r ′ 1 )ψ † ↓ (r ′ 2 )ψ ↓ (r 2 )ψ ↑ (r 1 ) ,
has a normalised eigenvector |ϕ 0 with an eigenvalue N0 of order N/2, and this is the only macroscopically populated two-particle mode. N0 is the mean number of condensed pairs and r 1 , r 2 |ϕ 0 = ϕ 0 (r 1 , r 2 ) is the corresponding pair condensate wavefunction. In this framework, the pairing field ψ ↓ (r 2 )ψ ↑ (r 1 ) is replaced by the pair-condensed field N 1/2 0 ϕ 0 (r 1 , r 2 ) so that the order parameter is replaced by

∆(r) = g 0 N 1/2 0 ϕ 0 (r, r) We then introduce the splitting ρ2 = N0 |ϕ 0 ϕ 0 | + δ ρ2
where both ρ2 and δ ρ2 are hermitian nonnegative, hence the chain

(ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ )(r) = r, r|ρ 2 |r, r ≥ N0 |ϕ 0 (r, r)| 2 = |∆(r)| 2 g 2 0
The summation over r as in the symmetry-breaking case leads to the announced minorization.

Our minorization extends to the continuous-space limit b → 0 where, in particular, our definition of the order parameter in 3D reconnects to the one (28) of reference [5], as can be shown from the normalisation condition (9) of φ(r) and from its value (15) at r = 0.

VI. RELATIONS FOR A FINITE-RANGE INTERACTION IN CONTINUOUS SPACE

In this Section VI, we restrict for simplicity to the case of a stationary state. It is then convenient to define C by [Tab. IV, Eqs. (1a,1b)].

A. Interaction energy

As for the lattice model, we find that the interaction energy is proportional to C, see [Tab. IV, Eqs. (2a,2b)]. It was shown in [99] that the 3D relation is asymptotically valid in the zero-range limit. Here we show that it remains exact for any finite value of the range and we generalize it to 2D.

For the derivation, we set

V (r) = g 0 W (r) (69)
where g 0 is a dimensionless coupling constant which allows to tune a. The Hellmann-Feynman theorem then gives

E int = g 0 dE/dg 0 .
g 0 d(-1/a) dg 0 = m 4π 2 d 3 r V (r)|φ(r)| 2 in 3D (70) g 0 d(ln a) dg 0 = m 2π 2 d 2 r V (r)|φ(r)| 2 in 2D. (71)
To derive these lemmas, we consider two values of the scattering length a i , i = 1, 2, and the corresponding scattering states φ i and coupling constants g 0,i . The corresponding two-particle relative-motion Hamiltonians are

H i = -( 2 /m) ∆ r + g 0,i W (r). Since H i φ i = 0, we have lim R→∞ r<R d d r (φ 1 H 2 φ 2 -φ 2 H 1 φ 1 ) = 0. ( 72
)
The contribution of the kinetic energies can be computed from the divergence theorem and the large-distance form of φ [229]. The contribution of the potential energies is proportional to g 0,2 -g 0,1 . Taking the limit a 2 → a 1 gives the results (70,71). Lemma (70) was also used in [99] and the above derivation is essentially identical to the one of [99]. For this 3D lemma, there also exists an alternative derivation based on the two-body problem in a large box [230].

Three dimensions Two dimensions

C ≡ 4πm 2 dE d(-1/a) (1a) C ≡ 2πm 2 dE d(ln a) (1b) Eint = C (4π) 2 d 3 r V (r)|φ(r)| 2 (2a) Eint = C (2π) 2 d 2 r V (r)|φ(r)| 2 (2b) E -Etrap = 2 C 4πma E -Etrap = lim R→∞ 2 C 2πm ln R a + σ d 3 k (2π) 3 2 k 2 2m nσ(k) - C (4π) 2 | φ′ (k)| 2 (3a) + σ d 2 k (2π) 2 2 k 2 2m nσ(k) - C (2π) 2 | φ′ R (k)| 2 (3b)
In the zero-range regime ktypb ≪ 1

d 3 R g (2) ↑↓ R + r 2 , R - r 2 ≃ C (4π) 2 |φ(r)| 2 for r ≪ k -1 typ (4a) d 2 R g (2) ↑↓ R + r 2 , R - r 2 ≃ C (2π) 2 |φ(r)| 2 for r ≪ k -1 typ (4b) nσ(k) ≃ C (4π) 2 | φ(k)| 2 for k ≫ ktyp (5a) nσ(k) ≃ C (2π) 2 | φ(k)| 2 for k ≫ ktyp (5b)
TABLE IV: Relations for spin-1/2 fermions with a finite-range interaction potential V (r) in continuous space, for any stationary state. C is defined in line 1. All relations remain valid at thermal equilibrium in the canonical ensemble; the derivatives of the energy in line 1 then have to be taken at constant entropy. Equations (1a,2a,4a) are contained in [99] (for ktypb ≪ 1). The functions φ ′ (r) and φ ′ R (r) are given by Eqs. (73,78) and φ′ (k), φ′ R (k) are their Fourier transforms.

B. Relation between energy and momentum distribution

Three dimensions: The natural counterpart, for a finiterange interaction potential, of the zero-range-model expression of the energy as a functional of the momentum distribution [Tab. II, Eqs. (5a)] is given by [Tab. IV, Eq. (3a)], where φ′ (k) is the zero-energy scattering state in momentum space with the incident wave contribution

∝ δ(k) subtracted out: φ′ (k) = φ(k) + a -1 (2π) 3 δ(k) = d 3 r e -ik•r φ ′ (r) with φ ′ (r) = φ(r) + 1 a . ( 73
)
This is simply obtained by adding the kinetic energy to [Tab. IV, Eq. (2a)] and by using the lemma:

d 3 r V (r)|φ(r)| 2 = 4π 2 ma - d 3 k (2π) 3 2 k 2 m | φ′ (k)| 2 . ( 74
)
To derive this lemma, we start from Schrödinger's equation -( 2 /m)∆φ + V (r)φ = 0, which implies

d 3 r V (r)|φ(r)| 2 = 2 m d 3 r φ∆φ. (75)
Applying the divergence theorem over the sphere of radius R, using the asymptotic expression (9) of φ and taking the limit R → ∞ then yields

d 3 r φ∆φ = 4π a -d 3 r (∇φ) 2 . ( 76
)
We then replace ∇φ by ∇φ ′ . Applying the Parseval-Plancherel relation to ∂ i φ ′ , and using the fact that φ ′ (r) vanishes at infinity, we get:

d 3 r (∇φ ′ ) 2 = d 3 k (2π) 3 k 2 | φ′ (k)| 2 (77)
The desired result (74) follows.

Two dimensions: An additional regularisation procedure for small momenta is required in 2D, as was the case for the zero-range model [Tab. II, Eq. (5b)] and for the lattice model [Tab. III, Eq. (3b)]. One obtains [Tab. IV, Eq. (3b)], where φ′

R (k) = d 2 r e -ik•r φ ′ R (r) with φ ′ R (r) = [φ(r) -ln(R/a)] θ(R -r). (78)
This follows from [Tab. IV, Eq. (2b)] and from the lemma:

d 2 r V (r)|φ(r)| 2 = lim R→∞ 2π 2 m ln R a - d 2 k (2π) 2 2 k 2 m | φ′ R (k)| 2 . ( 79
)
The derivation of this lemma again starts with the 2D version of (75). The divergence theorem then gives [229]

d 2 r φ∆φ = lim R→∞ 2π ln R a - r<R d 2 r (∇φ) 2 . ( 80
) We can then replace r<R d 2 r (∇φ) 2 by d 2 r (∇φ ′ R ) 2 , since φ ′ R (r) is continuous at r = R [229] so that ∇φ ′ R
does not contain any delta distribution. The Parseval-Plancherel relation can be applied to ∂ i φ ′ R , since this function is square-integrable. Then, using the fact that φ ′ R (r) vanishes at infinity, we get

d 2 r (∇φ ′ R ) 2 = d 2 k (2π) 2 k 2 | φ′ R (k)| 2 , (81)
and the lemma (79) follows.

C. Pair distribution function at short distances

In the zero-range regime k typ b ≪ 1, the short-distance behavior of the pair distribution function is given by the same expressions [Tab. III, Eqs. (9a,9b)] as for the lattice model. Indeed, Eq. ( 65) is derived in the same way as for the lattice model; one can then use the zero-range model's expressions [Tab. II, Eqs. (2a,2b)] of (A, A) in terms of C, since the finite range model's quantities C and A tend to the zero-range model's ones in the zerorange limit. In 3D, the result [Tab. III, Eq. (9a)] is contained in [99].

D. Momentum distribution at large momenta

In the zero-range regime k typ b ≪ 1 the momentum distribution at large momenta k ≫ k typ is given by

n σ (k) ≃ C (4π) 2 | φ(k)| 2 in 3D (82) n σ (k) ≃ C (2π) 2 | φ(k)| 2 in 2D. ( 83
)
Indeed, Eq. ( 68) is derived as for the lattice model, and (A, A) can be expressed in terms of C as in the previous subsection VI C.

VII. DERIVATIVE OF THE ENERGY WITH RESPECT TO THE EFFECTIVE RANGE

Assuming that the zero-range model is solved, we first show that the first correction to the energy due to a finite range of the interaction potential V (r) can be explicitly obtained and only depends on the s-wave effective range of the interaction. We then enrich the discussion using the many-body diagrammatic point of view, where the central object is the full two-body T -matrix, to recall in particular that the situation is more subtle for lattice models [START_REF] Burovski | [END_REF]. Finally, we relate ∂E/∂r e to a subleading term of the short distance behavior of the pair distribution function in Sec.VII D and to the coefficient of the 1/k 6 subleading tail of n σ (k) in Sec. VII E.

A. Derivation of the explicit formulas

Three dimensions:

In 3D, the leading order finite-range correction to the zero-range model's spectrum depends on the interaction potential V (r) only via its effective range r e , and is given by the expression [Tab. V, Eq. (1a)], where the derivative is taken in r e = 0 for a fixed value of the scattering length, the function A is assumed to be real without loss of generality. As a first way to obtain this result we use a modified version of the zero-range model, where the boundary condition [Tab. I, Eq. (1a)] is replaced by (84) where

ψ(r 1 , . . . , r N ) = rij →0 1 r ij - 1 a + m 2 2 Er e × A ij (R ij , (r k ) k =i,j ) + O(r ij ),
E = E-2U (R ij )-   k =i,j U (r k )   + 1 A ij (R ij , (r k ) k =i,j ) ×   2 4m ∆ Rij + 2 2m k =i,j ∆ r k   A ij (R ij , (r k ) k =i,j ) . (85)
Equations (84,85) generalize the ones already used for 3 bosons in free space in [122,123] (the predictions of [122] and [123] have been confirmed using different approaches, see [124] and Refs. therein, and [125,126] respectively; moreover, a derivation of these equations was given in [122]). Such a model was also used in the two-body case, see e.g. [127][128][129], and the modified scalar product that makes it hermitian was constructed in [130].

For the derivation of [Tab. V, Eq. (1a)], we consider a stationary state ψ 1 of the zero-range model, satisfying the boundary condition [Tab. I, Eq. (1a)] with a scattering length a and a regular part A (1) , and the corresponding finite-range stationary state ψ 2 satisfying (84,85) with the same scattering length a and a regular part A (2) . As in Appendix B we get (B3), as well as (B6) with 1/a 1 -1/a 2 replaced by mEr e /(2 2 ). This yields [Tab. V, Eq. (1a)].

A deeper physical understanding and a more selfcontained derivation may be achieved by going back to the actual finite range model V (r; b) for the interaction potential, such that the scattering length remains fixed when the range b tends to zero. The Hellmann-Feynman theorem gives

dE db = N ↑ i=1 N j=N ↑ +1 d 3 r 1 . . . d 3 r N |ψ(r 1 , . . . , r N )| 2 ∂ b V (r ij ; b). (86)
We need to evaluate |ψ| 2 for a typical configuration with two atoms i and j within the potential range b; in the limit b → 0 one may then assume that the other atoms are separated by much more than b and are at distances from R ij = (r i + r j )/2 much larger than b. This motivates the factorized ansatz

ψ(r 1 , . . . , r N ) ≃ χ(r ij )A ij (R ij , (r k ) k =i,j ). ( 87
)
We take a rotationally invariant χ, because we assume the absence of scattering resonance in the partial waves other than s-wave [231]: The p-wave scattering amplitude, that vanishes quadratically with the relative wavenumber k, is then O(b 3 k 2 ), resulting in an energy contribution O(b 3 ) negligible at the present order.

Three dimensions Two dimensions ∂E ∂re a = 2π(A, (E -H)A) (1a) ∂E ∂(r 2 e ) a = π(A, (E -H)A) (1b) Hij ≡ - 2 4m ∆R ij - 2 2m k =i,j ∆r k + 2U (Rij ) + k =i,j U (r k ) (2) Ḡ(2) ↑↓ (r) = r→0 C (4π) 2 1 r - 1 a 2 - m 2π 2 ∂E ∂re + O(r) (3a) Ḡ(2) ↑↓ (r) = r→0 C (2π) 2 ln 2 (r/a) - m 2π 2 ∂E ∂(r 2 e ) r 2 ln 2 r + O(r 2 ln r) (3b) nσ(k) - C k 4 ∼ k→∞ 1 k 6 16πm 2 ∂E ∂re -8π 2 (A, ∆RA) (4a) nσ(k) - C k 4 ∼ k→∞ 1 k 6 8πm 2 ∂E ∂(r 2 e ) -4π 2 (A, ∆RA) (4b)
TABLE V: For spin-1/2 fermions, derivative of the energy with respect to the effective range re, or to its square in 2D, taken at re = 0 for a fixed value of the scattering length. The functions A (assumed to be real) are the ones of the zero-range regime.

The compact notations for the scalar products and the matrix elements are defined in Tab. I. nσ(k) is the average of nσ(k) over the direction of k. Ḡ( 2) ↑↓ (r) is the pair distribution function integrated over the center of mass of the pair and averaged over the direction of r.

Inserting the ansatz (87) into Schrödinger's equation Hψ = Eψ, and neglecting the trapping potential within the interaction range r ij ≤ b, as justified in the Appendix D, gives [232]

Eχ(r ij ) ≃ [- 2 m ∆ rij + V (r ij ; b)]χ(r ij ), ( 88
)
where E is given by (85). For E > 0, we set E = 2 k 2 /m with k > 0, and χ is a finite energy scattering state; to match the normalization of the zero energy scattering state φ used in this article, see [START_REF] Lee | [END_REF], we take for r out of the interaction potential

χ(r) = r→∞ 1 f k sin(kr) kr + e ikr r , (89)
where f k is the scattering amplitude. The optical theorem, implying that

f k = - 1 ik + u(k) , (90)
where u(k) ∈ R, ensures that χ is real [233].

Inserting the ansatz (87) into the Hellmann-Feynman expression (86) gives

dE db ≃ i<j d 3 R ij ( k =i,j d 3 r k ) A 2 ij (R ij , (r k ) k =i,j ) × d 3 r ij χ 2 (r ij )∂ b V (r ij ; b) (91)
To evaluate the integral of χ 2 ∂ b V , we use the following lemma (whose derivation is given in the next paragraph):

4π 2 m [u 2 (k)-u 1 (k)] = R 3 d 3 r χ 1 (r)χ 2 (r)[V (r; b 1 )-V (r; b 2 )] (92
u(k) = 1 a - 1 2 r e k 2 + O(b 3 k 4 ) (93)
r e being the effective range of the interaction potential of range b, we obtain [Tab. V, Eq. (1a)] [234].

As a side result of this physical approach, the modified contact conditions (84) may be rederived. One performs an analytical continuation of the out-of-potential wavefunction (89) to the interval r ≤ b [103] and one takes the zero-r limit of that continuation [235]. In simple words, this amounts to expanding (89) in powers of r:

χ(r) = 1 r - 1 a + 1 2 k 2 r e + O(r). (94)
Inserting this expansion in (87) and using k 2 = mE/ 2 gives (84). The lemma ( 92) is obtained by multiplying Schrödinger's equations for χ 1 (respectively for χ 2 ) by χ 2 (respectively by χ 1 ), taking the difference of the two resulting equations, integrating this difference over the sphere r < R and using the divergence theorem to convert the volume integral of χ 2 ∆ r χ 1χ 1 ∆ r χ 2 into a surface integral, where the asymptotic forms (89) for r = R → +∞ may be used. When E < 0, we set E = -2 κ 2 /m with κ > 0 and we perform analytic continuation of the E > 0 case by replacing k with iκ. From (89) it appears that χ(r) now diverges exponentially at large distances, as e κr /r, if 1/f (iκ) = 0. If the interaction potential is a compact support potential, or simply tends to zero more rapidly than exp(-2κr), the lemma and the final conclusion [Tab. V, Eq. (1a)] still hold; the functions u 1 (iκ) and u 2 (iκ) remain real, since the series expansion of u(k) has only even powers of k.

Two dimensions:

The above physical reasoning may be directly generalized to 2D [236], giving [Tab. V, Eq. (1b)], where the derivative is taken for a fixed scattering length in r e = 0. The main difference with the 3D case [Tab. V, Eq. (1a)] is that the energy E now varies quadratically with the effective range r e , as already observed numerically for three-boson-bound states in [131]. In the derivation, the first significant difference with the 3D case occurs in the normalization of the two-body scattering state: (89) is replaced with

χ(r) = r→∞ π 2i 1 f k J 0 (kr) + H (1) 0 (kr) (95)
where H

(1) 0 = J 0 + iN 0 is a Hankel function, J 0 and N 0 are Bessel functions of the first and second kinds. The optical theorem implies |f k | 2 + Re f k = 0 so that

f k = -1 1 + iu(k) with u(k) ∈ R, ( 96
)
and χ is real. The low-k expansion for a potential of range b takes the form [132,133]

u(k) = 2 π ln (e γ ka/2) + 1 2 (kr e ) 2 + . . . , (97)
where γ = 0.577216 . . . is Euler's constant, the logarithmic term being obtained in the zero-range Bethe-Peierls model and the k 2 term corresponding to finite effective range corrections (with the sign convention of [132] such that r 2 e > 0 for a hard disk potential). The subsequent calculations are similar to the 3D case, also for the negative energy case where analytic continuation gives rises to the special functions I 0 (κr) and K 0 (κr). For example, at positive energy, the lemma (92) takes in 2D the form

π 2 2 m [u 1 (k)-u 2 (k)] = R 2 d 2 r χ 1 (r)χ 2 (r)[V (r; b 1 )-V (r; b 2 )] (98)
The fact that one can neglect the trapping potential within the interaction range is again justified in Appendix D. Finally, we note that the expansion of the asymptotic form (95) for r → 0, and for k → 0,

χ(r) = ln(r/a) - 1 2 (kr e ) 2 + O(r 2 ln r) (99)
allows to determine the 2D version of the modified zerorange model (84),

ψ(r 1 , . . . , r N ) = rij →0 ln(r ij /a) - m 2 2 Er 2 e × A ij (R ij , (r k ) k =i,j ) + O(r ij ) ( 100
)
where E is defined as in 3D by (85). To complete this 2D derivation, one has to check that the p-wave interaction brings a negligible contribution to the energy. The p-wave scattering amplitude at low relative wavenumber k vanishes as k 2 R 2 1 where R 2 1 is the p-wave scattering surface [134]. One could believe that r e ≈ R 1 ≈ b, one would then conclude that the p-wave contribution to the energy, scaling as R 2 1 , cannot be neglected as compared to the s-wave finite range correction, scaling as r 2 e . Fortunately, as shown in subsection VII B, this expectation is too naive, and [Tab. V, Eq. (1b)] is saved by a logarithm, r e being larger than R 1 by a factor ln(a/b) ≫ 1 in the zero range limit [237].

B. What we learn from diagrammatic formalism

In the many-body diagrammatic formalism [START_REF] Landau | Statistical Physics part[END_REF][START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF], the equation of state of the homogeneous gas (in the thermodynamic limit) is accessed from the single particle Green's function, which can be expanded in powers of the interaction potential, each term of the expansion being represented by a Feynman diagram. The internal momenta of the diagrams can however be as large as /b, where b is the interaction range. A standard approach to improve the convergence of the perturbative series for strong interaction potentials is to perform the so-called ladder resummation. The resulting Feynman diagrams then involve the two-body T -matrix of the interaction, rather than the bare interaction potential V . For the spin-1/2 Fermi gas, where there is a priori no Efimov effect, one then expects that the internal momenta of the Feynman diagrams are on the order of k typ only, where the typical wavenumber k typ was defined in subsection III B. As put forward in [START_REF] Burovski | [END_REF], the interaction parameters controlling the first deviation of the gas energy from its zero-range limit are then the ones appearing in the first deviations of the two-body T -matrix element k 1 , k 2 |T (E +i0 + )|k 3 , k 4 from its zero-range limit, where all the k i are on the order of k typ and E is on the order of 2 k 2 typ /m. The single particle Green's function is indeed a sum of integrals of products of T -matrix elements and of ideal-gas Green's functions.

We explore this idea in this subsection. For an interaction potential V (r), we confirm the results of subsection VII A. In addition to the effective range r e characterizing the on-shell T -matrix elements (that is the scattering amplitude), the diagrammatic point of view introduces a length ρ e characterizing the s-wave low-energy off-shell T -matrix elements, and a length R 1 characterizing the p-wave on-shell scattering; we will show that the contributions of ρ e and R 1 are negligible as compared to the one of the effective range r e . Moreover, in the case of lattice models, a length R e characterizing the breaking of the Galilean invariance appears [START_REF] Burovski | [END_REF]. Its contribution is in general of the same order as the one of r e . Both contributions can be zeroed for appropriately tuned matterwave dispersion relations on the lattice. Finally, in the case of a continuous space model with a delta interaction potential plus a spherical cut-off in momentum space, and in the case of a lattice model with a spherical momentum cut-off, we show that the breaking of Galilean invariance does not disappear in the infinite cut-off limit.

1. For the continuous space interaction V (r) When each pair of particles i and j interact in continuous space via the potential V (r ij ), one can use Galilean invariance to restrict the T -matrix to the center of mass frame, where k ′ ≡ k 1 = -k 2 and k ≡ k 3 = -k 4 . Further using rotational invariance, one can restrict this internal T -matrix to fixed total angular momentum l, with matrix elements characterized by the function t l (k ′ , k; E) whose low-energy behavior was extensively studied [134,[START_REF] Gibson | [END_REF]. This function is said to be on-shell iff k = k ′ = (mE) 1/2 / , in which case it is simply noted as t l (E), otherwise it is said to be off-shell.

Three dimensions:

We assume that the interaction potential, of compact support of range b, is everywhere non-positive (or infinite). We recall that we are here in the resonant regime, with a s wave scattering length a such that |a| ≫ b. The potential is assumed to have the minimal depth leading to the desired value of a, so as to exclude deeply bound dimers. In particular, at resonance (1/a = 0), there is no two-body bound state. To invalidate the usual variational argument [START_REF] Werner | [END_REF][138][START_REF] Blatt | Theoretical Nuclear Physics[END_REF][START_REF] Fregoso | [END_REF] (that shows, for a nonpositive interaction potential, that the spin-1/2 fermions have deep N -body bound states in the large N limit), we allow that V (r) has a hard core of range b hard < b. We directly restrict to the s-wave case (l = 0), since the non-resonant p-wave interaction bring a negligible O(b 3 ) contribution, as already discussed in subsection VII A.

The first deviation of the on-shell s-wave T -matrix from its zero-range limit is characterized by the effective range r e , already introduced in Eq. ( 93). The effective range is given by the well-known Smorodinski formula [133]:

1 2 r e = +∞ 0 dr (1 -r/a) 2 -u 2 0 (r) (101)
in terms of the zero-energy scattering state φ(r), with u 0 (r) = rφ(r) and φ is normalized as in Eq. ( 9). Note that u 0 (r) is zero for r ≤ b hard . As r e deviates from its resonant (|a| → ∞) value by terms O(b 2 /a), the discussion of its 1/a = 0 value is sufficient here. The function

u 0 then solves 0 = - 2 m u ′′ 0 (r) + V (r)u 0 (r) (102)
with the boundary conditions u 0 (b hard ) = 0 and u 0 (r) = 1 for r > b. Due to the absence of two-body bound states, u 0 is the ground two-body state and it has a constant sign, u 0 (r) ≥ 0 for all r. Since V ≤ 0, Eq. ( 102) implies that u ′′ 0 ≤ 0, the function u 0 is concave. Combined with the boundary conditions, this leads to 0 ≤ u 0 (r) ≤ 1, for all r. Then from Eq. ( 101):

2b hard ≤ r e ≤ 2b (103)
For the considered model, this proves that k typ r e → 0 in the zero-range limit b → 0, which is a key property for the present work. Note that the absence of two-body bound states at resonance is the crucial hypothesis ensuring that r e ≥ 0; it was not explicitly stated in the solution of problem 1 in Sec. 131 of [START_REF] Landau | Quantenmechanik[END_REF]. Without this hypothesis, r e at resonance can be arbitrarily large and negative even for V (r) ≤ 0 for all r, see an explicit example in [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF].

In the s-wave channel, the first deviations of the offshell T -matrix from its zero-range value introduces, in addition to r e , another length that we call ρ e , such that [START_REF] Gibson | [END_REF] [238]

t 0 (k, k ′ ; E) t 0 (E) -1 ∼ k,k ′ ,E→0 2mE 2 -k 2 -k ′2 1 2 ρ 2 e with 1 2 ρ 2 e = +∞ 0 dr r[(1 -r/a) -u 0 (r)]. ( 104
)
For our minimal-depth model at resonance, we conclude that 0 ≤ ρ 2 e ≤ b 2 , so it appears, in the finite-range correction to the energy, at a higher order than r e , and it cannot contribute to [Tab. V, Eq. (1a)].

Two dimensions:

The specific feature of the 2D case is that the minimaldepth attractive potential ensuring the desired scattering length a only weakly dephases the matter-wave over its range, when ln(a/b) ≫ 1. This is apparent e.g. if V (r) is a square-well potential of range b, V (r) = -

2 k 2 0 m θ(b -r): One has -k 0 bJ ′ 0 (k 0 b)/J 0 (k 0 b) = 1/ ln(a/b)
, where J 0 is a Bessel function, which shows that, for the minimaldepth solution, the matter-wave phase shift k 0 b vanishes as [2/ ln(a/b)] 1/2 in the zero-range limit. This property allows to treat the potential perturbatively.

There are three relevant parameters describing the lowenergy behavior of the T -matrix beyond the zero-range limit. The first one is the effective range r e for the swave on-shell T -matrix, see Eq. (97). It is given by the bidimensional Smorodinski formula [132,133]:

1 2 r 2 e = +∞ 0 dr r[ln 2 (r/a) -φ 2 (r)] (105)
where the zero-energy scattering state φ(r) is normalized as in Eq. (10). The second parameter is the length ρ e associated to the s-wave off-shell T -matrix: The 2D equivalent of Eq. ( 104) is [134]:

t 0 (k, k ′ ; E) t 0 (E) -1 ∼ k,k ′ ,E→0 2mE 2 -k 2 -k ′2 1 2 ρ 2 e with 1 2 ρ 2 e = +∞ 0 dr r[φ(r) -ln(r/a)]. ( 106
)
The third parameter is the length R 1 characterizing the low-energy p-wave scattering. For the l-wave scattering state of energy E = 2 k 2 /m, k > 0, we generalize Eq. ( 95) as

χ (l) (r) = r→∞ π 2i k l 1 f (l) k J l (kr) + H (1)
l (kr) .

The l-wave scattering amplitude then vanishes as

f (l) k ∼ k→0 i π 2 k 2l R 2l l ( 108
)
and the leading behavior of the off-shell l-wave T -matrix is characterized by the same length R l as the on-shell one [134].

The situation thus looks critical in 2D: Three lengths squared characterize the low-energy T -matrix, one may naively expect that they are of the same order ≈ b 2 and that they all three contribute to the finite-range correction to the gas energy at the same level, whereas [Tab. V, Eq. (1b)] singles out the effective range r e . By a perturbative treatment of the minimal-depth finite-range potential V (r) of fixed scattering length a, we however obtain in the zero-range limit the following hierarchy, see Appendix E:

r 2 e ∼ b→0 2ρ 2
e ln(a/b) (109)

ρ 2 e = b→0 1 2 R 2 d 2 r r 2 V (r) R 2 d 2 r V (r) 1 + O 1 ln(a/b) (110) R 2 1 ∼ b→0 ρ 2 e 2 ln(a/b) (111)
This validates [Tab. V, Eq. (1b)] when ln(a/b) ≫ 1.

Lattice models

We restrict here for simplicity to the 3D case. To obtain a non-zero T -matrix element k 1 , k 2 |T (E + i0 + )|k 3 , k 4 , due to the conservation of the total quasimomentum, we have to restrict to k 1 + k 2 = k 3 + k 4 ≡ K (modulo a vector of the reciprocal lattice). As the interactions in the lattice model are purely on-site, the matrix element only depends on the total quasi-momentum K and the energy E, and is noted as t(K, E) in what follows. We recall that the bare coupling constant g 0 is adjusted to have a fixed scattering length a on the lattice, see Eq. (11), which leads to

g 0 = 4π 2 a/m 1 -K 3 a/b (112)
where the numerical constant K 3 depends on the lattice dispersion relation ǫ k . One then gets [START_REF] Burovski | [END_REF] 1

t(K, E) = m 4π 2 a - D d 3 q (2π) 3 1 2ǫ q + 1 E + i0 + -ǫ 1 2 K+q -ǫ 1 2 K-q ( 113
)
where a is the s-wave scattering length and the dispersion relation ǫ q is extended by periodicity from the first Brillouin zone D to the whole space. The low-K and lowenergy limit of that expression was worked out in [START_REF] Burovski | [END_REF], it involves the effective range r e and an extra length R e quantifying the breaking of Galilean invariance:

1 t(K, E) = m 4π 2 1 a + ik - 1 2 r e k 2 - 1 2 R e K 2 + . . . (114
) where the relative wavenumber k such that E -

2 K 2 4m = 2 k 2
m is either real non-negative or purely imaginary with a positive imaginary part. The two lengths are given by

r e = R 3 \D d 3 q π 2 q 4 + D d 3 q π 2 1 q 4 - 2 2mǫ q 2 (115) R e = - • D d 3 q 4π 2 2 2mǫ q 2 1 - m 2 ∂ 2 ǫ q ∂q 2 x - π b -π b π b -π b dq y dq z 8π 2 2 mǫ 2 ( π b ,qy,qz ) ∂ǫ ( π b ,qy,qz) ∂q x ( 116
)
where the dispersion relation ǫ k was supposed to be twice differentiable on the interior • D of the first Brillouin zone and to be invariant under permutation of the coordinate axes. As compared to [START_REF] Burovski | [END_REF] we have added the second term (a surface term) in Eq. ( 116) to include the case where the dispersion relation has cusps at the border of the first Brillouin zone [239]. As mentioned in the introduction of the present section, we then expect that, in the lattice model, the first deviation of any manybody eigenenergy E from the zero-range limit is a linear function of the two parameters r e and R e with modelindependent coefficients:

E(b) = b→0 E(0) + ∂E ∂r e r e + ∂E ∂R e R e + o(b) (117)
This feature was overlooked in the early version [88] of this work. It invalidates the discussion of ∂T c /∂r e given in [88].

We illustrate this discussion with a few relevant examples. For a parabolic dispersion relation ǫ k = 2 k 2 2m , the constant K 3 = 2.442 749 607 806 335 . . . [15,143] and the effective range [START_REF] Castin | Ultracold Fermi gases[END_REF][START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF] were already calculated, first numerically then analytically; in the quantity R e , the first term vanishes but there is still breaking of Galilean invariance due to the non-zero surface term that can be deduced from Eq. (F6):

r e = b 12 √ 2 π 3 arcsin 1 √ 3 ≃ 0.337b and R e = - 1 12 r e (118)
In an attempt to reduce the dependence of the Monte Carlo results on the grid spacing b, a zero-effective-range dispersion relation was constructed [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF]144],

ǫ k = 2 k 2 2m [1 -C(kb/π) 2 ], (120)
with C ≃ 0.257 022, and used in real simulations [144].

The corresponding K 3 ≃ 2.899 952. Unfortunately this leads to a sizeable R e : R e ≃ -0.168b.

As envisioned in [START_REF] Burovski | [END_REF] one may look for dispersion relations with r e = R e = 0. We have found an example of such a magic dispersion relation:

ǫ k = ǫ Hub k [1 + αX + βX 2 ] with X = ǫ Hub k 6 2 /

The single-particle momentum cut-off model

A continuous space model used in particular in [57] takes a Dirac delta interaction potential g 0 δ(r ir j ) between particles i and j, and regularizes the theory by introducing a cut-off Λ on all the single-particle wavevectors. Due to the conservation of momentum one needs to evaluate the T -matrix only between states with the same total momentum K. Due to the contact interaction the resulting matrix element depends only on K and on E, and is noted as t(K, E). Expressing g 0 in terms of the s-wave scattering length as in [57] one gets

1 t(K, E) = m 4π 2 a - R 3 d 3 q (2π) 3 θ(Λ -q) 2ǫ q + θ(Λ -| 1 2 K + q|) θ(Λ -| 1 2 K -q|) E + i0 + -ǫ 1 2 K+q -ǫ 1 2 K-q (125)
where ǫ q = 2 q 2 /(2m) for all q. Introducing the relative wavenumber k such that E -

2 K 2 4m = 2 k 2 m , k ∈ R + or k ∈ iR + , we obtain the low wavenumbers expansion 1 t(K, E) = m 4π 2 1 a + ik - K 2π - 1 2 r e k 2 - 1 2 R e K 2 +. . . (126)
The effective range is given by r e = 4/(πΛ) and the length R e = r e /12 [240]. The unfortunate feature of this model is the occurrence of a term linear in K, that does not disappear even if Λ → +∞: The model thus does not reproduce the universal zero-range model in the large cut-off limit, as soon as pairs of particles have a non-zero total momentum. Note that here one cannot exchange the order of the integration over q and the Λ → ∞ limit. As a concrete illustration of the breaking of the Galilean invariance, for a > 0 and in the limit Λ → +∞, it is found (e.g. by calculating the pole of the T -matrix) that the total energy of a free-space dimer of total momentum K is

E model dim (K) = 2 K 2 4m - 2 m 1 a - K 2π 2 (127)
and that this dimer state exists only for K < 2π/a [241].

The single-particle momentum cut-off lattice model

A spherical momentum cut-off was also introduced for a lattice model in [53,56,146,147]. Our understanding is that this amounts to taking the following dispersion relation inside the first Brillouin zone: ǫ k = 2 k 2 /(2m) for k < π/b, ǫ k = +∞ otherwise. The T -matrix is then given by Eq. ( 113), where for K = 0 one extends ǫ k by periodicity out of the first Brillouin zone. By distinguishing three zones within the integration domain for q, similarly to the note [239], and restricting for simplicity to E = 2 K 2 /(4m), we find the same undesired term -K/(2π) as in Eq. ( 126), implying that the model does not reproduce the unitary gas even for b → 0.

C. The Juillet effect for lattice models

With the lattice dispersion relation ǫ k of (120), adjusted to have a zero effective range r e = 0, Olivier Juillet numerically observed, for two particles in the cubic box [0, L] 3 with periodic boundary conditions and zero total momentum, that the first energy correction to the zero-range limit b → 0 is linear in b [145], which seems to contradict [Tab. V, Eq. (1a)]. This is illustrated in Fig. 1. This cannot be explained by a non-zero R e [defined in Eq. ( 116)] because the two opposite-spin fermions have here a zero total momentum. This Juillet effect, as we shall see, is due to the fact that the integral of 1/ǫ k over k in the first Brillouin zone and the corresponding discrete sum for the finite size quantization box differ for b/L → 0 not only by a constant term but also by a term linear in b, when the dispersion relation has a cusp at the surface of the first Brillouin zone, such as Eq. ( 120). The Juillet effect thus disappears in the thermodynamic limit. This explains why it does not show up in the diagrammatic point of view of Sec. VII B, which was considered in the thermodynamic limit, so that only momentum integrals appeared. This also shows that the Juillet effect does not invalidate [Tab. V, Eq. (1a)] since it was derived for an interaction that is smooth in momentum space.

In [143] it was shown that the lattice model spectrally reproduces the zero-range model when the grid spacing b → 0. We now simply extend the reasoning of [143] for two particles to first order in b included. For an eigenenergy E which does not belong to the non-interacting spec- 120) with zero effective range re = 0 (in blue, lower set), and the magic one (122) with re = Re = 0 with the parameters of Eq. ( 123) (in black, upper set) and of Eq. ( 124) (in red, middle set). The fact that the energy varies linearly in b for the re = 0 quartic dispersion relation at zero total momentum is the Juillet effect explained in Sec. VII C, and the corresponding dashed line is the analytical result (134). At non-zero total momentum the quartic dispersion relation leads to an energy variation linear in b as expected e.g. from the fact that its has a non-zero Re [the dotted line is a linear fit for b/L ≤ 0.01]. The magic dispersion relations lead to a O(b 2 ) variation of the energy both at zero and non-zero total momentum [the dotted lines are purely quadratic fits performed for b/L ≤ 0.02].

trum, the exact implicit equation is

1 g 0 + 1 L 3 k∈D 1 2ǫ k -E = 0 (128)
where the notation with a discrete sum over k implicitly restricts k to 2π L Z 3 . By adding and subtracting terms, and using the expressions ( 11) and ( 115) for the bare coupling constant g 0 and the effective range r e , one obtains the useful form:

1 g - m 2 Er e 8π 4 + 1 L 3 - 1 E + k∈D * F (ǫ k )+ k∈R 3 * E ( 2 k 2 /m) 2 = R 1 + ER 2 -ER 3 (129) with g = 4π 2 a/m and F (ǫ) = (2ǫ -E) -1 -(2ǫ) -1 - E/(2ǫ) 2 . We have defined R 1 ≡ D d 3 k (2π) 3 1 2ǫ k - 1 L 3 k∈D * 1 2ǫ k , (130)
proportional to the function C(b) introduced in [143]. The quantities R 2 and R 3 have the same structure: R 2 is obtained by replacing in

R 1 the function 1/(2ǫ k ) by 1/(2ǫ k ) 2 -1/( 2 k 2 /m) 2
, in the integral and in the sum; R 3 is obtained by replacing in R 1 the function 1/(2ǫ k ) by 1/( 2 k 2 /m) 2 and the set D by R 3 \ D, both for the integration and for the summation. We now take b → 0 in Eq. ( 129), keeping terms up to O(b) included. Since F (ǫ) = O(1/ǫ 3 ) at large ǫ, we can replace F (ǫ k ) by its b → 0 limit F ( 2 k 2 /2m), and the summation set D * by its b → 0 limit [242]:

k∈D * F (ǫ k ) = b→0 k∈R 3 * F 2 k 2 2m + O(b 2 ) (131)
In the quantities R i , we perform the change of variables k = 2πq/b, and we write the dispersion relation as

ǫ k = (2π ) 2 mb 2 η kb/(2π) (132)
where the dimensionless η q does not depend on the lattice spacing b. We then find that bR 1 , R 2 /b and R 3 /b are differences between a converging integral and a threedimensional Riemann sum with a vanishing cell volume (b/L) 3 . As these differences vanish as O(b), we conclude that R 2 = O(b 2 ) and R 3 = O(b 2 ) can be neglected in Eq. ( 129). This however leads only to R 1 = O(1), so that more mathematical work needs to be done, as detailed in the Appendix F, to obtain

2 m LR 1 = b→0 C 4π 2 + πR surf e 2L + O(b/L) 2 . ( 133
)
The numerical constant C ≃ 8.91363 was calculated and called C(0) in [143]. R surf e remarkably is the surface contribution to the quantity R e in Eq. ( 116), it scales as b. It is non-zero only when the dispersion relation has a cusp at the surface of the first Brillouin zone. In this case, R 1 varies to first order in b, which comes in addition to the expected linear contribution of the Er e term in Eq. ( 129): This leads to the Juillet effect. More quantitatively, the first deviation of the eigenenergy from its zero-range limit E 0 , shown as a dashed line in Fig. 1a, is [243]:

E -E 0 ∼ b→0 m 2 E 0 re 8π 4 + mπR surf e 2 2 L 2 1 L 3 k∈R 3 1 2 k 2 m -E 0 2 (134)
D. Link between ∂E/∂re and the subleading short distance behavior of the pair distribution function

As shown by [Tab. II, Eqs. (3a,3b)] the short distance behavior of the pair distribution function (averaged over the center of mass position of the pair) diverges as 1/r 2 in 3D and as ln 2 r in 2D, with a coefficient proportional to C, that is related to the derivative of the energy with respect to the scattering length a. Here we show that a subleading term in this short distance behavior is related to the derivative of the energy with respect to the effective range r e . To this end, we explicitly write the next order term in the contact conditions [Tab. I, Eqs. (1a,1b)]. Three dimensions: Including the next order term in [Tab. I, Eq. (1a)] gives

ψ(r 1 , . . . , r N ) = rij →0 1 r ij - 1 a A ij (R ij , (r k ) k =i,j ) + r ij B ij (R ij , (r k ) k =i,j ) + 3 α=1 r ij,α L (α) ij (R ij , (r k ) k =i,j ) + O(r 2 ij ) ( 135
)
where we have distinguished between a singular part linear with the interparticle distance r ij and a regular part linear in the relative coordinates of i and j (r ij,α is the component along axis α of the vector r ij ). Injecting this form into Schrödinger's equation, keeping the resulting ∝ 1/r ij terms and using notation [Tab. V, Eq. ( 2)] gives

B ij (R ij , (r k ) k =i,j ) = - m 2 2 (E -H ij )A ij (R ij , (r k ) k =i,j ) (136)
[Tab. V, Eq. (1a)] thus becomes

∂E ∂r e = - 4π 2 m (A, B) (137)
We square [START_REF] Landau | Statistical Physics part[END_REF] and as in Sec. IV B we integrate over R ij , the r k 's and we sum over i < j. We further average G

↑↓ (r) over the direction of r to eliminate the contribution of the regular term L ij , defining Ḡ( 2)

↑↓ (r) = [G (2) ↑↓ (r) + G (2)
↑↓ (-r)]/2. We obtain [Tab. V, Eq. (3a)]. Two dimensions: Including next order terms in [Tab. I, Eq. (1b)] gives [244]:

ψ(r 1 , . . . , r N ) = rij →0 ln(r ij /a) A ij (R ij , (r k ) k =i,j ) +r 2 ij ln r ij B ij (R ij , (r k ) k =i,j )+ 2 α=1 r ij,α L (α) ij (R ij , (r k ) k =i,j ) + O(r 2 ij ) (138)
Proceeding as in 3D we obtain

B ij (R ij , (r k ) k =i,j ) = - m 4 2 (E -H ij )A ij (R ij , (r k ) k =i,j ) (139)
[Tab. V, Eq. (1b)] thus becomes

∂E ∂(r 2 e ) = - 4π 2 m (A, B) (140)
These equations finally leads to [Tab. V, Eq. (3b)].

E. Link between ∂E/∂re and the 1/k 6 subleading tail of the momentum distribution

A general idea given in [89] is that singular terms in the dependence of ψ on the interparticle distance r ij (at short distances) reflect into power-law tails in the momentum distribution n σ (k) given by Eq. ( 23). In Sec. IV A, we restricted to the leading order. Here we include the subleading term and we perform the same reasoning as in Sec. IV A to obtain

[245] [246] nσ (k) = k→∞ C k 4 + D k 6 + . . . (141)
where nσ (k) = 1 d d i=1 n σ (ku i ) and D is the linear combination of ∂E/∂r e and (A, ∆ R A) given in [Tab. V, Eqs. (4a,4b)]. Physically, the extra term (A, ∆ R A) results from the fact that the wavevector k 1 of a particle in an ↑↓ colliding pair is a linear combination of the relative wavevector k rel and of the center-of-mass wavevector K of the pair, so that, even if the probability distribution of k rel was exactly scaling as 1/k 4 rel , a non-zero K would generate a subleading 1/k 6 1 contribution in the single particle momentum distribution. This is apparent for the simple case of a free space dimer: When the dimer is at rest, ψ(r 1 , r 2 ) = φ dim (r 12 ), A 12 (R 12 ) is uniform and the extra term vanishes. When it has a momentum K, ψ(r 1 , r 2 ) = e iK•R12 φ dim (r 12 ), which shifts the single particle momentum distribution,

n mov ↑ (k) = n rest ↑ (k -K/2
). Applying this shift to the momentum tail C/k 4 gives, after continuous average over the direction of k, a subleading δD mov /k 6 contribution, with δD mov = CK 2 /2 in 3D and δD mov = CK 2 in 2D. Remarkably, the ratio of the extra term to C is proportional to the pair-center-of-mass kinetic energy.

In the N -body case, one can generalize this property by defining the mean center-of-mass kinetic energy of a ↑↓ pair at vanishing pair diameter, which is allowed in quantum mechanics since the center-of-mass operators and the relative-particle operators commute [247]. By a direct generalisation of the pair distribution function of Sec. IV B, one has for the opposite-spin pair density operator r ↑ , r ↓ |ρ

(2) ↑↓ |r ′ ↑ , r ′ ↓ = ψ † ↑ (r ′ ↑ ) ψ † ↓ (r ′ ↓ ) ψ↓ (r ↓ ) ψ↑ (r ↑ ) .
Whereas the usual pair-center-of-mass density operator is obtained by taking the trace over the relative coordinates r = r ↑r ↓ , we rather define it here by taking the limit of vanishing relative coordinates,

R|ρ (2) CoM |R ′ = N lim r→0 R + r 2 , R -r 2 |ρ
(2)

↑↓ |R ′ + r 2 , R ′ -r 2 φ 2 (r) ( 142
)
where the factor N is such that ρ(2)

CoM has a unit trace and φ(r) is the zero-energy scattering state of Eqs. [START_REF] Lee | [END_REF]10). Proceeding as in Sec. IV B we obtain

R|ρ (2) CoM |R ′ = N i<j ( k =i,j d d r k )A * ij (R ′ , (r k ) k =i,j ) × A ij (R, (r k ) k =i,j ) (143)
By taking the expectation value of -( 2 /4m)∆ R within ρ(2) CoM , we finally obtain for the mean pair-center-of-mass kinetic energy at vanishing diameter:

E r ↑↓ →0 kin pair-CoM = - 2 4m (A, ∆ R A) (A, A) (144)
where the denominator is ∝ C, see [Tab. II, Eqs. (2a,2b)].

VIII. GENERALIZATION TO ARBITRARY STATISTICAL MIXTURES

In this section, we generalize some of the relations derived in the previous sections for pure states to the case of arbitrary statistical mixtures. Let us first discuss zerorange interactions. We consider a statistical mixture of pure states ψ n with occupation probabilities p n , which is arbitrary, but non-pathological in the following sense [95]: Each ψ n satisfies the contact condition [Tab. I, Eqs. (1a,1b)]; moreover, p n decays sufficiently quickly at large n so that we have C = II, which were derived in Sec. IV for any pure state satisfying the contact conditions, obviously generalize to such a statistical mixture. The relations for the time derivative of E (Tab. II line 12) hold for any time-evolving pure state satisfying the contact conditions for a timedependent a(t), and thus also for any statistical mixture of such time-evolving pure states.

For lattice models, one can obviously take an average of the definition of Ĉ [Tab. III, Eqs. (1a,1b)] to define C = Ĉ for in any statistical mixture; taking averages of the relations between operators [Tab.III, lines 2,3,8] then gives relations valid for any statistical mixture.

IX. THERMODYNAMIC EQUILIBRIUM IN THE CANONICAL ENSEMBLE

We turn to the case of thermal equilibrium in the canonical ensemble. We shall use the notation

λ ≡ -1/a in 3D 1 2 ln a in 2D. ( 145
)
A. First order derivative of E

The thermal average in the canonical ensemble dE/dλ can be rewritten in the following more familiar way, as detailed in Appendix G:

dE dλ = dF dλ T = d Ē dλ S (146)
where (. . . ) is the canonical thermal average, F is the free energy and S is the entropy. Taking the thermal average of [Tab. II, Eqs. (4a,4b)] (which was shown above for any stationary state) thus gives [Tab. II, Eqs. (9a,9b)].

B. Second order derivative of E

Taking a thermal average of the line 8 in Tab. II we get after a simple manipulation:

1 2 d 2 E dλ 2 = 4π 2 m 2 1 2Z n,n ′ ;En =E n ′ e -βEn -e -βE n ′ E n -E n ′ × |(A (n ′ ) , A (n) )| 2 (147)
where Z = n exp(-βE n ). This implies

d 2 E dλ 2 < 0. ( 148
)
Moreover one can check that

d 2 F dλ 2 T - d 2 E dλ 2 = -β dE dλ 2 - dE dλ 2 < 0, (149)
which implies [Tab. II, Eqs. (10a,10b)]. In usual cold atom experiments, however, there is no thermal reservoir imposing a fixed temperature to the gas, one rather can achieve adiabatic transformations by a slow variation of the scattering length of the gas [148][149][150] where the entropy is fixed [START_REF] Landau | Physique Théorique: Physique Statistique[END_REF][START_REF] Carr | [END_REF][153]. One also more directly accesses the mean energy Ē of the gas rather than its free energy, even if the entropy is also measurable [38,39]. The second order derivative of Ē with respect to λ for a fixed entropy is thus the relevant quantity to consider. As shown in Appendix G one has in the canonical ensemble: For lattice models, the inequalities [Tab. III, Eq. ( 7)] are derived in the same way, by taking λ now equal to g 0 , and starting from the expression [Tab. III, Eq. ( 6)] of d 2 E n /dg 2 0 . For the case of a finite-range interaction potential V (r) in continuous space, the relations [Tab. IV, lines 1-3] which were derived for an arbitrary stationary state are generalized to the thermal equilibrium case in the same way. Finally, the relations which asymptotically hold in the zero-range regime, [Tab. III lines 9-10] for lattice models and [Tab. IV lines [4][5] for finite-range interaction potential models, which were justified for any eigenstate in the zero-range regime k typ b << 1 where the typical relative wavevector k typ is defined in terms of the considered eigenstate as described in Section III, remain true at thermal equilibrium with k typ now defined as the typical density-and temperature-dependent wavevector described in Section III, since all the eigenstates which are thermally populated with a non-negligible weight are expected to have a typical wavevector smaller or on the order of the thermal-equilibrium typical wavevector.

d 2 Ē dλ 2 S = d 2 E dλ 2 + Cov E, dE

C. Quantum-mechanical adiabaticity

To be complete, we also consider the process where λ is varied so slowly that there is adiabaticity in the many-body quantum mechanical sense: The adiabatic theorem of quantum mechanics [154] implies that in the limit where λ is changed infinitely slowly, the occupation probabilities of each eigenspace of the many-body Hamiltonian do not change with time, even in presence of level crossings [155]. We note that this may require macroscopically long evolution times for a large system. For an initial equilibrium state in the canonical ensemble, the mean energy then varies with λ as

E quant adiab (λ) = n e -β0En(λ0) Z 0 E n (λ) ( 151
)
where the subscript 0 refers to the initial state. Taking the second order derivative of ( 151) with respect to λ in λ = λ 0 gives

d 2 E quant adiab dλ 2 = d 2 E dλ 2 < 0. ( 152
)
Note that the sign of the second order derivative of E quant adiab remains negative at all λ provided one assumes that there is no level crossing in the many-body spectrum when λ is varied:

E n (λ) -E n ′ (λ) has the same sign as E n (λ 0 ) - E n ′ (λ 0
) for all indices n, n ′ , which allows to conclude on the sign with the same manipulation as the one having led to Eq. ( 147).

Thermodynamic vs quantum adiabaticity: The result of the isentropic transformation (150) and the one of the adiabatic transformation in the quantum sense [START_REF] Carr | [END_REF] differ by the second term in the right hand side of (150). A priori this term is extensive, and thus not negligible as compared to the first term. We have explicitly checked this expectation for the Bogoliubov model Hamiltonian of a weakly interacting Bose gas, which is however not really relevant since this Bogoliubov model corresponds to the peculiar case of an integrable dynamics. For a quantum ergodic system we now show that the second term in the right hand side of ( 150) is negligible in the thermodynamic limit, as a consequence of the Eigenstate Thermalization Hypothesis [156][157][158][159]. This Hypothesis was tested numerically for several interacting quantum systems [160][161][162]. It states that, for a large system, the expectation value ψ n | Ô|ψ n of a few-body observable Ô in a single eigenstate |ψ n of energy E n can be identified with the microcanonical average O mc (E n ) of Ô at that energy. Here the relevant operator Ô is the two-body observable (the so-called contact operator) such that d dλ E n = ψ n | Ô|ψ n . In the canonical ensemble, the energy fluctuations scale as V 1/2 where V is the system volume. We can thus expand the microcanonical average around the mean energy Ē:

O mc (E) = O mc ( Ē) + (E -Ē)O ′ mc ( Ē) + O(1) (153)
To leading order, we then find that Cov E,

dE dλ ∼ O ′ mc ( Ē)Var E and Var dE dλ ∼ [O ′ mc ( Ē)] 2
Var E, so that the second term in the right hand side of ( 150) is O(V 1/2 ) which is negligible as compared to the first term in that right hand side. For the considered quantity, this shows the equivalence of the thermodynamic adiabaticity and of the quantum adiabaticity for a large system. A microcanonical detour: We now argue that the quantum adiabatic expression [START_REF] Landau | Physique Théorique: Physique Statistique[END_REF] for the mean energy as a function of the slowly varying parameter λ can be obtained by a purely thermodynamic reasoning. This implies that the exponentially long evolution times a priori required to reach the quantum adiabatic regime for a large system are actually not necessary to obtain [START_REF] Landau | Physique Théorique: Physique Statistique[END_REF]. The first step is to realize that the initial canonical ensemble (for λ = λ 0 ) can be viewed as a statistical mixture of microcanonical ensembles [163]. These microcanonical ensembles correspond to non-overlapping energy intervals of width ∆, each interval contains many eigenstates, but ∆ is much smaller than the width of the probability distribution of the system energy in the canonical ensemble. For further convenience, we take ∆ ≪ k B T . One can label each energy interval by its central energy value, or more conveniently by its entropy S. If the eigenenergies E n (λ) are numbered in ascending order, the initial microcanonical ensemble of entropy S contains the eigenenergies with n 1 (S) ≤ n < n 2 (S) and

S = k B ln[n 2 (S) -n 1 (S)].
When λ is slowly varied, the entropy is conserved for our isolated system, and the microcanonical ensemble simply follows the evolution of the initial n 2 (S)n 1 (S) eigenstates, which cannot cross for an ergodic system and remain bunched in energy space. Furthermore, according to the Eigenstate Thermalization Hypothesis, the energy width E n2 -E n1 remains close to its initial value ∆: Each eigenenergy varies with a macroscopically large slope dE n /dλ but all the eigenenergies in the microcanonical ensemble have essentially the same slope [248]. The mean microcanonical energy for this isentropic evolution is thus

E mc (S, λ) = 1 n 2 (S) -n 1 (S) n2(S)-1 n=n1(S) E n (λ) (154)
Finally, we take the appropriate statistical mixture of the microcanonical ensembles (so as to reconstruct the initial λ = λ 0 canonical ensemble): The microcanonical ensemble of entropy S has an initial central energy E mc (S, λ 0 ), it is weighted in the statistical mixture by the usual expression P (S) = e S/kB e -βEmc(S,λ0) . Since ∆ ≪ k B T , one can identify e -βEmc(S,λ0) with e -βEn(λ0) , for n 1 (S) ≤ n < n 2 (S). The corresponding statistical average of ( 154) with the weight P (S) gives [START_REF] Landau | Physique Théorique: Physique Statistique[END_REF].

X. APPLICATIONS

In this Section, we apply some of the above relations in three dimensions, first to the two-body and three-body problems and then to the many-body problem. Except for the two-body case, we restrict to the infinite scattering length case a = ∞ in three dimensions.

A. Two-body problem in a harmonic trap: Finite range corrections

Two particles interact with the compact-support potential V (r 12 ; b) of range b and scattering length a in an isotropic harmonic potential U (r i ) = 1 2 mω 2 r 2 i . One separates out the center of mass, in an eigenstate of energy E cm . The relative motion is taken with zero angular momentum; its wavefunction ψ(r) is an eigenstate of energy E rel = E -E cm for a particle of mass µ = m/2 in the potential V (r; b) + µω 2 r 2 /2. We take in this subsection ω as the unit of energy and [ /(µω)] 1/2 as the unit of length. For r ≥ b the solution may be expressed in terms of the Whittaker function W , or equivalently of the Kummer function U , see §13 in [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]:

ψ(r) C 3 3D = W E rel 2 , 1 4 (r 2 ) r 3/2 = e -r 2 2 U ( 3 4 - E rel 2 , 3 2 , r 2 )(155) ψ(r) C 2 2D = W E rel 2 ,0 (r 2 ) r = e -r 2 2 U ( 1 -E rel 2 , 1, r 2 ) ( 156
)
where the factors C 2 and C 3 ensure that ψ is normalized to unity. The zero-range limit, where V (r; b) is replaced by the Bethe-Peierls contact conditions at the origin, is exactly solvable; it gives eigenenergies E 0 rel . We give here the finite range corrections to the energy in terms of r e . Three dimensions: Imposing the contact condition ψ(r) = A[r -1a -1 ] + O(r) to Eq. ( 155) gives an implicit equation for the spectrum in the zero-range limit, obtained in [165] with a different technique:

f (E 0 rel ) = - 1 a with f (E) ≡ - 2Γ( 3 4 -E 2 ) Γ( 1 4 -E 2 ) (157)
We have calculated the finite range corrections up to order two in b included, they remarkably involve only the effective range:

E rel = E 0 rel + E 0 rel r e f ′ + E 0 rel r e f ′ 2 1 E 0 rel - f ′′ 2f ′ + O(b 3 ) (158)
where the first and second order derivatives f ′ and f ′′ of f (E) are taken in E = E 0 rel . To obtain this expansion, we have used the result of Appendix D that one can neglect, at this order, the effect of the trapping potential for r ≤ b, so that the wavefunction is proportional to the free space scattering state at energy E rel = 2 k 2 /(2µ), ψ(r) = Aχ(r). Such an approximation was already proposed in [129,166,167], without analytical control on the resulting spectral error [249]. We have checked that the term of Eq. ( 158) linear in r e coincides with the prediction of [Tab. V, Eq. (1a)], due to the fact that, from relation 7.611(4) in [START_REF] Gradshteyn | Tables of integrals, series, and products[END_REF], the normalization factor in the zerorange limit obeys

(C 0 3 ) 2 2π 2 f ′ (E 0 rel )/Γ 2 ( 3 4 - E 0 rel 2 ) = 1.
The term in Eq. ( 158) linear in r e was already written explicitly in [START_REF] Werner | [END_REF]. This corresponds to the first order perturbative use of the modified version of the zero-range model, as put forward in [122]. It can also be obtained by solving to first order in r e the self-consistent equation considered in [127] obtained by replacing a 0 by a E [see Eq. ( 5) of [127]] into Eq. ( 6) of [127]. This self-consistent equation was also introduced in [166], and in [167] [see Eqs. (11,12,30) of that reference] with more elaborate forms for a E . With our notations and units this selfconsistent equation is simply

f (E) = -u(k = √ 2E) (159)
where u(k) is related to the s-wave scattering amplitude by Eq. ( 90). The self-consistent equation of [127] corresponds to the choice u(k) = 1 a -1 2 k 2 r e in Eq. ( 159). We have checked that solving that equation to second order in r e then exactly gives the term of Eq. ( 158) that is quadratic in r e . Our result of Appendix D shows that going to order three in r e with the self-consistent equation should not give the correct result, since one can then no longer neglect the effect of harmonic trapping within the interaction range. This clarifies the status of that self-consistent equation.

To ascertain this statement, we have calculated the ground state relative energy up to third order included in b, restricting for simplicity to an infinite scattering length, 1/a = 0 [250]. We find

E rel = 1 2 + r e 2π 1/2 + 2 -ln 2 4π r 2 e + (1 -ln 2)(2 -ln 2) 4π 3/2 - π 2 + 12 ln 2 2 192π 3/2 r 3 e - λ 2 + Λ 2 π 1/2 + O(b 4 ) (160)
Here

λ 2 is the coefficient of k 4 in the low-k expansion of u(k), u(k) = 1 a -1 2 k 2 r e + λ 2 k 4 + O(k 6
), it can be evaluated by a generalized Smorodinski relation [169]. On the contrary, Λ 2 is a new coefficient containing the effect of the trapping potential within the interaction range. It can be expressed in terms of the zero-energy free space scattering state φ(r), normalized as in Eq. ( 9):

Λ 2 = +∞ 0 dr r 2 [1 -u 2 0 (r)] (161)
with u 0 (r) = rφ(r). Although our derivation is for a compact support potential, we expect that our result is applicable as long as λ 2 and Λ 2 are finite. For both quantities, this requires (for 1/a = 0) that the interaction potential drops faster than 1/r 5 [169]. Interestingly, if one expands the self-consistent Eq. ( 159) up to order b 3 included, one exactly recovers Eq. ( 160), except for the term Λ 2 . This was expected from the fact that the derivation of ( 159) in [167] indeed neglects the trapping potential within the interaction range. This discussion is illustrated for the particular case of the square-well potential (182) in Fig. 2, with the exact spectrum obtained by matching the logarithmic derivative of a Whittaker M function for r = b -with the logarithmic derivative of a Whittaker W function for r = b + as in Eqs. (6.16,6.17,6.18) of [113] [251]. In this case, one finds r e = b [START_REF] Castin | Ultracold Fermi gases[END_REF] and, remarkably, Λ 2 = -2λ 2 so that the difference between the ground state energy of ( 159) and the exact ground state energy obeys

E self rel -E rel = Λ 2 π 1/2 + O(b 4 ) = 1 6 1 π 2 b 3 π 1/2 + O(b 4 ). ( 162
) Note that the case of two fermions with a square-well interaction in a harmonic trap was numerically studied in [170], for the s-wave and also for the p-wave case, with the exact spectrum compared to the self-consistent equation (159) or to its p-wave equivalent. No conclusion was given on the scaling with b of the difference between the exact and the approximate spectrum. Two dimensions: Imposing the contact condition ψ(r) = A ln(r/a) + O(r) to Eq. ( 156) gives an implicit equation for the spectrum in the zero-range limit [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF]165]:

ψ 1 -E 0 rel 2 -2ψ(1) = -2 ln a ( 163
)
where ψ is the digamma function. We have obtained the finite range correction

E rel = E 0 rel + 4r 2 e E 0 rel ψ ′ ( 1-E 0 rel 2 ) + O(b 4 ln 4 b) (164)
by neglecting the trapping potential for r ≤ b as justified by Appendix D, and by matching in r = b the scattering state Aχ(r) to Eq. ( 156). The bound on the error results in particular from the statement that . . . in Eq. ( 97) are FIG. 2: For two opposite spin fermions interacting in 3D via a potential of short range b in an isotropic harmonic trap, the self-consistent equation (159), derived e.g. in [167], gives the eigenenergies with an error of order b 3 , due to the fact that it neglects the effect of the harmonic trap within the interaction range, see Appendix D. This is illustrated with the ground state relative energy for a square-well potential of infinite scattering length: The deviation (solid line) between the approximate energy E self rel [solving Eq. ( 159)] and the exact one E rel (calculated as in [START_REF] Werner | [END_REF]) vanishes as b 3 , with a coefficient given by Eq. ( 162) (dotted line). µ is the reduced mass, ω is the angular oscillation frequency in the trap and

a ho = [ /(µω)] 1/2 .
O[(kb) 4 ln(a/b)], that one can e.g. check for the squarewell potential. As expected, the value of ∂E rel /∂(r 2 e ) in r e = 0 obtained from Eq. ( 164) coincides with [Tab, V, Eq. (1b)], knowing that the normalization factor in the zero-range limit, according to relation 7.611 (5) in [START_REF] Gradshteyn | Tables of integrals, series, and products[END_REF], obeys (C 0 2 ) 2 πψ ′ (

1-E 0 rel 2 )/Γ 2 ( 1-E 0 rel 2 ) = 1.

B. Three-body problem: corrections to exactly solvable cases and comparison with numerics

In this Subsection, we use the known analytical expressions for the three-body wavefunctions to compute the corrections to the spectrum to first order in the inverse scattering length 1/a and in the effective range r e . We shall consider not only spin-1/2 fermions, but also spinless bosons restricting to the universal stationary states [171,172] which do not depend on the three-body parameter.

The problem of three identical spinless bosons [171,172] or spin-1/2 fermions (say N ↑ = 2 and N ↓ = 1) [171,173] is exactly solvable in the unitary limit in an isotropic harmonic trap U (r) = 1 2 mω 2 r 2 . Here we restrict to zero total angular momentum (see however the last line of Appendix H) with a center of mass in its ground state, so that the normalization constants of the wavefunctions are also known analytically [START_REF] Werner | [END_REF]. Moreover we restrict to universal eigenstates [252]. The spectrum is then

E = E cm + (s + 1 + 2q) ω (165)
where E cm is the energy of the center of mass motion, s belongs to the infinite set of real positive solutions of

-s cos s π 2 + η 4 √ 3 sin s π 6 = 0 (166)
with η = +2 for bosons and -1 for fermions, and q is a non-negative integer quantum number describing the degree of excitation of an exactly decoupled bosonic breathing mode [121,174]. We restrict to states with q = 0. The case of a non-zero q is treated in subsection X C. a. Derivative of the energy with respect to 1/a. Injecting the expression of the regular part A of the normalized wavefunction [START_REF] Werner | [END_REF] into [Tab. II, Eqs. (2a,4a)] or its bosonic version (Tab. V, line 1 in [START_REF] Werner | Sections III and IV of[END_REF]) we obtain

∂E ∂(-1/a) a=∞ = 3 ω m Γ(s + 1 2 ) √ 2s sin s π 2 /Γ(s + 1) -cos s π 2 + s π 2 sin s π 2 + η 2π 3 √ 3 cos s π 6 (167)
For the lowest fermionic state, this gives (∂E/∂(1/a)) a=∞ ≃ -1.1980

3 ω/m, in agreement with the value -1. 19(2) which we extracted from the numerical solution of a finite-range model presented in Fig. 4a of [118], where the error bar comes from our simple way of extracting the derivative from the numerical data of [118].

b. Derivative of the energy with respect to the effective range. Using relation [Tab. V, Eq. (1a)], which holds not only for fermions but also for bosonic universal states, we obtain

∂E ∂r e a=∞ = mω 3 8 Γ(s -1 2 )s(s 2 -1 2 ) sin(s π 2 )/Γ(s + 1) -cos(s π 2 ) + s π 2 sin(s π 2 ) + η 2π 3 √ 3 cos(s π 6 ) (168)
For bosons, this result was derived previously using the method of [122] and found to agree with the numerical solution of a finite-range separable potential model for the lowest state [START_REF] Werner | [END_REF]. For fermions, [START_REF] Gradshteyn | Tables of integrals, series, and products[END_REF] agrees with the numerical data from Fig. 3 of [118] to ∼ 0.3% for the two lowest states and 5% for the third lowest state [253]; (168) also agrees to 3% with the numerical data from p. 21 of [START_REF] Werner | [END_REF] for the lowest state of a finite-range separable potential model. All these deviations are compatible with the estimated numerical accuracy.

C. N -body problem in an isotropic trap: Non-zero 1/a and re corrections

We now generalize subsection X B to the case of an arbitrary number N of spin-1/2 fermions (with an arbitrary spin configuration) at the unitary limit in an isotropic harmonic trap. Although one cannot calculate ∂E/∂(1/a) and ∂E/∂r e , some useful information can be obtained from the following remarkable property: For any initial stationary state, and after an arbitrary change of the isotropic trap curvature, the system experiences an undamped breathing at frequency 2ω, ω being the single atom oscillation frequency in the final trapping potential [121]. From this one can conclude that, in the case of a time independent trap, the system exhibits a SO(2, 1) dynamical symmetry [174]: The spectrum is a collection of semi-infinite ladders indexed by the natural integer q. Another crucial consequence is that the eigenstate wavefunctions are separable in N -body hyperspherical coordinates, with a know expression for the dependence with the hyperradius [174]. This implies that the functions A ij are also separable in (N -1)-body hyperspherical coordinates and that their hyperradial dependence is also known. As the eigenstates within a ladder have exactly the same hyperangular part, one can relate the energy derivatives (with respect to 1/a or r e ) for step q of a ladder to the derivative for the ground step of the same ladder, as detailed in Appendix H:

∂E ∂(1/a) q = ∂E ∂(1/a) 0 Γ(s + 1) Γ(s + q + 1) × q k=0 Γ(k + 1 2 ) Γ(k + 1)Γ( 1 2 ) 2 Γ(s + q -k + 1 2 )Γ(q + 1) Γ(s + 1 2 )Γ(q -k + 1) (169)
with the eigenenergy of step q is written as Eq. ( 165), s being now unknown for the general N -body problem. We have checked that this explicit result is consistent with the recursion relations derived in [175]. A similar type of result holds for the derivative with respect to r e :

∂E ∂r e q = ∂E ∂r e 0 Γ(s + 1) Γ(s + q + 1)

× q k=0 Γ(k + 3 2 ) Γ(k + 1)Γ( 3 2 ) 2 Γ(s + q -k -1 2 )Γ(q + 1) Γ(s -1 2 )Γ(q -k + 1) (170)
For non-zero 1/a or r e , the level spacing is not constant within a ladder, the system will not respond to a trap change by a monochromatic breathing mode. In a small system, a Fourier transform of the system response can give access to the Bohr frequencies (E q -E q-1 )/ , which would allow an experimental test of Eqs. (169,170). In the large N limit, for a system prepared in its ground state, we now show that the main effects of non-zero 1/a or r e on the breathing mode are a frequency change and a collapse. Let us take the macroscopic limit of Eqs. (169,170) for a fixed q: Using Stirling's formula for s → +∞ we obtain (173) For a non-polarized gas (with the same number N/2 of particles in each spin state) the local density approximation gives 4s ∼ (3N ) 4/3 ξ 1/2 [117,173] and it allows to obtain the derivative of the energy with respect to 1/a [100] or to r e in terms of ξ, ζ and ζ e , defined in Eqs. (179,186), so that

[∂E/∂(1/a)] q [∂E/∂(1/a)] 0 = 1 - q 4s + q(9q +
δω breath = 256ω 525πξ 5/4 ξ 1/2 ζ k F a + 2ζ e k F r e (174)
where we have introduced the Fermi momentum k F of the unpolarized trapped ideal gas with the same atom number N as the unitary gas, with 2 k 2 F /(2m) = (3N ) 1/3 ω. For r e = 0, we recover the superfluid hydrodynamic prediction of [176][177][178]. We have checked that the change of the mode frequency due to finite range effects can also be obtained from hydrodynamics [254]; this change in typical experiments is of the order of 0.1% for lithium and 0.5% for potassium, see subsection X E.

Furthermore, due to the presence of q 2 terms in Eqs. (171,172), the Bohr frequencies (E q -E q-1 )/ depend on the excitation degree q of the mode: If many steps of the ground state ladder are coherently populated, this can lead to a collapse of the breathing mode, which constitutes a mechanism for zero-temperature damping [179,180]. To coherently excite the breathing mode, we start with a ground state gas, with wavefunction ψ old , and we abruptly change at t = 0 the trap frequency from ω old to ω = λ 2 ω old . For the unitary gas, ψ old is deduced from the t = 0 + ground state ψ 0 by a dilation with scaling factor λ,

|ψ old = e -i D ln λ |ψ 0 ( 175
)
where D is the generator of the dilations [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF]174]. Using the representation of D in terms of the bosonic operator b [174], that annihilates an elementary excitation of the breathing mode ( b|q = q 1/2 |q -1 ), and restricting to |ǫ| ≪ 1, where ǫ = ln λ, one has

D ≃ -is 1/2 ( b † -b) (176)
so that the trap change prepares the breathing mode in a Glauber coherent state with mean occupation number q = ǫ 2 s and standard deviation ∆q = q1/2 . Similarly, the fluctuations of the squared radius of the gas i r 2 i /N , that can be measured, are given by -s 1/2 mω ( b + b † ) for small ǫ. In the large system limit, one can have q ≫ 1 so that 1 ≪ ∆q ≪ q. At times much shorter than the revival time 2π /|∂ 2 q E q |, one then replaces the discrete sum over q by an integral to obtain b (t) b (0) = e -t 2 /(2t 2 c ) with t c = ∆q ∂ 2 q E q q=q (177)

0 0.2 0.4 k F r 4×10 -4 8×10 -4 g (2)
(r) r

↑↓ (r) = ψ † ↑ (r) ψ † ↓ (0) ψ↓ (0) ψ↑ (r) of the homogeneous non-polarized unitary gas at zero temperature. Circles: fixed-node Monte Carlo results from Ref. [107]. Solid line: analytic expression (181), where the value ζ = 0.95 was taken to fit the Monte Carlo results. The arrow indicates the range b of the square-well interaction potential. Dashed line: analytic expression (184), with ζe = 0.12 [144].

For an unpolarized gas, using Eqs. (171,172) and the local density approximation, we obtain the inverse collapse time due to non-zero 1/a or r e :

(ωt c ) -1 = 64|ǫ| 35πξ(3N ) 2/3 3ζ 5k F a + 2ζ e k F r e 3ξ 1/2 (178)
For lithium experiments, t c is more than thousands of mode oscillation periods. To conclude with an exotic note, we recall that the q 2 terms in Eqs. (171,172) lead to the formation of a Schrödinger-cat-like state for the breathing mode at half the revival time [181].

D. Unitary Fermi gas: comparison with fixed-node Monte Carlo

For the homogeneous non-polarized unitary gas (i.e. the spin-1/2 Fermi gas in 3D with a = ∞ and N ↑ = N ↓ ) at zero temperature, we can compare our analytical expressions for the short-distance behavior of the one-body density matrix g (1) σσ and the pair distribution function g

(2) ↑↓ to the fixed-node Monte Carlo results in [START_REF] Astrakharchik | [END_REF][106][107]. In this case, g

σσ (R-r/2, R+r/2) and g ↑↓ (R-r/2, R+r/2) depend only on r and not on σ, R and the direction of r. Expanding the energy to first order in 1/(k F a) around the unitary limit yields:

E = E ideal ξ - ζ k F a + . . . ( 179
)
where E ideal is the ground state energy of the ideal gas, ξ and ζ are universal dimensionless numbers, and the Fermi wavevector is related to the density through

σσ (r) = ψ † σ (r) ψσ(0) of the homogeneous non-polarized unitary gas at zero temperature: comparison between the fixed-node Monte Carlo results from Ref. [106] (black solid line) and the analytic expression (180) for the small-kF r expansion of g ↑↓ , see Fig. 3.

k F = (3π 2 n) 1/3 . Expressing C in terms of ζ thanks to [Tab. II, Eqs. (2a,4a)] and Eq. ( 179), and inserting this into [Tab. II, Eq. (7a)], we get g (1) σσ (r) ≃

n 2 1 - 3ζ 10 k F r - ξ 10 (k F r) 2 + . . . . (180)
For a finite interaction range b, this expression is valid for b ≪ r ≪ k -1 F [255]. [Tab. IV, Eq. (4a)] yields

g (2) ↑↓ (r) ≃ kF r≪1 ζ 40π 3 k 4 F |φ(r)| 2 . ( 181
)
The interaction potential used in the Monte Carlo simulations [START_REF] Astrakharchik | [END_REF][106][107] is a square-well:

V (r) = - π 2 2 2 mb 2 θ(b -r) (182)
The corresponding zero-energy scattering state is

φ(r) = sin πr 2b r for r < b, φ(r) = 1 r for r > b (183)
and the range b was taken such that nb 3 = 10 -6 i.e. k F b = 0.0309367 . . . . Thus we can assume that we are in the zero-range limit k F b ≪ 1, so that (180,181) are applicable.

Figure 3 shows that the expression (181) for g

↑↓ fits well the Monte Carlo data of [107] if one adjusts the value of ζ to 0.95. This value is close to the value ζ ≃ 1.0 extracted from (179) and the E(1/a)-data of [START_REF] Astrakharchik | [END_REF].

Using ζ = 0.95 we can compare the expression (180) for g (1) σσ with Monte Carlo data of [106] without adjustable parameters. Figure 4 shows that the first order derivatives agree, while the second order derivatives are compatible within the statistical noise. This provides an interesting check of the numerical results, even though any wavefunction satisfying the contact condition [Tab. I, Eq. (1a)] leads to g A more interesting check is provided by our expression [Tab. V, Eq.(3a)] for the subleading term in the short range behavior of g

↑↓ (r) = ζ 40π 3 k 4 F r 2 - ζ e 20π 3 k 6 F + O(r) (184)
where ζ e is defined in Eq. (186). Remarkably, this expression is consistent with the fixed node Monte Carlo results of [107] if one uses the value of ζ e of [144], see Fig. 3.

E. Finite-range correction in simulations and experiments

We recall that, as we have seen in Section VII, the finite-range corrections to eigenenergies are, to leading order, of the form (∂E/∂r e ) r e for continuous-space models or (117) for lattice models, where the coefficients ∂E/∂r e , and ∂E/∂R e for lattice models, are modelindependent. This can be used in practice by extracting the values of these coefficients from numerical simulations, done with some convenient continuous-space or lattice models (usually a dramatic simplification of the atomic physics reality); then, knowing the value of r e in an experiment, one can compute the finite-range corrections present in the measurements, assuming that the universality of finite range corrections, derived in section VII for compact support potentials, also applies for multichannel O(1/r 6 ) models. The value of r e is predicted in Ref. [182] to be

r e = -2 R * 1 - a bg a 2 + 4π b 3 Γ 2 (1/4) Γ 2 (1/4) 2π - b a 2 + b 2 a 2 ( 185
)
where b is the van der Waals length b = (mC 6 / 2 ) 1/4 , a bg is the background scattering length and R * is the so-called Feshbach length [123]. We recall that the magnetic-field dependence of a close to a Feshbach resonance reads a(B) = a bg [1-∆B/(B-B 0 )] where B 0 is the resonance location and ∆B is the resonance width, and that R * = 2 /(ma bg µ b ∆B) where µ b is the effective magnetic moment of the closed-channel molecule. We note that the a-dependent terms in the second term of (185) are O(b 2 ) and thus do not contribute to the leading-order correction in b. In contrast, the a-dependence of the first term of ( 185) can be significant since a bg can be much larger than b (this is indeed the case for 6 Li) [256]. A key assumption of Ref. [182] is that the open-channel interaction potential is well approximated by -C 6 /r 6 down to interatomic distances r ≪ b. This assumption is well satisfied for alkali atoms [182,183]. Although we have not calculated the off-shell length ρ e explicitly, we have checked that it is finite for a -C 6 /r 6 potential [169].

As an illustration, we estimate the finite-range corrections to the non-polarized unitary gas energy in typical experiments. Similarly to (179), we have the expansion

E = E ideal (ξ + ζ e k F r e + . . . ) (186)
where E and E ideal are the ground state energies of the homogeneous Fermi gas (of fixed density n = k 3 F /(3π 2 )) for 1/a = 0 and a = 0 respectively. The value of ζ e was estimated both from fixed-node Monte Carlo and Auxiliary Field Quantum Monte Carlo to be ζ e = 0.12(3) [144] [257]. The value of r e as given by Eq. ( 185) is 4.7 nm for the B 0 ≃ 834 G resonance of 6 Li (in accordance with [184]) and 6.7 nm for the B 0 ≃ 202.1 G resonance of 40 K. The typical value of 1/k F is ≃ 400 nm in [48], while 1/k F at the trap center is ≃ 250 nm in [START_REF] Partridge | [END_REF] and ≃ 100 nm in [185], which respectively leads to a finite range correction to the homogeneous gas energy: δE E ≃ 0.4%, 0.6% and 2%.

In the case of lithium, this type of analysis was used in [48] to estimate the resulting experimental uncertainty on ξ.

XI. CONCLUSION

We derived relations between various observables for N spin-1/2 fermions in an external potential with zerorange or short-range interactions, in continuous space or on a lattice, in two or three dimensions. Some of our results generalize the ones of [89,95,96,99,102,103]: Large-momentum behavior of the momentum distribution, short-distance behavior of the pair distribution function and of the one-body density matrix, derivative of the energy with respect to the scattering length or to time, norm of the regular part of the wavefunction (defined through the behavior of the wavefunction when two particles approach each other), and, in the case of finite-range interactions, interaction energy, are all related to the same quantity C; and the difference between the total energy and the trapping potential energy is related to C and to a functional of the momentum distribution (which is also equal to the second order term in the short-distance expansion of the one-body density matrix). We also obtained entirely new relations: The second order derivative of the energy with respect to the inverse scattering length (or to the logarithm of the scattering length in two dimensions) is related to the regular part of the wavefunctions, and is negative at fixed entropy; and the derivative of the energy with respect to the effective range r e of the interaction potential (or to r 2 e in 2D) is also related to the regular part, to the subleading short distance behavior of the pair distribution function, and to the subleading 1/k 6 tail of the momentum distribution. We have found unexpected subtleties in the validity condition of the derived expression of this derivative in 2D: Our expression for ∂E/∂(r 2 e ) applies because, for the class of interaction potentials that we have specified, the effective range squared r 2 e is much larger than the true range squared b 2 , than the length squared ρ 2 e characterizing the low-energy s-wave off-shell T -matrix, and than the length squared R 2 1 characterizing the low energy p-wave scattering amplitude, by logarithmic factors that diverge in the zero-range limit. In 3D, for lattice models, our expression for ∂E/∂r e applies only for magic dispersion relations where an extra parameter R e quantifying the breaking of Galilean invariance (as predicted in [START_REF] Burovski | [END_REF]) vanishes; also, the magic dispersion relation should not have cusps at the border of the first Brillouin zone otherwise the so-called Juillet effect compromises the validity of our ∂E/∂r e expression for finite size systems. We have explicitly constructed such a magic relation, that may be useful to reduce lattice discretization effects in Quantum Monte Carlo simulations. We also considered models with a momentum cut-off used in Quantum Monte Carlo calculations, either in continuous space [57] or on a lattice [53,56,146,147]: Surprisingly, in the infinite cut-off limit, the breaking of Galilean invariance survives and one does not exactly recover the unitary gas.

Applications of general relations were presented in three dimensions. For two particles in an isotropic harmonic trap, finite-interaction-range corrections were obtained, and were found to be universal up to order r 2 e included in 3D; in particular, this clarifies analytically the validity of some approximation and self-consistent equation introduced in [127,129,166,167] that neglect the effect of the trapping potential within the interaction range. For the universal states of three particles with an infinite scattering length in an isotropic harmonic trap, the derivatives of the energy with respect to the inverse scattering length and with respect to the effective range were computed analytically and found to agree with available numerics. For the unitary gas in an isotropic harmonic trap, which has a SO(2, 1) dynamical symmetry and an undamped breathing mode of frequency 2ω, we have determined the relative finite-1/a and finite range energy corrections within each SO(2, 1) ladder, which allows in the large-N limit to obtain the frequency shift and the collapse time of the breathing mode. For the bulk unitary Fermi gas, existing fixednode Monte Carlo data were checked to satisfy exact relations. Also, the finite-interaction-range correction to the unitary gas energy expected from our results to be (to leading order) model-independent and thus extractable from Quantum Monte Carlo results, was estimated for typical experiments: This quantifies one of the experimental uncertainties on the Bertsch parameter ξ.

The relations obtained here may be used in various other contexts. For example, the result [Tab. II, Eqs. (11a,11b)] on the sign of the second order derivative of E at constant entropy is relevant to adiabatic ramp experiments [38,39,150,[START_REF] Carr | [END_REF]185], and the relation [Tab. III, Eq. (8a)] allows to directly compute C using determinantal diagrammatic Monte Carlo [186] and bold diagrammatic Monte Carlo [59,[START_REF] Talk | INT Program 11-1 "Fermions from cold atoms to neutron stars: Benchmarking the manybody problem[END_REF][START_REF] Van Houcke | [END_REF]. C is directly related to the closed-channel fraction in a two-channel model [98,100], which allowed to extract it [100] from the experimental photoassociation measurements in [START_REF] Partridge | [END_REF]. C was measured from the tail of the momentum distribution [52]. For the homogeneous gas C was extracted from measurements of the equation of state [45]. C also plays an important role in the theory of radiofrequency spectra [99,[189][190][191][192][193] and in finite-a virial theorems [97,194,195], as verified experimentally [52]. C was also extracted from the momentum tail of the static structure factor S(k), which is the Fourier transform of the spin-independent pair distribution function n(r)n(0) and was measured by Bragg spectroscopy [50,51]. In principle one can also measure via S(k) the parameter ζ e quantifying the finite range correction to the unitary gas energy, from the relation

∂E ∂r e = - π 2 m d 3 k (2π) 3 S(k) - C 4k (188)
resulting from [Tab. V, Eq. (3a)]. This procedure is not hampered by the small value of k F r e in present experiments, contrarily to the extraction of ζ e from a direct measurement of the gas relative energy correction ∝ ζ e k F r e 10 -2 . We can think of several generalizations of the relations presented here. All relations can be extended to the case of periodic boundary conditions. The techniques used here can be applied to the one-dimensional case to generalize the relations of [89]. For two-channel or multichannel models one may derive relations other than the ones of [98][99][100]. Generalization of the present relations to arbitrary mixtures of atomic species, and to situations (such as indistinguishable bosons) where the Efimov effect takes place, was given in [START_REF] Werner | Sections III and IV of[END_REF].

taken for fixed r j and fixed (r k ) k =i,j . Using R ij = r j + ǫu/2 with u ≡ (r ir j )/r ij we get

ψ n = ǫ→0 1 ǫ - 1 a n A (n) ij + 1 2 u • ∇ Rij A (n) ij + O(ǫ) (B4) ∇ ri ψ n = ǫ→0 - u ǫ 2 A (n) ij + 1 2ǫ ∇ Rij A (n) ij -u(u •∇ Rij A (n) ij ) +O(1) (B5)
where n equals 1 or 2, and the functions A

(n) ij and ∇ Rij A (n)
ij are taken at (r j , (r k ) k =i,j ). This simply gives

Sǫ(rj ) [ψ * 1 ∇ ri ψ 2 -ψ 2 ∇ ri ψ * 1 ] • dS = ǫ→0 4π 1 a 1 - 1 a 2 × A (1) * ij A (2)
ij + O(ǫ) (B6) because the leading order term cancels and most angular integrals vanish. Inserting this into (B3) gives the desired lemma (33).

Two dimensions:

The derivation is analogous to the 3D case. In (B3), the double integral on the sphere of course has to be replaced by a simple integral on the circle. Instead of (B4,B5), we now obtain, from the 2D contact condition [Tab. I, Eq. (1b)],

ψ n = ǫ→0 ln(ǫ/a n ) A (n) ij + O(ǫ ln ǫ) (B7) ∇ ri ψ n = ǫ→0 u ǫ A (n) ij + O(ln ǫ), (B8)
which gives

Sǫ(rj ) [ψ * 1 ∇ ri ψ 2 -ψ 2 ∇ ri ψ * 1 ] • dS = ǫ→0 2π ln(a 2 /a 1 ) × A (1) * ij A (2) ij + O(ǫ ln 2 ǫ) (B9)
and yields the lemma [START_REF] Partridge | [END_REF].

Appendix C: Zero-range limit of the lattice model's contact

In this appendix, we show that our definition [Tab. III, Eqs. (1a,1b)] of the contact operator Ĉ within the lattice model agrees in the zero-range limit b → 0 with the way [Tab. II, Eq. ( 1)] C is usually defined within the zero-range model. [defined by (19)]. Linearity immediately gives

A = n c (n) A (n) and A b = n c (n) A (n) b , as well as C b = n,m (c (n) b ) * c (m) b ψ (n) b | Ĉ|ψ (m) b
. Expressing Ĉ in terms of H int thanks to [Tab. III, Eq. ( 2)], and using the lemma (56) as well as (15), we get ψ (19) The motivation of this Appendix is to justify the fact that, in Eq. ( 88) and in its equivalent form in 2D for a N -body problem, we have neglected the effect of the trapping potential within the interaction range. In the case of an isotropic harmonic trap, the exact form of Eq. ( 88) contains the external potential term 1 4 mω 2 r 2 ij . This issue is thus mappable to the two-body problem in a trap with a finite range interaction, which was the object of numerous studies in 3D [127,129,166,167] that have however not analytically quantified the effect of the trapping potential within the interaction range. After elimination of the center of mass motion and restriction to a zero angular momentum, one faces the 3D or 2D eigenvalue problem The starting point is the Hellmann-Feynman theorem, with ψ real and normalized to unity:

(n) b | Ĉ|ψ (m) b = (4π) 2 (A (n) b , A (m) b ). When b → 0, we expect that this last quantity tends to (4π) 2 (A (n) , A (m) ) because A (n) b → A (n) [see
Eψ(r) = - 2 m ∆ψ(r) + 1 4 mω 2 r 2 + V (r; b) ψ(r) (
dE db = d d r ψ 2 (r) ∂ b V (r; b). (D2)
To reexpress this integral in a more operational way, we introduce the solution ψ(r) of Schrödinger's equation with the same eigenvalue E but for the interaction potential V (r; b) of a different range b. This solution ψ(r) remains finite in r = 0 but it diverges at infinity and cannot be L 2 -normalized. In what follows we take a convenient normalization of ψ such that lim b→b ψ = ψ.

We multiply Schrödinger's equation for ψ (respectively ψ) by ψ (respectively ψ) and we integrate the difference of the two resulting equations over the domain r < R. Using the divergence theorem, the Wronskian W (R) appears,

W (r) ≡ ψ(r)ψ ′ (r) -ψ(r) ψ′ (r). (D3)
For r > b, b, the Wronskian satisfies the differential equation

W ′ (r) = -d-1 r d-1 W (r), so that, for large R, W (R) = w R d-1 and w = m 2 +∞ 0 dr r d-1 [V (r; b) -V (r; b)] ψ(r)ψ(r). (D4)
Turning back to the Hellmann-Feynman formula (D2), we obtain the exact relation

dE db = 2(d -1)π 2 m lim b→b w b - b (D5)
It remains to calculate w treating perturbatively the trapping potential within the interaction range.

To zeroth order, one neglects the trapping potential for r ≤ b [or r ≤ b for ψ], so that ψ (0) (r) = Aχ(r), where χ is the scattering state of energy E for V (r; b). Taking for simplicity E > 0, we set E = 2 k 2 /m, k > 0, and χ is normalized as in Eqs. (89,95). Note that A is then fully specified by the continuous matching of ψ (0) in r = b to the outer solution in the trapping potential (that can be expressed in terms of Whittaker functions, see subsection X A) and by the fact that ψ is normalized to unity. We also have ψ(0) (r) = A χ(r) for r ≤ b, where χ is the scattering state of energy E for V (r; b) and the same prefactor A was taken for convenience. The zerothorder Wronskian W (0) can then be calculated explicitly, in particular using relations 8.477(1) and 8.473(4,5) of [START_REF] Gradshteyn | Tables of integrals, series, and products[END_REF]. We use Eqs. (93,97) We have checked that the b → 0 limit of these relations coincide with [Tab. V, Eqs. (1a,1b)].

To first order, we treat the trapping potential perturbatively within the interaction range. We rescale the distance by b, so that ψ (1) (r) = f (x), and χ(r) = N u(x), where x = r/b and the function u(x) is normalized by the condition u(0) = 1. The function f solves the inhomogeneous Schrödinger equation:

f ′′ (x)+ d -1 x f ′ (x)+ k 2 b 2 - mb 2 2 V (bx; b) f (x) = F x 2 u(x) with F = 1 4 AN m 2 ω 2 2 b 4 (D8)
The function u(x) is a solution of the corresponding homogeneous equation. A second solution v(x) can be constructed, that diverges for x → 0. It is of the form This can be seen with the method of variation of constants, where one sets (f (x), f ′ (x)) = α(x)(u(x), u ′ (x))+ β(x)(v(x), v ′ (x)), with the boundary conditions α(0) = 0 (so that ψ (1) does not duplicates the zeroth order solution) and β(0) = 0 (so that ψ (1) does not diverge in r = 0). This leads to where the factor A converges to a finite, energydependent value for b → 0.

v(x) = -u(x) x + Z 3 (x) with Z 3 (x) = O(x) for x → 0 in 3D, and v(x) = u(x) ln x + Z 2 (x) with Z 2 (x) = O(x 2 ) for x → 0 in 2D. More precisely, one has Z d (x) = u(x) x 0 dy y 1-d [-1 + 1/u 2 (
Here, in the context of the Juillet effect for lattice models, we justify the expansion (133). The quantity R 1 defined in Eq. ( 130) may be expressed in terms of the difference between an integral and a 3D Riemann sum. We are then guided by the following type of results: If f (x) is a C ∞ function inside the cube B = [-1/2, 1/2] 3 , then for ε = 1/(2N + 1), with the integer N → +∞:

B d 3 x f (x) -ε 3 n f (εn) = ε 2 24 B d 3 x ∆f (x) + O(ε 4 ) (F1)
where ∆f is the Laplacian of f and the sum over n ranges over {-N, . . . , N } 3 . To show this lemma, we introduce the short-hand notation S[f ] for the left-hand side of (F1) and we pave B with little cubes of volume ε 3 and of centers εn:

S[f ] = n ε 3 B d 3 x[f (εn + εx) -f (εn)]. (F2)
Then we use the fourth-order Taylor-Lagrange formula for f restricted to the line connecting εn to εn + εx: f (εn+εx)-f (εn) = ε 2 2 i,j x i x j ∂ i ∂ j f (εn)+odd+O(ε 4 ) where "odd" stands for terms that are linear and cubic in the components of x, and O(ε 4 ) results from the fact that the fourth-order derivatives of f are uniformly bounded on B. Integration over x inside the cube B eliminates the odd terms, and the i = j quadratic terms, so that

S[f ] = ε 5 24 n ∆f (εn) + O(ε 2 ) . ( F3
)
A Riemann sum thus deviates from the integral by O(ε 2 ), for a C ∞ integrand. Applying this conclusion to Eq. (F3), where ∆f is C ∞ , we obtain the desired Eq. (F1). This result is however not immediate to apply to the quantity R 1 because the integrand of R 1 is singular in k = 0. We thus use several steps. We first consider the quantity R 1 for a quadratic dispersion relation that is cut in a smooth way: One twice replaces 1/(2ǫ k ) in Eq. ( 130) by φ(kb/2π)/( 2 k 2 /m) where φ(x) is a C ∞ rotationally invariant function, equal to 1 in x = 0, and of compact support included inside B ≡ [-1/2, 1/2] 3 (which allows to replace the set D by R 3 in the integration and in the summation). After the change of variable k = 2πx/L, we decompose R 3 as a collection of cubes of size unity (as in [143]), to obtain

(F4)

with h = 2π is Planck's constant and ε ≡ b/L is the small parameter. As shown in [143], the right-hand side of Eq. (F4) has a finite limit when ε → 0, here called C ≃ 8.91363, that one can obtain by taking ε to zero inside the sum and the integral, which amounts to replacing φ by unity. The deviation of Eq. (F4) from its ε → 0 limit can thus be exactly written as {S[f ] + ε 3 f (0)}/ε, with S[f ] = ε 3 n∈Z 3 B d 3 x[f (εn+εx)-f (εn)], that we treat as we did for Eq. (F2). Here f (x) = [φ(x) -1]/x 2 (extended by continuity to x = 0) is a C ∞ function since φ is rotationally invariant. In the fourth-order Taylor-Lagrange formula, O(ε 4 ) is replaced with the more accurate O( ε 4 (1+ε 2 n 2 ) 3 ), due to the fact that the fourth order derivatives of f (x) are uniformly bounded and decrease as 1/x 6 at infinity. The integral of the Laplacian of f appears as in Eq. (F1), except that is in integrated over the whole R 3 space, which gives zero. We finally obtain

h 2 L m R φ 1 = b→0 C + b L 2 lim x→0 φ(x) -1 x 2 + O(b/L) 3 . (F5)
Turning back to the lattice model, we now evaluate how R 1 deviates from its b → 0 limit for the uncut parabolic dispersion relation k → 2 k 2 /(2m). The difference between the smoothly-cut R φ 1 and the uncut R parab 1 (times h 2 L/m) is now of the form ε 2 f (0) plus 1 ε times the difference S[f ] between an integral and a Riemann sum, with f (x) = [φ(x) -1]/x 2 as before is C ∞ . We then use the result (F1), the key point being that the integration domain is B (rather than the whole space), so that the integral of the Laplacian of f over B gives a non-zero surface contribution, equals to the flux of the gradient of f through the surface of B. This leads to Eq. ( 133) for the particular case of the parabolic dispersion relation. The surface term can be evaluated explicitly, as in subsection VII B, from the integral evaluated in polar coordinates:

[-1,1] 2 dxdy (1 + x 2 + y 2 ) 2 = √ 8 arcsin 1 √ 3 (F6)
Finally, we consider a general dispersion relation (132), with η x = 1 2 x 2 + O(x 4 ) for x → 0. One can consider the difference between the corresponding R 1 and R parab The corresponding function f (x) = 1/(2η x ) -1/x 2 is then not C ∞ in x = 0. E.g. for the Hubbard model, η x = [3i cos 2πx i ]/(2π) 2 is not rotationally invariant and f (x) behaves as i x 4

i /x 4 at low x, its x → 0 limit depends on the direction of x. This limiting behavior is however scaling invariant, a feature that holds for a general dispersion relation. The n th order derivatives of f are then O(1/x n ) for x → 0. For this class of functions, we introduce S * [f ] defined as S[f ] in Eq. (F1) except that one excludes the term n = 0 in the sum. This implies that in the equivalent of Eq. (F2), there is an isolated contribution, the integral of f over εB, which is O(ε 3 ) and negligible. Then reproducing the analysis with the fourth-order Taylor-Lagrange formula we obtain

S * [f ] = ε 2 24 B d 3 x ∆f (x) + O(ε 3 ). (F7) As h 2 L m (R 1 -R parab 1
) = S * [f ]/ε, we obtain Eq. ( 133).

Appendix G: Isentropic derivatives of the mean energy in the canonical ensemble One considers a system with a Hamiltonian H(λ) depending on some parameter λ, and at thermal equilibrium in the canonical ensemble at temperature T , with a density operator ρ = exp(-βH)/Z. In terms of the partition function Z(T, λ) = Tr e -βH(λ) , with β = 1/(k B T ), one has the usual relations for the free -E dE dλ . Insertion of (G6) into (G4) gives (146). Insertion of (G7,G8,G9) into (G5) gives (150).

Appendix H: Non-zero 1/a and re corrections within a ladder of the trapped unitary gas

For N spin-1/2 fermions at the unitary limit in an isotropic harmonic trap, there is separability of the wavefunction in internal hyperspherical coordinates [174]: ψ(r 1 , . . . , r N ) = ψ cm (C)R -(3N -5)/2 F (R)Φ(Ω) (H1) where C is the center-of-mass location of the N fermions, R is the hyperradius and Ω is a set of 3N -4 hyperangles constructed from the Jacobi coordinates (see e.g. [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF]). One has the general formulas C = N i=1 m i r i /M and R 2 = N i=1 m i (r i -C) 2 / m, where M = N i=1 m i is the total mass, m an arbitrary mass unit, and m i is the mass of particle i (here equal to m). We shall not need the expression of the hyperangles. Eq. (H1) is due to the separability of the non-interacting Hamiltonian in a harmonic trap, and to the fact that the Bethe-Peierls contact condition do not break this separability for 1/a = 0. One finds that Φ(Ω) is an eigenstate of the Laplacian on the unit sphere of dimension 3N -4, with contact conditions. Corresponding eigenvalues are conveniently written as ( 3N -5

2 ) 2s 2 , s > 0. In the N -body case, s is not known analytically. On the contrary, F (R) solves a simple 2D Schrödinger-like equation

(E -E cm )F (R) = - 2 2 m F ′′ (R) + 1 R F ′ (R) + 2 s 2 2 mR 2 + 1 2 mω 2 R 2 F (R) (H2)
This leads to a spectrum of the form (165), with eigenfunctions expressed in terms of generalized Laguerre polynomials multiplied by a Gaussian [174].

To derive Eqs. (169,170), one uses the fact that this separability extends to the functions A ij (R ij , (r k ) k =i,j ). One takes the limit r ij → 0 for a fixed R ij in (H1): Φ(Ω) diverges as R/r ij (since it depends on the hyperangles only), C and R respectively tend to the center-of-mass position Č and the hyperradius Ř of a fictitious system of N -1 particles of total mass M = N m, composed of a particle of position R ij and mass 2m, and N -2 fermions of positions r k , k = i, j and mass m [261]. We thus obtain the form A ij (R ij , (r k ) k =i,j ) = ψ cm ( Č) Ř-(3N-7)/2 F ( Ř) Φ( Ω) (H3) It remains to express the Hamiltonian [Tab. V, Eq. ( 2)] of the fictitious system in terms of its center-of-mass Č and hyperspherical coordinates ( Ř, Ω):

H ij = - 2 2M ∆ Č + 1 2 M ω 2 Č2 - 2 2 m [∂ 2 Ř + 3N -7 Ř ∂ Ř + 1 Ř2 ∆Ω] + 1 2 mω 2 Ř2 .
In the integral over Ř, we use the fact that F solves (H2) and we integrate by parts to obtain for s > 1/2 [262]:

(A, A) = ∞ 0 d Ř F 2 ( Ř) d Ω Φ2 ( Ω) (H4) (A, (H-E)A) = ∞ 0 d Ř 2 F 2 ( Ř) 2 mR 2 d Ω Φ( Ω)[Λ-∆Ω] Φ( Ω)(H5)
with Λ = ( 3N -8 2 ) 2 + 1 4s 2 . Within a given SO(2, 1) energy ladder, Φ is fixed, only F depends on the quantum number q. The normalization of ψ to unity imposes that +∞ 0 dR RF 2 (R) is also fixed within a ladder. From known integrals involving the Laguerre polynomials, see e.g. Eq. (F7) in [START_REF] Castin | [END_REF], one gets Eqs. (169,170). Another byproduct is for N = 3, where Φ( Ω) is a spherical harmonic of spin l: This leads to ∂ re E/∂ (-1/a) E = parison with a coupled-channel calculation, provided that the separable-potential range in [100] was adjusted to reproduce the correct value of re at resonance. [257] As discussed around Eq. ( 117), one has to take into account not only re but also Re for lattice models, which was not done in [144].

[258] In the case of an infinite scattering length, one has to take a finite a so that this expression makes sense, and only then take the limit |a| → ∞ (this comes from the fact that the scattering amplitude at zero energy is infinite in this case).

[259] In other words, the Dirac distributions originating from the action of the Laplacian onto the 1/rij divergences can be ignored.

[260] This also results from the fact that u(1) is not particularly close to zero: For 1/a = 0 in 3D, u(1)/u ′ (1) = -1.

[261] If the first Jacobi coordinates of the N particles are chosen to be ∝ rij , the other ones tend to the Jacobi coordinates of the fictitious system.

[262] Note that F (R) scales as R s for R → 0. Also, each term of the sum over i < j gives the same contribution, due to the fermionic antisymmetry, and we have dropped this sum and the ij indices for simplicity.

2 2m

 2 both in the first Brillouin zone [217]. The operator ∆ in (3) is the lattice version of the Laplacian defined by -r|∆ r |k ≡ ǫ k r|k . The simplest choice for the dispersion relation is

  The result then follows by writing dE/dg 0 = dE/d(-1/a) • d(-1/a)/dg 0 in 3D and dE/dg 0 = dE/d(ln a) • d(ln a)/dg 0 in 2D, and by using the definition [Tab. IV, Eqs. (1a,1b)] of C as well as the following lemmas:

  ) where χ 1 and χ 2 are the same energy E scattering states for two different values b 1 and b 2 of the potential range. Then dividing this expression by b 1b 2 , taking the limit b 1 → b 2 , and afterwards the limit b 2 → 0 for which the low-k expansion holds:

2 mb 2 [ 3 -

 223 A popular model for Quantum Monte Carlo simulations is the Hubbard model, that leads to the dispersion relation ǫ Hub k = cos(k x b)cos(k y b)cos(k z b)] (as already mentioned in subsection III B). This leads to K 3 ≃ 3.175 911 6. Again, both r e and R e differ from zero: r e ≃ -0.305 718b and R e ≃ -0.264 659b

FIG. 1 :

 1 FIG.1:(Color online) Illustration of the Juillet effect for the lattice model: In the cubic box [0, L] 3 with periodic boundary conditions, ground state energy of two opposite spin fermions as a function of the grid spacing b, for an infinite scattering length (1/a = 0), for a total momentum equal to 0 in (a) and equal to 2πL ez in (b). Three dispersion relations ǫ k are considered, the quartic one of Eq. (120) with zero effective range re = 0 (in blue, lower set), and the magic one(122) with re = Re = 0 with the parameters of Eq. (123) (in black, upper set) and of Eq. (124) (in red, middle set). The fact that the energy varies linearly in b for the re = 0 quartic dispersion relation at zero total momentum is the Juillet effect explained in Sec. VII C, and the corresponding dashed line is the analytical result(134). At non-zero total momentum the quartic dispersion relation leads to an energy variation linear in b as expected e.g. from the fact that its has a non-zero Re [the dotted line is a linear fit for b/L ≤ 0.01]. The magic dispersion relations lead to a O(b 2 ) variation of the energy both at zero and non-zero total momentum [the dotted lines are purely quadratic fits performed for b/L ≤ 0.02].

n

  p n C n , where C n (resp. C) is defined by [Tab. II, Eq. (1)] with n σ (k) = c † σ (k)c σ (k) and . = ψ n | . |ψ n (resp. . = n p n ψ n | . |ψ n ). Then, the relations in lines 3, 5, 6 and 7 of Table

dλ 2 -

 2 Var(E)Var dE dλ k B T Var(E) (150) where Var(X) and Cov(X, Y ) stand for the variance of the quantity X and the covariance of the quantities X and Y in the canonical ensemble, respectively. From the Cauchy-Schwarz inequality [Cov(X, Y )] 2 ≤ Var(X)Var(Y ), and from the inequality (148), we thus obtain [Tab. II, Eqs. (11a,11b)].

  FIG. 3:Pair distribution function g

FIG. 4 :

 4 FIG. 4: (Color online) One-body density matrix g

  σσ up to first order (red dashed straight line) and second order (blue dotted parabola) where we took the value ζ = 0.95 extracted from the Monte Carlo data for g(2)

  ↑↓ functions satisfying [Tab. II, Eqs. (3a,6a)] with values of C compatible with each other.

1 .

 1 Stationary state Let us first consider an eigenstate |ψ of the zero-range model with an energy E. Let |ψ b denote the eigenstate of the lattice model which tends to |ψ when b → 0, and E b the corresponding eigenenergy. Then, C b ≡ ψ b | Ĉ|ψ b tends to the contact C of the state ψ [defined in Tab. II, Eq. (1)] when b → 0. Indeed, C is related to dE/d(-1/a) by [Tab. II, Eq. (4a)]; C b is related to dE b /d(-1/a) by [Tab. II, Eq. (4a)]; and the function E b (1/a) should tend smoothly to E(1/a) when b → 0.

2 .

 2 Arbitrary pure state We now consider any pure state |ψ satisfying the contact condition [Tab. I, Eq. (1a)]. We will show that C b ≡ ψ b | Ĉ|ψ b tends to the contact C of the state |ψ [defined in Tab. II, Eq. (1)] when b → 0, where |ψ b is defined as follows: Writing |ψ as a linear combination n c (n) |ψ (n) of the zero-range model's eigenstates |ψ (n) , we define the linear combination |ψ b ≡ n c (n) |ψ (n) b of the lattice-model's eigenstates |ψ (n) b. We consider only the 3D case, the derivation being almost identical in 2D. Let A and A (n) denote the regular parts of ψ and ψ(n) [defined by the contact condition Tab. I, Eq. (1a)], and A b and A (n) b denote the regular parts of ψ b and ψ (n) b

  and the discussion thereafter]. Thus C b indeed tends to C. Appendix D: Spectral effect of the trapping potential within the interaction range

  D1)with the conditions that ψ diverges neither in r = 0 nor at infinity. The rotationally invariant compact support potential V (r; b) of range b is of the minimal depth ensuring a fixed scattering length a (as discussed in subsection VII B). In the limit b → 0, where E converges to a finite value, we show that neglecting the effect of the trapping potential within the interaction range r ≤ b, as done in subsection VII A, introduces on the eigenenergy E an error O(b 3 ) in 3D and O[b 4 ln 2 (a/b)] in 2D, which thus does not affect the results [Tab. V, Eqs. (1a,1b)].

  , with . . . = O[(kb) 4 ln(a/b)] in (97) [as we have checked for the square well], to obtain dE db O[b 3 ln(a/b)] (D7)

  y)]. Since the expression between square brackets in the left-hand side of Eq. (D8) is O(1), u(x) and Z d (x) are O(1) for x ≤ 1. A first consequence is that the factor N scales as 1/b in 3D and as ln(a/b) in 2D [260]. A second consequence is that, both in two and three dimensions, ψ (1) (b) and b ψ (1) ′ (b) = O(F ). (D9)

  for ψ(1) . From Eq. (D9) and its counterpart for ψ(1) ( b), ψ(1) ′ ( b), we can estimate the variation of the Wronskian W (R) for R close to b, b, and thus the variation w(1) of w due to the trapping potential. Dividing by b -b and taking the limit b → b as in Eq. (D5) amounts to taking a derivative with respect to b, which gives an additional factor O(1/b). Finally the error δE introduced on the eigenenergy by the neglection of the trapping potential within the interaction range is bounded in the zero range limit b → 0 as δE 3D = O(mω 2 b 3 A 2 ) (D12) δE 2D = O[mω 2 b 4 A 2 ln 2 (a/b)] (D13)

1

 1

.

2 dλ 2 [ 2 S = ∂ 2 λ 2 ∂ 2 T∂

 222222 energy F , the mean energy Ē = Tr(ρH) and the entropy S = -k B Tr(ρ ln ρ):F (T, λ) = -k B T ln Z(T, λ) (G1) F (T, λ) = Ē(T, λ) -T S(T, λ) (G2) ∂ T F (T, λ) = -S(T, λ). (G3)One now varies λ for a fixed entropy S. The temperature is thus a function T (λ) of λ such that S(T (λ), λ) = constant. The derivatives of the mean energy for fixed entropy are (d Ē dλ ) S ≡ d dλ [ Ē(T (λ), λ)] and ( d 2 Ē dλ 2 ) S ≡ d Ē(T (λ), λ)].Writing (G2) for T = T (λ) and taking the first order and the second order derivatives of the resulting equation with respect to λ, one findsd Ē dλ S = ∂ λ F (T (λ), λ) F (T (λ), λ) -[∂ T ∂ λ F (T (λ), λ)] F (T (λ), λ) (G5)It remains to use (G1) to obtain a microscopic expression of the above partial derivatives of F , from the partition function expressed as a sum Z = n e -βEn over the eigenstates n of the Hamiltonian:∂ λ F (T, λ) T ∂ λ F (T, λ) = Cov(E, dE/dλ) k B T 2 . (G9)Here the expectation value (. . .) stands for a sum over the eigenenergies with the canonical probability weights, and Var and Cov are the corresponding variance and covariance, e.g. Cov(E, dE/dλ)≡ n E n dEn dλ e -βEn Z

  /2) [s 2 -1 2l(l + 1)].

TABLE II

 II

: Relations for spin-1/2 fermions with zero-range interactions. The definition (1) of C, as well as the relations in lines 3, 5, 6 and 7, concern any (non-pathological) statistical mixture of states which satisfy the contact conditions [Tab. I, line 1]

we use the state |ψ b defined in Appendix C 2. As shown in that appendix, the expectation value of Ĉ taken in this state

  |ψ b tends to the contact C of |ψ [defined in Tab. II, Eq. (1)]. Moreover the expectation values of H -H trap and of nσ (k), taken in this state |ψ b , should tend to the corresponding expectation values taken in the state |ψ . This yields the desired relation.

					35)
	as shown in Appendix B. As in 3D, we deduce that
	dE d(ln a)	=	2π 2 m	(A, A),	(36)
	which gives the desired [Tab. II, Eq. (4b)] by us-
	ing [Tab. II, Eq. (2b)].
	D. Expression of the energy in terms of the
	momentum distribution
	Three dimensions:
	As shown in [95], the mean total energy E minus the
	mean trapping-potential energy E trap ≡ H trap , has the simple expression in terms of the momentum distribution
	given in [Tab. II, Eq. (5a)], for any pure state |ψ satis-fying the contact condition [Tab. I, Eq. (1a)]. We give
	a simple rederivation of this result by using the lattice
	model (defined in Sec. III B).
	We first treat the case where |ψ is an eigenstate of the zero-range model. Let |ψ b be the eigenstate of the lattice model that tends to |ψ for b → 0. We first note that C b ≡ ψ b | Ĉ|ψ b , where Ĉ is defined by [Tab. III, Eqs. (1a,1b)], tends to the contact C of the state ψ
	[defined in Tab. II, Eq. (1)] when b → 0, as shown in Appendix C. Then, the key step is to use [Tab. III,
	Eq. (3a)], which, after taking the expectation value in
	the state |ψ b , yields the desired [Tab. II, Eq. (5a)] in the zero-range limit since D → R 3 and ǫ k → 2 k 2 /(2m) for b → 0. To generalize [Tab. II, Eq. (5a)] to any pure state |ψ satisfying the contact condition [Tab. I, Eq. (1a)],

  mb 2 .(122)

	Two sets of parameters are possible. The first choice is
	α ≃ 1.470 885 and β ≃ -2.450 725,	(123)
	which leads to K 3 ≃ 3.137 788. The second choice is
	α ≃ -1.728 219 and β ≃ 12.838 540,	(124)
	which leads to K 3 ≃ 1.949 671. Other examples of magic dispersion relation can be found [145].

  The first deviations from unity are thus linear in q, and correspond to a shift of the breathing mode frequency ω breath to the new value 2ω + δω breath , that can be obtained to leading order in 1/a and r e from

	∂ω breath ∂(1/a)	= -	ω 4E 0	∂E 0 ∂(1/a)	and	∂ω breath ∂r e	=	3ω 4E 0	∂E 0 ∂r e
										7) 64s 2 + . . . (171)
									[∂E/∂r e ] q [∂E/∂r e ] 0	= 1 +	3q 4s	-	3q(5q + 11) 64s 2	+ . . . (172)

Acknowledgments

We thank E. Burovski, J. Dalibard, B. Derrida, V. Efimov, O. Goulko, R. Ignat, O. Juillet, D. Lee, S. Nascimbène, M. Olshanii, N. Prokof'ev, B. Svistunov, S. Tan, for useful discussions, as well as S. Giorgini and J. von Stecher for sending numerical data from [106,107,118,196]. The idea of using the range corrections present in Monte Carlo calculations to estimate the range corrections in cold atom experiments was given by W. Ketterle during the Enrico Fermi School in 2006. The question of a possible link between ∂E/∂r e and the 1/k 6 subleading tail of n σ (k) was asked by L. Platter during a talk given by F.W. at the INT workshop 10-46W in Seattle. The work of F.W. at UMass was supported by NSF under Grants No. PHY-0653183 and No. PHY-1005543. Our group at ENS is a member of IFRAF. We acknowledge support from ERC Project FERLODIM N.228177.

Note

[Tab. II, Eq. (4b)], as well as [Tab. II, Eq. (12b)], were obtained independently by Tan [197] using the formalism of [102]. After our preprint [88] appeared, some of our 2D relations were tested in [65] and some of them were rederived in [198].

Appendix A: Two-body scattering for the lattice model

For the lattice model defined in Sec. III B, we recall that φ(r) denotes the zero-energy two-body scattering state with the normalization [START_REF] Lee | [END_REF]10). In this Appendix we derive the relation (11,12) between the coupling constant g 0 and the scattering length, as well as the expressions (15,16,17,18) of φ(0). Some of the calculation resemble the ones in [15,199].

We consider a low-energy scattering state Φ q (r) of wavevector q ≪ b -1 and energy E = 2ǫ q ≃ 2 q 2 /m, i.e. the solution of the two-body Schrödinger equation (with the center of mass at rest):

where H 0 = D d d k/(2π) d 2ǫ k |k k| and V = g 0 |r = 0 r = 0|, with the asymptotic behavior Φ q (r) = r→∞ e iq•r + f q e iqr r + . . . in 3D (A2) Φ q (r) = r→∞ e iq•r + f q 2 iπqr e iqr + . . . in 2D.(A3)

Here f q is the scattering amplitude , which in the present case is independent of the direction of r as we will see. Note that, in 2D, the present definition corresponds to the convention (96), it differs e.g. from [200] by a factor 1/(4i). Also √ i ≡ e iπ/4 . We then have the well-known expression

where G ≡ (E + i0 + -H) -1 . Since G = G 0 + G 0 V G, with G 0 ≡ (E + i0 + -H 0 ) -1 , Eq. (A4) is equivalent to

where the T -matrix is T = V + V GV . Indeed, (A4) clearly solves (A1), and one can check [using the fact that r|G 0 |r = 0 behaves for r → ∞ as -m/(4π 2 ) e iqr /r in 3D and -(m/ 2 ) i/(8πqr)e iqr in 2D] that (A5) satisfies (A2,A3) with

In 3D the scattering length in defined by f q → q→0 -a, which gives the relation (11) between a and g 0 . In 2D,

where a is by definition the 2D scattering length. Identifying the inverse of the right-hand-sides of Eqs. (A7) and (A9) and taking the real part gives the desired (12). We note that Eqs. (A9,12) remain true if q → 0 is replaced by the limit b → 0 taken for fixed a.

To derive (15,16) we start from V |Φ q = T |q , which directly follows from (A4). Applying r = 0| on the left and using (A6,A7) yields

In 3D, we simply have φ = -a -1 lim q→0 Φ q [258], and the result (15) follows. In 2D, the situation is a bit more tricky because lim q→0 Φ q (0) = 0. We thus start with q > 0, and we will take the limit q → 0 later on. At finite q, we define φ q (r) as being proportional to Φ q (r), and normalized by imposing the same condition (10) than at zero energy, but only for b ≪ r ≪ q -1 . Inserting (A9) into (A11) gives an expression for Φ q (0). To deduce the value of φ(0), it remains to calculate the r-independent ratio φ q (r)/Φ q (r). But for r ≫ b we can replace φ q (r) and Φ q (r) by their values within the zero-range model (since we also have b ≪ q -1 ) which we denote by φ ZR q (r) and Φ ZR q (r). The two-body Schrödinger equation

where N is a constant and H

(1) 0 is an outgoing Hankel function. The contact condition

together with the known short-r expansion of the Hankel function [START_REF] Lebedev | Special functions and their applications[END_REF] then gives

Of course we also have Φ ZR q /φ ZR q = A, which gives (16). Finally, Eqs. (17,18) are obtained from (15,16) using the relations d(m/(4π 2 a))/d(1/g 0 ) = 1 in 3D and d(1/g 0 )/d(ln a) = -m/(2π 2 ) in 2D, which are direct consequences of the relations (11,12) between g 0 and a.

Appendix B: Derivation of a lemma

In this Appendix, we derive the lemma (33) in three dimensions, as well as its two-dimensional version [START_REF] Partridge | [END_REF]. Three dimensions: By definition we have

Here the notation ′ means that the integral is restricted to the set where none of the particle positions coincide [259]. We rewrite this as:

We note that this step is not trivial to justify mathematically. The order of integration has been changed and the limit ǫ → 0 has been exchanged with the integral over r i . We expect that this is valid in the presently considered case of equal mass fermions, and more generally provided the wavefunctions are sufficiently regular in the limit where several particles tend to each other. Since the integrand is the divergence of

1 , the divergence theorem gives

where the surface integral is for r i belonging to the sphere S ǫ (r j ) of center r j and radius ǫ, and the vector area dS points out of the sphere. We then expand the integrand by using the contact condition, in the limit

We derive the hierarchy [START_REF] Castin | Coherent atomic matter waves[END_REF][START_REF] Petrov | [END_REF]111) for a 2D nonpositive minimal-depth potential of finite range b, V (r) = 2 k 2 0 m v(r/b), for b → 0 and k 0 adjusted to have a constant s-wave scattering length a. The key point is then that k 0 b → 0 (differently from 3D).

In the s-wave channel, we write the zero-energy scattering wavefunction as

. This is integrated with the method of variation of constants,

Expressing that f 1 (x) ≃ β(+∞) ln(r/a) at infinity gives

and further using Eq. ( 106) leads to 109). Reproducing this perturbative expansion with the same v(x) in the l-wave, one gets

This relation for l = 1, combined with Eq. ( 110), gives Eq. (111). ber of equal-mass fermions, the energy is unbounded from below, is actually incorrect, since the fermionic antisymmetry was not properly taken into account. A theorem was published without proof in [START_REF] Minlos | On pointlike interaction between n fermions and another particle[END_REF] implying that the spectrum of the Hamiltonian of N ↑ same-spin-state fermions of mass m ↑ interacting with a distinguishable particle of mass m ↓ is unbounded below, not only for m ↑ = m ↓ and large enough N ↑ , but also for N ↑ = 3 and m ↑ /m ↓ larger than the critical mass ratio 5.29 . . . . No proof was found yet for this theorem; it was only proven that no Efimov effect occurs for N ↑ = 3, N ↓ = 1 provided m ↑ /m ↓ is sufficiently small [START_REF] Minlos | [END_REF]. It was recently shown that a four-body Efimov effect occurs in this 3+1 body problem (for an angular momentum l = 1 and not for any other l ≤ 10) and makes the spectrum unbounded below, however for a widely different critical mass ratio m ↑ /m ↓ ≃ 13.384 [78], which sheds some doubts on [START_REF] Minlos | On pointlike interaction between n fermions and another particle[END_REF].

[216] The corresponding state vector is

where there are N ↑ spins ↑ and N ↓ spins ↓, and the operator  antisymmetrizes with respect to all particles. The wavefunction ψ(r1, . . . , rN ) is then proportional to ( ↑, . . . , ↑, ↓, . . . , ↓ | ⊗ r1, . . . , rN |) |Ψ , with the proportionality factor

[218] For purely attractive interaction potentials such as the square-well potential, above a critical particle number, the ground state is a collapsed state and the zero-range regime can only be reached for certain excited states (see e.g. [START_REF] Castin | The BCS-BEC Crossover and the Unitary Fermi Gas[END_REF] and refs. therein).

[219] The existence of the 1/k 4 tail had already been observed within a self-consistent approximate theory [205].

[220] E.g. for n ↓ (k) in the trapped three-body case, with particles 1 and 2 in state ↑ and particle 3 in state ↓, one has i = 3 and j, j ′ = 1 or 2. Then the crossed term A31(r1, r2)A32(r2, r1) has to all orders finite derivatives with respect to r1 and r2, except if r1 = r2 where it vanishes as |r1 -r2| 2s-2 , s > 0 not integer, see e.g. Eq. (H3) and below that equation. By a power counting argument, its Fourier transform with respect to r1 -r2 contributes to the momentum distribution tail as 1/k 2s+5 = o(1/k 4 ); one recovers the "three-closeparticle" contribution mentioned in a note of [96].

[221] For simplicity, we refrain here from expressing C as the integral of a "contact density" C(R) related to the small-r behavior of the local pair distribution function g

↑↓ (R + r/2, Rr/2) as was done for the 3D case in [95][96][97]; this C(R) is then also related to the large-k tail of the Wigner distribution [i.e. the Fourier transform with respect to r of the one-body density matrix ψ † σ (Rr/2)ψσ(R + r/2) ], see Eq. ( 30) of [95].

[222] Our derivation is similar to the one given in the twobody case and sketched in the many-body case in Section 3 of [96].

[223] This relation was written in [102] in a form containing a generalised function η(k) (i.e. a distribution). We have checked that this form is equivalent to our Eq. ( 38), using Eq. (16b) of [102],

. This last property is implied in Eq. (16a) in [102].

[224] These last relations also hold if one averages over all directions of r uniformly on the unit sphere or unit circle.

[225] Our result does not follow from the well-known fact that, for a finite-range interaction potential in continuous space, -

σσ (r = 0) equals the kinetic energy; indeed, the Laplacian does not commute with the zero-range limit in that case [cf. also the comment below Eq. ( 180)].

[226] As suggested by a referee, [Tab. II, Eq. (7b)] can be tested for the dimer wavefunction ψ(r1, r2) = φ dim (r12) = -κK0(κr)/π 1/2 [111], which has the energy E = -2 κ 2 /m and the momentum distribution nσ(k) = 4πκ 2 /(k 2 + κ 2 ) 2 , where κ = 2/(ae γ ) and K0 is a Bessel function. From Eq. ( 42) we find G

σσ (r) = κrK1(κr). From C/(4π) = -mE/ 2 = κ 2 and the known expansion of K1 around zero, we get the same low-r expansion as in [Tab. II, Eq. (7b)]. To calculate G

(1) σσ (r), we used the fact that K0(κr) is the 2D Fourier transform of 2π/(k 2 + κ 2 ): it remains to take the derivative with respect to κ and to realize that K ′ 0 = -K1.

[227] In the lattice model in 3D, the coupling constant g0 is always negative in the zero-range limit |a| ≫ b, and is an increasing function of -1/a, as seen from ( 11).

[228] Eq. (32) shows that dE0/d(ln a) is proportional to the probability of finding two particles very close to each other, and it is natural that this probability decreases when one goes from the BEC limit (ln a → -∞) to the BCS limit (ln a → +∞), i.e. when the interactions become less attractive [in the lattice model in 2D, the coupling constant g0 is always negative in the zero-range limit a ≫ b, and is an increasing function of ln a, as can be seen from ( 12)].

[229] We assume, to facilitate the derivation, that V (r) = 0 for r > b, but the result is expected to hold for any V (r) which vanishes quickly enough at infinity. [230] We consider two particles of opposite spin in a cubic box of side L with periodic boundary conditions, and we work in the limit where L is much larger than |a| and b. In this limit, there exists a "weakly interacting" stationary state ψ whose energy is given by the "mean-field" shift E = g/L 3 with g = 4π 2 a/m. The Hellmann-Feynman theorem gives g0 dE/dg0 = Eint[ψ].

But the wavefunction ψ(r1, r2) ≃ Φ(r12)/L 3 where Φ is the zero-energy scattering state normalized by Φ → 1 at infinity. Thus Eint = d 3 r V (r)|Φ(r)| 2 /L 3 . The desired Eq. ( 70) then follows, since Φ = -aφ.

[231] More precisely, one first takes a general, nonrotationally invariant function χ(r), that one then expands in partial waves of angular momentum l, that is in spherical harmonics. Performing the reasoning to come for each l, one finds at the end that the l = 0 channel finite range correction dominates for small b, in the absence of l-wave resonance for l = 0. [232] Since E depends on Rij and the (r k ) k =i,j , χ actually depends on these variables and not only on rij. This dependence however rapidly vanishes in the limit b → 0, if one restricts to the distances rij b, for the normalization ( 89):

[234] In general, when N ↑ ≥ 2 and N ↓ ≥ 2, the functions Aij have 1/r kl divergences when r kl → 0. This is apparent in the dimer-dimer scattering problem [206]. As a consequence, in the integral of [Tab. V, Eq. ( 1a)], one has to exclude the manifold where at least two particles are at the same location. The same exclusion has to be performed in 2D [235] The wavefunction is not an analytic function of r for a compact support interaction potential, since a non-zero compact support function is not analytic.

[236] We consider here a truly 2D gas. In experiments, quasi-2D gases are produced by freezing the z motion in a harmonic oscillator ground state of size az = [ /(mωz)] 1/2 : At zero temperature, a 2D character appears for 2 k 2 F /(2m) ≪ ωz. From the quasi-2D scattering amplitude given in [207] (see also [208]) we find the effective range squared, r 2 e = -(ln 2) a 2 z . Anticipating on subsection VII B we also find ρe = R1 = 0. It would be interesting to see if the finite range energy corrections dominate over the corrections due to the 3D nature of the gas, both effects being controlled by the same small parameter (kF re) 2 .

[237] As in 3D one may also be worried by the dependence of χ with Rij and the (r k ) k =i,j via its dependence with the energy E . We reach the estimate ∂E χ(b)/χ(b) ≈ mr 2 e /[ 2 ln(a/b)] that vanishes more rapidly than r 2 e in the zero-range limit.

[238] We have checked that the hypothesis of a non-resonant interaction in [START_REF] Gibson | [END_REF] is actually not necessary to obtain (C16) and (C18) of that reference, that lead to (104).

[239] This term is obtained by distinguishing three integration zones before taking the limit Kx → 0, so as to fold back the vectors q ± 1 2 K inside the first Brillouin zone: the left zone

where ǫ q+ 1 2 K is written as ǫ q-2π b ex+ 1 2 K , and the central zone. The surface term can also be obtained by interpreting ∂ 2 qx in the sense of distributions, after having shifted the integration domain D by π b ex for mathematical convenience. The second order derivative in the first term of Eq. ( 116) is of course taken in the sense of functions.

[240] The integration can be performed in spherical coordinates of polar axis the direction of K.

[241] This problem does not show up in recent studies of the fermionic polaron problem [209,210] since the momentum cut-off is introduced only for the majority atoms and not for the impurity, see [211].

For the finite number low energy terms, we directly use this fact. For the other terms, such that ǫ k ≫ |E| and ≫ (2π ) 2 /(mL 2 ), we use

which is integrable at large k in 3D and leads to a total error O(b 2 ).

[243] The contribution proportional to re in Eq. ( 134) can also be obtained from [Tab. V, Eq. (1a)] and from the fact that k =0 e ik•r /k 2 ∼ L 3 /(4πr) for r → 0.

[244] From Schrödinger's equation, ∆r ij ψ diverges at most as ψ itself, that is as ln rij, for rij → 0. The particular solution f (r) = 1 4 r 2 (ln r -1) of ∆rf (r) = ln r fixes the form of the subleading term in ψ.

[245] In 3D we used the identity d 3 re ik•r /r = 4π/k 2 and its derivatives with respect to kα; e.g. taking the Laplacian with respect to k gives d 3 r e ik•r r = -8π/k 4 . Equivalently, one can use the relation

) and its derivatives with respect to kα; this relations holds for any u(r) which has a series expansion in r = 0 and rapidly decreases at ∞. In 2D for k > 0 we used the identity d 2 re ik•r ln r = -2π/k 2 and its derivatives with respect to kα. The regular terms involving L (α) ij have (as expected) a negligible contribution to the tail of nσ(k).

[246] The configurations with three close particles contribute to the tail of nσ(k) as 1/k 5+2s , see a note of [96], with s defined in Sec. X B, which is negligible for s > 1/2. [247] Similarly, a "contact current" was recently introduced in [212], whose spatial integral is proportional to (A, ∇RA). [251] In [213] a similar calculation was performed, except that the harmonic trap was neglected within the interaction range.

[252] For Efimovian eigenstates, computing the derivative of the energy with respect to the effective range would require to use a regularisation procedure similar to the one employed in free space in [122,124]. However the derivative with respect to 1/a can be computed [START_REF] Werner | [END_REF].

[253] Here we used the value of the effective range re = 1.435 r0 [196] for the Gaussian interaction potential V (r) = -V0e -r 2 /r 2 0 with V0 equal to the value where the first two-body bound state appears.

[254] The hydrodynamic frequencies Ω are given by the eigenvalue problem -mΩ 2 δρ = div [ρ0 ∇ (µ ′ hom [ρ0] δρ)] where δρ(r) is the infinitesimal deviation from the stationary density profile ρ0(r), µ hom [ρ] is the ground state chemical potential of the homogeneous gas of density ρ and the appex ′ indicates derivation. For the equation of state µ hom [ρ] = Aρ 2/3 + Bρ γ , where B is arbitrarily small, we treat the term in B to first order in perturbation theory around the breathing mode to obtain Ω = 2ω + ω 96 π (γ -2 3 ) B µ µ A 3γ/2 1 0 duu 2 (1 -2u 2 )(1 -u 2 ) (3γ+1)/2 where µ = ωN 1/3 (2mA/π 4/3 ) 1/2 is the unperturbed chemical potential of the trapped gas. To zeroth order in B, scaling invariance gives δρ (0) (r) = d dλ [ρ0(r/λ)/λ 3 ] λ=1 . To use perturbation theory, we made the differential operator Hermitian with the change of function δf (r) = (µ ′ hom [ρ0(r)]) 1/2 δρ(r). Hermiticity of the perturbation is guaranteed (i.e. surface terms coming from the divergence theorem vanish) for γ larger than 1/3. For finite-re corrections, γ = 1.

[255] For a finite-range potential one has g (1) σσ (r) = n/2r 2 mE kin /(3 2 V)+. . . where V is the volume; the kinetic energy diverges in the zero-range limit as E kin ∼ -Eint, thus E kin ∼ -C/(4π) 2 d 3 r V (r)|φ(r)| 2 from [Tab. IV, Eq. (2a)], so that E kin ∼ Cπ 2 /(32mb) for the squarewell interaction. This behavior of g (1) (r) only holds at very short distance r ≪ b and is below the resolution of the Monte Carlo data.

[256] The general structure of Eq. ( 185) already appeared for a simple separable two-channel model [100] with exactly the same expression for the first term, which explains why the a-dependence is correctly reproduced by the simple expression of [100], as observed in [214] by com-