F Werner

L Tarruell

Y Castin

Number of closed-channel molecules in the BEC-BCS crossover

Keywords: 03, 75, Ss Degenerate Fermi gases -67, 90, +z Other topics in quantum fluids and solids

Using a two-channel model, we show that the number of closed-channel molecules in a twocomponent Fermi gas close to a Feshbach resonance is directly related to the derivative of the energy of the gas with respect to the inverse scattering length. We extract this quantity from the fixed-node Monte Carlo equation of state and we compare to the number of closed-channel molecules measured in the Rice experiment with lithium [Partridge et al., Phys. Rev. Lett. 95, 020404 (2005)]. We also discuss the effect of a difference between the trapping potentials seen by a closed-channel molecule and by an open-channel pair of atoms in terms of an effective position-dependent scattering length.

Introduction

It is now experimentally possible to prepare two-component atomic Fermi gases at low temperature with a fully controlable interaction strength, in a regime where the interaction range is negligible and the atomic interactions are characterized by a single parameter, the s-wave scattering length a between two opposite "spin" atoms. The value of a can be adjusted at will thanks to a magnetically induced Feshbach resonance. The weakly interacting limits k F a → 0 + and k F a → 0 -, where k F is the Fermi wavevector of the gas, correspond respectively to the BEC limit (a Bose-Einstein condensate of dimers, observed in [1,2,3,4,5]) and the BCS limit (a "condensate" of Cooper pairs approximately described by the BCS theory). Furthermore, the experiments can access the so-called crossover regime between BEC and BCS, where the gas is strongly interacting (k F |a| 1), which includes the celebrated unitary limit k F |a| = ∞, where the gas acquires fully universal properties [6,7,8,9,10,11,12,13,14,15,16,17].

In this context, early theoretical studies [18,19,20] put forward a many-body Hamiltonian which accurately models the microscopic two-body physics of the Feshbach resonance. The interaction potential is described by two channels, an open channel and a closed channel. The atoms exist in the form of fermionic particles in the open channel and in the form of bosonic short range molecules in the closed channel. Two atoms may be converted into a short range closed-channel molecule and vice versa due to the interchannel coupling. The corresponding Hamiltonian is

H 2 = d 3 r    σ=↑,↓ ψ † σ (r) - 2 2m ∆ + U (r) ψ σ (
δ r 1 + r 2 2 - r 3 + r 4 2 ψ † ↑ (r 1 )ψ † ↓ (r 2 )ψ ↓ (r 3 )ψ ↑ (r 4 ) (1)
where the atomic fields ψ σ (r) obey the usual fermionic anticommutation relations and the field ψ b (r) describing the closed-channel molecules obeys the usual bosonic commutation relations. Λ gives the coupling between the closed and open channels, with a short range cut-off function χ(r), of range b, and such that d 3 r χ(r) = 1. A typical value for b is in the nanometer range 1 . The coupling constant g 0 represents the background atomic interaction in the open channel, which is conveniently modeled by a separable potential with the same cut-off function χ. The energy E b is the energy of a closed-channel molecule, in the absence of the coupling Λ, counted with respect to the dissociation limit of the open channel. This energy is adjusted with a magnetic field B, using the fact that the different magnetic moments in the open and closed channels lead to a differential Zeeman shift. The resulting effective magnetic moment is

µ b ≡ dE b dB . (2)
U (r) and U b (r) are the trapping potentials experienced by the atoms and the closed-channel molecules respectively. On the contrary, recent Monte Carlo simulations of the many-body problem in the BEC-BCS crossover simply use single channel model Hamiltonians [22,23,24,25,26,27,28]. In most unbiased Quantum Monte Carlo simulations, a lattice model is used:

H 1 = r b 3    σ=↑,↓ ψ † σ (r) - 2 2m ∆ + U (r) ψ σ (r) +g 0 ψ † ↑ (r)ψ † ↓ (r)ψ ↓ (r)ψ ↑ (r) . ( 3
)
Here b is the lattice spacing (which coincides with the "range" of the discrete delta interaction potential), ∆ is a discrete representation of the Laplacian, and g 0 is a coupling constant adjusted to have the desired scattering length: As was shown e.g. in [24,25,[START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School Quantum Gases in Low Dimensions[END_REF][START_REF] Pricoupenko | [END_REF],

1 g - 1 g 0 = [-π/b,π/b] 3 d 3 k (2π) 3 1 2ǫ k (4)
where g = 4π2 a/m is the effective coupling constant and ǫ k the energy of a single particle with wavevector k on the lattice. In fixed-node Monte Carlo calculations, one rather uses an interaction potential in continuous real space, e.g. a square well potential [31,32,33,34,35,36].

There is now a consensus about the fact that, in the zero-range limit, there is a convergence of predictions for the two models H 1 and H 2 , in particular for the equation of state of the gas. The zero-range limit means, for single channel models, that k F b ≪ 1, b ≪ λ, and b ≪ |a|, where k F is the Fermi wavevector and λ is the thermal de Broglie wavelength. 2 For two-channel models, one has to impose the same conditions, not only for the interaction range b but also for the effective range r e (one indeed has |r e | ≫ b on narrow Feshbach resonances [37,[START_REF] Castin | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF]); we shall see in appendix A that one also needs to impose a condition on the background scattering length, and that all conditions are satisfied in the experiment [5]. 3 One can however consider observables which exist in reality and are included in the two-channel model, but which are absent in a single channel model. The most natural example is the number of closed-channel molecules, which was recently measured at Rice using laser molecular excitation techniques [5].

In this paper, we show that the equilibrium number of closed-channel molecules N b is directly related to the equation of state of the gas, more precisely to the derivative of the gas energy (or free energy at non-zero temperature) with respect to 1/a. Paradoxically, we can thus use the equation of state calculated in [32] within a single channel model, in order to predict N b in the crossover and compare to Rice measurements, see section 2. Since the derivative of the energy with respect to 1/a is related to other observables involving atomic properties only [START_REF] Tan | [END_REF]40], it is also possible in principle to access N b by a pure atomic measurement, see section 3. We conclude in section 4.

2 Prediction for the number of closed-channel molecules

General result

We first suppose that the system is in an arbitrary eigenstate |ψ of the two-channel Hamiltonian H 2 , with eigenergy E. From Hellmann-Feynman theorem we have

dE dB = ψ| dH 2 dB |ψ . (5)
The only magnetic field dependent part of the Hamiltonian H 2 is the bare closed-channel molecule energy E b (B) so that

dE dB = µ b N b (6)
where

N b = d 3 r ψ † b (r)ψ b (r) (7)
is the mean number of closed-channel molecules.

We then eliminate the magnetic field by using as a parameter the scattering length a experienced by the trapped atoms rather than B. In the general case where the trapping potential U b (r) for the closed-channel molecule differs from twice the trapping potential U (r) experienced by an open-channel atom, the scattering length of two atoms around point r depends on r, see appendix B. To lift the ambiguity we thus take for a the scattering length in the trap center taken as the origin of coordinates, r = 0. We then rewrite (6) as

N b = dE d(1/a) d(1/a) dB 1 µ b . (8)
Assuming that µ b is magnetic field independent, a as a function of B can be calculated explicitly for H 2 by solving the zero-energy two-body scattering problem in free space (see appendix A), and by including the energy shifts due to the trapping potentials U b (0) and U (0), as explained in the appendix B. The function takes the form

a = a fs (B -δB 0 ) (9)
where we have introduced the function giving the free space scattering length as a function of B,

a fs (B) = a bg 1 - ∆B B -B 0 . ( 10
)
Here a bg is the background scattering length, B 0 is the location of the Feshbach resonance in the absence of trapping potential, δB 0 = [2U (0) -U b (0)]/µ b the shift of the resonance location due to the trapping potential, and ∆B the resonance width. Equation ( 8) is directly applicable at zero temperature. At non-zero temperature, one has to take a thermal average of (8), keeping in mind that the mean value of a derivative is not necessarily equal to the derivative of the mean value. In the canonical ensemble one can check the exact relation [41]

dE d(1/a) = dF d(1/a) T = d E d(1/a) S (11)
where . . . is the thermal average, the derivative of the free energy F is taken for a fixed temperature T , and the derivative of E is taken for a fixed entropy S. We thus have in the canonical ensemble

N b = d E d(1/a) S d(1/a) dB 1 µ b . (12)

Analytical results in a trap in limiting cases

We now restrict to a spin balanced gas where the number of particles is equal to N/2 in each spin component. We also assume that the gas is at zero temperature (or at temperatures much smaller than the Fermi temperature T F ), in a harmonic trap

U (r) = U (0) + 1 2 m α=x,y,z ω 2 α r 2 α , (13)
in the macroscopic regime where k B T F is much larger than the oscillation quanta ω α . The usual definition of the Fermi temperature T F and of the Fermi wavevector k F is then

k B T F = 2 k 2 F 2m = (3N ) 1/3 ω ( 14
)
where ω is the geometric mean of the three oscillation frequencies ω α . In this macroscopic regime, we rely on the local density approximation to calculate dE/d(1/a) for the trapped gas: This is expected to be exact in the large N limit, and to already hold for the Rice experiment [5] since

k B T F /( ω α ) 9.
The key ingredient of the local density approximation is the equation of state of the interacting homogeneous gas. At zero temperature, and in the zero-range limit detailed in the introduction, the equation of state is expected to be a universal function of the gas density and of the scattering length, independent of the microscopic details of the interaction. Furthermore, the spatial dependence of the scattering length in presence of the trapping potential, due to U b (r) = 2 U (r), is expected to be negligible in the zero-range limit, as argued in the appendix B. Under these assumptions, the local-density approximation for the trapped gas leads to the universal form

- dE d(1/a) = 2 k F m N F 1 k F a , ( 15
)
where F is a dimensionless function, expressed in the appendix C in terms of the dimensionless universal function f giving the energy per particle ǫ hom of the homogeneous interacting gas:

ǫ hom = ǫ hom 0 f 1 k hom F a . ( 16
)
Here ǫ hom 0 and k hom F are the energy per particle and the Fermi wavevector of the ideal Fermi gas with the same total density n as the interacting gas:

k hom F = (3π 2 n) 1/3 ( 17
)
ǫ hom 0 = 3 10
2 k hom F 2 m . (18)
In the BEC regime 0 < k F a ≪ 1 we have [42,43,44]

f 1 k hom F a = - 5 3(k hom F a) 2 (19)
+ 5 18π k hom F a d 1 + 128 15 √ 6π 3 (k hom F a d ) 3/2 + . . .
where the first term comes from the dimer binding energy 2 /(ma 2 ), the second term is the mean field interaction energy among the dimers, proportional to the dimer-dimer scattering length a d , and the last term is the (bosonic) Lee-Huang-Yang correction. The solution of the 4-fermion problem gives [45,46]

a d ≃ 0.60a. ( 20
)
The local-density approximation then leads to

F 1 k F a = 1 k F a + 5 2/5 2 12/5 • 7 a d a 2/5
(k F a) 7/5 + . . . (21) where we have omitted the Lee-Huang-Yang correction.

In the BCS limit 0 < -k F a ≪ 1, we have [43,44,47,48]

f 1 k hom F a = 1 + 10 9π k hom F a + 4(11 -2 ln 2) 21π 2 (k hom F a) 2 + . . . ( 22
)
where the first term is the ideal gas result, the second term is the Hartree-Fock mean field term, and the last term is the (fermionic) Lee-Huang-Yang correction. This leads to 4

F 1 k F a = 512 945π 2 (k F a) 2 1+ 256 35π 2 - 63 + 189 ln 2 1024 k F a + . . . . ( 23
)
Around the unitary limit 1/(k F a) → 0 we use the expansion

f 1 k hom F a = ξ - ζ k hom F a + . . . , (24)
which gives for the trapped gas:

F (0) = 2 7 5 2 • 7 • π ζ ξ 1/4 ≃ 0.28 (25)
where we took the estimates

ξ ≃ 0.41 (26) ζ ≃ 0.95 . ( 27
)
The estimate for ξ is obtained by averaging the predictions of [32] (ξ = 0.42(1)), of [26] (ξ = 0.449(9)), of [33] (ξ = 0.42(1)) and of [27] (ξ = 0.37 (5)). The estimate for ζ is obtained from the calculation of the short range pair correlation of [34], see subsection 3.1, and is close to the value ζ ≃ 1.0 that may be extracted directly from the values of f in [32] by approximating the derivative by a two-point formula. 5

Comparison to Rice experiment

At Rice [5] the quantity Z = N b /(N/2) was measured for a lithium gas, by resonantly exciting the closed-channel molecules with a laser that transfers them to another, short-lived molecular state; the molecule depletion reflects into a reduction of the atom number that can be measured. We now compare the Rice results to our prediction (8). We first insert (9,10) into (8) which gives

N b = N k F R * F 1 k F a 1 - a bg a 2 (28)
4

For a small positive scattering length 0 < kF a ≪ 1, Eqs. (22,23) also apply to the atomic gas state; since the interaction potentials considered here have a two-body bound state for a > 0, this atomic gas state is only metastable [49], the ground state being the BEC of dimers considered previously [cf. Eq.( 21)]. 5 The values (26,27) of ξ and ζ can be inserted into the expressions [50,51] for the derivatives of the frequencies of the hydrodynamic collective modes with respect to 1/(kF a) taken at unitarity. The resulting slope is compatible with the experimental data for the radial mode [12,52], while the agreement with the data for the axial mode [53] is only marginal. The data from [12] agree very well for all positive values of a with the theoretical result [54] obtained from the fixed-node Monte Carlo equation of state [32], while in [52,53] finite-temperature effects may play a role [12,54,55]. (19,22) in the weakly interacting regimes with an interpolation of the fixed-node Monte Carlo results of [32,34] in the strongly interacting regime (see appendix D). Dashed lines: Analytical predictions in the weakly interacting regimes, on the BEC side (21) and on the BCS side (23). Cross: Result (25) at unitarity. Circles (with error bars): Experimental measurement of the number N b of closed-channel molecules in a lithium gas at Rice [5], combined with the present theory linking dE/d(1/a) to N b , see (28). where the function F is defined in (15). The length R * depends on the interchannel coupling and is related to the width of the Feshbach resonance as:

R * = 2 ma bg µ b ∆B . ( 29
)
For the B 0 ≃ 834G Feshbach resonance in 6 Li, we have µ b ≃ 2µ B where µ B is the Bohr magneton [21], ∆B = -300G, a bg = -1405a 0 where a 0 = 0.0529177 nm is the Bohr radius [56]; this gives R * = 0.0269 nm. From the Rice data for N b /N , using (28), we thus get values of F , i.e. of dE/d(1/a), that we compare in Fig. 1 to theoretical predictions, given in the local-density approximation either by the expansions (21,23) in the weakly interacting regime, or by (25) at the unitary limit, or by an interpolation of the Monte Carlo results of [32,34] in the strongly interacting regime (see appendix D for details) 6 . The agreement is satisfactory, the experimental points being scattered around the theoretical prediction. In particular, a possible shift δB 0 of the resonance location due to the fact that U b (r) = 2U (r), see ( 9), an effect not explicitly included in [5] for the calculation of 1/(k F a) as a function of B, and which would manifest itself in horizontal shifts of the experimental data with respect to theory, does not show up in Fig. 1. For the considered broad resonance of lithium, this is not surprising, see appendix B.

Our approach even allows to predict the time evolution of the particle number, in presence of the depleting laser, under the assumption that the molecular depletion rate is so weak that the equilibrium relation ( 28) holds at all times. Neglecting the slow off-resonant free-bound photoassociation process [5], we have

dN dt = -2N b (t) Ω 2 γ ( 30
)
where Ω 2 /γ is the effective decay rate of the closed-channel molecule, the Rabi frequency Ω and the spontaneous emission rate γ being defined in [5]. 7 Simple expressions can be obtained in limiting cases: In the BEC limit, the dominant term of the expansion (21) gives

N (t) = N (0) e -Γ0t ; (31)
in the unitary limit, (25) gives

N (t) = N (0) (1 + Γ 0 t/6) 6 ; (32)
and in the BCS limit the dominant term of (23) gives

N (t) = N (0) (1 + Γ 0 t/2) 2 . (33)
Here, Γ 0 = -[dN (0)/dt]/N (0) is the inital atom loss rate, which can be expressed in terms of the initial atom number N (0) using Eqs. (28,[START_REF] Pricoupenko | [END_REF]. We see that in general N (t) is 6 The error bars in Fig. 1 include the systematic theoretical uncertainty resulting from the leading order correction to (10) given in [56], but are largely dominated by the experimental error bars from [5] because we did not include the three experimental points from [5] which are deep in the BEC regime [1/(kF a) > 5]; in this regime, two-body theories which go beyond our two-channel model give good agreement with the Rice measurements [5,21,57]. 7 Eq. ( 30) can be obtained from the following master equation for the density matrix ρ of the gas: [START_REF] Pricoupenko | [END_REF] stands for N f + 2N b , where N f is the number of fermions; thus N ≃ N f in the regime considered in this paper. not an exponential function of time. For the last point on the BCS side in [5], the (unpublished) data for N (t) are better fitted by Eq.( 33) than by an exponential; however the value of Γ 0 resulting from this new fitting procedure does not differ significantly from the published value [58].

dρ/dt = (Ω 2 /γ) R d 3 r h ψ b (r)ρψ † b (r) -1 2 {ψ † b (r)ψ b (r), ρ} i + [H2, ρ]/(i ). N in

Related observables

In this section, contrarily to subsections 2.2 and 2.3, we return to the general case and we do not assume that the temperature is ≪ T F and that the numbers of atoms in each spin component are equal. For concreteness we take the lattice model ( 3), even if the reasonings may be generalized to a continuous space model and to a twochannel model.

Short range pair correlations

As was shown in [START_REF] Tan | [END_REF]40] and recently rederived in [59], in the zero interaction range limit, the derivative of the gas energy with respect to the inverse scattering length can be expressed in terms of the short range behavior of the pair distribution function of opposite spin particles:

- d E d(1/a) S = 4π 2 m d 3 R lim r→0 r 2 g ↑↓ (R + r/2, R -r/2) (34)
where the pair distribution function is given by

g ↑↓ (r 1 , r 2 ) = ψ † ↑ (r 1 )ψ † ↓ (r 2 )ψ ↓ (r 2 )ψ ↑ (r 1 ) . ( 35
)
Since we are dealing here with a finite range interaction model, we have taken a zero-range limit in (34), introducing the notation lim

r→0 = lim r→0 lim b→0 . (36)
Note that the order of the limits is important. The limit in ( 34) is reached for a distance r much smaller than |a|, the mean distance between particles, and the thermal de Broglie wavelength λ, but still much larger than the interaction range b.

We thus see that a measurement of the pair distribution function of the atoms gives access, using (8), to the number of closed-channel molecules. Experimentally, the pair distribution function was measured via the photoassociation rate in a 1D Bose gas [60,61], or simply by dropping the gas on a detector with high enough spatiotemporal resolution [62,63,64]. The possibility of measuring the pair distribution function in the BEC-BCS crossover was studied in [34,65]. In the unitary limit, by inserting into Eq.( 34) the value of lim r 2 g ↑↓ calculated in [34] with the fixed-node Monte Carlo technique, we get the value (27) for the parameter ζ defined in (24).

The relation (34) is the three-dimensional version of the relation obtained in [66] for a one-dimensional Bose gas with contact interaction, using the Hellmann-Feynman theorem. We now show that the Hellmann-Feynman technique also provides a simple derivation of (34). Taking the derivative of the gas eigenergy E with respect to the effective coupling constant g = 4π 2 a/m, for a fixed lattice spacing b, the Hellmann-Feynman theorem gives

d E dg S = dg 0 dg r b 3 ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ (37) = dg 0 dg r b 3 g ↑↓ (r, r), (38)
where g 0 is a function of g as given in ( 4), and where we used the relation dE/dg = (d E /dg) S [cp. (11)]. In a gas, with an interaction range much smaller than the mean distance ∼ 1/k F between particles and than the thermal de Broglie wavelength λ, the pair distribution function is dominated by two-body physics if |r 1 -r 2 | is much smaller than 1/k F and λ (and also much smaller than a in the BEC regime), a property contained in the Yvon [67] and Waldmann-Snider [68] ansatz of kinetic theory, which also appeared in the context of quantum gases [34,[START_REF] Huang | Statistical Mechanics[END_REF][START_REF] Holzmann | [END_REF][START_REF] Leggett | Lecture 2 at the international quantum gases workshop[END_REF]:

g ↑↓ (r 1 , r 2 ) ≃ C r 1 + r 2 2 |φ(r 1 -r 2 )| 2 ( 39
)
where φ is the two-body zero-energy free space scattering state for two atoms, here in the spin singlet state. With the lattice model we find

8 |φ(0)| 2 = dg dg 0 , (40)
so that Eq.( 38) simplifies to

d E dg S ≃ r b 3 C(r). (41)
On the other hand, at interparticle distances much larger than b, the zero-energy scattering wavefunction φ behaves as 1

-a/|r 1 -r 2 | so that lim r→0 r 2 g ↑↓ (R + r/2, R -r/2) = a 2 C(R). (42)
Eq.( 41) thus leads to (34).

Tail of the momentum distribution

As was shown in [40] and recently rederived in [59], in the zero-range limit, the momentum distribution of the gas has a large momentum tail scaling as 1/k 4 , with a

8 An elegant derivation is to enclose two atoms in a fictitious cubic box with periodic boundary conditions of volume L 3 , with L ≫ |a|, b. To leading order in a/L, the wavefunction is φ box (r1 -r2) = φ(r1 -r2)/L 3/2 with an energy given by the "mean-field" shift E box = g/L 3 . The Hellmann-Feynman theorem (37) then gives the result (40). An alternative derivation is to directly calculate φ in Fourier space [see (47)] which gives

φ(r) = 1 -g0φ(0) R [-π/b,π/b] 3 [d 3 k/(2π) 3 ] exp(ik • r)/(2ǫ k )
; this, combined with (4), gives (40).

coefficient common to the two spin states and proportional to -dE/d(1/a):

lim k→+∞ k 4 n σ (k) = - 4πm 2 d E d(1/a) S (43)
where lim

k→+∞ = lim k→+∞ lim b→0 , (44)
σ =↑ or ↓ and the momentum distribution is normalized as

d 3 k (2π) 3 n σ (k) = N σ , (45)
where N σ is the mean number of particle in the spin state σ. When combined with (8), this relation reveals that the number of closed-channel molecules may be deduced from a measurement of the momentum distribution of the atoms.

The momentum distribution in a two-component Fermi gas was actually already measured [START_REF] Regal | [END_REF][START_REF] Tarruell | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF][START_REF] Chen | [END_REF], simply by switching off the interaction and measuring the cloud density after a ballistic expansion. Two practical problems however arise, the effect of the non-zero switch-off time of the interaction [START_REF] Regal | [END_REF]75] and the signal to noise ratio in the far tails of the distribution; we have therefore not been able to extract the value of dE/d(1/a) from the published data for n σ (k). The value of dE/d(1/a) has not been extracted either from the fixed-node Monte Carlo calculations of n σ (k) performed in [76].

The relation ( 43) is a three-dimensional version of the relation derived for a one-dimensional Bose gas in [77]. We now present a very simple rederivation of (43).

In a gas with short range interactions, i.e. b ≪ 1/k F and b ≪ λ (and b ≪ a on the BEC side), the tail of the momentum distribution is dominated by binary collisions in the spin singlet channel between a particle with a large momentum k and a particle with a momentum ≃ -k, so that at large k n

σ (k) ≃ B| φ(k)| 2 (46)
where B is momentum-independent and φ(k) is the Fourier transform r b 3 e -ik•r φ(r) of the zero-energy two-body free space scattering state φ(r). For the contact interaction in the lattice model the two-particle Schrödinger's equation takes the form [START_REF] Pricoupenko | [END_REF]

2ǫ k φ(k) + g 0 φ(r = 0) = 0, ( 47
)
where ǫ k is the free wave dispersion relation in the lattice model, so that φ(k) drops as 1/ǫ k for large k, and the associated kinetic energy diverges for b → 0. Eq.( 46) implies that the kinetic energy of the gas is proportional to the kinetic energy of the zero-energy two-body scattering state in the zero-range limit:

E kin [gas] ≃ B • E kin [φ]. ( 48
)
Since the total energies have a finite limit for b → 0, we can replace the kinetic energies by the opposite of the interaction energies in the above formula:

E int [gas] ≃ B • E int [φ]. ( 49
)
The interaction energy of the gas is g 0 r b 3 ψ † ↑ ψ † ↓ ψ ↓ ψ ↑ , equal to (d E /dg) S g 0 (dg/dg 0 ) according to (37). The interaction energy of φ is simply g 0 |φ(0)| 2 , equal to g 0 (dg/dg 0 ) according to (40). The common factor g 0 (dg/dg 0 ) simplifies and

d E dg S = B. (50)
It remains to relate the coefficient B to lim k→+∞ k 4 n σ (k), using (46). In the zero-range limit, the scattering state φ(r) tends to 1 -a/r, so that, by Fourier transform, φ(k) tends to (2π) 3 δ(k) -4πa/k 2 . We then get (43).

Conclusion

We have calculated the number N b of closed-channel molecules in the whole BEC-BCS crossover for a two-component Fermi gas. Our result is in satisfactory agreement with the experimental results from Rice [5]. The expression (8) for N b is proportional to the quantity -dE/d(1/a) where E is the energy of the gas; this quantity is universal in the zero interaction range limit. At unitarity we find (25) and the length R * is related to the resonance parameters by [START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School Quantum Gases in Low Dimensions[END_REF]; in the zero range and zero effective range limit where k F b ≪ 1 and k F |r e | ≪ 1 we have k F R * ≪ 1 [see Eq.( 65)], we thus confirm that N b ≪ N i.e. the gas mainly populates the open channel. In the BCS limit, -dE/d(1/a) and thus N b are determined by the Hartree-Fock mean field energy, and do not tend exponentially to zero with 1/(k F |a|), contrarily to the prediction of usual BCS theory [5]. This is also related to the fact that the BCS theory does not predict the correct short-distance pair correlations [78].

N b /N = k F R * F (0) where F (0) is a universal constant given in
The quantity -dE/d(1/a) is related to the short-distance behavior of the pair correlation function and to the large-k tail of the momentum distribution, as discovered by Tan [START_REF] Tan | [END_REF]40] and rederived in [59] and in this paper. The quantity -dE/d(1/a) is also proportional to the average radiofrequency shift [79,80]. The general idea that various short-range two-body quantities such as the number of closed-channel molecules, the interaction energy of the gas and the average radiofrequency shift are the product of the value of the considered observable in the zero-energy twobody scattering state and a single universal many-body quantity was given by Leggett [START_REF] Leggett | Lecture 2 at the international quantum gases workshop[END_REF]. Our work shows that the Rice experiment [5] is the first direct measurement of this fundamental quantity in the BEC-BCS crossover.

It would be interesting to study the effect of a finite temperature. Experimentally, this could be achieved by measuring the number of closed-channel molecules together with the density profile of the trapped gas, the released energy [6,4], or the entropy [14,15]. Theoretically, Eq.( 12) can be used when more Monte Carlo results at finite temperature will become available. A useful tool in this context is the virial theorem: at unitarity it directly gives the total energy in terms of the density profile [11,15,81,82] while for a finite a it also involves [d E /d(1/a)] S [83,84].

We thank A. Bulgac, F. Chevy, M. Colomé-Tatché, R. Hulet, M. Jona-Lasinio, F. Laloë, P. Naidon, S. Nascimbène, G. Partridge, L. Pricoupenko and C. Salomon for very useful discussions and comments, and G. Partridge et al. for providing us with the data of [5] in numerical form. Our group is a member of IFRAF.

Note: We became aware of the related work [85] while completing this paper, and of [86] while revising it.

A Two-body analysis of the two-channel model in free space

A.1 Scattering amplitude

In this appendix we briefly discuss the two-body scattering properties for our model Hamiltonian H 2 in free space, that is in the absence of trapping potential (U b (r) ≡ 0, U (r) ≡ 0). The scattering problem may be considered in the center of mass frame so that we take the two-body scattering state

|ψ = βb † 0 |0 + d 3 k (2π) 3 A(k)a † k↑ a † k↓ |0 , (51)
(E b -E)β + Λγ = 0 (52) 2 k 2 m -E A(k) + χ(k)(βΛ + g 0 γ) = 0 ( 53
)
where we used the fact that the function χ(r) is real and even, and where we introduced the amplitude

γ = d 3 k (2π) 3 χ(k)A(k). ( 54
)
We now specialize to a state corresponding to the scattering of a spin ↑ atom with incoming wavevector k 0 onto a spin ↓ atom with incoming wavevector -k 0 , so that E = 2 k 2 0 /m is the energy of the incoming state. Eliminating β in terms of γ with (52) and using (53) leads to

A(k) = (2π) 3 δ(k -k 0 ) + γ χ(k)[g 0 + Λ 2 /(E -E b )] E + i0 + -2 k 2 /m (55)
where we have singled out the incoming state, and the scattered component is guaranteed to correspond to a pure outgoing wave in real space by the introduction of +i0 + in the denominator. Injecting ( 55) in (54) gives a closed linear equation for γ. From the known large r expression for the kinetic energy Green's function we extract the scattering amplitude

f k0 = -m χ(k 0 ) 2 /(4π 2 ) [g 0 + Λ 2 2 k 2 0 m -E b ] -1 -d 3 k (2π) 3 χ2 (k) 2 k 2 0 m +i0 + - 2 k 2 m . (56)
Since we are considering a s-wave Feshbach resonance, the function χ(r) is rotationally invariant, a property that we have used in (56). We now perform explicit calculations of various quantities characterizing the low energy scattering, such as the scattering length and the effective range. To this end, it is convenient as in [START_REF] Castin | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF] to perform the specific choice

χ(k) = exp(-k 2 b 2 /2) (57)
and to extend (56) to negative energies, setting k 0 = iq, q > 0 and E = -2 q 2 /m. One then obtains an explicit expression

- 1 f iq = e -q 2 b 2 4π 2 /m g 0 - Λ 2 E b + 2 q 2 /m + 1 b √ π + q erf (qb) -q (58
) where the last term -q is imposed by the optical theorem and erf is the error function.

The inverse of the scattering length is obtained by taking q → 0 in (58):

1 a = 4π 2 /m g 0 -Λ 2 /E b + 1 b √ π . ( 59
)
The large E b limit of this expression is the inverse of the background scattering length:

1 a bg = 4π 2 /m g 0 + 1 b √ π . ( 60
)
The location of the Feshbach resonance corresponds to an energy E 0 b of the closed-channel molecule such that 1/a = 0, so that

E 0 b = Λ 2 g 0 + 4π 3/2 2 b/m . ( 61
)
Supposing that the energy E b of the closed-channel molecule varies linearly with the magnetic field B,

E b (B) = E 0 b + µ b (B -B 0 ), (62)
and introducing the (possibly negative) resonance width ∆B such that

µ b ∆B = Λ 2 g 0 -E 0 b , (63)
we get from (59) exactly the B-field dependence (10) for the scattering length a in free space.

We now consider the effective range r e defined by the low-q expansion -1/f iq = a -1 -q + 1 2 q 2 r e + O(q 4 ). Its general expression for our model is

r e = -8π 4 Λ 2 /m 2 (g 0 E b -Λ 2 ) 2 + 4b √ π - 2b 2 a . ( 64
)
Of particular interest is the value of the effective range right on resonance,

r res e = -2R * + 4b √ π , (65)
where the length R * is always non-negative,

R * = Λ 2πbE 0 b 2 . ( 66
)
This definition of R * agrees with Eq. ( 29).

A very convenient rewriting of the scattering amplitude f iq is obtained by taking as independent parameters, in addition to the van der Waals length b, the two lengths a bg , R * rather than the bare parameters g 0 , Λ of the Hamiltonian. The quantity E b -Λ 2 /g 0 appearing in ( 58) is conveniently rewritten using the quantity

Q 2 ≡ m 2 g 0 g 0 E b -Λ 2 (67) = -1 a bg R * 1 1 -a bg /a (68)
which can be positive or negative. The ratio g 0 /Λ is also easily eliminated with the identity

g 2 0 4πΛ 2 = R * a 2 bg . (69)
Our final expression for the scattering amplitude is thus

- 1 f iq = e -q 2 b 2 a 1 -(1 -a/a bg ) q 2 q 2 + Q 2 + q erf (qb) -q (70) i.e. - 1 f k = e k 2 b 2 a 1 -(1 -a/a bg ) k 2 k 2 -Q 2 -ik erf (-ikb)+ik.
(71) We deduce the effective range in terms of a, a bg , R * and b:

r e = -2R * (1 -a bg /a) 2 + 4b √ π - 2b 2 a . (72)
A.2 Conditions for reaching the zero-range limit

In a degenerate gas, we expect that the zero-range limit is reached and the equation of state E(1/a) becomes model independent provided the conditions kb ≪ 1 ( 73)

- 1 f k ≃ 1 a + ik (74)
hold for k ∼ k F and for k ∼ 1/a on the BEC side [START_REF] Castin | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF]. We now assume that (73) holds, and we derive validity conditions for [START_REF] Chen | [END_REF] in the case of our two-channel model. Following the usual procedure, we use the second order expansion

- 1 f k = 1 a + ik - 1 2 r e k 2 + . . . , (75)
which reduces to [START_REF] Chen | [END_REF] if the effective range is small enough to have 1 2

|r e |k 2 ≪ 1 a + ik . ( 76
)
Taking into account the condition (73) as well as the expression (72) of r e for the two-channel model, and using |i + 1/(ka)| ≥ 1 and |1 + ika| ≥ 1, we find that the condition (76) becomes equivalent to

k 2 R * 1 - a bg a 2 1 a + ik -1 ≪ 1. (77)
A second condition is obtained by imposing that the higher order terms . . . in ( 75) are indeed negligible compared to k 2 |r e |/2. For this to hold we need not only ( 73), but also

k 2 ≪ |Q 2 | i.e. k 2 |a bg |R * 1 - a bg a ≪ 1. (78)
We conclude that to reach the universal zero-range limit, it is sufficient to fulfill the three conditions [START_REF] Tarruell | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF]77,78). All three conditions are satisfied for the experimental data on 6 Li from [5] reproduced in Fig. 1, where 1/k F > 250 nm and |a| > 200 nm: using the values |a bg | = 74 nm and R * = 0.027 nm (see section 2.3) and taking the estimate b ≃ 2.1 nm obtained by imposing that the effective range on resonance (65) should be equal to the value 4.7 nm obtained in the multi-channel calculation of [87], we find that the left hand sides of Eqs. [START_REF] Tarruell | Proceedings of the International School of Physics "Enrico Fermi" on Ultra-Cold Fermi gases[END_REF]77,78) are respectively smaller than 10 -2 , 3 • 10 -4 and 7 • 10 -5 .

To be complete, let us point out that the conditions (77,78) are not stricly necessary since it is possible that the . . . in (75) are not negligible compared to k 2 |r e |/2 and (74) nevertheless holds. We have chosen these conditions because they simplify the discussion, in particular they automatically ensure that the analytic continuation f iq has a pole at q ≃ 1/a so that the dimer energy is ≃ -2 /(ma 2 ). For this last property, in the regime b ≪ a, one directly gets from (70) the truly necessary condition

(1 -a/a bg ) -1 + (1 -a bg /a) -2 a/R * ≫ 1. ( 79
)
Close to the Feshbach resonance, a ≫ |a bg |, this reduces to the simpler condition

a ≫ R * . ( 80
)

B Effect of a difference between the trapping potentials seen by a closed-channel molecule and by an open-channel pair of atoms

In this appendix, we consider the physical consequences of the fact that, in general, the trapping potential U b (r) seen by a closed-channel molecule is not simply the trapping potential 2 U (r) seen by two separated atoms in the open channel:

U b (r) = 2 U (r), (81)
and we derive conditions for the gas to behave approximately as if we had U b (r) = 2 U (r). Note that a significant deviation in ( 81) is not a priori excluded for the laser traps used e.g. in [5,56] for lithium, with a laser wavelength λ L ≃ 1 µm: Contrarily to the atomic lightshift, the lightshift of the closed-channel molecule may be sensitive to the laser detuning from transitions to bound states of the 1 1 Σ + u (A) Born-Oppenheimer interaction potential between a ground-state atom and an excited atom, considering the fact that the minimum of this excited-ground interaction potential is separated from the dissociation limit of the ground-ground interaction potential by only ≃ 0.6 µm -1 < 1/λ L in spectroscopic units [88].

B.1 Position dependent scattering length and effective magnetic field

At the two-body level, (81) leads to an effective position dependent scattering length for two opposite spin atoms. To define such a local scattering length around point r 0 , we consider the "tangential" free space problem where the trapping potentials are replaced by the constant values U b (r 0 ) and U (r 0 ) and we call a(r 0 ) the resulting scattering length. We simply have to revisit the calculations of appendix A: In ( 52 m + 2U (r 0 ) to include the trapping energy shift of the two atoms; finally, for atoms with incoming wavevectors ±k 0 , the scattering state eigenenergy is now E = 2 k 2 0 m + 2U (r 0 ). The net effect on the scattering amplitude f k0 in ( 56) is to apply the substitution

E b -→ E b + U b (r 0 ) -2U (r 0 ), (82)
so that a(r 0 ) is simply obtained by applying ( 82) to (59).

A powerful interpretation of this substitution procedure may be given for a magnetic field independent µ b : We can rewrite the term

E b (B) + U b (r) appearing in (1) as E b [B eff (r)] + 2 U (r) provided we set B eff (r) = B + U b (r) -2 U (r) µ b . ( 83
)
Thus we can replace U b (r) by 2 U (r) in (1) provided we also replace the homogeneous magnetic field B in (1) by the inhomogeneous effective magnetic field B eff (r). One then immediately obtains the position dependence of the scattering length,

a(r) = a fs [B eff (r)], (84)
the function a fs (B) being given in (10). This effective magnetic field picture is also very efficient at the many-body level. In the grand canonical point of view, one adds the term -µ N to the Hamiltonian H 2 , where µ is the chemical potential, and the particle number operator N is defined as N = Nf + 2 Nb , Nf giving the number of fermions and Nb the number of bosons, since this is the quantity which commutes with H 2 . In the local density approximation, in the general case U b = 2U , we thus see that the gas behaves locally like a free space, homogeneous gas of chemical potential µ -U (r), in presence of the locally homogeneous magnetic field B eff (r).

B.2 Shift of the resonance location

We now consider the case of an optical trap such as the one at Rice [5]: U (r) and U b (r) are the lightshift potentials experienced by an atom and by a closed-channel molecule. Being both proportional to the laser intensity, they are mutually proportional:

U b (r) = 2(1 + η)U (r) ( 85
)
where the dimensionless parameter η is a convenient measure of the relative deviation between U b and 2U . At large r the laser intensity vanishes, so do U b and U . In practice, the gas is trapped around the intensity maximum of the laser field, located in r = 0, and the extension of the cloud along the axial (resp. radial) direction is small compared to the Rayleigh length (resp. the waist) of the laser beam; thus, to leading order, U (r) ≃ U (0), we can temporarily forget about the spatial variation of B eff and take B eff ≃ B -δB 0 , where

δB 0 = -U b (0) + 2 U (0) µ b = - 2ηU (0) µ b . (86)
Since the free space scattering length a fs (B) only depends on the difference B -B 0 , replacing B by B eff is equivalent to shifting the position B 0 of the Feshbach resonance by the quantity δB 0 . Introducing the length R * defined in [START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School Quantum Gases in Low Dimensions[END_REF], the ratio of this shift to the width ∆B of the Feshbach resonance may be written as

δB 0 ∆B = - 2mU (0) 2 a bg R * η. (87)
In the Rice experiment, the trap depth -U (0) is at most k B × 2.71 µK [89], so that we have

|δB 0 | |∆B| < |η| 7500 ( 88
)
where we have used the values of a bg and R * given below [START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School Quantum Gases in Low Dimensions[END_REF] and the mass of lithium. Except if η takes truly large values for the laser wavelength λ L = 1064 nm used in [5], we expect only a small shift of the resonance location due to the laser trap, relative to the resonance width.

In [5] the value of the scattering length as a function of B was calculated essentially from the free space formula a fs (B) ignoring the shift δB 0 due to the trap, and taking for the resonance location B 0 the value measured in [56] in a laser trap with different laser wavelength λ L = 1030 nm and intensity. On the contrary, our definition of a as a function of B includes the shift δB 0 , a = a fs (B -δB 0 ), as already defined in (9). This may introduce a systematic discrepancy between our theory curves and the experimental data in Fig. 1, in the form of shifts of the abscissas of the experimental points. We thus define the horizontal shift

δ 1 k F a ≡ 1 k F a fs (B -δB 0 ) - 1 k F a fs (B) (89)
= δB 0 ∆B 1 k F a bg (1 -a bg /a) 2 1 -(1 -a bg /a)δB 0 /∆B .( 90
)
To leading order in η (that is for |δB 0 /∆B| ≪ 1), and using the fact that the data from [5] reproduced in Fig. 1 is sufficiently close to resonance to have |1 -a bg /a| < 1.4, this simplifies to:

δ 1 k F a ≃ δB 0 ∆B 1 k F a bg (1 -a bg /a) 2 . ( 91
)
Using the parameters of the Rice experiment [89] and Eq.( 87) this gives

δ 1 k F a < |η| 1000 . (92)
Even without knowing η, we can thus resonably expect that this shift is negligible on the scale of Fig. 1, i.e. that

δ 1 k F a ≪ 1. (93)
This is also consistent with the fact that no systematic horizontal shift is apparent in Fig. 1 between the experimental data and the theory curves.

B.3 Effect of the position dependence of B eff (r)

We now discuss the effect on the many-body properties of the gas of the inhomogeneity of B eff (r), that is of a spatial dependence of the scattering length a(r), see (84), now taking into account the r dependence of U (r). We rewrite (83) as

B eff (r) = B+ 2ηU (r) µ b = B-δB 0 + 2η[U (r) -U (0)] µ b . (94)
When inserted in (84) this leads to the following position dependence for the scattering length:

1 a(r) -a bg = 1 a -a bg -η 2mR * 2 [U (r) -U (0)], (95)
where we recall that a is the scattering length in the trap center. Right on resonance, |a| = +∞, R * is one of the contributions to the interaction effective range, see (65), so that the zero-range limit implies k F |R * | ≪ 1. Since the typical value of U (r) -U (0) over the cloud size is of the order of 2 k 2 F /2m, we see that the small parameter k F R * appears in (95) so that we may expect that the position dependence of a(r) is negligible in the zero-range limit (for a fixed value of η).

Let us investigate this expectation more quantitatively and beyond the unitary limit, however still in a limiting case with simplifying assumptions. First, we assume that we are sufficiently close to resonance to have |a(r)| ≫ |a bg | over the atomic cloud size. This implies that the excursion of B eff (r) over the cloud size is small as compared to |∆B|; since the excursion of U (r) over the laser trap size is at most |U (0)|, we note that this condition is automatically satisfied for |δB 0 | ≪ |∆B|. Second, we shall assume that the gas may be treated in the zero temperature local density approximation. Third, we use the quadratic approximation (13) for U (r), further supposing for simplicity that the harmonic trap is isotropic, a case to which one can always reduce within the local density approximation by a rescaling of the coordinates, see appendix C. Equation (95) then reduces to

1 a(r) ≃ 1 a -η mR * 2 mω 2 r 2 , (96)
and the gas density profile n(r) is given by

µ -U (0) = 1 2 mω 2 r 2 + µ hom [n(r), 1 a(r) ] (97)
where µ hom [n, 1/a] is the chemical potential of a free space homogeneous gas of density n and scattering length a.

Let us derive simple conditions to have a weak effect of the scattering length inhomogeneity on the density profile. Some of these conditions will be expressed in terms of k F defined by (14).

In the BCS regime, a < 0 and k F |a| < 1, we have µ -U (0) ≃ k B T F = 2 k 2 F /2m. Over the Thomas-Fermi radius R of the cloud, such that 1 2 mω 2 R 2 ≃ k B T F , we see that the relative change of 1/a(r) is negligible if

|1 -a/a(R)| ≃ (k F |a|)(|η|k F R * ) ≪ 1. ( 98
)
This is satisfied in the zero-range limit, for fixed η, according to the condition (77) written for k = k F . In the unitary limit k F |a| ≫ 1, the previous reasoning fails, because requiring that the relative change of a(r) is small is too stringent. Taking for simplicity the limiting case |a| = +∞, it is sufficient to check that the gas remains locally unitary except in a small region near the edge of the cloud. First, to zeroth order in η, we calculate the local Fermi wavevector k hom F (r) using the unitary equation of state µ hom [n, 1/a = 0] = ξ 2 (k hom F ) 2 /2m, where k hom F is defined in (17), and using the value of the chemical potential µ -U (0) = √ ξk B T F . We obtain

k hom F (r) = k F ξ 1/4 1 - r 2 R 2 1/2 ( 99
)
where R is the unperturbed Thomas-Fermi radius such that mω 2 R 2 /2 = √ ξk B T F . Then, we include the η position dependent term in 1/a(r) to calculate at which distance l from the edge of the cloud the gas starts leaving the unitary regime, that is

k hom F (R -l)|a(R -l)| = 1.
(100)

Assuming a priori l ≪ R we obtain from (99)

l R ≃ ξ 3/2 2 (ηk F R * ) 2 . ( 101
)
We conclude a posteriori that l/R ≪ 1 and that the spatial dependence of a(r) is negligible when

|η|k F R * ≪ 1. ( 102
)
As in the BCS case, this is equivalent for a fixed |η| to the zero range condition (77) written here with k = k F and |a| = +∞. Finally we turn to the BEC regime a > 0 and k F |a| < 1. We then use the approximate equation of state

µ hom [n(r), 1/a(r)] ≃ - 2 2ma 2 (r) + π 2 2m αa(r)n(r) (103)
where α ≃ 0.60 is the proportionality factor in (20) between the dimer-dimer scattering length a d and a. In particular, this requires a ≫ R * , see (80). Replacing 1/a(r) by (96) in the above expression, we find that the leading term in 1/a 2 (r) gives rise to an external potential that sums up with the trapping potential to give

U eff (r) ≡ U (r) - 2 2ma 2 (r) (104) = U eff (0) + 1 2 ω 2 eff r 2 - m(ηR * ) 2 2 2 (mω 2 r 2 ) 2 . ( 105
)
We see that the atomic oscillation frequency is renormalized into

ω eff = 1 + 2ηR * a 1/2 ω. (106)
To have a small effect of the scattering length inhomogeneity on the trapped gas, we first require

ω eff ≃ ω so that |η|R * /a ≪ 1, (107)
a condition with no counterpart in the BCS regime and more stringent than (102) since k F a < 1 here. Again, we find that we recover, for a fixed |η|, the zero range condition (77) written here for k = 1/a. Second we require that the r 4 term in U eff (r) remains small as compared to the r 2 one over the Thomas-Fermi radius R of the condensate of dimers, calculated here to zeroth order in η; we reach the second condition

m(ηR * ) 2 2 mω 2 R 2 ≃ (ηk F R * ) 2 5 64 k F a d 2/5 ≪ 1. ( 108
)
Since k F a < 1 here, this second condition is automatically ensured by the first condition (107). Third, we require that the relative variation of 1/a(r) is small over the Thomas-Fermi radius R of the condensate, to ensure that the spatially inhomogeneity of a(r) in the mean field term [π 2 /(2m)]αa(r)n(r) has a small impact on the density profile:

|1-a/a(R)| ≃ (k F a)(|η|k F R * ) 5 64 k F a d 2/5 ≪ 1, (109)
which is also implied by (107).

To summarize, discussing perturbatively the effect of the position dependence of B eff (r) on the gas density profile in the local density approximation, we find close to the Feshbach resonance (|a| ≫ |a bg |, and also a ≫ R * in the BEC regime) that the effect is small under the condition (98) in the BCS regime, under the condition (102) in the unitary regime, and under the condition (107) in the BEC regime: In all three cases, this is the zero range condition (77), written respectively for k = k F , k = k F and k = 1/a, and multiplied by the parameter |η|. If |η| is not too large, neglecting the r dependence of the scattering length is thus automatically justified in the zero-range limit.

Interestingly, in the unitary limit, it is possible to relate ηk F R * to the quantity δ(1/k F a) introduced in the previous subsection:

|η|k F R * ξ 1/2 ≃ δ 1 k F a µ -U (0) U (0) . (110)
Experimentally, the gas is confined at the bottom of the laser trap so that [µ -U (0)]/|U (0)| ≪ 1 (its value near unitarity at Rice is 0.15 [89]). Together with the assumption (93) this implies that the right-hand-side of ( 110) is ≪ 1, and that the position dependence of a(r) has a negligible effect on the gas density.

A quantitative extension of the present qualitative discussion is to evaluate the density change δn(r) of the gas to first order in η. To this end, one simply expands (97) to first order in η around the solution n 0 (r) with η = 0, writing to first order n(r) = n 0 (r)+δn(r) and µ = µ 0 +δµ:

δn(r) = δµ -δ 1 a (r)∂ 1/a µ hom [n 0 (r), 1/a] ∂ n µ hom [n 0 (r), 1/a] (111)
where we have introduced the variation of 1/a(r) around η = 0 to first order in η, not supposing |a(r)| ≫ |a bg | in (95):

δ 1 a (r) = -1 - a bg a 2 η mR * 2 mω 2 r 2 . ( 112
)
δµ is then determined by the condition of a fixed number of particles r<R0 d 3 r δn(r) = 0, where R 0 is the Thomas-Fermi radius of the gas for η = 0. In this way, taking the variation of (130) to first order in η around η = 0, one can calculate the effect of a non-uniformity of a(r) on dE/d(1/a) to first order in η. We give here the result for the unitary gas:

δ [dE/d(1/a)] dE/d(1/a)| η=0 = ηk F R * 64 25π - 63π 256 ξ -1/4 ζ + 315π 2048 ξ 3/4 ζ -1 f (2) (0) + 2 3 ξ 1/2 k F a bg , (113)
where we recall that ξ = f (0), ζ = -f ′ (0) , see (24), and f (2) (0) is the second order derivative of f (x) with respect to x in x = 0. Note that, in this first order variation with respect to η, the trapping potential U (r) is kept fixed so that k F is also fixed. The calculation (113) shows that, at unitarity, a measurement of the variation of the number of closed-channel molecules N b as a function of U b (r)/U (r) may give access to f (2) (0). It also shows that, in the unitary limit, the condition to neglect a bg with respect to a(r) as done in the qualitative reasoning around (102), is k F |a bg | ≪ 1.

B.4 Metastability issues for a position dependent scattering length

We now discuss whether it is energetically possible for a pair of atoms to form a dimer which could escape from the trap by tunneling. The energy variation in such a process is

∆E = - 2 ma 2 (∞) -2µ ( 114
)
where µ is the chemical potential of the trapped gas and a(∞) is the effective scattering length at infinity, which we assume to be positive. We discuss under which conditions this process is energetically forbidden, i.e. ∆E > 0.

We assume that the gas has the usual Thomas-Fermi density profile and is weakly affected by the position dependence of B eff (r), see subection B.3. We distinguish between two regimes for the scattering length a in the trap center: the negative-a regime where 1/a ≤ 0, and the BEC regime where 0 < k F a < 1. Under the usual assumption that the gas is confined in a small region around the bottom of the laser trap we have µ ≃ U (0) in the negative-a regime and µ ≃ U (0) -2 /(2ma 2 ) in the BEC regime. We also make the simplifying assumption a(∞) ≫ |a bg | so that (95) reduces to

1 a(∞) ≃ 1 a + 2η m 2 R * U (0). (115)
In the BEC regime one easily deduces

∆E 2|U (0)| ≃ 1 + A -B (116)
where

A ≡ 2η R * a (117) B ≡ 2η 2 mR 2 * 2 |U (0)|; (118)
since we already assumed (107) we automatically have |A| ≪ 1. In the negative-a regime, Eq.( 115) immediately gives an upper bound on 1/a(∞), from which we get

∆E 2|U (0)| 1 -B. (119)
In both regimes we conclude that under the condition |B| ≪ 1 (120)

we have ∆E 2|U (0)| so that the process is energetically forbidden.

Finally we note that under the assumption |a| ≫ |a bg | we have from (90,87) Restricting to a spin balanced gas, we first consider the general case of the two-channel model (1), in which U b (r) = 2U (r). As explained in the appendix B, where an effective position dependent magnetic field B eff (r) is introduced, this leads to a position dependent scattering length a(r).

δ 1 k F a ≃ ηk F R * |U (0)| k B T F (121) so that B ≃ δ 1 k F a 2 k B T F |U (0)| . ( 122
As a consequence, in the local density approximation, here at zero temperature, the gas density n(r) is given by

µ = U (r) + µ hom [n(r), 1/a(r)] (123)
where µ hom [n, 1/a] is the chemical potential of the free space homogeneous gas of density n and scattering length a. The mean energy and particle numbers in the gas are given in the same approximation by

E = d 3 r {U (r) + ǫ hom [n(r), 1/a(r)]} n(r) (124) N = d 3 r n(r) (125)
where ǫ hom [n, 1/a] is the mean energy per particle of the free space gas. We wish to calculate the derivative of E with respect to 1/a, where a is the scattering length in the trap center r = 0, with the constraint that N should be fixed. Let us perform an infinitesimal variation δ 1 a of 1/a. This leads to a variation δµ of the chemical potential and δn(r) of the density, that may be calculated from (123) but this is not required here. This also leads to a variation of 1/a(r), according to (95): Using this expression and (123) one finds that a simplification occurs in (127). Furthermore, the fact that the particle number is fixed imposes [n, 1/a] for a free space gas of density n and scattering length a, a quantity that can be related to the derivative of ǫ hom [n, 1/a] with respect to 1/a simply by applying the general formula (8) to a free space system. Then setting for the trapped system in the local density approximation N b = d 3 r n hom b [n(r), 1/a(r)] and using (8) for the trapped system leads to (130).

In practice, in the present paper, we perform explicit calculations in the limiting cases presented in §2.2: One takes the zero-range limit, which also allows to neglect the spatial variation of a(r) under conditions discussed in the appendix B. Then ǫ hom is given by the relation (16) involving some universal function f (x). For the free space chemical potential we then introduce the convenient parametrization

µ hom [n, 1/a] = 2 2ma 2 u 1 k hom F a (131)
where k hom F = (3π 2 n) 1/3 is the ideal gas Fermi wavector, and the function u is

u(x) = f (x) -1 5 xf ′ (x) x 2 . ( 132
)
One also restricts to a harmonic trap for U (r), see (13).

The surface n(r) = 0 is then an ellipsoid with Thomas-Fermi radii R α along the directions α.

For a < 0 we set µ = 2 q 2 2m (133) with q > 0. Then ω α R α = q/m, and the density depends only on X ≥ 0 such that X 2 = α r 2 α /R 2 α . Equation (123) then reduces to the implicit equation

q 2 a 2 (1 -X 2 ) = u 1 k hom F (X)a ( 134
)
that one can solve numerically, using the interpolation for f (x) given in appendix D. For a > 0, µ hom [n, 1/a] tends to -2 /(2ma 2 ) for n → 0 so that we set

µ = 2 q 2 2m - 2 2ma 2 (135)
with q > 0. Then again ω α R α = q/m, and (123) reduces to

q 2 a 2 (1 -X 2 ) = 1 + u 1 k hom F (X)a ( 136
)
to be solved numerically in general.

One can then calculate dE/d(1/a) and k F defined in ( 14) from the expressions valid on both sides of the resonance: 3 (137)

dE d(1/a) = 3N 2 10m 1 0 dX X 2 k hom F (X) 4 f ′ 1 k hom F (X)a 1 0 dX X 2 k hom F (X)
k F = √ 2q π 1/3 4π 1 0 dX X 2 k hom F (X) 3 1/6 (138)
The analytical calculations performed in the weakly interacting regimes are also obtained from the above formulas.

D Interpolation formula for the homogeneous gas energy

We explain how an interpolation formula was constructed for the function f defined in (16), which was then used to produce the solid line in Fig. 1a and Fig. 1b.

On the BEC side of the resonance, we applied a cubic spline interpolation to the points in Table I of [32] with values of x = 1/(k hom F a) equal to 0, 1, 2, 4 and 6, using the last column of Table I (the one said to give E/N -ǫ b /2) and adding to the corresponding data the quantity -2 /(2ma 2 ). We added a home made point at x = 8 where the value of f (x) and its first order derivative f ′ (x) are obtained from the bosonic Lee-Huang-Yang expansion (19).

On the BCS side of the resonance we used the empirical interpolation formula of [90], f S (x) = α 1 -α 2 arctan[α 3 x(β 1 -x)/(β 2 -x)]. (139) Constraints on α 1 , α 2 , α 3 and β 1 , β 2 are obtained from the values f (0), f ′ (0) and from the large |x| expansion f (x) = 1 + c 1 /x + c 2 /x 2 + O(1/x 3 ), where the coefficients c 1 and c 2 are given in (22). This differs from [90] where the coefficient c 2 was not used and β 2 was deduced from a fit to the Monte Carlo data of [32] on the BCS side. Another difference with [90] is that we took the value f ′ (0) = -0.95, that we extracted from [34] using the relation (34). Also, we set the value of f (0) to 0.41, which is within the error bars of [32], in order to reduce the spurious discontinuity of the second order derivative of f (x) in x = 0 between the interpolation formula for x < 0 and the spline for x > 0. All this leads to α 1 ≃ 0.41, α 2 ≃ 0.3756, α 3 ≃ 1.062, β 1 ≃ 0.9043, β 2 ≃ 0.3797.

A test of this interpolation formulation is to compare its prediction for f ′ (-1) with the value ≃ -0.13 that one can extract from [34] using (34); one gets the satisfactory result f ′ S (-1) ≃ -0.14.

  r) +ψ † b (r) E b (B) -2 4m ∆ + U b (r) ψ b (r) +Λ d 3 r 1 d 3 r 2 χ(r 1 -r 2 ) ψ † b [(r 1 + r 2 )/2]ψ ↓ (r 1 )ψ ↑ (r 2 ) + h.c. + g 0 d 3 r 1 d 3 r 2 d 3 r 3 d 3 r 4 χ(r 1r 2 )χ(r 3r 4 )

Fig. 1 .

 1 Fig. 1. Derivative -dE/d(1/a) of the energy E of a trapped two-component Fermi gas, where a is the s-wave scattering length. Solid line: Theoretical prediction combining the Lee-Huang-Yang formulas(19,22) in the weakly interacting regimes with an interpolation of the fixed-node Monte Carlo results of[32,34] in the strongly interacting regime (see appendix D). Dashed lines: Analytical predictions in the weakly interacting regimes, on the BEC side(21) and on the BCS side(23). Cross: Result (25) at unitarity. Circles (with error bars): Experimental measurement of the number N b of closed-channel molecules in a lithium gas at Rice[5], combined with the present theory linking dE/d(1/a) to N b , see(28).(a): Linear scale on the vertical axis. (b): Logarithmic scale. -dE/d(1/a) is expressed in units of N 2 kF /m. All the theoretical predictions are obtained in the local-density approximation.

  Fig. 1. Derivative -dE/d(1/a) of the energy E of a trapped two-component Fermi gas, where a is the s-wave scattering length. Solid line: Theoretical prediction combining the Lee-Huang-Yang formulas(19,22) in the weakly interacting regimes with an interpolation of the fixed-node Monte Carlo results of[32,34] in the strongly interacting regime (see appendix D). Dashed lines: Analytical predictions in the weakly interacting regimes, on the BEC side(21) and on the BCS side(23). Cross: Result (25) at unitarity. Circles (with error bars): Experimental measurement of the number N b of closed-channel molecules in a lithium gas at Rice[5], combined with the present theory linking dE/d(1/a) to N b , see(28).(a): Linear scale on the vertical axis. (b): Logarithmic scale. -dE/d(1/a) is expressed in units of N 2 kF /m. All the theoretical predictions are obtained in the local-density approximation.

  with the unknown amplitude β in the closed channel and the unknown momentum dependent amplitude A(k) in the open channel. The annihilation operators for the bosonic molecules, b k , and for the fermionic atoms a kσ , σ =↑ or ↓, are the Fourier transforms of the corresponding field operators ψ b (r) and ψ σ (r), where the Fourier transform of a function f (r) is f (k) = d 3 r exp(-ik • r)f (r). Inserting the state (51) in Schrödinger's equation (H 2 -E)|ψ = 0 and projecting respectively on the molecular subspace and the atomic subspace leads to

  ) one has to replace E b by E b + U b (r 0 ) to include the trapping energy shift of the closed-channel molecule; similarly in (53) one has to replace

= 1 -

 1 a bg /a(r) 1 -a bg /a the gas energy is thenδE = d 3 r δn(r) U (r) + ǫ hom [n(r), 1/a(r)] r), 1/a(r)]. (127) Writing ǫ hom [n, 1/a] = E hom [N, V, 1/a]/N ,where V is the volume and N the particle number of the homogeneous gas, we find the relation n ∂ǫ hom ∂n [n, 1/a] = µ hom [n, 1/a] -ǫ hom [n, 1/a]. (128)

