Asymptotic expansion of the mean-field approximation

We established and estimate the full asymptotic expansion in integer powers of 1 N of the [ $\sqrt$ N ] first marginals of N-body evolutions lying in a general paradigm containing Kac models and non-relativistic quantum evolution. We prove that the coefficients of the expansion are, at any time, explicitly computable given the knowledge of the linearization on the one-body meanfield kinetic limit equation. Instead of working directly with the corresponding BBGKY-type hierarchy, we follows a method developed in [22] for the meanfield limit, dealing with error terms analogue to the v-functions used in previous works. As a by-product we get that the rate of convergence to the meanfield limit in 1 N is optimal.

1 N is optimal.

Introduction: motivation and main results
Mean field limit concerns systems of interacting (classical or quantum) particles whose number diverges in a way linked with a rescaling of the interaction insuring an equilibrium between interaction and residual kinetic energies. In the case of an additive one-body kinetic energy part and a twobody interaction, and without taking in consideration quantum statistics, this equilibrium is reached by putting in front of the interaction a coupling constant proportional to the inverse of the number of particles.
The system is then described by isolating the evolution of one (or j) particle(s) and averaging over all the other. This leads to a partial information on the system driven by the so-called j-marginals. The mean field theory insures that the j-marginals tend, as the number of particles diverges, to the j-tensor power of the solution of a non-linear one-body meanfield equation (Vlasov, Hartree,...) issued from the 1-marginal on the initial N -body state. This program has be achieved in many different situations, and the literature concerning the mean field approach is protuberant. We refer to the review article [27] for a reasonable bibliography.
Much less is known about the fluctuations around this limit, namely the correction to be added to the factorized limit in order to get better approximations of the true evolution of the j-marginals.
The identification of the leading order of these fluctuations with a Gaussian stochastic process has been established in the quantum context first in [15] and in the classical one in [5]. For the classical dynamics of hard spheres, the fluctuations around the Boltzmann equation have been computed at leading order in [26], generalizing to non-equilibrium states the results of [3]. More recently, for the quantum case, fluctuations near the Hartree dynamics has been derived in [21] ( after [20]) and in [2] also for the grand canonical ensamble formalism (number of particles non fixed), using in both cases the methods of second quantization (Fock space).
Recently, we developed (together with S. Simonella) in [24] a method to derive mean field limit, alternative to the ones using empirical measures or direct estimates on the "BBGKY-type" hierarchies (systems of coupled equations satisfied by the set of the j-marginals). This method rather uses the hierarchy followed by the "kinetic errors" E j−k (defined below), already used (under the name "v-functions") to deal with kinetic limits of stochastic models [10,7,4,11,12,6,8,13] and recently investigated in the more singular low density limit of hard spheres [25]. These quantities are, roughly speaking, the coefficient of the decomposition of the j-marginal as a linear combination of the kth tensor powers, k = 1, . . . , j, of the solution of the mean-field equation issued from of the 1-marginal of the initial full state. We developed in [24] a strategy suitable in particular for Kac models (homogeneous original one [16,17] and non-homogeneous [9]) and quantum mean field theory. This strategy allowed us to derive the limiting factorization property of the j-marginals up to, roughly speaking, j ≲ √ N . This threshold is, on the other side, the one obtained by heuristic arguments as shown in [24].
Here and in all this article, N denotes the number of particles of the system under consideration.
In the present note we provide and estimate a full asymptotic expansion in powers of 1 N of the difference between the evolution of j-marginals and its factorized leading order form (Corollary 1.5), following a similar result for the kinetic errors E j (t) (Theorem 1.4). Our results are valid for j ≲ √ N in an abstract paradigm, generalization of the abstract formalism developed in [24] and described in Appendix A, and applies of course to the different Kac models and quantum mean field theory treated in [24]. Moreover, we show that the additional knowledge of the linearization of the mean field flow, around the meanfield solution issued from the 1-marginal of the initial data, gives an explicit construction of the full asymptotic expansion of the j-marginals in powers of 1 N . We will state in this section the quantum results and postpone in Section 5 and in the Appendix A the Kac's type and the abstract results, respectively. Sections 2 and 3 contain the algebraic and anlytical proofs of our results in the quantum case, immediately transposable to the Kac and abstract situations as shown in Section 5 and Appendix A. In Section 4 one compute explicitly the first terms of the asymptotic expansions obtained in the quantum case and rely them to previous works.
1.1. Quantum mean-field. Let L 1 (L 2 (R d ) be the space of trace class operators on L 2 (R d ), with their associated norms.
We consider the evolution of a system of N quantum particles interacting through a (real-valued) two-body, even potential V , described for any value of the Planck constant ̵ h > 0 by the Schrödinger equation We will suppose in the whole presnt paper that the N -body Hamiltonian H N is essentially self-adjoint. Instead of the Schrödinger equation written in terms of wave functions, we shall rather consider the quantum evolution of density matrices. An N -body density matrix is an operator F N such that The evolution of the density matrix F N ↦ F N (t) of a N -particle system is governed for any value of the Planck constant ̵ h > 0 by the von Neumann equation equivalent to the Schrödinger equation when F N (0) is a rank one projector, modulo a global phase.
Positivity, norm and trace are obviously preserved by (1) since H N is taken essentially self-adjoint. For each j = 1, . . . , N , the j-particle marginal F N j (t) of F N (t) is the unique operator on H j such that . for all A 1 , . . . , A j bounded operators on H. Alternatively and equivalently, the F N j can be defined by the partial trace of F N on the N − j last "particles": defining F N through its integral kernel It will be convenient for the sequel to rewrite (1) in the following operator form As mentioned before the (essantial) self-adjointness of H N implies that We will denote (6) L ∶= L 1 (L 2 (R d )) so that L ⊗n = L 1 (L 2 (R nd )), n = 1, . . . , N, and, with a slight abuse of notation, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ⋅ 1 the trace norm on any L ⊗j , ⋅ the operator norm on any L(L ⊗i , L ⊗j ) A density matrix F n ∈ L ⊗n is called symmetric if its integral kernel . . , n. Note that the symmetry of F N is preserved by the equation (1) due to the particular form of the potential.
We define, for n = 1, . . . , N , Note that F N j ∈ L ⊗j (F N 0 = 1 ∈ L ⊗0 ∶= C) and F N j > 0, F N j 1 = F N 1 , and obviously F N j is symmetric as F N . That is to say: The family of j-marginals, j = 1, . . . , N , are solutions of the BBGKY hierarchy of equations (see [1]) where: where Tr j+1 is the partial trace with respect to the (j + 1)th variable, as in (2).
The Hartree equation is Note that (15) reads also (16) Since V is bounded, (15) has for all time a unique solution F (t) > 0 and F (t) = 1 (see again [1]).
In order to define the correlation error in an easy way, we need a bit more of notations concerning the variables of integral kernels.
. In [24] was shown that (18) is inverted by the following equality.

Asymptotic expansion and main result of the present article.
Two questions arise naturally: (1) are the estimates (23) sharp?
(2) could (24) be improved with a r.h.s. of any order we wish? 1 We will see below that, indeed, not only the estimates (23) are true, but N j 2 E N j (t) has a full asymptotic expansion in positive powers of ( 1 N ) (actually we will show that this expansion contains only powers ( 1 N ) k 2 with k + j even) if N j 2 E N j (0) do posses such an expansion in half powers of 1 N .
More precisely we will show that, under the hypothesis (22) on the initial data and for all time t and all j = 1, . . . , N , when the same is true at t = 0.
Moreover we will see that all the E ℓ j can be explicitly recursively computed after the knowledge of the linearization of the mean field equation (15) around the solution of (15) with initial condition F (0) = (F N (0)) 1 . Indeed the proof will involve the "j-kinteic linear mean field flow" defined by the linear kinetic mean field equation of order j: (27) is solved by the two parameter semigroup U 9 U 0 j (s, s) = I. Note that U 0 j exists since K j generates a unitary flow and ∆ j is bounded. The reason of the terminology comes form the fact that, as shown by (37), ∆ 1 = Q(F, ⋅)+Q(⋅, F ) so that, for j = 1, (27) is the linearization of the mean field equation (15) around its solution F (t). Note moreover that, for G 1 , G 2 ∈ L, and therefore . More generally, if P j ∶ L j → L ⊗j is any homogeneous polynomial invariant by permutations, That is: U 0 J drives each G j along the linearized mean field flow "factor by factor". Denoting by L ⊗j sym the subspace of symmetric (by permutations) vectors, we just prove the following result.
Note also that U 0 j (t, s) is given by a convergent Dyson expansion and that, by the isommetry of the flow generated by K j and the bound (54) below, we have by Gronwall Lemma that U o j (t, s) ≤ e j t−s . We will also need in the sequel the semigroup defined by U j (t) exists by the same argument as for U 0 j (t). Moreover, in the regime j 2 N small, U j can be also computed out of U 0 j by a convergent perturbation expansion (in j 2 N ). Indeed (37), (92) and (54) show clearly that Therefore, for any t, U j (t) can be approximated, up to any power of j 2 N , by a finite Dyson expansion.
Finally, we will extend (24) as we will show that F N j has an asymptotic expansion in positive powers of 1 N whose partial sums up to any order The main results of the present note are the following.
(that is (101) truncated at order n, same slight abuse of notation). F N,n j is therefore a polynomial of order n in 1 N .
the solution of the Hartree equation (15) with initial datum F . Then, for all n ≥ 1 and N ≥ 4(eA 2n t j) 2 , Remark 1.6. If one is interested only to the expansion up to order n < j, we can change the sum in the l.h.s. of the inequality in Corollary 1.5 by a sum up to ℓ = n.
Proof. The proof is similar to the one of Corollary 2.2 in [24].
Corollary 1.7. The rate of convergence to the meanfield limit in 1 N is optimal.
, depend on N as well: first by the dependence of ∆ + j = (1 − j N )C j+1 and also by U j (t, s) defined by (42). As mentioned already the latter can be expressed as a (convergent) series in 1 N out of the linearization of the meanfield equation so that obtaining a full asymptotic expansion of E j (t) with the only knwoledge of the linearization of the meanfield equation is (tedious but) elementary. Let us note also that E j (0) is allowed to depend on N , without any restriction as soon as it satisfies (22).

The recursive construction
Let us recall from [24] that the hierarchy of error terms saisfy the following equation: Here the four operator D j , D 1 j , D −1 j , D −2 j , j = 0, . . . , N , are defined as follows (here again, J = {1, . . . , j}): for any operator G ∈ L ⊗n , n = 1, . . . , N , we denote by G(Z n ) its integral kernel, and for any function F (Z n ), n = 1, . . . , N, we define F (Z n ) as being the operator on L ⊗n of integral kernel F (Z n ) then where, by convention, In (37) We define H j (t) = K j + T j N + D j (t) and recall the definition of the (two parameters) semigroup U j (t, s) satisfying, for all s, t ∈ R, Let us perform the following rescaling We find easily that, again with the convention (41), (45) and (46) are two closed equations whose solutions are given by Iterating till j, we get explicitly the solution of (45)-(47) given by Therefore, for j = 1, . . . , N, t ∈ R, the knowledge of U j (t, s), s ≤ t ,, and E 0 j ′ (0), j ′ = 1, . . . , j guarantees the knowledge of E 0 j ′ (t), t ∈ R, j ′ ≤ j. We write this fact as Making now the ansatz E j (t) ∼ ∞ ∑ k=0 E k j (t)N −k 2 we find that the family (E k j (t)) j=1,...,N,k=0,... must satisfy with again the conventions (41) and E l j = 0 for l < 0, solved by Since E k −1 (t) = 0 by convention and E k 0 (t) = 0 for k ≥ 1 since E 0 (t) ∶= 1, we find after (51) that E 1 1 (t) and E 1 2 (t) are determined by E 1 1 (0) and E 1 2 (0). Therefore, by (53), E 1 j (t), j = 1, . . . , N are determined by (E 1 j (0)) j=1,...,N , and determine E 2 1 (t) and E 2 2 (t). These ones determine in turn all the E 2 j (t), j = 1, . . . , N and so on.
Formula (53) will give easily the following result.
has an asymptotic expansion in powers of N −1 with 3. Estimates and proof of Theorem 1.4 In order to simplify the expressions, we will first suppose that V L ∞ ̵ h = 1.
Note that one has therefore the following estimates: Let us first notice that (21) expressed on the E j s reads (22) and (23) can be rephrased as for some explicit constants A ′ , C. We get Let us define the mapping In other words, the family (U N j (t, s)) j=1,...,N solves the equation: s)), U N j (s, s) = I. Hence, the solution of (58) reads with again the same convention on negative indices.
By hypothese, R n j (0) = 0 since E n j (0) = δ n,0 E 0 j (0). Let us suppose now that D n (t) = tC n (t) and D ′ n (t) = C ′ C ′ n (t)e C t . It remains to prove an estimate like (60). We will obtain such an estimate by iterating (53). We first remark that, since e K j +T j N is unitary and D j ≤ j, the Gronwall Lemma gives that (62) U j (t, s) ≤ e j t−s .
We conclude by (64): by the same argument as in [24], Section 3, namely a rescaling of the time and the kinetic part of the Hamiltonian, (70) Therefore (60) is satisfied and Theorem 1.4 is proven.
The values of the two constants D n (t), D ′ n (t) in (61) can be expressed out of (70) by taking, by Theorem 1. [24].

