Søren Fournais
email: fournais@math.au.dk

Mathieu Lewin
email: mathieu.lewin@math.cnrs.fr

Jan Philip Solovej
email: solovej@math.ku.dk

THE SEMI-CLASSICAL LIMIT OF LARGE FERMIONIC SYSTEMS

We study a system of N fermions in the regime where the intensity of the interaction scales as 1/N and with an effective semiclassical parameter = N -1/d where d is the space dimension. For a large class of interaction potentials and of external electromagnetic fields, we prove the convergence to the Thomas-Fermi minimizers in the limit N → ∞. The limit is expressed using many-particle coherent states and Wigner functions. The method of proof is based on a fermionic de Finetti-Hewitt-Savage theorem in phase space and on a careful analysis of the possible lack of compactness at infinity.

Large interacting quantum systems are extremely difficult to describe accurately, due to the high complexity of the many-particle Schrödinger equation. It is therefore useful to derive approximate models that are much easier to handle and, at the same time, have a proper physical behavior. In the mean-field model, the real interaction between the particles is replaced by a Date: May 11, 2019. self-consistent effective potential depending on the average density of particles in the system. It is believed that this model becomes a good approximation to the true many-particle problem when the system is sufficiently dense (such that the particles meet very often) and the interactions are weak (for a law of large numbers to hold). The purpose of this paper is to discuss the validity of the mean-field model for fermions in a semi-classical-type limit.

The mean-field Hamiltonian for bosons and fermions. The Hamiltonian that will be the object of our attention describes N spinless quantum particles in R d , where d 1 is arbitrary (of course physically d = 1, 2, 3). In units such that m = 1/2, it takes the form

H ,N = N j=1 -i ∇ j + A(x j ) 2 + V (x j ) + 1 N 1 k<ℓ N w(x k -x ℓ ). (1.1)
For bosons this operator must be restricted to the subspace N s L 2 (R d ) ⊂ L 2 (R dN ) containing all the functions that are symmetric under exchange of variables, that is, Ψ(x σ(1) , ..., x σ(N ) ) = Ψ(x 1 , ..., x N ) for every permutation σ. For fermions, it must be restricted to the subspace N L 2 (R d ) ⊂ L 2 (R dN ) containing all the functions that are anti symmetric, Ψ(x σ(1) , ..., x σ(N ) ) = sgn(σ)Ψ(x 1 , ..., x N ).

In (1.1), A plays the role of a magnetic vector potential but can have a different physical origin (e.g. for rotating gases), V is an external potential that traps the particles (possibly only in a finite region of space if V → 0 at infinity), and w is the two-particle interaction potential (that could in principle depend on N as well, but that will be kept fixed here). The factor 1/N in front of the interaction term is typical of the mean-field scaling and its role is to make the two sums of the same order N in the Hamiltonian. In the limit of large N , the expectation is that the particles behave independently. In that case, the law of large numbers tells us that the jth particle experiences the mean-field potential 1 N k =j w(x j -x k ) ≃ ˆRd w(x j -y)ρ(y) dy where ρ is the density of particles in the system. The validity of the mean-field approximation can be studied in the timedependent or in the time-independent case. For the time-dependent equation, one usually assumes that the particles are independent at the initial time and then proves the propagation of chaos (that is, of the independence) for all times in the limit of large N . The situation is slightly different for the stationary problem where one has to prove that (say, in the ground state) the particles are independent. This second situation will be the object of our study.

For bosons, the validity of the mean-field approximation has been studied in the limit N → ∞ with fixed, first in some specific situations [START_REF] Lieb | Exact analysis of an interacting Bose gas. I. The general solution and the ground state[END_REF][START_REF] Benguria | Proof of the Stability of Highly Negative Ions in the Absence of the Pauli Principle[END_REF][START_REF]Gravitational collapse in quantum mechanics with relativistic kinetic energy[END_REF][START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF][START_REF] Solovej | Asymptotics for bosonic atoms[END_REF][START_REF] Bach | Ionization energies of bosonic Coulomb systems[END_REF][START_REF] Bach | On the number of bound states of a bosonic N -particle Coulomb system[END_REF][START_REF] Fannes | Equilibrium states for mean field models[END_REF][START_REF] Van Den Berg | The large deviation principle and some models of an interacting boson gas[END_REF][START_REF] Raggio | Quantum statistical mechanics of general mean field systems[END_REF][START_REF] Werner | Large deviations and mean-field quantum systems[END_REF][START_REF] Kiessling | The Hartree limit of Born's ensemble for the ground state of a bosonic atom or ion[END_REF][START_REF] Lieb | The mathematics of the Bose gas and its condensation[END_REF][START_REF] Lieb | Derivation of the Gross-Pitaevskii equation for rotating Bose gases[END_REF][START_REF] Seiringer | The excitation spectrum for weakly interacting bosons[END_REF][START_REF] Grech | The excitation spectrum for weakly interacting bosons in a trap[END_REF][START_REF] Seiringer | Disordered Bose-Einstein condensates with interaction in one dimension[END_REF]. The most general case has been addressed in a recent series of works [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF][START_REF]Remarks on the quantum de Finetti theorem for bosonic systems[END_REF][START_REF]The mean-field approximation and the non-linear Schrödinger functional for trapped Bose gases[END_REF][START_REF] Lewin | A note on 2D focusing many-boson systems[END_REF][START_REF] Rougerie | De Finetti theorems, mean-field limits and Bose-Einstein condensation[END_REF] by Nam, Rougerie and the second author of the present article. Under rather general assumptions on A, V and w, it is possible to prove that the bosonic ground state energy per particle converges to that of the nonlinear Hartree functional In mathematical terms, the bottom of the spectrum of H ,N restricted to the bosonic subspace N s L 2 (R d ), divided by N , converges to the infimum of E V,A Hartree :

lim

N →∞ inf σ N s L 2 (R d ) (H ,N ) N = inf ´Rd |u| 2 =1 E V,A Hartree (u).
For the convergence of states and the link with the minimizers of E V,A Hartree , we refer to [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF] and to the discussion below. The link between the timedependent Schrödinger equation i Ψ = H ,N Ψ and the nonlinear time-dependent Hartree equation

i ∂ ∂t u = -i ∇ + A 2 u + V u + (|u| 2 * w)u
has stimulated many works as well [START_REF] Hepp | The classical limit for quantum mechanical correlation functions[END_REF][START_REF] Ginibre | The classical field limit of scattering theory for nonrelativistic many-boson systems. I[END_REF][START_REF] Spohn | Kinetic equations from Hamiltonian dynamics: Markovian limits[END_REF][START_REF] Bardos | Weak coupling limit of the Nparticle Schrödinger equation[END_REF][START_REF] Elgart | Gross-Pitaevskii equation as the mean field limit of weakly coupled bosons[END_REF][START_REF] Elgart | Mean field dynamics of boson stars[END_REF][START_REF] Ammari | Mean field limit for bosons and infinite dimensional phasespace analysis[END_REF][START_REF] Erdös | Rigorous derivation of the Gross-Pitaevskii equation with a large interaction potential[END_REF][START_REF] Fröhlich | On the mean-field limit of bosons with Coulomb two-body interaction[END_REF][START_REF] Rodnianski | Quantum fluctuations and rate of convergence towards mean field dynamics[END_REF][START_REF] Knowles | Mean-field dynamics: singular potentials and rate of convergence[END_REF][START_REF] Pickl | A simple derivation of mean-field limits for quantum systems[END_REF][START_REF] Lewin | Fluctuations around Hartree states in the mean-field regime[END_REF]. The coupled limit N → ∞ and → 0 has been investigated for the time-dependent problem in [START_REF] Lions | Sur les mesures de Wigner[END_REF][START_REF] Graffi | Mean-field approximation of quantum systems and classical limit[END_REF][START_REF] Fröhlich | Mean-field and classical limit of manybody schrödinger dynamics for bosons[END_REF]. The purpose of this work is to address the case of fermions. The antisymmetry (called the Pauli principle) implies that two particles cannot be at the same place and this usually makes the kinetic energy grow much faster than N , by the Lieb-Thirring inequality [START_REF] Lieb | Bound on kinetic energy of fermions which proves stability of matter[END_REF][START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Stability of Matter in Quantum Mechanics[END_REF]. More precisely, if Ψ N is an N -particle normalized antisymmetric function with support in a bounded domain Ω N ⊂ R dN , then we have

N j=1 ˆΩN |∇ j Ψ| 2 C|Ω| -2 d N 1+ 2 d .
In order to ensure that all the terms in the Hamiltonian are of the same order, it is therefore necessary for fermions to take = 1

N 1 d
, which means that the mean-field limit is coupled to a semi-classical limit.

We therefore end up with the Hamiltonian

H N := N j=1 -i∇ j N 1 d + A(x j ) 2 + V (x j ) + 1 N 1 k<ℓ N w(x k -x ℓ ), (1.3)
that will be our main object of interest in this work. Its fermionic ground state energy (that is, the bottom of the spectrum in the fermionic subspace N L 2 (R d )) will be denoted by

E(N ) = inf σ N L 2 (R d ) (H N ). (1.4)
A given physical Hamiltonian is not necessarily in the form (1.1) but it can sometimes be recast in that form after an appropriate scaling. This is for instance the case for atoms [START_REF]The Thomas-Fermi theory of atoms, molecules and solids[END_REF][START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF] (in units of length of the order Z -1/3 ) and non-relativistic fermion stars [START_REF]Gravitational collapse in quantum mechanics with relativistic kinetic energy[END_REF][START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF], see (1.11) and (1.12) below.

Vlasov and Thomas-Fermi theories. The Hamiltonian H N has been studied a lot in the time-dependent setting [START_REF] Narnhofer | Vlasov hydrodynamics of a quantum mechanical model[END_REF][START_REF] Spohn | On the Vlasov hierarchy[END_REF][START_REF] Bardos | Mean field dynamics of fermions and the time-dependent Hartree-Fock equation[END_REF][START_REF] Elgart | Nonlinear Hartree equation as the mean field limit of weakly coupled fermions[END_REF][START_REF] Fröhlich | A microscopic derivation of the time-dependent Hartree-Fock equation with Coulomb two-body interaction[END_REF][START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF][START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF][START_REF] Bach | Kinetic Energy Estimates for the Accuracy of the Time-Dependent Hartree-Fock Approximation with Coulomb Interaction[END_REF][START_REF] Petrat | A New Method and a New Scaling For Deriving Fermionic Mean-field Dynamics[END_REF][START_REF] Benedikter | From the Hartree dynamics to the Vlasov equation[END_REF], where its dynamics is known to converge to the Vlasov time-dependent equation in the limit N → ∞. This is the Hamiltonian dynamics associated with the Vlasov energy, which is the semi-classical equivalent to the Hartree functional (1.2):

E V,A Vla (m) = 1 (2π) d ˆRd ˆRd |p + A(x)| 2 m(x, p) dx dp + ˆRd V (x)ρ m (x) dx + 1 2 ¨Rd ×R d w(x -y)ρ m (x) ρ m (y) dx dy. (1.5)
Here

ρ m (x) = 1 (2π) d ˆRd m(x, p) dp.
and m(x, p) is a probability measure on the phase space R d × R d that must satisfy the additional constraint 0 m(x, p) 1 a.e. (1.6) This condition says that one cannot put more than one particle at x with a momentum p and it is inherited from the Pauli principle.

We will look at the ground state problem, that is, we study the minimization of the Vlasov energy (1.5). If we minimize in the variable p for any fixed x without the constraint (1.6), the ground state measures all have the trivial momentum distribution δ 0 (p + A(x)) and the first term disappears. With the fermionic constraint (1.6), the Vlasov energy is minimized for measures of the form

m ρ (x, p) = 1 |p + A(x)| 2 c TF ρ(x) 2/d (1.7)
where ρ now minimizes the Thomas-Fermi energy

E V TF (ρ) := E V,A Vla (m ρ ) = d d + 2 c TF ˆRd ρ(x) 1+ 2 d dx + ˆRd V (x)ρ(x) dx + 1 2 ¨Rd ×R d w(x -y)ρ(x) ρ(y) dx dy (1.8)
and

c TF = 4π 2 d |S d-1 | 2 d .
The functional is similar to the Hartree energy (1.2) with ρ = |u| 2 except that the quantum term ´Rd |∇u| 2 has been replaced by the fermionic classical kinetic energy proportional to ´Rd ρ 1+2/d and that the magnetic field has been discarded (but it appears in the formula (1.7) of the minimizers on phase space).

We are interested in the link between the many-particle ground state energy E(N ) and that of the Vlasov and Thomas-Fermi energies in the limit N → ∞ with = N -1/d . To this end we introduce the Thomas-Fermi ground state energy

e V TF (λ) := inf E V TF (ρ) : 0 ρ ∈ L 1 (R d ) ∩ L 1+2/d (R d ), ˆRd ρ = λ = inf E V,A Vlas (m) : 0 m 1, (2π) -d ˆR2d m = λ ,
where λ = 1 in our case. The equality of the two infima is obtained by inserting the formula (1.7) in the Vlasov energy (1.5).

Convergence of E(N )/N . Our first result is that the leading order of E(N ) is given by the Thomas-Fermi ground state energy e V TF (1). Theorem 1.1 (Convergence of the ground state energy). Assume that w is even, that w,

|A| 2 ∈ L 1+d/2 (R d ) + L ∞ ε (R d ) and that either V ∈ L 1+d/2 (R d ) + L ∞ ε (R d ) or V -∈ L 1+d/2 (R d ) + L ∞ ε (R d ), V + ∈ L 1 loc (R d ), lim |x|→∞ V + (x) = +∞.
Then we have

lim N →∞ E(N ) N = e V TF (1). (1.9) We recall that f ∈ L p (R d ) + L ∞ ε (R d ) means that for every ε > 0, we can write f = f p + f ∞ with f p ∈ L p (R d ) and ||f ∞ || ∞ ε.
The assumptions on the potentials cover most physically interesting systems in this regime. Our result can easily be generalized in many directions but we do not state precisely the corresponding theorems. Our two assumptions on V cover either systems which are locally confined (the assumption

V ∈ L 1+d/2 (R d ) + L ∞ ε (R d
) essentially means that V → 0 at infinity) or confined (when V + → +∞ at infinity). Our proof works the same if V has a different limit at infinity or if V = +∞ outside of a domain Ω, which corresponds to Dirichlet boundary conditions. We could also weaken the assumptions on w + but we refrain from doing it in order to simplify the next statement. It is also possible to use a pseudo-relativistic kinetic energy

-i∇ j N 1 d + A(x j ) 2 + m 2 (1.10)
instead of the non-relativistic kinetic energy, but this generalization (and the appropriate modifications on E TF and on the assumptions on the potential) will not be discussed further. Finally, we notice that the assumption on the magnetic potential A is probably far from optimal, but it already covers the physical case of a potential in 3D satisfying ∇ • A = 0 and whose magnetic field B = ∇ × A is square-integrable (since then A ∈ L 6 (R 3 ) by the Sobolev inequality).

Results of the form of (1.9) have been proved for particular models. For atoms as in the work of Lieb and Simon [START_REF] Lieb | Thomas-Fermi theory revisited[END_REF][START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF], we have

d = 3, A(x) = 0, V (x) = - 1 t|x| , w(x) = 1 |x| , (1.11)
where t = lim(N/Z) is the limiting proportion of electrons and protons.

In [START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF], Lieb and Yau studied pseudo-relativistic stars but their results also apply to the simpler non-relativistic model which corresponds to

d = 3, A(x) = 0, V (x) = 0, w(x) = - 1 |x| .
(1.12)

In Section 3, we provide an elementary proof of Theorem 1.1 which is very much inspired of [START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF] and was recently written for bosons in [START_REF]Mean-field limit of Bose systems: rigorous results[END_REF].

Coherent states and Wigner functions. The main results of the paper concern the convergence of states, which requires more subtle tools. We will express it using coherent states and Wigner functions, but other choices are possible.

Let f be any fixed normalized real-valued function in L 2 (R d ). For every fixed (x, p) in the phase space R d × R d , we introduce the coherent state

f x,p (y) = -d 4 f y -x √ e i p•y , (1.13)
where we recall that = N -1/d . Then we have the resolution of the identity

1 (2π ) d ˆRd ˆRd |f x,p f x,p | dx dp = 1 (1.14) in L 2 (R d ).
A typical choice is f a Gaussian, but any real-valued f can indeed be used. For any such f and a fermionic N -particle state Ψ N , we introduce the corresponding k-particle Husimi function [START_REF] Husimi | Some Formal Properties of the Density Matrix[END_REF][START_REF] Takahashi | Wigner and husimi functions in quantum mechanics[END_REF][START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF] m

(k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) := Ψ N , a * (f x 1 ,p 1 ) • • • a * (f x k ,p k )a(f x k ,p k ) • • • a(f x 1 ,p 1 )Ψ N , (1.15)
for k = 1, ..., N , where a and a * are the fermionic annihilation and creation operators. In practice f is a very well localized function and the measure m

(k)
f,Ψ N describes how many particles are in the k semi-classical boxes with length scale √ , centered at (x 1 , p 1 ),...,(x k , p k ) in the phase space

R d × R d . Two equivalent formulas are m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) = N ! (N -k)! Ψ N , P x 1 ,p 1 ⊗ • • • ⊗ P x k ,p k ⊗ 1 N -k Ψ N L 2 (R dN ) (1.16) = N ! (N -k)! ˆRd(N-k) f x 1 ,p 1 ⊗ • • • ⊗ f x k ,p k , Ψ N (•, y) L 2 (R dk ) 2
dy (1.17)