d 3 r 2 ×

 32 δn(r) = 0 (129) so that (127) leads todE d(1/a) N = d 3 r n(r) 1 -a bg /a(r) 1 -a bg /a ∂ǫ hom ∂(1/a) [n(r), 1/a(r)]. (130)We note that an alternative derivation of this result is possible: One can introduce the density of closed-channel molecules n hom b

)

  In a typical experiment such as the one at Rice we have k B T F < |U (0)|; thus the condition |B| ≪ 1 is satisfied as soon as |δ(1/k F a)| ≪ 1, which is likely to hold at Rice, see (92).

C Calculation of dE/d(1/a) in the local density approximation

b can be estimated by the van der Waals length, that is the length that one forms with , m and the C6 coefficient of the van der Waals interaction[21].

The assumption b ≪ |a| is not necessary on the BCS side if one restricts to macroscopic observables such as the total energy, but it is necessary for microscopic observables such as the ones considered in section

3.1.3 A simple formal definition of the zero-range limit for the homogeneous gas is to take the limit of vanishing density kF → 0 while tuning the magnetic field B in such a way that 1/(kF a) remains constant, and keeping constant all parameters of the Hamiltonian other than B ; one then also has to take the limit T → 0 so that T /TF remains constant.


Mickaël Péchaud Di