Explicit computations of first orders
We have

The Kac and "soft spheres" models
In this section we consider the two following classes of mean field models (see [24] for details).
• Kac model. In this model, the N -particle system evolves according to a stochastic process. To each particle i, we associate a velocity v i ∈ R 3 . The vector V N = {v 1 , ⋯, v N } changes by means of two-body collisions at random times, with random scattering angle. The probability density F N (V N , t) evolves according to the forward Kolmogorov equation gives the outgoing velocities after a collision with scattering (unit) vector ω and incoming velocities v i , v j .
is the differential cross-section of the two-body process. The resulting mean-field kinetic equation is the homogeneous Boltzmann equation • 'Soft spheres' model. A slightly more realistic variant, taking into account the positions of particles X N = {x 1 , ⋯, x N } ∈ R 3N and relative transport, was introduced by Cercignani [9] and further investigated in [18]. The probability density F N (X N , V N , t) evolves according to the equation Here h ∶ R + → R + is a positive function with compact support. Now a pair of particles collides at a random distance with rate modulated by h. The associated mean-field kinetic equation is the Povzner equation v 1 )}, which can be seen as an h−mollification of the inhomogeneous Boltzmann equation (formally obtained when h converges to a Dirac mass at the origin). Both classes have be treated in [24] and Theorem 1.2 apply to them, in the following sense.
The underlying space L is now L 1 (R d , dv) (resp. L 1 (R 2d , dxdv))) for the Kac model (resp. soft spheres) both endowed with the L 1 norms ⋅ 1 . For for the Kac (resp. soft spheres) model. In both cases E j (t) is defined by (18), inverted by (20), and it was proven in [24] that Theorem 1.2 holds true verbatim in both cases.
Stating now the dynamics driven by (72) and (74) under the form (3) Kac (resp. soft spheres) model and V N given by the right hand-sides of (72),(74) respectively, one sees immediately that the proofs contained in Sections 2,3 remain valid after an elementary redefinition of the operators D j , D −1 j , D −2 j in (37)-(40) consisting in removing the bottom and overhead straight lines in the right hand sides and, by a slight abuse of notation, identifying functions with their evaluations. The convention (41) remains verbatim the same, together with the estimates Therefore, by Remark 3.1, the statements contained in Theorem 1.4 and Corollary 1.7 hold true, in both cases, verbatim. Moreover defining F N,n j by (35) in both cases, Corollary 1.5 reads now as follows f (x, v)dxdv = 1), and F (t) the solution of the homogeneous Boltzmann equation (73) (resp. the Povzner equation (75)) with initial datum F . Then, in both cases, for all n ≥ 1 and N ≥ 4(eA 2n Appendix A. The asbtract model A.1. The model. We will show in this section that the main results of [24] and of Section 1 of the present paper remain true in the "abstract'" mean field formalism for a dynamics of N particles that we will describe now. The present formalism contains the abstract formalism developed in [24], without requiring a space of states endowed with a multiplicative structure. States of the particle system and evolution equations. Let L be a vector space on the complex numbers. We suppose the family of (algebraic) tensor products {L ⊗n , n = 1, . . . , N } equipped with a family of norms ⋅ n satisfying assumption (A) below. the N -body dynamics will be driven on L ⊗N by a one-and two-body interaction satisfying assumption (B) and the mean field limit equation will be supposed to satisfy assumption (C). Assumptions (A) − (C) below will be followed by their incarnations in the K(ac), S(oft spheres) and Q(uantum) models.
By convention we denote L ⊗0 ∶= C, z 0 = z and we denote by L⊗ n the completion of L ⊗n with respect to the norm ⋅ n .
For the K, S and Q models, L ⊗n is L 1 (R d , dv), L 1 (R 2d , dxdv) and L 1 (L 2 (R d ), the space of trace class operators on L 2 (R d ), with their associated norms.
(A) There exists a family of subsets L⊗ n + of L⊗ n , n = 1, . . . , N , of positive elements F denoted by F > 0 stable by addition, multiplication by positive reals and tensor product and there exists a linear function Tr ∶ L → C, called trace. For every 1 ≤ k, n ≤ N and 1 ≤ i ≤ j ≤ n ≤ N , let Tr k n and σ n i,j be the two mapping defined by 3 3 The fact that the second and fourth lines of (76) define a mapping on the whole tensor space L ⊗n results easily from the definition of tensors products through the so-called universal property [19]. Indeed, let ϕn be the natural embedding L ×n → L ⊗n , (v 1 , . . . , vn) ↦ v 1 ⊗ ⋅ ⋅ ⋅ ⊗ vn, and let h be any mapping L ×n → L ×n ′ , then the universal property of tensor products says that there is a unique maph ∶ L ⊗n → L ⊗n ′ such thath ○ ϕn = ϕ n ′ ○ h. Taking  give the desired extensions.
We will suppose that Tr k N and σ n i,j , i, j, k ≤ n ≤ N , satisfy, for any F ∈ L ⊗n , Tr k N (F ), σ n i,j (F ) > 0, Tr k n (F ) n−1 = F n when F > 0 σ n i,j (F ) n = F n Tr k n (F ) n−1 ≤ F n In particular one has that F n = Tr n . . . Tr 1 F when F > 0 and Tr n . . . Tr 1 F ≤ F n in general. Note that (77) allows to extend Tr k n and σ n i,j to L⊗ n by continuity. We will use the same notation for these extensions.
For the K, S and Q models, Tr k is ∫ R d ⋅dv k , ∫ R 2d ⋅dx k dv k as indicated in Section 5, and the partial traces defined in Section 1.1. The action of σ n i,j consists obviously in exchanging the variables v i and v j , (x i , v i ) and (x j , v j ) and (x i , x ′ i ) and (x j , x ′ j ), (in the integral kernel), respectively. Finally (77) is satisfied in the three cases.
From now on and when no confusion is possible, we will identify L ⊗n with its completion L⊗ and we will denote Tr k N = Tr k (note also that Tr = Tr 1 1 = Tr 1 ), σ N i,j = σ i,j and Tr(= Tr n ) = Tr n n Tr n−1 n . . . Tr 1 n . Moreover, with a slight abuse of notation, we will denote We call symmetric any element of L ⊗n invariant by the action of σ n i,j , i, j ≤ n. We call state of the N −particle system an element of (79) D N = {F ∈ L ⊗n F > 0, F = 1 and F is symmetric}.
For j = 0, . . . , N , the j-particle marginal of F N ∈ (L ⊗N ) + 1 is defined as the the partial trace of order N − j of F N , that is Note that F N j ∈ L ⊗j (F N 0 = 1 ∈ L ⊗0 ∶= C) and F N j > 0, F N j j = F N N since Tr is positivity and norm preserving, and obviously F N j is symmetric as F N . That is to say: (B) The evolution of a state F N in L ⊗N is supposed to be given by the N −particle dynamics associated to a two-body interaction: where the operators on the right hand side are constructed as follows.
for a (possibly unbounded) operator K acting on L and a bounded two-body (potential) operator V acting on L ⊗2 . We assume furthermore that K is the generator of a strongly continuous, isometric, positivity preserving semigroup (in L) and K N + V N is the generator of a strongly continuous, isometric, positivity preserving semigroup (in L ⊗N ) Finally, for any F ∈ L, F N ∈ L ⊗N and i, r > j, we assume (86) Tr(KF ) = 0 and Tr j, This last property is necessary to deduce the forthcoming hierarchy. Hierarchies. The family of j-marginals, j = 1, . . . , N , are solutions of the BBGKY hierarchy of equations where: , Indeed, thanks to (86) we get easily by applying Tr j,N on (81) that By symmetry of F N and V i,k we get Tr j,N (V i,k F N ) = Tr j+1 (V i,j=1 F N j+1 ) for all k > j and (87) follows.
We introduce the non-linear mapping Q(F, F ), Q ∶ L × L → L by the formula Eq. (94) is the Boltzmann, Povzner or Hartree equation according to the specifications established in the table above. In full generality we will assume (C) (94) has for all time a unique solution F (t) > 0 and F (t) = 1.
For the K, S and Q models, (C) is true by standard perturbations methods.
Correlation error. To introduce the correlation errors, we need to extend slightly the above structure.
For any subset J ⊂ {1, . . . , N } we first define where χ J is the characteristic function of J and L ⊗0 = C.
Then we introduce L ⊗J , the subspace of L ⊗J N formed by vectors of the form We define a norm on L ⊗J by For F ∈ L and K ⊂ J ⊂ {1, . . . , N } we introduce the linear operator [F ] ⊗K J , defined through its action on factorized elements as where Note that, for K, K ′ ⊂ J, K ∩ K ′ = ∅, we have the composition where F solves (94), the operator [F ] ⊗K J is defined by (96) and F N L ∈ L ⊗L is defined through its decomposition on factorized states. Namely if The link between the definition of F N L and the definition of the marginals F N j given in (80) is the following: The formula inverse to (99) reads Note that the contribution in the right hand side of (101) corresponding to K = J and K = ∅ are F ⊗ J and E J respectively. To prove (101), we plug (99) in the r.h.s. of (101) and we use (97): One notices that since F N j is the marginal of some F N which decomposes on elements of the form v 1 ⊗ ⋯ ⊗ v N , F N j decomposes on elements of the form ( Trv k )v 1 ⊗ ⋅ ⋅ ⋅ ⊗ v j . Since one knows that F N j is symmetric, it is enough to choose one bijection i J ∶ {1, . . . , j} → J, J = j, and consider the mapping where a s = 1 if i ∉ J and a i J (j) = v j . Φ i J is obviously one-to-one since i J is so, and, though (102) depends on the embedding chosen, (103) does not: Φ i J restricted to the space L ⊗ J S of symmetric-by-permutation elements of L ⊗ J , depends only on J and not on i J . We will call Φ J this restriction, The same argument is also valid for E J which enjoys the same symmetry property than F N J and we define (105) E J = Φ −1 J E J . Φ J is obviously isometric and we have that Therefore, considering the one-to-one correspondence Φ J , it is enough to compute/estimate the quantities E j , j = 1, . . . , N . E j and F N j are linked by For the K, S and Q models, the corrsponding expression are given in Sections 5 and 1.1.
We get the following result. To compute T 2 we make use of the inverse definition (101): Distinguishing among the belonging or not to K ′ of i and j + 1 in the r.h.s. of (113), we decompose (114) We have = α(j, N ) i∈J C i,j+1 [F ] ⊗{i,j+1} D j ∶ L ⊗j → L ⊗j , j = 1, . . . , N, where, by convention, . Note that one has the following estimates: