A New Method and a New Scaling for Deriving Fermionic Mean-Field Dynamics

We introduce a new method for deriving the time-dependent Hartree or Hartree-Fock equations as an effective mean-field dynamics from the microscopic Schrödinger equation for fermionic many-particle systems in quantum mechanics. The method is an adaption of the method used in Pickl (Lett. Math. Phys. 97 (2) 151–164 2011) for bosonic systems to fermionic systems. It is based on a Gronwall type estimate for a suitable measure of distance between the microscopic solution and an antisymmetrized product state. We use this method to treat a new mean-field limit for fermions with long-range interactions in a large volume. Some of our results hold for singular attractive or repulsive interactions. We can also treat Coulomb interaction assuming either a mild singularity cutoff or certain regularity conditions on the solutions to the Hartree(-Fock) equations. In the considered limit, the kinetic and interaction energy are of the same order, while the average force is subleading. For some interactions, we prove that the Hartree(-Fock) dynamics is a more accurate approximation than a simpler dynamics that one would expect from the subleading force. With our method we also treat the mean-field limit coupled to a semiclassical limit, which was discussed in the literature before, and we recover some of the previous results. All results hold for initial data close (but not necessarily equal) to antisymmetrized product states and we always provide explicit rates of convergence.


Introduction
The behavior of an interacting many-body system in classical or quantum mechanics can be very complicated and the microscopic equations governing its behavior are usually practically impossible to solve for more than three or four particles. For large interacting systems it is therefore essential to use a statistical description in order to make statements about the typical behavior of a large number of particles. One type of such a statistical description is to approximate the microscopic dynamics by an effective one-body dynamics, i.e., replace the microscopic evolution equation with very many degrees of freedom by a simpler, usually non-linear equation with very few degrees of freedom. Famous examples are the Boltzmann, Navier-Stokes and Vlasov equations in classical mechanics, and the Hartree, Hartree-Fock and Gross-Pitaevskii equations in quantum mechanics; different regimes can lead to different effective evolution equations, e.g., kinetic equations or mean-field equations. Apart from the study of these equations and their interesting properties and consequences, it is an ongoing project of mathematical physics to derive them from the microscopic dynamics, i.e., to prove rigorously that the solutions to the effective equations approximate the solutions to the microscopic equations well in certain situations. Only some cases of such a rigorous derivation are known. In classical mechanics, this could for example be shown for the Vlasov equation [17], however, for the Boltzmann equation it has been shown only for very short times [40], and for the Navier-Stokes equation it is still an open problem (see [58] for an excellent overview). In quantum mechanics, starting in the 70s and by a series of recent works, the derivation of effective dynamics for bosons near a condensate is well understood, see [28,31,37,39,51,54,56] for the case of the Hartree equation and [10, 23-27, 50, 52] for the Gross-Pitaevskii equation. Derivations of mean-field dynamics for fermions have been discussed in [8,13,14,22,30,47,57]; in particular [14] treats an interesting scaling limit very comprehensively. In general, the Fermi statistics makes the situation more complicated, see the discussion in Sections 1.1 and 1.2.
For fermionic many-particle systems in quantum mechanics, the microscopic evolution is given by the Schrödinger equation (we set = 1 = 2m throughout this work) i∂ t ψ t = H ψ t (1) for antisymmetric complex-valued N-particle wave functions ψ t ∈ L 2 (R 3N ) (for simplicity, we neglect spin and always work in three dimensions). Antisymmetry means that ψ t (. . . , x j , . . . , x k , . . .) = −ψ t (. . . , x k , . . . , x j , . . .)∀j = k. We consider Hamiltonians where H 0 j acts only on x j and v (N) is a real-valued pair-interaction potential (the superscript (N) denotes a possible scaling with the particle number N which is explained below). According to (1), the unitary time evolution of an initial wave function ψ 0 is given by ψ t = e −iH t ψ 0 if H is self-adjoint which we henceforth assume. Note that for antisymmetric initial conditions ψ 0 , the wave function ψ t remains antisymmetric under the Schrödinger evolution (1) with Hamiltonian (2) for all times. For the desired mean-field description, consider N orthonormal one-particle wave functions (also called orbitals) ϕ t 1 , . . . , ϕ t N ∈ L 2 (R 3 ) which are solutions to the fermionic Hartree equations (sometimes called reduced Hartree-Fock equations). These are the coupled system of non-linear differential equations for j = 1, . . . , N, where * denotes convolution, and ρ t N = N i=1 |ϕ t i | 2 is the spatial density. Note that for orthonormal initial conditions ϕ 0 1 , . . . , ϕ 0 N , (3) preserves the orthonormality for all times. The term is called the mean-field. It can be viewed as the average value of the interaction potential at point x, created by particles distributed according to the density ρ t N . Note that closely related effective equations for fermions are the Hartree-Fock equations, where an additional exchange term is present on the right-hand side of (3). In general, the Hartree-Fock equations are expected to be a better approximation than the fermionic Hartree equations. However, the exchange term is always smaller than the direct term (4), and, as we see later, is negligibly small for the scaling limits that we consider in the sense that including the exchange term does not improve our error estimates. It would be interesting to study other scaling limits where taking into account the exchange term gives a more accurate approximation. Now suppose that some initial ϕ 0 1 , . . . , ϕ 0 N are given. Let the initial N-particle wave function be ψ 0 is the antisymmetrized product of ϕ 1 , . . . , ϕ N , with S N the symmetric group and (−1) σ the sign of the permutation σ . Then, under the Schrödinger evolution (1), this initial wave function evolves to ψ t = e −iH t ψ 0 . We want to compare this ψ t to the wave function N j =1 ϕ t j , where the ϕ t j are the solutions to the fermionic Hartree equations (3). In other words, if still ψ t ≈ N j =1 ϕ t j at some time t, then the Schrödinger dynamics is approximated well by the Hartree dynamics. To show such a statement is the goal of this article.
Note that in the presence of an interaction potential v (N) it is in general never true that e −iH t N j =1 ϕ 0 j = N j =1 ϕ t j , since the interaction leads to correlations between the particles, i.e., the wave function evolves into a superposition of many antisymmetrized product states. These correlations are caused by deviations from the mean-field behavior, i.e., by fluctuations around the mean-field. For a rigorous derivation we thus need to consider a regime where the fluctuations become small for large N, on the relevant time scales. Such an appropriate limit N → ∞ is then called a mean-field limit.

Density ∝ N Regime
One such regime describes a system of fermions with very high density proportional to N. The microscopic wave function ψ t (x 1 , . . . , x N ) is a solution to the Schrödinger equation (here for non-relativistic particles) and ϕ t 1 , . . . , ϕ t N are solutions to the corresponding fermionic Hartree equations for j = 1, . . . , N (recall ρ t N = N i=1 |ϕ t i | 2 ), where w (N) is a real-valued external field. The evolution (7) is considered for initial data in a volume of order 1 1 which can be realized if the particles are confined by a nice external trapping potential w (N) . Note, that due to the antisymmetry of ψ t , the total kinetic energy for particles in a volume of order 1 is always bigger or equal CN 5/3 (which can, e.g., be read off from the kinetic energy inequality in Section 8.1). For states close to the ground state and nice external fields, the kinetic energy is of order N 5/3 (e.g., for free particles in a box of constant side length, this can easily be checked). Then, the kinetic term on the right-hand side of (7) is of order N −2/3 N 5/3 = N. The interaction term is of the same order, since N 2 terms are scaled down with N −1 .
Let us describe heuristically two ways of how to arrive at the scaling limit of (7) (and let us for simplicity disregard external fields here). One way is to start with initial conditions in a volume of order 1 and kinetic energy of order N 5/3 and then add a scaling for the interaction, such that the scaled interaction and the kinetic energy are of the same order. The Schrödinger equation is then Since the average momentum per particle is of order N 1/3 , but the average scaled interaction or force per particle is of order N −1/3 N = N 2/3 , one can, for large N, expect a nice limiting equation only for short timest of order N −1/3 (change in momentum = force × time). By introducing t = N 1/3t (and then dividing both sides of the equation by N 2/3 ), we get (7). Let us describe another, more physical way to arrive at (7), which, however, only works for Coulomb interaction. Consider initial data in a small volume of order N −1 and kinetic energy of order N 7/3 . This is realized, e.g., for electrons in the "core" region of a large-Z atom, see [41]. The unscaled Schrödinger equation for Coulomb interacting particles is A calculation for the free ground state in a box shows that the average force per particle is of order N 5/3 . Since the average momentum is of order N 2/3 , from "change in momentum = force × time" we can expect nice behavior for large N for times of order N −1 . If we now introduce rescaled variables t = Nt and x = N 1/3x (such that in the variables t, x we consider the same kind of initial conditions as in (7)), we arrive at (7).
Note that (7) has a semiclassical structure due to the N −1/3 factors, which play the role of a very small parameter, like the in the Schrödinger equation with units. As a consequence, the solutions to the Schrödinger equation (7) are close to solutions to the classical Vlasov equation, e.g., in the sense that the Wigner transform of a solution to (7) is close to a classical phase space density ρ t (x, p) that solves the Vlasov equation. Such a derivation of the Vlasov equation from the microscopic Schrödinger equation has been considered in [47] and improved in [57]. Several other works [1,2,4,12,33,45,46,49] also study the derivation of the Vlasov equation starting from the Hartree equation (8) or the Hartree-Fock equation. However, the mean-field dynamics (8) is a better approximation to the dynamics (7) than the Vlasov dynamics. A derivation of the mean-field dynamics (8) from the microscopic dynamics (7) has first been given in [22] for bounded analytic interactions and short times. In [13,14], the derivation was crucially improved in the sense that it was shown to hold for all times, for initial data with a certain semiclassical structure. Furthermore, the result holds for pseudo-relativistic free Hamiltonians and with fewer regularity assumptions on the interaction. The interaction is, however, still assumed to be bounded. Since many applications concern electrons and in light of the discussion around (10), a derivation for Coulomb interaction would be desirable. Let us also mention [11] where the analysis of [14] is extended to fermionic mixed states.

Density ∝ 1 Regime
Let us now introduce a new scaling limit for fermions. The inspiration for this new limit is twofold. On the one hand, we want to study a scaling limit where the solution to the Schrödinger equation is not close to a solution to the Vlasov equation. In such a limit, one would see more quantum effects like interference of wave packets. On the other hand, we are inspired by the application of the Hartree-Fock equations to large molecules, where the size of the system is of order N. This is the case for the scaling limit we now introduce.
We consider initial conditions in a volume of order N, i.e., with average density of order 1. The idea is that for such initial data a mean-field approximation is valid for long-range interactions like Coulomb interaction. The microscopic evolution for the wave function ψ t (x 1 , . . . , x N ) is given by the Schrödinger equation (for nonrelativistic particles) and the corresponding fermionic Hartree equations are a real-valued long-range interaction. The Schrödinger equation (11) is interesting for initial data with total kinetic energy of order N. (Note that by the kinetic energy inequality (see Section 8.1), this implies that the system volume is at least of order N.) Now the crucial observation is that the total unscaled interaction energy for longrange v is not of order N 2 , as one might expect from the order N 2 terms in the double sum i<j in (11). It is in fact smaller, namely of order N 1+β , with some 0 ≤β ≤ 1 depending on the long-range behavior of v. This can be seen heuristically by considering the mean-field for constant density and interactions |x| −s , since (V N denotes a volume of order N) for appropriate s > 0 (see Lemma 8.1). Thus, if we set β =β in (11), the kinetic term and the interaction term are both of order N.
Let us now consider the average force per particle. For simplicity, let us discuss it heuristically in the mean-field approximation (a similar consideration can be done for the microscopic wave function ψ t ). There, the average force, including the scaling withβ = 1 − s/3, is given by the gradient of the mean-field, i.e., it is of the order Thus, the mean-field is almost constant on scales of order 1 and only varies over the whole system size. One would therefore heuristically expect free evolution to leading order. This classically inspired heuristics is actually only almost right, since orbitals which are spread over the whole system should feel an effect coming from the mean-field. We expect, however, that this could be easily taken care of by adding a time and space dependent phase to the freely evolving orbitals. Those orbitals can be expected to approximate the Schrödinger equation (11). However, similar to before, one would expect that the fermionic Hartree equations (12) provide a more accurate approximation to the dynamics. In fact, we prove this in Section 2.1 for 0 < s < 3/5, where we also discuss relevant observables.
Note that one can repeat the heuristic calculation (14) and ask for the variation of the force. One finds, e.g., that for s < 1, and (16) for Coulomb interaction (where one would actually expect an ln N correction to the variation of the force). Thus, while the leading order in this scaling limit is free dynamics, the next to leading order is a constant force (at least when taking appropriate phase factors into account, see the discussion after (14)). For interactions |x| −s , we can compare the new scaling limit to the one from Section 1.1. Let us choose s = 1 here. The new scaling limit for initial conditionsψ 0 in a volume of order N is then If we rescale the spatial variables such that the initial conditions are in a volume of order 1, i.e., we introduce x = N −1/3x , then the rescaled wave function ψt solves This is similar to (7), but on a different time scale. The time scales are related bỹ t = N 1/3 t. In other words, if we formulate the new scaling limit for initial conditions in a volume of order 1, then the time scales are much shorter than those considered in Section 1.1. The important difference between the two scaling limits lies in the kind of initial conditions one would naturally consider. In (7), the initial conditions would vary over spatial scales of order 1, i.e., the whole system size, and not over scales of order, say N −1/3 . This is also assumed in [14,22]. On the other hand, the initial conditions for (11) would naturally vary over scales of order 1. If we would rescale those initial conditions to a volume of size 1, then they would have a small scale structure on spatial scales of order N −1/3 . In fact, in our main results about the new scaling limit, we have no restrictions for the initial conditions, except that the kinetic energy is appropriately bounded by CN.
The described regime, to our knowledge, has not been considered in the literature before for a derivation of mean-field dynamics. Note, however, that [7][8][9]30] consider the case β = 1 which, for Coulomb interaction, leads to an interaction term of order N −1/3 . Compared to that, the results in this article are an improvement by a factor N 1/3 in the interaction strength. However, in [30] Coulomb interaction without any cutoff is considered, while we have to introduce a very mild cutoff on scales much shorter than the average particle distance, in order to treat the Coulomb singularity. Note that in [6], by using the method introduced in this article and additional techniques like the Fefferman-de la Llave decomposition of the Coulomb potential, also Coulomb interaction without cutoff can be dealt with.