where P x,p := |f x,p f x,p | is the orthogonal projection onto f x,p . As will be proved below in Lemma 2.2, we have 0 m

(k) f,Ψ N 1 and 1 (2π) dk ˆR2dk m (k) f,Ψ N = N (N -1) • • • (N -k + 1) N k -→ N →∞ 1,
for all k 1. Next we turn to the k-particle Wigner function which is defined as in [START_REF] Lions | Sur les mesures de Wigner[END_REF][START_REF] Fröhlich | Mean-field and classical limit of manybody schrödinger dynamics for bosons[END_REF] by

W (k) Ψ N (x 1 , p 1 , ..., x k , p k ) := ˆRdk ˆRd(N-k) e -i k ℓ=1 p ℓ •y ℓ × × Ψ N (x 1 + y 1 /2, • • • , x k + y k /2, x k+1 , ..., x N )× × Ψ N (x 1 -y 1 /2, • • • , x k -y k /2, x k+1 , ..., x N ) dy 1 • • • dy k dx k+1 • • • dx N . (1.18) Contrary to m (k) f,Ψ N , the Wigner function W (k)
Ψ N is not necessarily positive, but it will have the same positive limit as m

(k) f,Ψ N in the semi-classical regime.
Convergence of states: confined case. In our main results about the convergence of states, we for simplicity distinguish the confined and unconfined situations, which we state in two separate theorems. We start with the former.

Theorem 1.2 (Convergence of states, confined case). Assume that w is even, that w,

|A| 2 , V -∈ L 1+d/2 (R d ) + L ∞ ε (R d ) and that V + ∈ L 1 loc (R d ), lim |x|→∞ V + (x) = +∞.
Let {Ψ N } ⊂ N L 2 (R d ) be any sequence such that Ψ N = 1 and

Ψ N , H N Ψ N = E(N ) + o(N ). (1.19)
Then there exists a subsequence {N j } and a probability measure P on the set

M = 0 ρ ∈ L 1 (R d ) ∩ L 1+2/d (R d ) : ˆRd ρ = 1, E V TF (ρ) = e V TF (1)
of all the minimizers of the Thomas-Fermi functional, such that the following limit holds:

m (k) f,Ψ N j (x 1 , p 1 , ..., x k , p k ) → ˆM k ℓ=1 1 |p ℓ + A(x ℓ )| 2 c TF ρ(x ℓ ) 2/d =mρ(x ℓ ,p ℓ ) dP(ρ)
(1.20) weakly in L 1 (R 2dk ) and weakly- * in L ∞ (R 2dk ), for all k 1 and every real-

valued normalized f ∈ L 2 (R d ) or, in other words, ˆR2dk m (k) f,Ψ N j ϕ → ˆM ˆR2dk (m ρ ) ⊗k ϕ dP(ρ) (1.21) for every test function ϕ ∈ L 1 (R 2dk ) + L ∞ (R 2dk ). For the Wigner function W (k) Ψ N j
defined in (1.18), we have the same limit

ˆR2dk W (k) Ψ N j ϕ → ˆM ˆR2dk (m ρ ) ⊗k ϕ dP(ρ), (1.22)
this time for every test function ϕ such that

∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k ϕ ∈ L ∞ (R 2dk ), max(α j , β j ) 1. (1.23)
Furthermore, we have the convergence of the k-particle probability density

ˆRd • • • ˆRd |Ψ N j (x 1 , ..., x N j )| 2 dx k+1 • • • dx N j → ˆM k j=1 ρ(x j ) dP(ρ) (1.24) weakly in L 1 (R d )∩L 1+ 2 d (R d ) for k = 1
, and weakly- * in the sense of measures for k 2. Finally, we have the convergence of the k-particle kinetic energy density

ˆRd • • • ˆRd F [Ψ N j ](p 1 , ..., p N j ) 2 dp k+1 • • • dp N j → ˆM k ℓ=1 ρ |p ℓ + A| d c -d/2 TF dP(ρ), (1.25)
weakly- * in the sense of measures for k 1.

In the statement,

F [f ](p) := 1 (2π ) d/2 ˆRd f (x)e -i p•x dx (1.26)
is the -dependent Fourier transform. Later we also use the notation f := F 1 [f ] for the unscaled Fourier transform. The condition max(α j , β j ) 1 in (1.23) means that all the components of the d-dimensional multi-indices α j and β j are 1. The condition on ϕ is taken from [START_REF] Hwang | The L 2 -Boundedness of Pseudodifferential Operators[END_REF] but it can be replaced by any other for which the Calderon-Vaillancourt theorem holds true.

The result says that, in the limit N → ∞, the many-body approximate minimizers Ψ N become purely semi-classical to leading order and that the corresponding semi-classical measures are a convex combination of factorized states involving the Vlasov minimizers m ρ with ρ ∈ M. Note that if the Thomas-Fermi energy has a unique minimizer ρ 0 , then there is no need to extract subsequences and the probability measure P has to be a delta measure at ρ 0 . In particular, all the limits are now factorized, which corresponds to independent probabilities. The convergence of the k-particle densities (1.24) and kinetic energy densities (1.25) follows easily from the convergence (1.20) of the Husimi measures and the fact that the system is confined (and from the Lieb-Thirring inequality [START_REF] Lieb | Bound on kinetic energy of fermions which proves stability of matter[END_REF][START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Stability of Matter in Quantum Mechanics[END_REF] for the one-particle density).

The simplest fermionic trial states Ψ N are the Slater determinants Ψ N = (N !) -1/2 det(ψ i (x j )) (a.k.a. Hartree-Fock states) where ψ 1 , ..., ψ N form an orthonormal system [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF][START_REF]Solutions of Hartree-Fock equations for Coulomb systems[END_REF][START_REF] Solovej | The ionization conjecture in Hartree-Fock theory[END_REF]. The one-particle Husimi measure of a Slater determinant is given by

m (1) f,Ψ N (x, p) = N i=1 ψ i , f x,p 2 .
It can be proved that

Ψ N satisfies (1.19) if m (1)
f,Ψ N (x, p) is converging to a Vlasov minimizer m ρ with ρ ∈ M (see Lemma 3.2 below). For a Slater determinant, the k-particle semi-classical measures are always factorized in the semi-classical limit and we then end up with a Dirac delta P = δ ρ . Our proof of the convergence of the energy (Theorem 1.1) will indeed rely on the Hartree-Fock problem.

An equivalent way to formulate the limits in Theorem 1.2 is in terms of quantization of observables. In other words, for every test function ϕ we associate an operator Op N, (ϕ) acting on N L 2 (R d ) and we look at the limit of Ψ N , Op N, (ϕ)Ψ N . The quantization associated with the coherent states f x,p is defined by

Op N, f (ϕ) := ˆ(R d ×R d ) k ϕ(x 1 , ..., p k )a * (f x 1 ,p 1 ) • • • a * (f x k ,p k )× × a(f x k ,p k ) • • • a(f x 1 ,p 1 ) dx 1 • • • dp k .
Similarly, associated with the Wigner function there is an operator defined in terms of its integral kernel by Op N, Weyl (ϕ)(x 1 , ..., x N , y 1 , ..., y N )

:= N k -1 1 i 1 <•••<i k N Op k, Weyl (ϕ)(x i 1 , ..., x i k , y i 1 , ..., y i k )
where

Op k, Weyl (ϕ)(x 1 , ..., x k , y 1 , ..., y k ) := -dk ˆRdk ϕ x 1 + y 1 2 , p 1 , ..., x k + y k 2 , p k e i k j=1 p j •(x j -y j ) dp 1 • • • dp k .
Then the theorem gives the limit

Ψ N j , Op N j ,N -1/d j f /Weyl (ϕ)Ψ N j -→ j→∞ ˆM ˆR2dk m ⊗k ρ ϕ dP(ρ)
for any fixed ϕ in the function spaces mentioned in the statement.

Convergence of states: unconfined case. In the unconfined case we have a similar result, except that the limits are a priori local. Since some of the particles can escape to infinity, our result will involve the minimizers of the problems e V TF (λ) for a mass 0 λ 1. Theorem 1.3 (Convergence of states, unconfined case). Assume that w is even and that w,

|A| 2 , V ∈ L 1+d/2 (R d ) + L ∞ ε (R d ). Let {Ψ N } ⊂ N L 2 (R d ) be any sequence such that Ψ N = 1 and Ψ N , H N Ψ N = E(N ) + o(N ).
Then there exists a subsequence {N j } and a probability measure P on the set

M = 0 ρ ∈ L 1 (R d ) ∩ L 1+2/d (R d ) : ˆRd ρ 1, E V TF (ρ) = e V TF ˆRd ρ = e V TF (1) -e 0 TF 1 - ˆRd ρ
which coincides with the set of all the possible weak limits of minimizing sequences for the Thomas-Fermi problem, such that

• the limit (1.21) for m (k) f,Ψ N j holds for every ϕ ∈ L 1 (R 2dk )+L ∞ ε (R 2dk ); • the limit (1.22) for W (k) Ψ N j
holds for every ϕ satisfying (1.23) and tending to zero at infinity;

• the limit (1.24) for ρ (k) Ψ N j holds weakly in L 1 loc (R d ) ∩ L 1+ 2 d (R d ) for k = 1, and weakly- * on C 0 0 (R dk ) instead of C 0 (R dk ) for k 2. If ´Rd ρ = 1 for P-almost-every ρ ∈ M, then these limits hold as in Theo- rem 1.2.
In the unconfined case some particles may be lost at infinity (if not all), and the limiting minimizing densities ρ might not be probability measures. Nevertheless, the result says that the remaining particles must solve the minimization problem e V TF ( ´ρ), corresponding to the fraction ´Rd ρ of the N particles which have not escaped to infinity. Furthermore, if no particle is lost ( ´Rd ρ = 1 on M), then the convergence is the same as in the confined case.

Note that the weak limit of the k-particle kinetic energy density (1.25) is unknown in the unconfined case, due to the lack of control in the x-variable. However, it follows from (1.21) that the amount of kinetic energy in the limit cannot exceed the Vlasov one

lim inf j→∞ Ψ N j , (-∆)Ψ N j N 1+ 2 d j ˆB d d + 2 c TF ˆRd ρ 1+ 2 d dP(ρ).
Discussion and strategy of proof. Our theorems 1.2 and 1.3 seem to be the first dealing with the convergence to the Thomas-Fermi problem for general potentials A, V, w. Given that the Thomas-Fermi model is at the core of many approximate models in chemistry and material science [START_REF] Elliott | Semiclassical origins of density functionals[END_REF], this result is important for applications. For atoms and stars, our result provides the following information.

Example 1.4 (Large atoms). For the choice (1.11), the Thomas-Fermi energy is

ˆR3 3 5 c TF ρ(x) 5/3 - ρ(x) t|x| dx + 1 2 ˆR3 ˆR3 ρ(x) ρ(y) |x -y| dx dy
and the set M consists of a unique ρ t which has the mass ´R3 ρ t = min(1, 1/t). Theorem 1.3 therefore provides the strong convergence to m ρt for 0 < t 1 and the weak convergence to m ρ 1 for t > 1. This was proved in [START_REF] Lieb | Thomas-Fermi theory revisited[END_REF][START_REF]The Thomas-Fermi theory of atoms, molecules and solids[END_REF][START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF].

Example 1.5 (Non-relativistic stars). For stars (1.12), the corresponding Thomas-Fermi energy is

3 5 c TF ˆR3 ρ(x) 5/3 dx - 1 2 ˆR3 ˆR3 ρ(x) ρ(y) |x -y| dx dy
and we have , 27, 61, 62, 60]. Therefore, if we can find a translation of the system for which the limiting semi-classical measures are non trivial, the limit can only involve the Thomas-Fermi minimizing density ρ 0 . The convergence of states seems to be new. It was only discussed for a locally perturbed model in [START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF]Thm. 6]. As mentioned previously, we have a similar convergence result for pseudo-relativistic stars.

M = {0} ∪ ρ 0 (• -τ ), τ ∈ R 3 for some unique ρ 0 with ´R3 ρ 0 = 1 [3, 2
Our main new tool in this paper is a fermionic (weak) version of the de Finetti-Hewitt-Savage theorem for classical measures (Theorem 2.7 below), which implies that the weak limits of the semi-classical measures must always be a combination of factorized probabilities (Theorem 2.1 below). We follow here ideas of [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF], where a weak version of the quantum de Finetti theorem was introduced for the mean-field limit of bosons. The fact that we do not need a quantum version for fermions is of course due to the semi-classical feature of our model. The Pauli principle only persists in the constraint that 0 m 1. We refer to [START_REF] Rougerie | De Finetti theorems, mean-field limits and Bose-Einstein condensation[END_REF] for a general presentation of classical and quantum de Finetti theorems, with applications to the mean-field limit for bosons, and to [START_REF] Golse | On the Dynamics of Large Particle Systems in the Mean Field Limit[END_REF][START_REF] Ammari | Mean field limit for bosons and infinite dimensional phasespace analysis[END_REF] for the time-dependent problem.

Organization of the paper. The rest of the paper is devoted to the proofs of our results. In Section 2 we give the main properties of the Husimi and Wigner measures. The main result there is Theorem 2.1 in which we explain that a sequence {Ψ N } has always convergent semi-classical measures (up to extraction of a subsequence), with a limit that is a convex combination of factorized 'fermionic" measures on the phase space:

m (k) f,Ψ N j ⇀ ˆ0 µ 1 ´R2d µ (2π) d µ ⊗k dP(µ), ∀k 1.
Sections 3 and 4 are then devoted to the proof of Theorems

Fermionic semi-classical measures

This section is devoted to the study of the general properties of the Husimi measures m (k) f,Ψ N which we have defined in (2.9), and to the relation with the Wigner measures. The main result of the section is a general theorem about the properties of the measures obtained in the limit N → ∞, whose proof will be given later in Section 2.5.

Theorem 2.1 (Convergence to factorized fermionic measures on phase space). Let Ψ N be a sequence of normalized fermionic functions, = N -1/d and m

(k) f,Ψ N , W (k)
Ψ N be defined as in (1.15) and (1.18). Then, there exists a subsequence N j and a probability measure P on the set

B = µ ∈ L 1 (R 2d ) : 0 µ 1, (2π) -d ˆR2d µ 1 such that, for every k 1, ˆR2dk m (k) f,Ψ N j ϕ → ˆB ˆR2dk µ ⊗k ϕ dP(µ), (2.1
)

for every normalized real-valued function f ∈ L 2 (R d ) and every ϕ ∈ L 1 (R 2dk )+ L ∞ ε (R 2dk ), and ˆR2dk W (k) Ψ N j ϕ → ˆM ˆR2dk µ ⊗k ϕ dP(µ) (2.2)
for every function ϕ tending to zero at infinity, satisfying (1.23). Furthermore, if Ψ N satifies the kinetic energy bound

Ψ N ,   N j=1 -∆ j   Ψ N CN 1+2/d , (2.3
)

then ˆB ˆR2d |p| 2 µ(x, p) dx dp dP(µ) C(2π) d
(hence ´R2d |p| 2 µ(x, p) dx dp is finite P-almost surely) and the k-particle densities converge weakly ˆRdN

U (x 1 , ..., x k )|Ψ n (x 1 , ..., x N )| 2 dx 1 • • • dx N → ˆB ˆRdk ρ ⊗k µ U dP(µ) (2.4) for every k 1 and every U ∈ L 1+d/2 (R d ) + L ∞ ε (R d ) for k = 1 and U ∈ C 0 0 (R dk ) for k 2, where ρ µ (x) = 1 (2π) d ˆRd µ(x, p) dp.
The result says that, whatever the sequence {Ψ N }, the semi-classical measures always get factorized in the limit N → ∞ or, more precisely, they are a convex combination of factorized phase-space measures uniformly bounded by 1 and with a mass 1. In the rest of this section, we first derive some elementary properties of the measures m (k) f,Ψ N , then we state a general convergence theorem in the spirit of the classical de Finetti-Hewitt-Savage theorem, before we provide the detailed proof of Theorem 2.1 in Section 2.5.

2.1.

Measures on phase space for N -body states. In this section, we recall the elementary properties of the measures m (k) f,Ψ N . Some of the results of this section are well-known, but we gather them all here for the convenience of the reader.

Coherent states. As before, we fix a real-valued normalized function f ∈ L 2 (R d ) and define

f x,p (y) = -d 4 f y -x √ e i p•y , (2.5)
as well as

f := f 0,0 = -d/4 f (•/ √ ) and g := F [f ] = -d/4 f (•/ √ ).
(2.6)

We will later have to take = N -1/d but for the moment could be arbitrary. We recall the resolution of the identity

1 (2π ) d ˆRd ˆRd |f x,p f x,p | dx dp = 1. (2.7)
The latter follows from the useful formulas

f x,p , u = ˆRd f (y -x)u(y)e -i p•y dy = (2π ) d/2 F f (• -x)u (p) (2.8) = ˆRd g (k -p)F [u](k)e -i k•x dk = (2π ) d/2 F g (• -p)F [u] (x)
which imply immediately that

1 (2π ) d ˆRd ˆRd | f x,p , u | 2 dx dp = ˆRd ˆRd F f (• -x) u (p) 2 dp dx = ˆRd ˆRd f (y -x)u(y) 2 dy dx = f 2 L 2 (R d ) ||u|| 2 L 2 (R d ) = ||u|| 2 L 2 (R d ) .
Semi-classical measures on phase space. Next, we recall the definition (1.15) of the k-particle semi-classical (Husimi) measure

m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) := Ψ N , a * (f x 1 ,p 1 ) • • • a * (f x k ,p k )a(f x k ,p k ) • • • a(f x 1 ,p 1 )Ψ N , (2.9)
for k = 1, ..., N . Here a * (f ) and a(f ) are, respectively, the fermionic creation and annihilation operators, satisfying the Canonical Anticommutation Relations [START_REF] Dereziński | Mathematics of quantization and quantum fields[END_REF] a