Renaud Keriven

Gabriel Peyré

Extraction of Tubular Structures over an Orientation Domain

Figure 1. Centerlines and radii extraction of vessels in a retinal image, vessels in a cortical image and roads in a satellite image.

Introduction

The problem of tubular structures extraction has received considerable attention in the computer vision and medical imaging communities -see for instance the recent review by Kirbas and Quek [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF]. They generally rely on the use of a local detector, postprocessed by a method that links locally detected tubular structure.

Various methods have been developed to locally detect tubular structures, for instance thresholding, ridge or crest detection [1,[START_REF] Prinet | Thin network extraction in 3d images: Application to medical angiograms[END_REF], wavelets [9], and matching filters [8,2].

Many processes allow to link or post-process the locally detected points. Besides classical fusion processes [3,25] and region growing algorithm [START_REF] Fang | Reconstruction of vascular structures in retinal images[END_REF]21], tracking-based approaches such as [18,[START_REF] Canny | A computational approach to edge detection[END_REF][START_REF] Schaap | Bayesian tracking of tubular structures and its application to carotid arteries in cta[END_REF] are an important attempt to to address this issue. They are based on a propagation inside the tubular structure from a known starting point by using a local search.

The notion of shortest path has proved to be efficient in order to perform global segmentation of tubular structures, in images, see for instance [START_REF] Cohen | Fast marching the global minimum of active contours[END_REF][START_REF] Poon | Live-vessel: Extending livewire for simultaneous extraction of optimal medial and boundary paths in vascular images[END_REF]. These geodesic curves can also be used to extract tubular structure centerlines, as proposed by Deschamps and Cohen [START_REF] Deschamps | Fast extraction of minimal paths in 3D images and applications to virtual endoscopy[END_REF] and by Santamaria et al. [19].

In order to enhance geodesic-based methods, Li and Yezzi [14] have proposed to extend the shortest path computation to a higher dimensional domain. They include the local radius of the vessels as an additional scale dimension in order to stabilize the computations and to select the centerline without any post-processing. This method is however facing difficulties on 2D images because of the design of their local vessel detector and of the junctions that are not handled in this 2D+1D space and scale domain.

This paper shows how to robustly detect tubular structures in 2D images by the use of shortest paths. It intro-1 duces an additional local orientation dimension that disambiguates the crossing singularities using pattern-recognition inspired techniques. Centerline detection and radius estimation are completely embedded in the proposed segmentation method and does not rely on additional properties of the images (such as centerlines being the darker part of vessels), nor on a post-processing step. Numerical results on synthetic and real-life examples shows the advantage of using this orientation information. The database will available online so that future method can compare their results to the proposed approach.

Shortest paths in orientation domain

This paper considers a gray scale image I : [0, 1] 2 → [0, 1] and in numerical computation, this image is sampled regularly on a grid of n × n pixels. In the following derivative operators in the continuous domain are used although the computations are done on the discrete grid using finite differences.

Local Tubular Structure Model

At the core of the proposed approach is the local detection of tubular structures. This paper uses a normalized cross-correlation with a tubular model as a local feature indicator.

The local geometry of a tubular structure is captured with a model

M (x) ∈ R for x = (x 1 , x 2 ) ∈ Λ = [-Λ 1 , Λ 1 ] × [-Λ 2 , Λ 2 ].
The models considered is this article are of the form M (x 1 , x 2 ) = m(x 2 ) and therefore only depend on a 1D profile m. A model is thus a small template which represents a typical horizontal tubular structure of normalized width one wishes to detect. This normalized pattern is then rotated and scaled in order to define warped models M r,θ (x) for x ∈ Λ(r, θ) = rR θ (Λ)

∀ x ∈ Λ(r, θ), M r,θ (x) def. = M (R -θ (x/r)) (1)
where R θ is the planar rotation of angle θ. Each pattern M r,θ looks like a typical tubular structure oriented along direction θ ∈ [0, π) and of width r > 0.

The model M can be tuned for each targeted application. Figure 2 shows some examples of these patterns for various r and θ.

The choice of Λ is a trade-off between robustness and accuracy of local structure detector. Experiments where carried out to determine an "optimal" Λ, which will not be detailled in this article due to the lack of space. Λ was chosen as a vertical rectangle of dimension 1 × 2 (cf. figure 2). The experiments show low sensitivity of the proposed method w.r.t. this choice.

Model for road extraction. A typical road in satellite imaging has a slow variation of intensity along a curve. It is thus efficiently captured by a binary model defined as

m(x 2 ) def. = 0 for |x 2 | > Λ 2 /2, 1 otherwise. ( 2
)
Model for blood vessels. Two important applications of vessels extraction are cortical optical imaging [12] and retinal imaging [START_REF] Niemeijer | Comparative study of retinal vessel segmentation methods on a new publicly available database[END_REF]. In order to capture efficiently blood vessels in these medical images, one designs a specific local pattern derived from a physical model of these structures.

The major difference with road tracking is that a vessel is not of constant intensity across its section.

The model is defined as

m(x2) def. =  1 for |x2| > Λ2/2, exp(-α p (1/2) 2 -(x2/Λ2) 2 ) otherwise. ( 3
)
This model assumes that the image has some basal intensity outside the vessels -e.g. light is reflected without absorption. For a point inside a vessel of circular section, the intensity is assumed to result from a light absorption (with coefficient α) proportional to the vessel width at this point.

Orientation Domain Lifting

In order to extract tubular structures that might exhibit self-crossings, the 2D image I is lifted to a 4D representation

F F : Ω def. = [0, 1] 2 × [r min , r max ] × [0, π) -→ | (4)
defined as the normalized cross-correlation between the image and the local model introduced in equation ( 1)

∀ (x, r, θ) ∈ Ω, F (x, r, θ) def. = N CC Λ(r,θ) (M r,θ (•), I(x+•))
(5) where I(x + •) is the image translated by x, N CC A (f, g) is the normalized cross-correlation between f and g over the domain A, defined by:

N CC A (f, g) def. = A (f -f )(g -ḡ) A (f -f ) 2 A (g -ḡ) 2 (6)
where h = R A h µ(A) , µ(A) being the area of A. Function F ranges from -1 to 1, is invariant under local intensity changes and is an indicator of how likely a tubular structure of radius r and orientation θ is present at location x.

The main interest of this lifting I → F , as already noticed in [START_REF] Jonasson | Representing Diffusion MRI in 5D Simplifies Regularization and Segmentation of White Matter Tracts[END_REF], is that it disambiguates situations where two tubular structures are crossing each other. Figure 3 shows an example of a medical image where orientation lifting is crucial to distinguish locally between two orientations. This lifting is computed for n r radii evenly spaced in [r min , r max ] and n θ orientations evenly spaced in [0, π), with n r = 12 and n θ = 12 in the experiments. This requires O((r max n) 2 n 2 n r n θ ) operations with r max ≪ 1 and n r , n θ ≪ n.

The 4D lifting (5) defines an isotropic metric ρ over the 4D domain Ω

∀ ω ∈ Ω, ρ(ω) def. = max(1 -F (ω), ε). (7)
The parameter ε prevents the metric to vanish and is set to ε = 10 -3 in the numerical tests.

Distance Map and Geodesic Computation

The length of a lifted curve γ : [0, 1] → Ω over the lifted domain is defined as

L F (γ) def. = 1 0 ρ(F (γ(t)))||γ ′ (t)||dt, (8)
The length of the speed vector

v = γ ′ (t) = (v x , v r , v θ ) is ||v|| 2 def. = v 2 x + λv 2 r + µv 2 θ . (9)
where (λ, µ) are normalizing constants set to λ = 0.5 and µ = 0.1 in numerical experiments. In practice one experiments strong robustness with respect to the choice of (λ, µ). Given a set A ⊂ Ω of seeds points and a set B ⊂ Ω of ending points, a shortest lifted curve γ * (t) ⊂ Ω joining A to B is defined as a shortest path for the metric L F γ * (A, B) def.

= argmin γ∈π(A,B)

L F (γ), [START_REF] Jonasson | Representing Diffusion MRI in 5D Simplifies Regularization and Segmentation of White Matter Tracts[END_REF] where π(A, B) is the set of curves γ such that γ(0) ∈ B and γ(1) ∈ A. The corresponding geodesic distance is

d F (A, B) = L F (γ * ).
In practice, γ * is estimated as follows [23]: the distance to the seeds A is the geodesic action map U A (ω) = d F (A, ω) is the unique viscosity solution of the Eikonal equation

||∇U A (ω)|| = ρ(ω), with ∀ ω ∈ A, U A (ω) = 0 ( 11
)
where

∇U A = ∂U A ∂x , λ ∂U A ∂θ , µ ∂U A ∂r T
, and where one has to be careful about computing the derivatives with respect to θ modulo π.

Calling ω 1 ∈ B the point in B with lowest distance to A, the geodesic curve γ * between A and B is then obtained by a gradient descent of

U A dγ * dt (t) = -∇U A (γ * (t)) with γ(0) = ω 1 (12)
The numerical computation of U A on a discrete grid of N

4D Curves Extraction

Shortest Paths and 4D curves

In order to extract a 4D curve c(x, x ′ ) between two points x, x ′ ∈ [0, 1] 2 , one needs to compute a shortest lifted tubular structure γ * between the lifted points in Ω that correspond to the x i . Let A(x) and A(x ′ ) be defined by

A(x ( ′ ) ) def. = (x ( ′ ) , r, θ) \ r ∈ [r min , r max ], θ ∈ [0; π) .
(13) This 4D curve is then defined

c x,x ′ def. = argmin γ∈π(ω,ω ′ ),ω ( ′ ) ∈A(x ( ′ ) ) L(γ). (14)
This 4D curve contains three components c x1,x2 (t) = (x(t), r(t), θ(t)). The path x(t) ⊂ [0, 1] 2 is the actual centerline over the image plane, whereas r(t) and θ(t) give the local width and orientation of the tubular structure (figure 4), which can be important in medical application, for instance to compute the local blood flow or to detect tissue diseases.

Numerical Experiments

Phantom experiments.

Experiments were carried out on several phantoms images, for which the centerlines positions and the radii have known analytical forms with sub-pixelic accuracy. The cross section of these phantoms corresponds to the model (3) with parameter α = 0.01. An additive Gaussian white noise with various amplitudes are added to the phantoms. Ten phantoms are generated for each condition, and each noise level, see Figure 5. In each case, the true starting and ending points of each phantom are used, as well as the true starting and ending radii for the [14] method. The resulting benchmark database will be available online in order to allow comparison with future works. For each phantom and level of noise, three algorithms were tested: the classical 2D method Fast Marching [START_REF] Cohen | Minimal Paths and Fast Marching Methods for Image Analysis[END_REF] with pre-smoothing of the image (pre-smoothing is classically used in Fast Marching based methods to obtain a centered path), Li/Yezzi method[14] and the method proposed in this article. The smoothing parameter of the classical 2D method, the smoothing parameter as well the three more parameters of the 3D method were optimized for this benchmark. In order not to take advantage of the exact knowledge of the profile, the proposed method was applied with a vol-untary wrong model parameter α = 0.1. In each case, the true starting and ending points of each phantom are used, as well as the true starting and ending radii for the [14] method.

For each method, the error of the vessels centerline position and radius w.r.t the ground truth was computed using the following formula.

Error C (c) 2 = 1 0 ||x(t) -x * (t * )|| 2 dt ( 15
)
where t * is such that x * (t) is the ground truth centerline point closest to the computed centerline x(t).

The error of the radii estimation ([14] and proposed method only) was computed using the following formula:

Error R (c) 2 = 1 0 |r(t) -r * (t * )| 2 dt ( 16
)
where r is the radius computed by the method, and r * the ground truth radius.

Figure 6 shows Error C (c) and Error R (c) curves for several synthetic images as a function of the noise level. Using the 3D space+scale lifting [14] produces results of varying quality, and requires a careful tuning for each noise level of the parameter to achieve the optimal error rate. [START_REF] Cohen | Minimal Paths and Fast Marching Methods for Image Analysis[END_REF] with an optimal smoothing provides a precise evaluation of the centerline locations, but without any evaluation of the local radius. The proposed method provides both positions and radii with more robustness and accuracy.

Real images experiments.

Figure 1 shows results of tubular structure extraction for the three modalities considered in this paper:

Left: vessels extraction on a retinal image from the DRIVE database [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF][START_REF] Niemeijer | Comparative study of retinal vessel segmentation methods on a new publicly available database[END_REF]. Center: vessels extraction for a complex optical imaging of the cortex with several branches and intersections.

Right: road extracted from a satellite image over an urban area. The starting area A is shown with a white square and several end points ω are shown with black squares.

The crossings in the retinal image (left) show the interest of the 4D lifting. Note also the large overlap of the centerlines computed from several different points, the subpixelic precision of the centerline and the correct handling of intersections on the cortical image (center). Another interest of this method is that it automatically computes the radius r(t) and orientation θ(t) parameters.

The precision of the 4D lifting method is evaluated on the DRIVE database [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF][START_REF] Niemeijer | Comparative study of retinal vessel segmentation methods on a new publicly available database[END_REF]. Approximate ground truth centerlines positions and radii are computed from the binary masks available with the database. Figure 7, shows the binary segmented vessels together with the ground trust centerline position and boundaries (top and middle), as well as the result of the 4D lifting method (bottom). The three geodesic extraction algorithms are applied to these three images between the indicated starting and ending points. Table 1 report the centerlines position and radii errors Error C (γ) and Error R (γ) for each method.

For the centerline extraction, due to the lack of precision of the ground truth, there is no significant difference between the proposed 4D lifting method and the space only geodesic extraction with smoothing of the metric. The 3D space+scale lifting [14] method showed unstable behavior with respect to its initialization and parameters, which had to be chosen carefully -for the second image, no parameters giving a correct result were found. The 4D lifting method is also more precise for the radii estimation than the 3D lifting.

Sensibility of the proposed method with respect

to the tubular structure model.

In order to show that the proposed method in not sensitive to the choice of the tubular structure model, some tests were performed with different values of α, showing that the proposed method is able to recover with good accuracy the centerlines and radii, even if the vessel model is unknown. Table 2 shows the results for one of these tests: for phantoms generated with α = 0.01, estimation with different fixed values of α are conducted for different noise levels. The obtained average errors remains comparable and in any case better than with [14] method.

Advantage of 4D modelling.

As already noted, figure 1 (left) shows a crossing situation that cannot be handled by methods that do not take orientation into account. In the case of vessel extraction, going straightforward is much more probable when encountering a crossing. Figure 8 shows a synthetic example with a such simple crossing configuration. Both classical 2D Fast Marching and Li/Yezzi methods give an undesired path, while the proposed method finds a correct result.

Figure 9 shows a real-life example with a simple road crossing. In this experiment, two starting points x 1 and x 2 were provided -e.g. the set of seed points is A(x 1 )∪A(x 2 ). Shortest paths from two different ending points are shown for the three tested algoritms. The proposed method again finds the correct paths.

Classical FM

Li/Yezzi Proposed method

Left blank during the reviewing process

Conclusion

This article introduces a novel framework for the automatic extraction of tubular structures, estimating both centerlines and radii. Examples on synthetic, medical and satellite images, with potentially intersecting curves, show the necessity to incorporate the orientation information in the model. The proposed method was compared to state of the art methods, proving to be more accurate and more stable. It shows robust results on different types of images. Its independence with respect to the local tubular structures detector makes it adaptable to any modality. A synthetic dataset will be available online in order to establish a benchmark for future contributions on centerlines/radii extraction.

Figure 2 .

 2 Figure 2. Vessel model (left) and road model (right) for different orientations and scales.

Figure 3 .

 3 Figure 3. left: original 2D image. right: 4D lifting (fixed radius), ranging from -1 (black) to 1 (white).

  def. = n 2 n r n θ pixels is done in O(N log(N )) operations with the Fast Marching algorithm of Sethian [22]. Furthermore, the Fast Marching is stopped as soon as the set B is reached, leading to important computation savings. During the Fast Marching computation, one should treat the orientation dimension θ ∈ [0, π) with periodic boundary conditions 0 ≃ π.

Figure 4 .

 4 Figure 4. Left: extraction of a vessel in a cortical image. Starting point: white square. Ending point: black square. Right: corresponding orientation θ(t) and radius r(t).

Figure 5 .

 5 Figure 5. Some of the sample phantoms used as benchmarks (basic intensities range from 0 to 1), shown here with a spatially independent Gaussian noise of variance 0.15.

Figure 6 .

 6 Figure 6. Centerlines and radii estimation for the first 2 phantoms presented fig. 5. Top row : error ErrorC (γ) (in pixel) between the computed and theoretical centerlines for the four methods versus noise level (100σ where σ is the independent Gaussian noise variance) (circles: proposed method, diamonds: Li/Yezzi[14], crosses: 2D Fast Marching with pre-smoothing[4] ) Bottom row : error ErrorR(γ) (in pixel) for the two first methods versus noise ratio.

Figure 7 .

 7 Figure 7. binary segmented images from the DRIVE database, together with the extracted ground trust. middle: corresponding images with ground-truth centerline and boundary and initial and ending points. bottom: centerlines positions and boundaries computed with the proposed method.

  C Err R Err C Err R Err C Err R 2D 033 1.67 3.13 3.31 0.53 1.90 4D 0.31 0.43 0.35 0.44 0.40 0.47

Figure 8 .

 8 Figure 8. Comparison of the different methods when encountering a crossing. Only the proposed method goes from the start to the end point following the natural path for vessel extraction.

Table 1 .

 1 Centerlines positions ErrorC and radii ErrorR estimation errors on retina images for the three different methods.

Table 2 .

 2 Centerlines positions ErrorC and radii ErrorR estimation errors for two different vessel models (gaussian noise with σ = 0.2)

		Phantom 1	Phantom 2
		Error C Error R Error C Error R
	α = 0.01	0.09	0.24	0.10	0.23
	α = 0.1	0.09	0.24	0.11	0.31
	α = 1	0.10	0.53	0.08	0.60

Acknowledgments