Outline
Next, we present the main results of this article. The results for the new scaling (11), which are the focus of this work, are presented and discussed in Section 2.1. In Section 2.2, we derive the mean-field equations (8) for a class of bounded interactions and initial states with a semiclassical structure. For this case, we reproduce some of the results from [14], with minor improvements on the initial conditions. We give this derivation in order to demonstrate the generality of our approach and since the proof is very short. In Section 3, we introduce our method for deriving fermionic mean-field dynamics. In Section 3.1, we explain the connection between the functional we use in the proof and reduced density matrices. In Section 3.2, we state a slightly more general version of our main result in terms of the newly defined functional and in Section 3.3 we sketch the proof of these theorems. All proofs are given in Sections 4 to 9. Finally, let us remark that the length of this article is due to the fact that we introduce our method in great detail and generality, in particular in Sections 4 to 7.2. The core parts of the proof are given in Sections 7.3 and 9.

Main Results
In order to state the main results of this article, we need a measure for the closeness of a many-body wave function ψ ∈ L 2 (R 3N ) to the antisymmetrized product state N j =1 ϕ j , with ϕ 1 , . . . , ϕ N ∈ L 2 (R 3 ). A natural choice is a trace norm estimate on the reduced density matrices. For any normalized antisymmetric ψ ∈ L 2 (R 3N ), the reduced one-particle density matrix is defined by its integral kernel If we want to have control over the statistics of one-body observables A : , we need to control the trace norm difference of the reduced one-particle density matrices of ψ and ϕ j , since, e.g., for bounded A, where ||·|| op denotes the operator norm and ||·|| tr the trace norm (see also Section 5). Thus, we express our main results in terms of the trace norm difference above. Note that we actually prove slightly stronger convergence statements, see Section 3.

Density ∝ 1 Regime with Interactions |x| −s
In this section we explicitly consider the non-relativistic Schrödinger equation (11) and the corresponding fermionic Hartree equations (12), as discussed in Section 1.2.
The results in this subsection are concerned with interaction potentials |x| −s with 0 < s < 6/5, sometimes with singularity weakened or cutoff, and the corresponding β = 1 − s/3 (see also Lemma 8.1). For the following results we assume the existence of solutions to the Eqs. (11) and (12) and that the total kinetic energy of the solutions to (12) (but not necessarily to (11)) is bounded by AN, i.e., N i=1 ||∇ϕ t i || 2 ≤ AN. There are several works about solution theory to the Hartree(-Fock) equations [15,16,18,19] which establish existence and uniqueness of solutions even for singular (attractive or repulsive) interactions and external fields like |x| −1 . A blowup of solutions is only expected for strong attractive interactions (e.g., for gravitating fermions) with semirelativistic free Hamiltonian, see, e.g., [32,35,36]. (Indeed, as is shown in [32], even for bounded interactions, v * ρ t N can become infinite in the limit of t → T , for some T < ∞.) Therefore, for non-relativistic Hamiltonians, due to the conservation of the total Hartree energy, N i=1 ||∇ϕ t i || 2 ≤ AN always holds if it holds for the initial states ϕ 0 1 , . . . , ϕ 0 N and if the external field is nice enough (e.g., scaled external Coulomb fields generated by nuclei with some N-independent distances to each other are ok).
Proof See Section 8.2.

Remark
1. Let us put the three different cases of Theorem 2.1 into perspective. In (a) we treat singular interactions, but with weaker singularity than Coulomb. Case (b) includes Coulomb interaction with the singularity cut off on a ball with radius much smaller than the average particle distance. However, in case (2.1) the convergence rate is only N −γ /2 , with γ < 1/3 in the Coulomb case, depending on the cutoff distance N −δ . In case (c) we treat Coulomb interaction without cutoff. We get the convergence rate N −1/2 . For this we need additional regularity assumptions on the solutions to the fermionic Hartree equations. We expect these conditions to hold for all times under suitable assumptions on the initial data, but we refrain from proving that in this paper. 2. Using Lemmas 3.2 and 3.3 we can also deduce bounds in Hilbert-Schmidt norm.
For example, in case (a) we have (Note that our choice of normalization is γ ϕ t j 1 tr = 1 and If we set the external field w (N) = 0, then, for the interactions in this theorem, the bound (21) holds for all times, if it holds for the initial conditions. Furthermore, we have existence and uniqueness for the solutions to (11) and (12) in this case. Thus, for w (N) = 0, all the results in this theorem hold under the sole assumption that (21) holds at t = 0. 4. Note that in the setup of Theorem 2.1, the kinetic energy per particle is of order 1, i.e., a particle can on average travel a distance of order 1, while the diameter of the system is of order N 1/3 . Heuristically speaking this means that the system looks static on very large scales. The results above are therefore relevant for observables that are sensitive for properties on small scales. Let us give a simple example of such an observable. We divide the volume of the system into 2N cells, where each cell has a volume of order 1. As initial state we choose smooth wave packets with disjoint supports, each with kinetic energy of order 1. Each packet occupies one of the cells, and we leave a neighboring cell unoccupied. Then a suitable observable A would be the number of particles in the unoccupied cells. After a short time of order 1, due to the diffusion of wave packets, there will be order N particles in the previously unoccupied cells. Furthermore, the wave packets will show interference effects. 5. Let us follow up on the example of the previous remark and discuss what result we would expect if we compare the Schrödinger evolution with the free evolution. As discussed in Section 1.2, the average force per particle is of order N −1/3 . This would be the expected relative error we make in estimating the expectation value of the observable A from above: after some time t the unoccupied cells contain order N ± N 2/3 particles. On the other hand, the estimates (23) and (27) would give a prediction of order N ±N 1/2 . Thus, in the cases (a) and (c), the prediction from the fermionic Hartree equations is better than what we could expect from the free evolution (or the free evolution with a time and space dependent phase).
Let us conclude this subsection with a remark about Coulomb interaction with N −1 scaling, and with N −2/3 scaling for times smaller than order 1. In the first case, we have free evolution for times of order N 1/3−ε , for any ε > 0. In the second case, we even have closeness to any state, without evolution, if, of course, the initial conditions are close. We summarize this in the following proposition which for simplicity we only state for Coulomb interaction without external fields.

Proposition 2.2 Let ψ t ∈ L 2 (R 3N ) be a solution to the Schrödinger equation
with antisymmetric initial condition Proof See Section 8.2.

Density ∝ N Regime
Let us now state the result for the derivation of the mean-field dynamics for the regime discussed in Section 1.1. Such a derivation has recently been given in [14] and here, we reproduce some of these results. We actually use estimates about the propagation of properties of the initial data from [14] (see Lemma 9.1). A slight improvement of our result is that our conditions on the closeness of the initial data are more transparent and general, see Remark 10. In contrast to [14], we express our result solely in terms of the one-particle reduced density matrices. We consider the Schrödinger equation (7) (for simplicity, without external fields) and the corresponding fermionic Hartree equations (8). Note that for the interactions we consider the exchange term can easily be shown to be subleading (compared to the error term in the bound (36)), so the theorem can be proven directly for both the fermionic Hartree and for the Hartree-Fock equations . Note that here we do not have to use the long-range behavior of the interaction, since we are interested in solutions in a volume of order 1. For technical reasons, we consider a class of interactions that is in particular bounded. Our theorem is analogous to [14,Thm. 2.1]. We denote the trace norm by ||·|| tr (see also Section 5) and the commutator by [A, B] = AB − BA.
be solutions to the fermionic Hartree equations (8) or to the Hartree-Fock equations (32), with orthonormal initial conditions ϕ 0 wherev is the Fourier transform of v, and that the initial conditions ϕ 0 1 , . . . , ϕ 0 N are such, that sup Proof See Section 9.

Remark
6. The condition (34) captures the semiclassical structure of the initial data, see [14] for a more detailed discussion. The condition (35) is used to propagate (34) in time. Note that for the initial state from Remark 4 the condition (35) does not hold. 7. The theorem also holds with external fields that are such that they preserve the bounds (34) and (35) for all t. 8. The double exponential in the bound (36) comes from using the Gronwall Lemma twice: once for propagating the conditions (34) and (35) in time, and once for the time derivative of a quantity similar to the left-hand side of (36) (see Section 3). 9. Using Lemma 3.2 we can also deduce a bound in Hilbert-Schmidt norm, i.e., ]. In particular, the convergence rates are the same. The bounds (36) and (37) are formulated for more general initial data, in the sense that convergence of the initial trace respectively ( √ N times) Hilbert-Schmidt norm is sufficient to control the left-hand sides of (36) and (37). In [14,Thm. 2.1], there are two additional bounds that improve the rate of convergence: in [14,Eq. (2.21)] this is achieved by additional assumptions on the closeness of the initial data to antisymmetric product states and in [14,Eq. (2.22)] by evaluating the trace norm difference only for certain observables. Furthermore, in [14,Thm. 2.2] bounds for the reduced k-particle density matrices are proven.