* (f )a(g) + a(g)a * (f ) = g, f , a * (f )a * (g) + a * (g)a * (f ) = 0. (2.10)
See (1.16) and (1.17) for two other formulas for m

(k) f,Ψ N .
The following is a simple consequence of these formulas and of the fact that Ψ N is an antisymmetric wave function.

Lemma 2.2 (Elementary properties of the phase space measures).

Let f ∈ L 2 (R d , R) and Ψ N ∈ N 1 L 2 (R d ) be two normalized functions. Then for every 1 k N , the function m (k)
f,Ψ N defined in (2.9) is symmetric and satisfies

0 m (k) f,Ψ N 1 a.e. on R 2dk , (2.11) 1 (2π) dk ˆR2dk m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) dx 1 • • • dp k = N (N -1) • • • (N -k + 1) dk , (2.12)
and

1 (2π) d ˆR2d m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) dx k dp k = d (N -k + 1)m (k-1) f,Ψ N (x 1 , p 1 , ..., x k-1 , p k-1 ). (2.13)
Proof. The symmetry of m

(k)
f,Ψ N follows from the fact that two creation (resp. annihilation) operators anti-commute due to (2.10) or, equivalently, from the formula (1.16) and the fact that Ψ N is anti-symmetric. Similarly, the uniform bound (2.11) is a consequence of the operator inequality 0 a * (f )a(f ) 1. The formulas (2.13) and (2.12) follow from the representation (1.16) and the resolution of the identity (2.7), which can be rewritten as (2π ) -d ˜R2d P x,p dx dp = 1.

From (2.12) and (2.13), the choice = N -1/d arises naturally, as it makes (2π) -dk m (k) f,Ψ N a probability measure in the limit N → ∞ for all k 1.

Remark 2.3. The measures m (k) f,Ψ N can be defined in a similar manner for bosons. Only the uniform bound (2.11) (which is the expression of Pauli's principle here) will be lost. Link with k-particle densities. For any fixed (normalized) fermionic function Ψ N , we denote by

ρ (k) Ψ N (x 1 , ..., x k ) := N k ˆRd • • • ˆRd |Ψ N (x 1 , ..., x N )| 2 dx k+1 • • • dx N (2.14) and t (k) Ψ N (p 1 , ..., p k ) := N k ˆRd • • • ˆRd |F [Ψ N ](p 1 , ..., p N )| 2 dp k+1 • • • dp N
(2.15) the position and momentum densities for k particles, with 1 k N . For later purposes, we also define the k-particle density matrix of Ψ N which is the operator γ

(k) Ψ N with integral kernel γ (k) Ψ N (x 1 , ..., x k ; x ′ 1 , ..., x ′ k ) := N k ˆRd • • • ˆRd Ψ N (x 1 , ..., x N )× × Ψ N (x ′ 1 , ..., x ′ k , x k+1 , ..., x N ) dx k+1 • • • dx N . (2.16
) All the physical observables can be expressed in terms of these objects only.

The following gives a link with the semi-classical measures we have defined in the previous section.

Lemma 2.4 (Particle densities and fermionic semi-classical measures). Let f be any normalized function in L 2 (R d , R) and f , g be defined as in (2.6). Then we have

1 (2π) dk ˆ(R d ) k m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) dp 1 • • • dp k = k! dk ρ (k) Ψ N * |f | 2 ⊗k (2.17) and 1 (2π) dk ˆ(R d ) k m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) dx 1 • • • dx k = k! dk t (k) Ψ N * |g | 2 ⊗k (2.18) for any (normalized) fermionic function Ψ N .
Proof. As in (2.8), we start by noticing that, for every fixed y

∈ (R d ) N -k , f x 1 ,p 1 ⊗ • • • ⊗ f x k ,p k , Ψ N (•, y) L 2 ((R d ) k ) = (2π ) dk F f x 1 ,0 ⊗ • • • ⊗ f x k ,0 Ψ N (•, y) (p 1 , ..., p k ) = (2π ) dk F g 0,p 1 ⊗ • • • ⊗ g 0,p k F [Ψ N ](•, y) (x 1 , ..., x k ).
Next we compute the integral of m (k) f,Ψ over the p j 's, using (1.16):

ˆ(R d ) k m (k) f,Ψ N (x 1 , p 1 , ..., x k , p k ) dp 1 • • • dp k = k! N k ˆ(R d ) k dp 1 • • • dp k ˆRd(N-k) dy f x 1 ,p 1 ⊗ • • • ⊗ f x k ,p k , Ψ N (•, y) 2 = (2π ) dk k! N k ˆ(R d ) k dp 1 • • • dp k ˆRd(N-k) dy × F f x 1 ,0 ⊗ • • • ⊗ f x k ,0 Ψ N (•, y) (p 1 , ..., p k ) 2 = (2π ) dk k! N k ˆRdN f (y 1 -x 1 ) • • • f (y k -x k )Ψ N (y) 2 dy 1 • • • dy N = (2π ) dk k! ρ (k) Ψ N * |f | 2 ⊗k (x 1 , ..., x k
). The proof of (2.18) is similar.

A simple consequence of (2.18) is the well-known formula for the kinetic energy and its generalization with magnetic field (see, e.g., [START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF]Eq. (5.20)]):

Corollary 2.5 (Kinetic energy). If f ∈ H 1 (R 3 , R) and |A| 2 ∈ L 1+d/2 (R d ) + L ∞ ε (R d ), we have Ψ N , N j=1
(-i ∇ j + A(x j )) 2 Ψ N = 1 (2π ) d ˆRd ˆRd |p + A(x)| 2 m (1) f,Ψ N (x, p) dx dp -N ˆRd |∇f | 2 + 2ℜ Ψ N , N j=1 (A -A * |f | 2 )(x j ) • (-i ∇ j ) Ψ N + Ψ N , N j=1 (|A| 2 -|A| 2 * |f | 2 )(x j ) Ψ N . (2.19)
Proof. By density we may assume that Ψ and A are regular enough, which allows us to make the following calculation. By definition, we have

Ψ N ,   N j=1 -2 ∆ j   Ψ N = ˆRd t (1)
Ψ N (p)|p| 2 dp = 1 (2π ) d ˆRd ˆRd m (1)
f,Ψ N (x, p)|p| 2 dx dp + ˆRd ˆRd t

(1)

Ψ N (p)|g (q -p)| 2 (|p| 2 -|q| 2 ) dp dq = 1 (2π ) d ˆRd ˆRd m (1) f,Ψ N (x, p)|p| 2 dx dp - ˆRd ˆRd t (1)
Ψ N (p)|g (q -p)| 2 |q -p| 2 dp dq -2 ˆRd pt (1)
Ψ N (p) dp • ˆRd p|g (p)| 2 dp .
The last term vanishes since f is real, hence g is even, and we obtain

Ψ N , N j=1 -2 ∆ j Ψ N = 1 (2π ) d ˆRd ˆRd |p| 2 m (1) f,Ψ N (x, p) dx dp -N ˆRd |∇f | 2 . (2.20)
With magnetic field, we expand the square

| -i ∇ + A(x)| 2 = -2 ∆ + A(x) • (-i ∇) + (-i ∇) • A(x) + |A(x)| 2
and apply (2.17) for the term with |A| 2 and (2.20) for the term with -∆.

The cross term is calculated using (2.8) as follows:

1 (2π ) d ¨Rd ×R d p • A(x)| u, f xp | 2 dx dp = ¨Rd ×R d p • A(x)|F (f x,0 u)(p)| 2 dx dp = ˆRd A(x) • ˆRd p|F (f x,0 u)(p)| 2 dp dx = ˆRd A(x) • ℑ ˆRd f x,0 (y)u(y) ∇(f x,0 u)(y) dy dx = ˆRd A(x) • ˆRd |f (y -x)| 2 ℑ(u(y) ∇u(y)) dy dx = ˆRd (A * |f | 2 )(y) • ℑ(u(y) ∇u(y)) dy.
2.2. Structure of the limiting measures: tight case. For an arbitrary sequence Ψ N with N → ∞, the functions (m

(k) f,Ψ N ) N k are bounded in L 1 (R 2dk ) ∩ L ∞ (R 2dk
), for every fixed k, by Lemma 2.2. It is therefore clear that we can find a subsequence such that m

(k) f,Ψ N ⇀ m (k) f weakly for every k 1, in the sense that ˆR2dk m (k) f,Ψ N ϕ → ˆR2dk m (k) f ϕ for every ϕ ∈ L 1 (R 2dk ) + L ∞ ε (R 2dk
). In the limit we obtain a hierarchy of symmetric functions (m

(k) f ) k 1 .
The purpose of this section and of the following is to study the properties of these limiting hierarchies. In general, the limiting functions m (k) f are not probability densities because some mass can be lost at infinity. However, if the sequence (m

(1) f,Ψ N ) is tight, that is, lim R→∞ lim sup N →∞ ˆ|x|+|p| R m (1) f,Ψ N (x, p) dx dp = 0, then the m (k)
f,Ψ N are also tight for k 2 and the limiting m

(k)
f are all probability measures. We consider this simpler case in this section. The tightness of the sequence also implies that the compatibility relation (2.13) is preserved in the limit:

1 (2π) d ˆR2d m (k) (x 1 , p 1 , ..., x k , p k ) dx k dp k = m (k-1) (x 1 , p 1 , ..., x k-1 , p k-1 )
(2.21) for all k 1. The famous de Finetti-Hewitt-Savage theorem deals with the structure of such infinite sequences of symmetric probability measures [START_REF] De Finetti | Funzione caratteristica di un fenomeno aleatorio[END_REF][START_REF]La prévision : ses lois logiques, ses sources subjectives[END_REF][START_REF] Dynkin | Classes of equivalent random quantities[END_REF][START_REF] Hewitt | Symmetric measures on Cartesian products[END_REF][START_REF] Diaconis | Finite exchangeable sequences[END_REF]. In our situation, the result can be stated as follows.

Theorem 2.6 (Fermionic semi-classical measures on phase space). Let m (k) be a family of symmetric positive densities in

L 1 (M k ), with M ⊂ R D , sat- isfying c ˆM m (k) (ξ 1 , ..., ξ k ) dξ k = m (k-1) (ξ 1 , ..., ξ k-1 ) (2.22)
with m (0) = 1 and 0 m (k) 1 for all k 1. Then there exists a Borel probability measure P on the set

S := µ ∈ L 1 (M ) : 0 µ 1, c ˆM µ = 1 such that m (k) = ˆS µ ⊗k dP(µ), (2.23)
for all k 1.

In our case M = R 2d and c = (2π) -d but the result is actually true for any Borel set M in R D and any c > 0. The theorem says that infinite exchangeable (i.e. symmetric) fermionic systems are always convex combination of independent ones.

Proof.

A very similar theorem is proved in [13, p. 510-511] (see also [START_REF] Messer | Statistical mechanics of the isothermal Lane-Emden equation[END_REF]Lemma 4]). The usual theorem from [START_REF] Hewitt | Symmetric measures on Cartesian products[END_REF] furnishes a probability measure P ∈ P(P(M)) such that (2.23) holds with S replaced by P(M ). The condition 0 µ 1 can be characterized by the property that µ(A) |A| for all A in a countable set of balls, showing that S is a measurable set in P(M ). We therefore only have to prove that this measure P is concentrated on S. Now, the assumption that 0 m (k) 1 implies m (k) (A k ) |A| k for any Borel set A ⊂ M , and this gives

ˆP(M) µ(A) |A| k dP(µ) 1
for every k 1. Taking k → ∞ proves that P is concentrated on the subset of P(M ) containing all the probability measures µ such that µ(A) |A| for all A. These measures are absolutely continuous with respect to the Lebesgue measure and the corresponding density is between 0 and 1.

2.3.