Exposition of the Method
We now present our method for deriving fermionic mean-field dynamics from the microscopic Schrödinger dynamics. This method does not rely on BBGKY hierarchies but on a Gronwall type estimate for a suitably defined functional, which leads to many technical advantages. The key of the method is to define the functional α f (ψ, ϕ 1 , . . . , ϕ N ) that measures the closeness of a many-body wave function ψ ∈ . For the definition of α f we need to introduce several projectors that play a crucial role in the proofs later.
i.e., its action on any ψ ∈ L 2 (R 3N ) is given by We define p m := N j =1 p ϕ j m and q m := 1 − p m . For any 0 ≤ k ≤ N we define with the set A k := a = (a 1 , . . . , a N ) ∈ {0, 1} N : N m=1 a m = k , i.e., P (N,k) is the symmetrized tensor product of q 1 , . . . , q k , p k+1 , . . . , p N . We define P (N,k) = 0 for all k < 0 and k > N. We call any f : For any weight function f we define the operator and the functional where ·, · denotes the scalar product on L 2 (R 3N ).
The functional α f has first been introduced by one of the authors (P.P.) for bosons [51], that is, with p m = |ϕ ϕ| m . The functional was used in [39,51] for the derivation of the bosonic Hartree equation, and in [50,52] for the derivation of the Gross-Pitaevskii equation. A variant of it has been used in [20,21] to derive the effective dynamics of a tracer particle in a Bose gas and in [3] to derive the Hartree-von Neumann limit.
Let us explain Definition 3.1 a little further. First, note that p m , q m and P (N,k) are indeed projectors, since ϕ 1 , . . . , ϕ N are assumed to be orthonormal. As operators on L 2 (R 3 ), p 1 projects on the subspace spanned by ϕ 1 , . . . , ϕ N , and q 1 on its complement, i.e., in particular, where the last equation is true since P (N,k) contains exactly k q-projectors in each summand and q m p m = 0. We can thus define the decomposition ψ = N k=0 P (N,k) ψ of the microscopic wave function, where now each contribution P (N,k) ψ has (N − k) particles in one of the orbitals ϕ 1 , . . . , ϕ N ("good particles") and k particles not in the orbitals ϕ 1 , . . . , ϕ N ("bad particles"). With the functional α f , each contribution P (N,k) ψ is given the weight f (k), i.e., f specifies how much weight is given to the number of particles outside the antisymmetrized product state N j =1 ϕ j . Since 0 ≤ α f ≤ 1, α f ≈ 0 means (for suitable weight functions f ) that the approximation of ψ by N j =1 ϕ j is very good, while α f ≈ 1 means that the approximation is not valid at all. By choosing an appropriate f we can fine tune what exactly is meant by the closeness of ψ to N j =1 ϕ j . One obvious and very simple weight is the relative number k/N. We always denote this weight function by n(k) = k/N (44) and the corresponding functional by α n . For this weight function, due to the antisymmetry of ψ and (43), the functional has the simple form (recall that each summand in P (N,k) contains exactly k projectors q) Let us note here that α n has been used before in [5,34] to measure deviation from the antisymmetrized product structure in the static setting; see also the remarks following (51). Furthermore, α n coincides with the measure introduced in [14] for states in Fock space with fixed particle number N. Another important weight is with some 0 < γ ≤ 1. The function m (γ ) (k) gives a much larger weight to already very few particles outside the antisymmetrized product structure. On the other hand, for k > N γ , i.e., very many particles outside the antisymmetrized product structure, m (γ ) (k) gives the same weight 1 for all k > N γ . These properties enable us to derive mean-field approximations for a much wider range of physical situations, e.g., singular or weakly scaled interaction potentials. Let us stress that the freedom in the choice of the weight function is a key feature and advantage of the described method.
Note that we could also define projectors p ϕ = N m=1 p ϕ m and q ϕ = 1 − p ϕ . For antisymmetric ψ, the α f functional defined with those projectors coincides with as can be seen from multiplying out the left-hand side. Note that the projectors p ϕ are related to fermionic creation and annihilation operators a(ϕ) and a † (ϕ) (see, e.g., [43]) by p ϕ = a † (ϕ)a(ϕ), i.e., p ϕ is the number operator for the state ϕ.
The goal of this work is to prove bounds on where ψ t is a solution to the Schrödinger equation and ϕ t 1 , . . . , ϕ t N are solutions to the fermionic Hartree equations. In more detail, we first look for a bound on the time derivative of α f (t) of the type which then, by Gronwall's Lemma, implies the bound where the function C(t) is independent of N, and δ > 0 is called the convergence rate. In the main theorems of Section 2, the weight function is either n from (44) or m (γ ) from (46). A bound of the form (48) implies that if initially (at time t = 0) α f is small, then it stays small for times t > 0 and N large enough. In the limit of N → ∞, we arrive at the statement that lim N→∞ α f (t = 0) = 0 implies lim N→∞ α f (t) = 0 for all t > 0. Note that for all f ≥ n, we have α f ≥ α n , i.e., if the bound (47) holds for such an f , we find Thus, by using a weight function f ≥ n (which can be advantageous in the estimates, see below) we can still control α n (t), but now under stronger conditions on the initial state, namely that α f (0) has to be small in N.

Connection to Density Matrices
The functional α f is closely related to the trace and Hilbert-Schmidt norms of the difference between reduced one-particle density matrices as defined in (19). Note that for an antisymmetrized product state N j =1 ϕ j we find Let us now consider α n , i.e., the α-functional with the weight n(k) = k/N. Note that (51) where tr(·) denotes the trace. It is the expression tr(γ ψ 1 q 1 ) that has been used before to measure deviation from the antisymmetrized product structure in the static setting, see [5,34]. The relations between α n and the differences between Schrödinger and Hartree reduced one-particle density matrices in trace norm ||·|| tr and Hilbert-Schmidt norm ||·|| HS are summarized in the following lemma (see Section 5 for the definition of the norms).
) be antisymmetric and normalized, and let Proof See Section 5.
Note that the extra N factors are due to the choice of normalization; indeed γ ϕ j 1 tr = 1 and γ ϕ j 1 HS Lemma 3.2 is the main result of this section. It implies in particular that (55) i.e., convergence of γ ψ 1 to γ ϕ j 1 in certain norms is equivalent to convergence of α n to zero. However, there is a difference in the convergence rates, e.g., if α n converges with rate N −1 , then the density matrices converge in trace norm only with rate N −1/2 , as can be seen from (52).
Let us now consider more general weight functions f with f (k) ≥ n(k) for all k and in particular the weight m (γ ) (k) from (46) which we use later. Note that for those weights α f ≥ α n . Lemma 3.3 Let ψ ∈ L 2 (R 3N ) be antisymmetric and normalized, and let Thus, convergence of α f to zero still implies convergence of γ ψ 1 to γ ϕ j 1 , but in general not the other way around. We use (57) to express some of our main result only in terms of trace norms.
Let us make a brief remark about convergence in operator norm. Note that ||γ ψ 1 || op ≤ N −1 for antisymmetric ψ, so a possible indicator of convergence would be the operator norm times N. This, however, is not a good type of convergence for our purpose, since the operator norm is given by the largest eigenvalue (which at most can be N −1 for fermionic density matrices; recall our choice of normalization). Thus, while convergence of N times the operator norm does imply convergence of α n , the opposite is not true. One orbital not in the antisymmetrized product of the ϕ 1 , . . . , ϕ N is enough to let the operator norm of N times the difference between the density matrices be equal to one, while α n converges to zero.

Main Theorems in Terms of α
Let us state and discuss two theorems with bounds on α n and α m (γ ) which hold for very general Hamiltonians. From these theorems we can then deduce Theorem 2.1.
Theorems 3.4 and 3.6 are of the form: Given certain properties of the solutions to the fermionic Hartree equations, the mean-field approximation for the dynamics is good, i.e., α f (t) ≤ C(t) α f (0) + N −δ for some δ > 0. We consider wave functions ψ t ∈ L 2 (R 3N ) that are solutions to the Schrödinger equation (1) with self-adjoint Hamiltonian (2) with real and possibly scaled interaction v (N) The most important example for the free part H 0 j is the non-relativistic free Hamiltonian with external field, H 0 j = − j + w (N) (x j ), but we could also replace the Laplacian by relativistic operators like √ − + m 2 − m (m > 0) or |∇|. The fermionic meanfield equations for the one-particle wave functions ϕ t 1 , . . . , ϕ t N ∈ L 2 (R 3 ) are given by (3). The first theorem gives a bound on α n (t) = ψ t , q t 1 ψ t .

11
. The corresponding bounds for the reduced one-particle density matrices follow from Lemma 3.2. 12. The condition (58) is only a condition on the solutions to the fermionic Hartree equations (3), and not on the solutions to the Schrödinger equation (1). Note that for solutions ϕ t 1 , . . . , ϕ t N ∈ H 1 (R 3 ) (the first Sobolev space) and interactions with at most a singularity |x| −1 , the left-hand side of (58) is always finite due to Hardy's inequality. Therefore, the challenge lies in the N-dependence on the right-hand side of (58). 13. Note that condition (58) implies that (by Cauchy-Schwarz and ρ t N = N) i.e., the scaled mean-field interaction is everywhere bounded. 14. There is an interesting connection between condition (58) and the fluctuations around the mean-field. At each point y ∈ R 3 , let us denote the "fluctuations" in the interaction at y in the state The condition (58) can be relaxed and replaced by other conditions if we use the weight function m (γ ) (k) from (46). This allows to treat more singular interactions and smaller scaling exponents. Let us first summarize the precise assumptions that we need on the scaled interaction v (N) and the density ρ t N and afterwards state the theorem.

Assumption 3.6 For all
sup sup Under this assumption we can conclude convergence of α m (γ ) (t).
Proof See Section 7.4.

Remark
16. The corresponding bounds for the reduced one-particle density matrices follow from Lemmas 3.2 and 3.3.