Structure of the limiting measures: general case. As we have said, in general the limiting functions m (k) f need not be probability measures, and they need not satisfy the compatibility condition (2.21). The idea that a de Finetti theorem nevertheless holds true in this case as well seems to have been advertized for the first time in the quantum case in [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF]. A similar result with a different interpretation had been published before by Ammari and Nier in [START_REF] Ammari | Mean field limit for bosons and infinite dimensional phasespace analysis[END_REF]. The following is inspired of [44, Thm 2.2]. N be its marginals defined recursively as in (2.22). We assume that m (0) = 1, that 0 m

m (k) = ˆB µ ⊗k dP(µ), (2.24)
for all k 1.

The theorem says that the natural object obtained in the semi-classical limit of large fermionic systems (with possible lack of compactness) is a Borel probability measure P on the set of all the functions 0 µ 1 on M such that ´M µ 1/c.

Proof. The proof follows step by step that of [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF]Thm 2.2] in the quantum case and we only outline it. We fix a ball B R = {ξ ∈ M : ||ξ|| R} in the space M , with R an arbitrary positive integer. We then look at the probability measures for the particles to be in B R . In particular, we will need to know the exact number of particles in B R , which leads us to introduce the functions

g N,R,0 := c N ˆ(M\B R ) N m (N ) N (ξ 1 , ..., ξ N )dξ 1 • • • dξ N and g N,R,n (ξ 1 , ..., ξ n ) := c N -n N n n j=1 1 B R (ξ j ) ˆ(M\B R ) N-n m (N ) N (ξ 1 , ..., ξ N )dξ n+1 • • • dξ N
for n = 1, ..., N , and where ξ ∈ M is the phase space variable. In words, g N,R,n is the probability density for n particles in B R , under the constraint that the other N -n are all outside of B R . Note that

N n=0 c n ˆ(B R ) n g N,R,n = c N ˆMN m (N ) N = 1 (2.25)
and that the k-particle density in B R can be written as

k j=1 1 B R (ξ j )m (k) N (ξ 1 , ..., ξ k ) = N n=k N -k n-k N n g (k) N,R,n (ξ 1 , ..., ξ k ) = N n=k n k N k g (k) N,R,n (ξ 1 , ..., ξ k )
where g (k)

N,R,n denotes the kth marginal of the function g N,R,n . It is useful to think of the fraction n/N of particles in B R as an additional variable 0 t 1. For fixed N this new parameter takes values in Z/N ∩ [0, 1]. In the limit N → ∞, the whole interval [0, 1] is filled and t will correspond to the radial integration in the ball B. We therefore introduce the probability measure

dg (k) N,R (t, ξ 1 , ..., ξ k ) := N n=k g (k) N,R,n (ξ 1 , ..., ξ k )δ n/N (t) + |B R | -k k-1 n=0 c k ˆ(B R ) k g N,R,n δ 0 (t) (2.26) on the compact set [0, 1] × (B R ) k .
Extracting subsequences, we may assume that g

(k)
N,R converges weakly to a probability measure g

(k)
R on [0, 1] × (B R ) k , for every k 1. Since R is an integer by assumption, the weak convergence can be assumed for every R as well. Now, from the inequality

0 n N k - n k N k (k -1) 2 N -k + 1
(see [44, Eq.(2.13)]) and the normalization (2.25), we deduce that

c k ˆ(B R ) k k j=1 1 B R (ξ j )m (k) N (ξ 1 , ..., ξ k ) - ˆ1 0 t k g (k) N,R (t, ξ 1 , ..., ξ k ) dt dξ 1 • • • dξ k (k -1) 2 N -k + 1 ,
from which we conclude that

1 ⊗k B R m (k) = ˆ1 0 t k dg (k) R (t, •) (2.27)
for every k, R 1.

Next we remark that

t c ˆBR dg (k+1) N,R (•, ξ k+1 ) -dg (k) N,R = tδ k/N (t) g (k) N,R = k N δ k/N (t) g (k) N,R
where the right side converges to 0 in the sense of measures. Passing to the weak limit N → ∞, we deduce that the limiting hierarchy is consistent for all 0 < t 1:

ct ˆBR dg (k+1) R (t, •, ξ k+1 ) = t dg (k) R (t, •).
At t = 0 this equation tells us nothing. However, we can always change the value of the limit dg

(k)
R as we want since the point t = 0 does not contribute in (2.27). Similarly to (2.26), we choose it to be the constant that makes it a probability measure and the consistency at t = 0 is then obvious. From the de Finetti-Hewitt-Savage theorem (or, rather, a simple one-parameter version of it which can be proved similarly as sketched in [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF]), we conclude that there exists a probability measure P R on [0, 1] × P(B R ) such that

dg (k) R (t, •) = ˆP(B R ) p ⊗k d P R (t, p).
Inserting in (2.27), we get the representation formula

(1 B R ) ⊗k m (k) = ˆ[0,1]×P(B R ) (tp) ⊗k d P R (t, p).
Using that 0 m (k) 1 we can show similarly as in the proof of Theorem 2.6, that P R concentrates on [0, 1] × S. The right side can then be uniquely interpreted as an integral for a measure P R on B R := {0 µ 1 B R : ´BR µ 1/c}, with t playing the role of the radial variable:

(1 B R ) ⊗k m (k) = ˆBR µ ⊗k dP R (µ).
Since the previous equality holds for every positive integer R, the measure P R is the cylindrical projection of P R ′ for every R ′ > R and the family (P R ) defines uniquely a measure P over B. We get the result by passing to the limit R → ∞.

2.4.

Link with Wigner functions. It is well known that the Wigner functions converge to the same limit as the Husimi measures and we quickly recall this here. It is indeed a general principle that all regular quantizations are equivalent in the semi-classical limit. First, we recall that when f = (π) -d/4 e -|x| 2 /2 , the Husimi and Wigner measures are related by a convolution:

m (k) Ψ,f = N (N -1) • • • (N -k + 1) N k W (k) Ψ * G , where G (x 1 , ..., p k ) = (π ) -dk exp(--1 k j=1 |x j | 2 + |p j | 2
) is a Gaussian in phase space (see, e.g. [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF]Prop. 21]). It is then clear that the two sequences must have the same weak limits. On the other hand, the Calderon-Vaillancourt theorem from [START_REF] Hwang | The L 2 -Boundedness of Pseudodifferential Operators[END_REF] states that

Op k, Weyl (ϕ) C sup max(α,β) 1 ∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k ϕ L ∞ (R 2dk )
, for any function ϕ for which the right side is finite. This bound is useful for showing that W

Ψ has weak limits. From this we obtain the estimate

Op N, Weyl (ϕ) -Op N, f (ϕ) C sup max(α,β) 1 ∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k (ϕ -ϕ * G ) L ∞ (R 2dk ) + C k 2 N sup max(α,β) 1 ∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k ϕ L ∞ (R 2dk )
.

The term ϕ -ϕ * G can be estimated uniformly by the gradient of ϕ, using the fundamental theorem of calculus, leading to

Op N, Weyl (ϕ) -Op N, f (ϕ) C k √ sup max(α,β) 2 ∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k ϕ L ∞ (R 2dk )
.

One more derivative of ϕ is necessary to prove that the limit of the Wigner function coincides with that of the Husimi measure, but the limit of W Ψ nevertheless holds with only max(α, β) 1, using the density of smooth test functions.

As a conclusion, these estimates give the convergence of the Wigner functions when tested against any ϕ satisfying

∂ α 1 x 1 • • • ∂ α k x k ∂ β 1 p 1 • • • ∂ β k p k ϕ ∈ L ∞ (R 2dk
) for max(α, β) 1, provided that the convergence of the Husimi measures has been established in the special case f = (π) -d/4 e -|x| 2 /2 . If the weak convergence for the Husimi measures only holds in L 1 loc (R dk ), then in addition ϕ has to tend to zero at infinity. 2.5. Proof of Theorem 2.1. We are now able to write the proof of Theorem 2.1.

Let f ∈ C ∞ c (R d ) be a real-valued L 2 -normalized function. As we have already explained in the beginning of Section 2.2, up to extraction of a subsequence, we have ˆMk m

(k) f,Ψ N ϕ → ˆMk m (k) f ϕ (2.28)
for every k 1 and every

ϕ ∈ L 1 (M k ) + L ∞ ε (M k )
. By Theorem 2.7, there exists a probability measure P on B such that

m (k) f = ˆB µ ⊗k dP(µ)
and therefore the convergence (2.1) follows for this f .

Next we establish the independence of the weak convergence with respect to f . This is the content of the following Lemma 2.8 (The limiting measures do not depend on f ). Let Ψ N be a sequence of normalized fermionic functions and let f ∈ L 2 (R d ) be a normalized function such that m

(k) f,Ψ N ⇀ m (k) f weakly as N → ∞ in the sense of (2.1), for every k 1. Then m (k) g,Ψ N ⇀ m (k) f weakly, for every normalized g ∈ L 2 (R d ).

Proof of Lemma 2.8. It clearly suffices to prove that lim

N →∞ ˆMk ϕ(ξ) m (k) f,Ψ N (ξ) -m (k) g,Ψ N (ξ) dξ = 0
for all functions ϕ in a dense set. It will be convenient to take ϕ of the form

ϕ(x 1 , p 1 , ..., x k , p k ) = u 1 (x 1 )v 1 (p 1 ) • • • u k (x k )v k (p k )
with u j , v j in the Schwartz class.

First we remark that we can assume that f and g are smooth with compact support. Indeed, using the definition of m

(k) f,Ψ N and that a(f -f ) = a * (f -f ) = f -f L 2 , we find that m (k) f,Ψ N -m (k) f ,Ψ N L ∞ (M k ) 2k f -f L 2 .
Therefore we can approximate f and g by smooth functions f and g in L 2 (R d ), at the expense of an error which is uniform in N . Then it suffices to prove the result for f and g. For simplicity, we just assume that f and g are smooth for the rest of the proof.

By definition of the phase space measure m

(k) f,Ψ N , we have 1 (2π) dk ˆMk ϕ(ξ) m (k) f,Ψ N (ξ) dξ = N (N -1) • • • (N -k + 1) Ψ N , Op k, f (ϕ) ⊗ 1 N -k Ψ N where Op k, f (ϕ) = 1 (2π) dk ˆMk ϕ(x 1 , p 1 , ..., x k , p k ) P x 1 ,p 1 ⊗ • • • ⊗ P x k ,p k dx 1 • • • dp k . (2.29) We now claim that Op k, f (ϕ) -Op k, g (ϕ) C kd √ = CN -k-1 2d (2.30)
where C depends on the regularity of f , g, ϕ and on k, but not on N . This will imply

ˆMk ϕ(ξ) m (k) f,Ψ N (ξ) -m (k) g,Ψ N (ξ) dξ C √ = CN -1 2d
and finish the argument. The proof of (2.30) is based on well-known techniques from semi-classical analysis. It consists in estimating the difference between the above quantization (2.29) of ϕ and another quantization which is independent of f . We provide the full proof here for the convenience of the reader. First we remark that for our specific function ϕ we can write

Op k, f (ϕ) = B f,1 ⊗ • • • ⊗ B f,k
where

B f,j = 1 (2π) d ˆM u j (x)v j (p)|f x,p f x,p | dx dp. Since by (2.7) we have B f,j d ||u j || L ∞ (R d ) ||v j || L ∞ (R d )
and the same estimate for B g,j , we deduce that

Op k, f (ϕ) -Op k, g (ϕ) d(k-1) k j=1 B f,j -B g,j ℓ =j ||u ℓ || L ∞ (R d ) ||v ℓ || L ∞ (R d ) .
Therefore we only have to prove that

B f,j -B g,j C d √
for all j = 1, ..., k. We will actually prove that

B f,j -d u j (x)v j (-i ∇) C d √ (2.31)
which will of course imply the previous bound. We ignore the index j for simplicity and compute the kernel

B f (y, y ′ ) = 1 (2π) d d 2 ˆRd ˆRd u(x)v(p)f y -x √ f y ′ -x √ e i p•(y-y ′ ) dx dp = 1 (2π) d 2 d 2 v y -y ′ ˆRd u(x)f y -x √ f y ′ -x √ dx
where v is the inverse Fourier transform of the function v. Now the idea is to replace u(x) by u(y) and f ((y -x)/ ) by f ((y ′ -x)/ ) and to then use that ´Rd f 2 = 1. Assuming that v and f have their support in a ball of radius R, we find

B f (y, y ′ ) -(2π) -d 2 u(y) v y -y ′ R √ (2π) d 2 v y -y ′ (||∇u| | L ∞ + ||∇f || L ∞ ||f || L 1 ||u|| L ∞ ) .
The operator with kernel |v((y-y ′ )/ )| has the operator norm

||v(•/ )|| L 1 (R d ) = d ||v|| L 1 (R d )
and this concludes the proof of (2.31), hence of the lemma.

Together with the content of Section 2.4 dealing with the link between the Husimi and Wigner measures, this concludes the proof of the first part of Theorem 2.1. If now the kinetic energy of Ψ N satisfies the semi-classical bound (2.3), then by (2.20) we have

(2π) -d ˆM |p| 2 m (1) f,Ψ N (x, p) dx dp C + N -1/d ˆRd |∇f | 2 .
Therefore the sequence is tight in the p variables and the weak-limit (2.1) remains valid for functions ϕ(x 1 , ..., x k , p 1 , ...,

p k ) = U (x 1 , ..., x k ) with U ∈ L ∞ c (R dk )
. By definition this is the weak convergence

ρ m (k) f,Ψ N ⇀ ˆB ρ ⊗k µ dP(µ), for every k 1. Now, ρ (k)
Ψ N /N k is bounded in L 1 (R dk ) and satisfies ρ m (k) f,Ψ N = k! N k ρ (k) Ψ N * (|f | 2 ) ⊗k by Lemma 2.4. Since |f | 2 ⇀ δ, this proves that k! N k ρ (k) Ψ N ⇀ ˆB ρ ⊗k µ dP(µ)
weakly- * on C 0 0 (R dk ). For k = 1, the Lieb-Thirring inequality [START_REF] Lieb | Bound on kinetic energy of fermions which proves stability of matter[END_REF][START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Stability of Matter in Quantum Mechanics[END_REF] gives us ˆRd ρ

(1)

Ψ N (x) 1+2/d dx C Ψ N N j=1 (-∆) x j Ψ N CN 1+2/d
and, therefore, ρ

Ψ N /N is bounded in L 1 ∩ L 1+2/d . It must thus converge weakly in that space. This concludes the proof of Theorem 2.1.

Convergence of the energy: proof of Theorem 1.1

This section is devoted to the proof of Theorem 1.1 and to the derivation of some a priori estimates that will be useful in the remaining of the article.

3.1. Upper bound on E(N ) using Hartree-Fock theory. Proposition 3.1 (Upper bound). Assume that w is even, that w,

|A| 2 , V -∈ L 1+d/2 (R d ) + L ∞ ε (R d ) and that V + ∈ L 1 loc (R d ). We have lim sup N →∞ E(N ) N e V TF (1)
.

Proof. Taking a Hartree-Fock state Ψ = (N !) -1/2 det(f i (x j )) and introducing the corresponding density matrix

γ = N j=1 |f j f j |, with ρ γ := N j=1 |f j | 2 , we see that E(N ) inf γ 2 =γ=γ * Tr γ=N Tr (-i ∇ + A) 2 γ + ˆRd V ρ γ + 1 2N ¨Rd ×R d w(x -y) ρ γ (x)ρ γ (y) -|γ(x, y)| 2 dx dy .
We first estimate the exchange term (the last term in the previous formula). We write

w = w 1 + w 2 with w 1 ∈ L 1+d/2 (R d ) and w 2 ∈ L ∞ (R d ). Then ¨Rd ×R d |w 2 (x -y)| |γ(x, y)| 2 dx dy ||w 2 || L ∞ Tr (γ 2 ) = ||w 2 || L ∞ N.
We have used here that Tr γ 2 = Tr γ = N . For the term involving w 1 , we use the Hölder and Gagliardo-Nirenberg-Sobolev inequalities to obtain

ˆRd |w 1 (x)| |f (x)| 2 dx ||w 1 || L 1+d/2 (R d ) ||f | | 2 L 2+4/d (R d ) ||w 1 || L 1+d/2 (R d ) ε ||∇f | | 2 L 2 (R d ) + C ε d/2 ||f || 2 L 2 (R d )
for all ε > 0. Applying this in the variable x with y fixed, this gives

¨Rd ×R d |w 1 (x -y)| |γ(x, y)| 2 dx dy ||w 1 || L 1+d/2 (R d ) εTr (-∆)γ + CN ε d/2 .
Combining the two inequalities, we have shown that

1 2N ¨Rd ×R d w(x -y) |γ(x, y)| 2 dx dy ε N Tr (-∆)γ + C(1 + ε -d/2 )
for all ε > 0. In dimensions d 3, we can simply take ε = 1. In dimensions d = 1 and d = 2 we take for instance ε = N d-2 (log N ) -1 . In all cases, we arrive at

1 2N ¨Rd ×R d w(x -y) |γ(x, y)| 2 dx dy ε N N 2/d Tr (-∆)γ + N ε N
for some ε N → 0 which is independent on γ. Finally, in order to relate Tr (-∆)γ to the magnetic kinetic energy, we may use that

| -i ∇ + A| 2 = -2 ∆ -i ∇ • A -i A • ∇ + |A| 2 - 2 2 ∆ -|A| 2 (3.1)
from which we deduce the final upper bound lim sup

N →∞ E(N ) N lim sup N →∞ inf γ 2 =γ=γ * Tr γ=N 1 + 2ε N N Tr | -i ∇ + A| 2 γ + 1 N ˆRd V + 2ε N |A| 2 ρ γ + 1 2N 2 ¨Rd ×R d w(x -y)ρ γ (x)ρ γ (y) dx dy . (3.2)
By using semi-classical analysis, we now construct an appropriate trial state for the variational problem on the right side, based on the classical probability density on phase space.

γ N := 1 (-i ∇ + A) 2 C R -c TF ρ(x) 2/d 0 where (-i ∇ + A) 2 C R is the magnetic Dirichlet Laplacian in the cube and c TF = 4π 2 (d/|S d-1 |) 2/d . Then we have lim N →∞ N -1 Tr (-i ∇ + A) 2 γ N = d d + 2 c TF ˆRd ρ(x) 1+2/d dx (3.3) and lim N →∞ N -1 Tr γ N = ˆRd ρ(x) dx. (3.4) Furthermore, ρ γ N /N ⇀ ρ weakly in L 1 (R d ) and weakly- * in L ∞ (R d ).
The same properties all hold true if γ N is replaced by the projection γ N onto the N lowest eigenvectors of

(-i ∇ + A) 2 C R -c TF ρ(x) 2/d . Remark 3.3.
The lemma is true under much weaker assumptions on A, but we have used the same hypothesis as in the rest of the paper for convenience.

The lemma is proved in Appendix A and the proof of Proposition 3.1 now follows immediately. Indeed, let ρ ∈ C ∞ c (R d ) and γ N be as in the lemma (extended by zero outside of C R ). Then, by (3.2) lim sup

N →∞ E(N ) N lim N →∞ 1 + 2ε N N Tr | -i ∇ + A| 2 γ N + 1 N ˆRd (V + 2ε N |A| 2 )ρ γ N + 1 2N 2 ¨Rd ×R d w(x -y)ρ γ N (x)ρ γ N (y) dx dy = d d + 2 4π 2 d |S d-1 | 2/d ˆRd ρ 1+2/d + ˆRd V ρ + 1 2 ¨Rd ×R d w(x -y)ρ(x)ρ(y) dx dy.
Here we have used that

N -1 ˆRd V ρ γ N = N -1 ˆCR V ρ γ N → ˆCR V ρ, that ´|A| 2 ρ γ N is uniformly bounded and that N -2 ¨CR ×C R w(x-y)ρ γ N (x)ρ γ N (y) dx dy → ¨CR ×C R w(x-y)ρ(x)ρ(y) dx dy since ρ γ N /N converges weakly in L 1 (R d
) and weakly- * in L ∞ (R d ), and V, w, |A| 2 ∈ L 1 loc (R d ) by assumption.

3.2.

A priori estimates. We now derive some a priori bounds on any (normalized) sequence {Ψ N } of N -particle states with Ψ N , H N Ψ N = O(N ), which will be useful for the lower bound.

Lemma 3.4 (Lieb-Thirring

). There exists a constant C such that

H N N j=1 - ∆ j 4N 2 d + V + (x j ) -CN -CN. (3.5)
In particular, any

{Ψ N } such that Ψ N , H N Ψ N = O(N ) must satisfy 1 N 1+ 2 d Ψ N ,   N j=1 -∆ j   Ψ N + ˆRd ρ (1)
Ψ N (x) N 1+ 2 d dx + ˆRd V + ρ (1)
Ψ N (x) N dx C, (3.6) with ρ (1)
Ψ N the one-particle densities defined in (2.14). In addition, for any potential

0 f = f 1 + f 2 ∈ L 1+d/2 (R d ) + L ∞ (R d ), we have 1 N ˆRd f (x)ρ (1)
Ψ N (x) dx + 1 N 2 ˆRd ˆRd f (x -y)ρ (2)
Ψ N (x, y) dx dy

C ||f 1 || L 1+d/2 (R d ) + ||f 2 || L ∞ (R d ) . (3.7)
Proof. First we get rid of the magnetic field at the expense of a factor in front of the kinetic energy, by using the simple bound

| -i ∇ + A| 2 - 2 2 ∆ -|A| 2
as in (3.1). Our assumption on A implies that

|A| 2 ∈ L 1+d/2 (R d ) + L ∞ ε (R d )
, just as for V -. We may therefore write

V -+ |A| 2 = V 1 + V 2 ∈ L 1+d/2 (R d ) + L ∞ (R d ).
The estimate on V 2 being obvious, we use the Lieb-Thirring inequality [START_REF] Lieb | Bound on kinetic energy of fermions which proves stability of matter[END_REF][START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Stability of Matter in Quantum Mechanics[END_REF] to obtain

N j=1 - ∆ j 8N 2 d -V 1 (x j ) - 1 4 Tr - ∆ j N 2 d -4V 1 - -CN ˆRd V 1+d/2 1
.

(3.8) In a similar manner, we note that |Ψ| 2 is symmetric and rewrite, following Lévy-Leblond [START_REF] Lévy-Leblond | Nonsaturation of Gravitational Forces[END_REF],

Ψ,   1 N 1 k<ℓ N w -(x k -x ℓ )   Ψ = N -1 2 Ψ, w -(x 1 -x 2 )Ψ = 1 2 Ψ,   2 j N w -(x 1 -x j )   Ψ .
We now use that Ψ(x 1 , •) is antisymmetric in the remaining N -1 variables and the above Lieb-Thirring estimate in those variables. Since

w -= w 1 + w 2 ∈ L 1+d/2 (R d ) + L ∞ (R d ), we obtain as before N j=1 - ∆ j 8N 2 d - 1 N 1 k<ℓ N w -(x k -x ℓ ) -CN ˆRd w 1+d/2 1 + ||w 2 || L ∞ .
(3.9) This concludes the proof of (3.5). The estimate (3.6) follows from the Lieb-Thirring inequality written in terms of the density

Ψ N ,   N j=1 -∆ j   Ψ N C ˆRd ρ (1) Ψ N 1+2/d , (3.10)
which is dual to (3.8). In (3.7), the one-body part follows directly from (3.6) and Hölder's inequality. For the two-body estimate, we use (3.9) with w - replaced by f /ε and obtain

ε N j=1 - ∆ j 8N 2 d - 1 N 1 k<ℓ N f (x k -x ℓ ) -CN ε -d/2 ˆRd f 1+d/2 1 + ||f 2 || L ∞ .
Taking the expectation value in the state Ψ N and using (3.6), we find (3.7) after optimizing over ε.

3.3.

Lower bound and end of the proof of Theorem 1.1. Here we establish the lower bound, which concludes the proof of Theorem 1.1. Proof. The first step is to regularize the interaction potential w. For every ε > 0, we may write w = f + g with f ∈ L 1+d/2 (R d ) and g L ∞ ε. We then consider the problem E ′ (N ) in which w has been replaced by another potential w ′ , and a state Ψ N such that Ψ N , H N Ψ N = E(N ) + o(N ). Then by (3.7), we obtain Ψ N ,

1 j<k N w(x j -x k )Ψ N Ψ N , 1 j<k N w ′ (x j -x k )Ψ N -CN f -w ′ L 1+d/2 (R d ) - N -1 2 ε
and, therefore,

E(N ) E ′ (N ) -CN f -w ′ L 1+d/2 (R d ) - N -1 2 ε. (3.12)
A similar bound holds for the Thomas-Fermi model and we conclude by an 'ε/2 argument' that it suffices to prove (3.11) for a smooth interaction potential w ′ which approximates f . For simplicity of notation, we will assume for the rest of the proof that w is itself smooth enough. The exact property we need is that w ∈ L 1 (R d ), which implies that w is a continuous bounded function that tends to zero at infinity. Now, we write w = w 1 -w 2 where w 1 = ( w) + and w 2 = ( w) -which are also in L 1 (R d ). Note that w 1 and w 2 are in addition both even since w is real. We will use the fact that Ψ N ,

1 j<k N w 1,2 (x j -x k )Ψ N = N (N -1) 2 Ψ N , w 1,2 (x ℓ -x m )Ψ N
(3.13) for any ℓ = m, by the symmetry of |Ψ| 2 . The idea is now to split the N particles into two groups of, respectively M and L particles, with L = ⌊ √ N ⌋. The first group will just be x 1 , ..., x M whereas the second (smaller) group will be denoted for convenience as y 1 = x M +1 , ..., y L = x N . Then, we express the repulsive part w 1 of the potential w using only the particles in the first group and write the attractive part w 2 as an interaction between the two groups. This means that we rewrite, using (3.13), Ψ N ,

1 j<k N w 1 (x j -x k )Ψ N = N (N -1) M (M -1) Ψ N , 1 m<m ′ M w 1 (x m -x m ′ )Ψ N
on one hand, and

-Ψ N , 1 j<k N w 2 (x j -x k )Ψ N = N (N -1) Ψ N , w 2 (y 1 -y 2 ) 2 -w 2 (x 1 -y 1 ) Ψ N = N (N -1) L(L -1) Ψ N , 1 ℓ<ℓ ′ L w 2 (y ℓ -y ℓ ′ )Ψ N - N (N -1) LM Ψ N , M m=1 L ℓ=1 w 2 (x m -y ℓ )Ψ N
on the other hand. Expressing the one-particle terms in a similar manner using only the M particles of the first group, we deduce that

Ψ N , H N Ψ N N = Ψ N , HΨ N M where H = M m=1 | -i ∇ m + A(x m )| 2 + V (x m ) + 1 -1/N M -1 1 m<m ′ M w 1 (x m -x m ′ ) + M (1 -1/N ) L(L -1) 1 ℓ<ℓ ′ L w 2 (y ℓ -y ℓ ′ ) - (1 -1/N ) L M m=1 L ℓ=1 w 2 (x m -y ℓ ).
This new Hamiltonian corresponds to a system with M ∼ N quantum particles which repel through the potential w 1 and L ∼ √ N classical particles that repel through the potential w 2 , with an additional attraction between the two species given by -w 2 . In other words, we have transformed the attractive part w 2 of w into an interaction with an auxiliary system of L particles that repel each other. A similar approach was used by Lieb-Thirring [START_REF]Gravitational collapse in quantum mechanics with relativistic kinetic energy[END_REF] and Lieb-Yau [START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF] and it was in turm inspired by arguments of Levy-Leblond [START_REF] Lévy-Leblond | Nonsaturation of Gravitational Forces[END_REF] and Dyson-Lenard [START_REF] Dyson | Stability of matter. I[END_REF]. For bosons the argument was explained in [START_REF]Mean-field limit of Bose systems: rigorous results[END_REF] (there one can take M and L of the order N , which simplifies a bit the argument). From the antisymmetry of Ψ N in the M first variables, we conclude that

Ψ N , H N Ψ N N inf y 1 ,...,y L inf σ M 1 L 2 (R d ) ( H) M
and it therefore remains to estimate the bottom of the spectrum of H for M fermions, uniformly with respect to the positions y 1 , ..., y L of the L classical particles. The rest of the argument will then be based on the following wellknown lemma, which allows to bound from below the two-body interaction by a one-particle term.

Lemma 3.6 (Estimating the two-body potential by a one-particle term).

Let f be a function such that f is in L 1 (R d ) and non-negative, and denote

D f (g, g) := ˆRd ˆRd g(x)g(y)f (x -y) dx dy = (2π) d 2 ˆRd | g(k)| 2 f (k) dk 0.
Then we have

1 k<k ′ K f (z k -z k ′ ) - f (0) 2 K + K k=1 f * η(z k ) - 1 2 D f (η, η) (3.14)
for every η ∈ L 1 (R d ) and every z 1 , ..., z K ∈ R d , and Ψ,

1 k<k ′ K f (z k -z k ′ ) Ψ - f (0) 2 K + 1 2 D f (ρ (1) Ψ , ρ (1) Ψ
).

(3.15)

for every K-particle state Ψ.

Proof of Lemma 3.6. Use that D f (g, g) 0 for g = K k=1 δ z k -η and then take η = ρ (1) Ψ . Now, let Ψ be any fermionic M -particle state. Using (3.15) with f = w 1 , we get

1 -1/N M -1 Ψ, 1 m<m ′ M w 1 (x m -x m ′ ) Ψ 1 -1/N 2(M -1) D w 1 (ρ (1) Ψ , ρ (1)
Ψ ) - M (1 -1/N ) 2(M -1) w 1 (0) 1 2M D w 1 (ρ (1) Ψ , ρ (1) Ψ
) -M 2(M -1) w 1 (0).

Using now (3.14) with f = w 2 and η = ρ (L -1)/M , we find

M (1 -1/N ) L(L -1) 1 ℓ<ℓ ′ L w 2 (y ℓ -y ℓ ′ ) - (1 -1/N ) L Ψ, M m=1 L ℓ=1 w 2 (x m -y ℓ ) Ψ - (L -1)(1 -1/N ) 2M L D w 2 (ρ (1) Ψ , ρ (1)
Ψ ) - M (1 -1/N ) 2(L -1) w 2 (0) - 1 2M D w 2 (ρ (1) Ψ , ρ (1)
Ψ ) - M 2(L -1)
w 2 (0).

We now fix a real-valued normalized function f ∈ L 2 (R d ) and use formulas (2.19) and (2.17) to obtain

Tr (| -i 2 ∇ + A| 2 + V )γ N = 1 (2π ) d ˆR2d |p + A(x)| 2 + V (x) m N (x, p) dx dp -N ˆRd |∇f | 2 -2Tr (A -A * |f | 2 ) • (-i ∇)γ N + ˆRd (|A| 2 -|A| 2 * |f | 2 + V -V * |f | 2 )ρ γ N and D w (ρ γ N , ρ γ N ) = -2d D w (ρ m N , ρ m N ) + D w-w * |f | 2 * |f | 2 (ρ γ N , ρ γ N )
where m N (x, p) = f x,p , γ N f x,p and ρ m N = (2π) -d ´Rd m N (x, p) dp. Note that 0 m N 1 and (2π) -d ˜mN (x, p) dx dp = 1, hence m N is a suitable trial function for the Vlasov problem. Using the bound (3.17) we have

Tr (A -A * |f | 2 ) • (-i ∇)γ -||(-i ∇) √ γ N || S 2 √ γ N (A -A * |f | 2 ) S 2 = - 2 Tr (-∆)γ N ˆRd (A -A * |f | 2 ) 2 ρ γ N -CN A -A * |f | 2 L 2+d (R d )
.

We can argue similarly for the other terms (when

V ∈ C ∞ (R d ) is equal to M outside of a large ball B R , we have to use that V -V * |f | 2 = (V -M ) -(V -M ) * |f | 2
has compact support and converges to 0 strongly in L 1+d/2 (R d )). We obtain Let {Ψ N } be as in the statement, that is, such that Ψ N , H N Ψ N = E(N ) + o(N ). By the estimates in Lemma 3.4, we know that Ψ N satisfies the semi-classical kinetic energy bound

1 N Tr (| -i 2 ∇ + A| 2 + V )γ N + 1 2N 2 D w (ρ γ N , ρ γ N ) 1 (2π) d N d ˆR2d |p + A(x)| 2 + V (x) m N (x, p) dx dp + 1 2N 2 2d D w (ρ m N , ρ m N ) + o(1). Since N d → 1,
Ψ N ,   N j=1 -∆ j   Ψ N CN 1+2/d . (4.1)
Finally, using the same argument as in the proof of Lemma 3.7 for the magnetic kinetic energy, we also deduce from Fatou's lemma that lim inf

N →∞ Ψ N , 1 N N j=1 | -i ∇ j + A(x j )| 2 Ψ N = lim inf N →∞ 1 (2π) d ¨R2d |p + A(x)| 2 m (1) f,Ψ N (x, p) dx dp ˆS 1 (2π) d ¨R2d |p + A(x)| 2 µ(x, p) dx dp dP(µ).
Therefore, we have proved that

e V TF (1) lim N →∞ E(N ) N = lim N →∞ Ψ N , H N Ψ N N ˆS E V,A Vla (µ) dP(µ).
The minimum of the Vlasov energy is precisely e V TF (1) and since P is a probability, we have

ˆS E V,A Vla (µ) dP(µ) = ˆS E V,A Vla (µ) -e V TF (1)
0 dP(µ) + e V TF (1)
and there must be equality everywhere. This implies that P is supported on the minimizers of E V,A Vla . These are all of the form m ρ (x, p) = 1(|p + A(x)| 2 c TF ρ(x) 2/d ) where ρ minimizes the Thomas-Fermi energy. Hence, P induces a probability density P(ρ) on the set of minimizers of the Thomas-Fermi energy and this concludes the proof of Theorem 1.2. 4.2. Unconfined case: proof of Theorem 1.3. The proof in the unconfined situation is much more subtle and we will follow here ideas from [START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF], based on localization methods in Fock space [START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF].

Let χ be a smooth function on

R + with 0 χ 1, χ |[0,1] ≡ 1 and χ |[2,∞) ≡ 0. Denote χ R (x) = χ(|x|/R) and η R = (1 -χ 2 R ) 1/2 . Let G - R,n = N n (χ R ) ⊗n Tr n+1,...,N (η R ) ⊗N -n |Ψ N Ψ N |(η R ) ⊗N -n (χ R ) ⊗n
be the projection onto the n-particle space of the many-particle state localized using χ R [START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF], and, in a similar fashion,

G + R,n = N n (η R ) ⊗n Tr n+1,...,N (χ R ) ⊗N -n |Ψ N Ψ N |(χ R ) ⊗N -n (η R ) ⊗n .
We remark that

Tr G + R,n = Tr G - R,N -n (4.4)
and that the k-particle density matrices are localized in the usual sense: We first get a uniform estimate on ρ γ N . Indeed, using that We have therefore shown that ρ γ N /N is bounded in L ∞ (C R ). Up to extraction of a subsequence, we may assume that ρ γ N /N ⇀ ρ weakly. Now we use that 0 -

(χ R ) ⊗k γ (k) Ψ N (χ R ) ⊗k = N n=k γ (k) G - n,R , (η R ) ⊗k γ (k) Ψ N (η R ) ⊗k = N n=k γ (k) G + n,R . (4
1 N Tr (-iN -1/d ∇ + A) 2 C R -c TF ρ - = 1 N Tr ((-iN -1/d ∇ + A) 2 -c TF ρ)γ N = 1 N Tr (-iN -1/d ∇ + A) 2 C R γ N -c TF 1 N ˆCR ρρ γ N from which we deduce that 1 N Tr (-iN -1/d ∇ + A) 2 C R γ N c TF 1 N ˆCR ρρ γ N C ||ρ|| L 1 ||ρ|| d/2 L ∞ ,
due to the uniform upper bound on ρ γ N /N . We then introduce the semiclassical measure m N (x, p) = f x,p , γ N f x,p and call its weak limit m. Arguing as before, we obtain lim inf The result (A.2) is standard for smooth vector potentials A. Let A ε be a smooth approximation of A in L 2 . Using the inequality

N →∞ 1 N Tr (-iN -1/d ∇ + A) 2 C R γ N -c TF 1 N ˆCR ρρ γ N
(1 -δ)(-iN -1/d ∇ + A ε ) 2 -δ -2 |A -A ε | 2 (-iN -1/d ∇ + A) 2 (1 + δ)(-iN -1/d ∇ + A ε ) 2 + δ -2 |A -A ε | 2 ,
and the uniform upper bound on ρ γ N /N , the result follows for general A. So we consider the Weyl asymptotics, Finally, the results are all the same for an orthogonal projection γ N on N first eigenfunctions of (-iN -1/d ∇ + A) 2 C R -c TF ρ since Tr |γ N -γ N | = o(N ) by (A.3) and therefore ρ γ N -ρ γ N L 1 = o(N ). For N large we have

1 (-iN -1/d ∇+A) 2 C R ρ(x)-ε γ N 1 (-iN -1/d ∇+A) 2 C R ρ(x)+ε
since, by the above arguments with ρ replaced by ρ ± ε,

Tr 1 (-iN -1/d ∇ + A) 2 C R ρ(x) ± ε ∼ (1 ± εR d )N.
From the above estimates we conclude that ρ γ N /N is bounded in L ∞ , and therefore ρ γ N /N has the same weak limit as ρ γ N in L 1 ∩ L ∞ . We also have

E 2

 2 V,A Hartree (u) = ˆRd | -i ∇u(x) + A(x)u(x)| 2 + V (x)|u(x)| 2 dx + 1 ˆRd ˆRd w(x -y)|u(x)| 2 |u(y)| 2 dx dy. (1.2)

Theorem 2 . 7 (

 27 Weak fermionic semi-classical measures on phase space). Let m (N ) N be a sequence of symmetric positive densities in L 1 (M N ), with M ⊂ R D , and let m (k)

(k) N 1

 1 for every 1 k N and that m (k) N ⇀ m (k) weakly in L 1 (M k ) and weakly- * in L ∞ (M k ) for every fixed k 1, as N → ∞. Then there exists a Borel probability measure P on the set B := {µ ∈ L 1 (M ) : 0 µ 1, c ´M µ 1} such that

Lemma 3 . 2 (

 32 Semi-classical analysis in a compact domain). Let ρ 0 be a fixed function in C ∞ c (R d ) with support in the cube C R = (-R/2, R/2) d , such that ρ 0 and ´CR ρ = 1. Let A ∈ L 2+d (C R ) be a magnetic vector potential. Define

Proposition 3 . 5 (

 35 Lower bound). Under the assumptions of Theorem 1

4 .

 4 we obtain the Vlasov energy on the right and the proof of Lemma 3.7 is complete. This now concludes the proof of the lower bound (3.11), hence of Theorem 1.1. Convergence of states: proof of Theorems 1.2 and 1.3

. 5 )Lemma 4 . 1 .χ 2 R

 5412 When n → ∞ and n/N → λ > 0, then a simple adaptation of Theorem 1w(x -y)ρ(x)ρ(y) dx dy = e 0 TF (λ). Arguing as in[44, p. 613-614], it is not difficult to see that lim - R,n ). Now the proof can be concluded based on the following Lemma, inspired of [44, Thm. 2.6]. For every continuous function f on [0, 1], we have lim ρ m dP(m)and therefore we have proved that limN →∞ Ψ N , H N Ψ N N ˆB E V,A Vla χ 2 R (x)m + e 0 TF 1 -ˆRd χ 2 R ρ m dP(m) + ε R .Taking finally the limit R → ∞ gives the final estimate

(4. 8 )Ψ

 8 Since E V,A Vla (m) e VTF ´Rd ρ m and since e V TF (1) e V TF (λ) + e V TF (1 -λ) for every 0 λ 1, we conclude that P must be supported on the set 0 m 1 :ˆRd ρ m 1, E V,A Vla (m) = e V TF ˆRd ρ m = e V TF (1) -e 0 TF 1 -ˆRd ρ mwhich can easily be seen to be the set of weak limits of minimizing sequences for the variational problem e V,A Vla (1). The probability measure P induces a probability measure P on the set M of the corresponding ρ's and the proof of Theorem 1.3 is finished.Proof of Lemma 4.1. The proof is the same as [44, Thm.2.6]. For k N (x 1 , ..., x k ) dx 1 • • • dx k -1) • • • (n -k + 1) N k Tr (G - R,n ) has the same limit as N n=1 (n/N ) k Tr (G - R,n). Therefore we have proved the lemma for f (x) = x k , for all k 1. The result for any continuous function f follows from the density of polynomials on [0, 1], together with the fact that N n=0 Tr (G - R,n ) = 1 for every N , by definition of the localized state. Appendix A. Proof of Lemma 3.2

γ N = 1 (.

 1 (-iN -1/d ∇ + A) 2 C R -c TF ρ 0) e -β (-iN -1/d ∇+A) 2 C R -c TF ρ ,for any β > 0, we deduce thatρ γ N (x) |e -β((-iN -1/d ∇+A) 2 C R -c TF ρ) (x, x)|.From the Feynman-Kac formula and the diamagnetic inequality, we have|e -β (-iN -1/d ∇+A) 2 C R -c TF ρ (x, y)| e c TF β||ρ|| L ∞ e N -2/d β∆ C R (x, y),where ∆ C R is the (non-magnetic) Dirichlet Laplacian on C R , and henceρ γ N (x) 2 R d e c TF β||ρ|| L ∞ k 1 ,...,k d ∈(π/R)N\{0} e -N -2/d β d j=1 |k j | 2 TF β||ρ|| L ∞ N π -d β -d/2 ˆR e -|k| 2 dk d

1 (→∞ - 1 N

 11 2π) d ˆRd ˆRd |p + A(x)| 2 -c TF ρ(x) m(x, p) dp dx.We now use the well-known Weyl asymptotics for the energy, i.e.lim N Tr (-iN -1/d ∇ + A) 2 C R -c TF ρ(x) - = -1 (2π) d ˆRd ˆRd |p + A(x)| 2 -c TF ρ(x) -dx dp. (A.2)

1 N 1 ˆCR ρ γ N → ˆCR ρ = 1 . (A. 3 )→∞ 1 N

 11131 Tr (-iN -1/d ∇ + A) 2 C R -c TF ρ(x) - = -1 (2π) d ˆRd ˆRd |p + A(x)| 2 -c TF ρ(x) -dx dp = inf 0 m ′ 1 1 (2π) d ˆRd ˆRd |p + A(x)| 2 -c TF ρ(x) m ′ (x, p) dx dp. with unique minimizer m ′ (x, p) = 1(|p + A(x)| 2 -c TF ρ(x) 0) in L ∞ (R d × R d ), we conclude that m = 1(|p + A(x)| 2 -c TF ρ(x) 0) a.e. This gives in particular thatN -1 ρ γ N (x) ⇀ ρ m (x) = 1 (2π) d ˆRd 1(|p + A(x)| 2 -c TF ρ(x) dp = ρ(x) weakly in L 1 ∩ L ∞ , hence that N -1 Tr (γ N ) = N -The latter says that (-iN-1/d ∇ + A) 2 C R -c TF ρ has N + o(N ) negative eigenvalues.From the above limits we also have as desiredlim N Tr (-iN -1/d ∇ + A) 2 C R γ N = 1 (2π) d ˆRd ˆRd |p| 2 1(|p| 2 -c TF ρ(x) dx dp = d d + 2 4π 2 d |S d-1 | 2/d ˆRd ρ(x) 1+2/d dx.

  Tr ((-iN-1/d ∇ + A) 2 C R -c TF ρ)(γ N -γ N ) = o(N ) which implies that Tr (-iN -1/d ∇ + A) 2 C R (γ N -γ N ) = o(N 1+2/d) and concludes the proof of Lemma 3.2.

  1.1, 1.2 and 1.3. The proof of Theorem 1.1 is based on ideas of Lieb and Yau[START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF], while that of Theorem 1.2 follows easily from our fermionic de Finetti-Hewitt-Savage theorem. The proof of Theorem 1.3 is more tedious and relies on the techniques introduced in[START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF][START_REF] Lewin | Derivation of Hartree's theory for generic mean-field Bose systems[END_REF].

	Acknowledgement. M.L. and J.P.S acknowledge financial support from the
	European Research Council (Grant Agreements MNIQS 258023 and MAS-
	TRUMAT 321029). S.F. acknowledges support from a Danish research coun-
	cil Sapere Aude grant. This work was started when the authors were at the

Centre Émile Borel of the Institut Henri Poincaré in Paris in 2013. Part of this work was done when S.F. was a visiting professor at the University Paris-Dauphine.

Recall that Ψ only depends on the x m 's. We have therefore proved the following lower bound, independent of the y ℓ 's,

where we recall that M = N -⌊ √ N ⌋ and that = N -1/d ≃ M -1/d . The term on the right is (up to the last constant) the reduced Hartree-Fock (rHF) energy of the M particles, that is well-known to converge to the Thomas-Fermi problem in the limit M ≃ N → ∞ [START_REF] Lieb | Thomas-Fermi theory revisited[END_REF][START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF][START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF]. Even if we only need the lower bound, we state and prove the full convergence. Lemma 3.7 (Reduced-Hartree-Fock model). Under the same assumptions on V , A and w as in Theorem 1.1, we have we deduce that any appropriate minimizing sequence γ N satisfies

, these estimates can be used to replace the potentials V , A and w by functions in C ∞ c (R d ). When V + → +∞ at infinity, we first bound from below V + by the truncated potential

and take M → ∞ at the very end of the argument, using that lim

The previous estimates then allow to replace A and w by functions in

), which is equal to M outside of a sufficiently large ball. This implies by Theorem 2.1 that, after extraction of a (not displayed) subsequence,

weakly for all k 1, where P is a probability measure on

Therefore, there only remains to show that P has its support in the set of minimizers of the Vlasov energy (in the confined case) and in the set of weak limits of minimizing sequences (in the unconfined case).

4.1. Confined case: proof of Theorem 1.2. We start with the much simpler confined case where V + → +∞ at infinity. In that case, the argument below actually gives directly the energy lower bound, and the arguments of Section 3.3 are not needed. From the bound (3.6) we have

Ψ N (x, y) dx dy C and infer that ρ

Ψ N /N and ρ

Ψ N /N 2 are tight at infinity. From this we infer that there is no loss of mass at infinity and that P has its support on

Using (4.3) and the tightness property at infinity, we are able to pass to the limit in the interaction term: lim

In order to pass from the second to the third line, one must first use the bound (3.7) to replace the interaction potential w by a bounded function, and then use the weak convergence of ρ

(2)

. By Fatou's lemma for V + and the weak convergence of ρ Ψ N for V -, we also deduce that lim inf

Next we split the N -body quantum energy using a partition of unity in the one-particle state as follows

Ψ N (x, y) dx dy

Ψ N (x, y) dx dy

Ψ N (x, y) dx dy + ε N,R , with the localization error

Ψ N (x, y) dx dy.

From the known weak convergence of the densities one has

(the term involving w is treated by first replacing w by a bounded function of compact support, using (3.7), and then by using the weak convergence of ρ

Ψ N ). On the other hand, the terms involving η R may be rewritten using the property (4.5) as