Sketch of Proof
Let us give a brief overview of how the bounds (59) and (65) can be obtained. We also aim to illustrate that the idea to "count the number of particles" not in the antisymmetrized product is very natural and has a clear physical interpretation which is reflected in the proof. We first consider the weight function n(k) = k/N. Note that q t 1 is a solution to the Heisenberg equation of motion i∂ t q t is the "mean-field Hamiltonian" from the right-hand side of (3), with V Using the antisymmetry of ψ t , we then find where, using p t 1 + q t 1 = 1 = p t 2 + q t 2 , Let us emphasize here that the kinetic and external field terms coming from the Schrödinger and the fermionic mean-field equations cancel, which is why the theorems in Section 3.2 hold for any H 0 j . Also, terms (I I ) n and (I I I ) n do not depend on the mean-field V (N) 1 ; it is only term (I ) n where a difference between Schrödinger interaction and the mean-field enters. Note that these three terms have an intuitive explanation. The change ∂ t α n (t) in the "number of bad particles" is given by three different kind of transitions due to the microscopic interaction: transitions between one good, one bad particle and two good particles (term (I ) n ); between two bad and two good particles (term (I I ) n ); and transitions between two bad and one good, one bad particle (term (I I I ) n ). Note, however, that for α n (0) = 0, i.e., ψ 0 = N j =1 ϕ 0 j , we find ∂ t α n (t) t=0 = 0. Thus, what causes the deviations from the antisymmetrized product structure in the first place, is a higher order effect, viz., the fluctuations around the mean-field. The exponential growth in time in (48) comes from the fact that the larger the number of bad particles, the more possibilities do good particles have to become correlated, i.e., as in (47), the change in the number of bad particles at some time t is expected to be proportional to the number of bad particles at this time, plus fluctuations.
Let us briefly discuss how the three terms in (66) can be bounded rigorously (Lemma 7.3). Let us first show why the terms can be bounded by C N α n (t) plus error terms of O(N −1 ), where C N denotes some possibly N-dependent constant, and afterwards discuss why C N in fact does not depend on N for the scalings discussed in Sections 1.1 and 1.2. Recall that q t 1 ψ t = ψ t , q t 1 ψ t = √ α n (t). Thus, by using Cauchy-Schwarz, term (I I I ) n is bounded by C N α n (t), assuming v 12 p t 1 op is bounded. In term (I I ) n we also have two q's available, but both q's are left from v (N) 12 , so we cannot directly apply Cauchy-Schwarz. However, there is a trick that uses the antisymmetry of ψ t to shift the q t 2 to the right side of the scalar product, on the expense of a boundary term of O(N −1 ). Then, by Cauchy-Schwarz, we can again conclude (I I ) n ≤ C N (α n (t) + N −1 ). Term (I ) n is the crucial term in the sense that there is only one q-projector available (which is not enough for the desired bound) and one thus has to use a cancellation between the microscopic interaction and the mean-field. From the proof of Lemma 7.3 we see that one can indeed extract an extra q-projector from the difference (N − 1)p t 2 v (N) 1 , again at the expense of an error term of O(N −1 ), which yields the desired bound. 2 For Coulomb interaction and the corresponding scaling v (N) = N −2/3 v, there is a prefactor N 1/3 in front of the scalar products in the three terms from (66). For the desired bound one thus needs to gain an extra N −1/3 from the scalar products, so that the constants C N from above are of O(1), i.e., N-independent. For the case discussed in Section 1.2, it is the long-range decay of the Coulomb interaction that gives us this additional factor. In the estimates, the p projectors are used to smooth out the singularity, which leads to terms like v * ρ t N or v 2 * ρ t N . From these terms we gain the additional N −1/3 by using kinetic energy inequalities (Lieb-Thirring inequalities, see Section 8) which take into account the orthonormality of ϕ t 1 , . . . , ϕ t N (fermi pressure), i.e., the fact that the system volume is of order N when the total kinetic energy is bounded by AN. This leads to the results in Section 2.1. In the semiclassical case, we have to gain the additional N −1/3 by other means, namely from condition (34). In the estimates, this is reflected by using N −1 p t 1 e ikx q t 1 tr ≤ c(t)N −1/3 , see Section 9. Finally, let us discuss where and how the use of the weight function (46) helps in the estimates. For strong singularities or weak scalings, term (I I I ) n can be problematic, since there is only one projector p available to smooth out the singularity and gain the additional N −1/3 by the kinetic energy inequalities. Here we would like to improve the estimate. For general weight functions, ∂ t α f (t) is calculated in Section 6. Each of the transitions between good and bad particles is now weighted with the corresponding change in the weight function, i.e., a derivative of the weight function appears. More exactly, ∂ t α f (t) is again given by three terms similar to (66), but with ψ t replaced by where f (d) is a discrete derivative of f . (This is consistent with (66), since the derivative of n(k) = k/N equals N −1 .) The key observation is now that the derivative of the weight function m (γ ) is zero for all k > N γ , i.e., all contributions coming from P (N,k) ψ t with k > N γ vanish. This helps us to make term (I I I ) m (γ ) small, which only gives contributions from large k due to the three projectors q. 2 Note that here lies the crucial difference compared to the bosonic Hartree equation, as treated in [39,51]. There, since there is just one orbital ϕ, i.e., p m = |ϕ ϕ| m , and since v (N) = N −1 v, one can directly calculate that V i.e., the corresponding term (I ) n can directly be seen to be bounded by C N α n (t) + N −1 . (Note that in [39,51], the term v * |ϕ| 2 (x 1 ) q 2 is usually regarded as being part of term (I I I ) n .)

Notation and Preliminaries
We summarize here some notation, list often used inequalities and establish some basic properties of the projectors p ϕ m that we use in the proofs throughout the following sections. The Hilbert space of complex square integrable functions on R d is denoted by We denote by ·, · a+1,...,N the scalar product only in the variables x a+1 , . . . , x N , i.e., it is a "partial trace" or "partial scalar product", formally defined for any χ, ψ ∈ L 2 (R 3N ) by which should be regarded as a vector in L 1 (R 3a ) (for χ = ψ, it is the diagonal of the reduced a-particle density matrix, a generalization of (19)). In the same style we denote by |· m a vector in L 2 (R 3 ) acting only on the m-th variable of L 2 (R 3 ) N , by ·| m its dual, and by ·, · m the scalar product only in the m-th variable. Given a function h : For the operator norm of the projectors p ϕ m it makes an important difference if it is calculated on all L 2 functions or only on antisymmetric functions in L 2 , as the following lemma shows.

Proof of Results About Density Matrices
In this section, we prove the results from Section 3.1 about the relation of α f to the reduced one-particle density matrices of ψ and of the antisymmetrized product state N j =1 ϕ j . For a bounded operator A on a Hilbert space H , the operator norm, trace norm and Hilbert-Schmidt norm are defined by where {φ i } i∈N is some orthornormal basis and A * denotes the adjoint of A. Note that (for proofs, see, e.g., [ Proof of Lemma 3.2 First, note that reduced one-particle density matrices have the well-known properties that they are non-negative, i.e., f, γ Recall that we are here concerned with α n = ψ, q 1 ψ , i.e., the α-functional with the weight n(k) = k/N. Also, recall that γ ϕ j 1 = N −1 p 1 . We first show that N −1 p 1 − p 1 γ ψ 1 p 1 is a non-negative operator with trace norm α n . Note that the operator p 1 γ ψ 1 p 1 maps the N-dimensional subspace span(ϕ 1 , . . . , ϕ N ) to itself, and is non-negative and self-adjoint. We can therefore diagonalize it, i.e., there is an orthonormal basis {χ 1 , . . . , χ N }, such that with 0 ≤ λ i ≤ N −1 ∀i = 1, . . . , N, since (see also Lemma 4.1) We now show γ Since 0 ≤ α n ≤ 1, it is indeed true that (91) Note that indeed ||p 1 − q 1 || op = 1, since for all f ∈ L 2 (R 3 ), We now show γ We now show α n ≤ √ N γ ϕ j 1 − γ ψ 1 HS . Using first (87) and then (84), ||p 1 || op = 1 and ||p 1 || HS = √ N , we find Proof of Lemma 3.3 The inequalities (56) follow directly from α n ≤ α f and Lemma 3.2. Let us denote the floor function by · , i.e., for any x ∈ R, x = max{m ∈ Z : m ≤ x}. Then

The Time Derivative of α f (t)
The expression for the time derivative of for arbitrary weight functions f (k) follows from direct calculation (recall that the projectors P (N,k) depend on t through their dependence on ϕ t 1 , . . . , ϕ t N ). We use that the time derivative of ψ t ∈ L 2 (R 3N ) is given by the Schrödinger equation (1) with self-adjoint Hamiltonian H given by (2) with real interaction v (N) (x) = v (N) (−x). As effective equations for the one-particle wave functions ϕ t 1 , . . . , ϕ t N ∈ L 2 (R 3 ) we consider the general equations where the operator V (N) can possibly depend on N, j and ϕ t 1 , . . . , ϕ t N . The two interesting cases are when there is only direct interaction, where ρ t N = N i=1 |ϕ t i | 2 , and when there is direct and exchange interaction, We calculate the time derivative of α f (t) in two steps. First, in Lemma 6.1 we directly calculate a straightforward expression for the time derivative. We then prove an auxiliary lemma which we use to bring the time derivative into a form (Lemma 6.5) that can nicely be estimated (later in Lemma 7.3). Lemma 6.1 Let ψ t ∈ L 2 (R 3N ) be an antisymmetric solution to the Schrödinger equation (1) and let ϕ t 1 , . . . , ϕ t N ∈ L 2 (R 3 ) be orthonormal solutions to the effective equations (97). Then where Proof Note that the operators p m , q m , P (N,k) all depend on t through the orbitals ϕ t 1 , . . . , ϕ t N . For ease of notation, we do not explicitly write out this t-dependence. Let us first prove that f fulfills the Heisenberg equation of motion Note that i∂ t p m = H mf m , p m and, using p m + q m = 1, i∂ t q m = H mf m , q m . Then, by applying the product rule and using H mf m , p n = 0∀m = n, it follows that which implies (101). Using this and the antisymmetry of ψ t , we find Let us now bring the expression (100) into a form we can use for the desired Gronwall estimate. For that, we need to define shifted operators and discrete derivatives of weight functions f . Definition 6.2 Let d ∈ Z and f be a weight function as in Definition 3.1.

We define the shifted operators
where for the last step we defined f (k) = 0 for all k < 0 and k > N. 2. For d > 0, we define discrete first derivatives of f by where 1 A (k) = 1 for k ∈ A and 0 otherwise. For Lemma 6.4 and the proof of Lemma 6.5 we need some more notation.
where the subscript sym means the symmetrized tensor product.
In accordance with Definition 3.1 we have P (k) 1...N = P (N,k) . For the proof of Lemma 6.5 we need an auxiliary lemma that shows how to shift an f to the other side of a two-particle operator. Lemma 6.4 As in Definition 6.3, we abbreviate P (0) 12 = p 1 p 2 , P (1) 12 = p 1 q 2 + q 1 p 2 and P (2) 12 = q 1 q 2 . Let h 12 be an operator that acts only on the first and second particle index. Then, for all a, b = 0, 1, 2, Proof We only prove (108) since (109) can be proven in just the same way. In the following calculation we use the splitting where P (d) 12 contains exactly d q-projectors, and P Note that Lemma 6.4 holds more generally for m-particle operators h 1...m , i.e., for all a, b ∈ {0, . . . , m}. As last preparatory remark, note that for all f ≥ 0, since P (N,k) P (N, ) = δ k P (N,k) . 3 With Lemma 6.4 we can simplify the expression (100) for the time derivative of α f (t) by splitting it into three parts, each of which will be estimated separately later in Lemma 7.3. The advantage of having the square root of the "derivatives" f (d) next to ψ t is made clear in Lemma 7.1 and the estimates in Lemma 7.3. with Remark 17. Note that the time derivative is formally the same as for bosons, where p 1 := |ϕ ϕ| 1 , see [50,51]. Note that in [50,51] the splitting into three summands is done slightly differently: compared to (114), an additional identity 1 = p 2 + q 2 is added in front of V is pulled over to the left side of the scalar product. 18. For the case f (k) = n(k) = k/N we find a simple expression for the time derivative of α n (t) = ψ t , q t 1 ψ t . It can easily be found by direct calculation as in (66) or from (114). Note that, in view of Definition 3.1 and the identity N k=0 P (N,k) = 1, and similarly Then (114) simplifies to 19. Note that the proof of Lemma 6.5 can easily be generalized to m-particle interactions W 1...m . With the proper definition of W 1...m coming from Lemma 6.1 and using (112), we find Proof of Lemma 6. 5 We calculate the time derivative of α f (t) starting from the expression (100) of Lemma 6.1. The idea of the proof is to insert two identities 1 = p 1 + q 1 and 1 = p 2 + q 2 in front of each ψ t (which leads to 16 summands) and then to use Lemma 6.4 in order to shift f . It turns out that a lot of terms drop out due to the commutator structure. Let us first note two auxiliary calculations. For all a, b = 0, 1, 2 with a > b, we have where the positivity is true since f was assumed to be monotone increasing, and Then, by inserting two identities 1 = 2 a=0 P (a) 12 (in the following sums, the indices a, b always run from 0 to 2), we find Equation (114) follows by inserting the definition 2 , using p m q m = 0 and using p 1 q 2 = p 1 − p 1 p 2 in the third summand.

Using the New Weight Function
In Section 6 we went through some trouble to handle general weight functions in the time derivative of α f (t). The next lemma shows what we gain by choosing the weight function m (γ ) (k) from (46). Morally, the lemma says that by choosing γ < 1 compared to γ = 1, we make the convergence rate worse when there is no q available (see (125)), we do not loose anything when there is one q available (see (126)) and we gain powers in N when there are two q's available (see (127)). and Remark 20. Note that more generally the estimate q 1 . . . q n m (±d) holds for all 1 < d < N and 1 < n < N, for some constant C that depends only on d and n.
Proof of Lemma 7.1 First, recall that for all antisymmetric φ ∈ L 2 (R 3N ), and in a similar way (130) Recall that we denote the floor function by · , i.e., for any x ∈ R, x = max{m ∈ Z : m ≤ x}. Then we find and, by using (129), and, by using (130), Careful consideration of the boundary terms around the k = N γ summands gives the values of the constants C(d) andC(d).

Diagonalization of p 2 h 12 p 2 and Related Lemmas
For handling the terms in the time derivative of α f (t) from (114), it is often useful to diagonalize operators p 2 h 12 p 2 . Let us briefly summarize what we mean by that and introduce some notation on the way. For any h : , the operator p 2 h 12 p 2 is a multiplication operator in x 1 and a projector onto the Ndimensional subspace span(ϕ 1 , . . . , ϕ N ) in the second variable. Therefore one can write it as an x 1 dependent non-negative self-adjoint (N × N)-matrix acting on the second variable, i.e., (recall our notation from Section 4) where the eigenvectors |χ x 1 i 2 are orthonormal and can be written as |χ . . . , ϕ N ) for all x 1 ∈ R 3 and that the projector p 2 is independent of the choice of basis, i.e., The eigenvalues λ i (x 1 ) have the properties that and furthermore, for all i = j , where this inequality remains true with (N(N −1)) −1 replaced by ((N −1)(N − 2)) −1 , when ψ is antisymmetric in all variables except x 3 .
Proof Recall our notation for the "partial scalar product" from (73). For proving (139) we use the diagonalization (134), Lemma 4.1 and (137). We find For proving (140) we diagonalize We call the eigenvalues μ j and the eigenvectorsφ j . With the diagonalization (134) and Lemma 4.1 we find If ψ is antisymmetric in all variables except x 3 , then by Lemma 4.1 one can only extract factors (N − 2) −1 instead of (N − 1) −1 , and (N − 1) −1 instead of N −1 from the antisymmetry of ψ.

Bounds on ∂ t α f (t)
We now give the rigorous bounds for the three terms in the time derivative of α f (t) given by (114). Here, we use the weight function m (γ ) (k) from (46). This also contains the case where γ = 1, i.e., where the weight function is n(k). The estimates are collected in the following lemma, which constitutes the heart of the proof of our main results.
We state this lemma only for positive v (N) . If v (N) contains both positive and negative parts, we later decompose v (N) = v (N) . Then, using the weight function m (γ ) (k) from (46), we find for all 0 < γ ≤ 1 for the three terms from (114), Proof Term (I ) m . In the following, we diagonalize the operator p m v (N) 1m p m and use the same notation as in (134). We split each eigenvalue into two parts, Note that by using Cauchy-Schwarz, and which remains true if |ϕ i m = |χ x 1 i m (m ≥ 2). In the following we abbreviate φ = m (1) 1/2 ψ andφ = m (−1) 1/2 ψ. Note that both φ andφ are still antisymmetric. Then we find For the first summand in (151), we find by Cauchy-Schwarz, For the second summand in (151), we find ψ. The idea of the bound for this term is to shift one q to the right side of the scalar product by using the antisymmetry of ψ. Using Cauchy-Schwarz and the antisymmetry of φ andφ, we find Term (I I I ) m . We abbreviate φ = m (1) 1/2 ψ andφ = m (−1) ψ. By Cauchy-Schwarz we find

Proof of Theorems 3.4 and 3.6
Proof of Theorems 3.4 and 3.6 First, we split v (N) = v (N) where Term ± refers to (I ) m (γ ) , (I I ) m (γ ) , (I I I ) m (γ ) from (115), (116) and (117) with interaction v (N) ± . We bound Term ± separately by using Lemma 7.3, which proves under the stated assumptions the bound Applying the Gronwall Lemma gives the desired bound The values of the constant C(t) in Theorems 3.4 and 3.6 can be obtained by using the respective assumptions, together with Lemma 7.3 and (156).

21.
Following up on Remark 15 after Theorem 3.4, let us consider the size of the error we make by neglecting the exchange term. Suppose that the exchange term is of O(N −δ ). It then gives an additional term CN −δ √ α n (t) ≤ C α n (t) + N −2δ in the time derivative of α n (t) (where the √ α n (t) comes from the q 1 in term (I ) n ).

Proof of Results for Density ∝ 1 Regime
In this section, C denotes a constant which can be different from line to line.

Kinetic Energy Inequalities
Let us first recall two well-known inequalities which we use in Section 8.2 to show that the conditions of Theorems 3.4 and 3.6 hold if the total kinetic energy is bounded by AN. A general version of the Lieb-Thirring or kinetic energy inequality [43,44,55] for orthonormal ϕ 1 , . . . , ϕ N ∈ L 2 (R 3 ) is for any a > 0, where ρ N = N i=1 |ϕ i | 2 . The Hardy-Littlewood-Sobolev inequality (see, e.g., [42,Thm. 4.3]) in three dimensions states that for f ∈ L p (R 3 ), h ∈ L r (R 3 ), p, r > 1 and 0 < λ < 3 with 1/p + λ/3 + 1/r = 2, there is a constant C = C(λ, p) such that Let us now show for interactions |x| −s , 0 < s < 6/5, that the mean-field term v (N) * ρ N is bounded independent of N, if it is scaled with β = 1 − s/3 and if the total kinetic energy is bounded by AN. We need this statement for the proofs in Section 8.2.
Proof of Theorem 2.4 The strategy of the proof is again to bound |∂ t α n (t)| ≤ Ce Ct α n (t) + N −1 (184) and use the Gronwall Lemma and then Lemmas 3.2 and 3.3 to conclude the desired bound (36). Recall that we use the weight function n(k) = k/N here, i.e., we can use the form (120) for the time derivative of α n (t). Using the scaling v (N) = N −1 v and noting the additional N −1/3 in front of the time derivatives in the Schrödinger and mean-field equations, we find that ∂ t α n (t) is given by the sum of the three terms (I ) n = 2N −2/3 Im ψ t , q t 1 (N − 1)p t 2 v 12 p t 2 − V 1 p t 1 ψ t , (I I ) n = 2N −2/3 Im ψ t , q t 1 q t 2 (N − 1)v 12 p t 1 p t 2 ψ t , (I I I ) n = 2N −2/3 Im ψ t , q t 1 q t 2 (N − 1)v 12 p t 1 q t 2 ψ t , with V 1 = V dir 1 in the case of the fermionic Hartree equations, and V 1 = V dir 1 + V exch 1 in the case of the Hartree-Fock equations. For ease of notation we write ψ t = ψ, ϕ t i = ϕ i in the following. The double exponential in (36) comes from one exponential in Lemma 9.1 and another exponential from the Gronwall Lemma applied to (184). In the following estimates, we decompose v(x) = d 3 kv(k)e ikx . Note that the assumption d 3 k (1 + |k| 2 ) |v(k)| < ∞ in particular implies that d 3 k |v(k)| < ∞ and d 3 k |k| |v(k)| < ∞.
Term (I ) n . Let us first bound the contribution from the exchange term. Using the Fourier decomposition of v and Cauchy-Schwarz we find Here we see explicitly that the contribution from the exchange term is of lower order in N. Let us now bound (I ) n only with direct interaction. Using the Fourier decomposition of v we find Similar to Lemma 7.4 we would like to diagonalize the operator p 2 e −ikx 2 p 2 . However, since it is not self-adjoint, we decompose e −ikx = cos(kx) − i sin(kx) and diagonalize the self-adjoint operators where the real eigenvalues λ j ,λ j and orthonormal eigenvectors χ j ,χ j depend on k. Note that λ j = χ j , cos(kx)χ j , so |λ j | ≤ χ j 2 = 1 and N j =1 λ j = N j =1 χ j , cos(kx)χ j = N j =1 ϕ j , cos(kx)ϕ j , and analogous forλ j . In the following, we use the projector q χ j =1 = 1 − N m=2 p χ j m introduced in (149) and the singular value decomposition q 1 e ikx 1 p 1 = μ |φ φ | 1 (see Lemma 9.2), with μ = q 1 e ikx 1 p 1 tr . Let us now decompose term (I ) n by using e −ikx = cos(kx) − i sin(kx). For the cos-term we find, using the antisymmetry of ψ and Cauchy-Schwarz, The same bound holds for the sin-term. Term (I I ) n . Similarly to Lemma 7.4, we use the antisymmetry of ψ to shift one q to the right side of the scalar product, i.e., |(I I ) n | = 2N −2/3 |Im ψ, q 1 q 2 (N − 1)v 12 p 1 p 2 ψ | = 2N −2/3 Im ψ, q 1 N m=2 q m v 1m p 1 p m ψ ≤ 2N −2/3 ||q 1 ψ|| q 1 N m=2 q m v 1m p 1 p m ψ ≤ 2N −2/3 ||q 1 ψ|| N 2 ψ, q 3 p 1 p 2 v 12 q 1 v 13 p 1 p 3 q 2 ψ + N ψ, p 1 p 2 v 12 q 1 q 2 v 12 p 1 p 2 ψ For the following estimate, note that for all trace class operators A 1 , B 2 , we find by the singular value decomposition and Lemma 4.1, ψ, q 3 A 1 B 2 q 3 ψ = j, μ j μ ψ, q 3 |φ j φ j | 1 |φ φ | 2 q 3 ψ ≤ C ||A|| tr ||B|| tr N −2 ψ, q 3 ψ .