Higher order corrections to the mean-field description of the dynamics of interacting bosons

In this paper, we introduce a novel method for deriving higher order corrections to the mean-field description of the dynamics of interacting bosons. More precisely, we consider the dynamics of $N$ $d$-dimensional bosons for large $N$. The bosons initially form a Bose-Einstein condensate and interact with each other via a pair potential of the form $(N-1)^{-1}N^{d\beta}v(N^\beta\cdot)$ for $\beta\in[0,\frac{1}{4d})$. We derive a sequence of $N$-body functions which approximate the true many-body dynamics in $L^2(\mathbb{R}^{d N})$-norm to arbitrary precision in powers of $N^{-1}$. The approximating functions are constructed as Duhamel expansions of finite order in terms of the first quantised analogue of a Bogoliubov time evolution.


Introduction
We consider a system of N bosons in R d , d ≥ 1, interacting with each other via pair interactions in the mean field scaling regime. The Hamiltonian of the system is given by Here, V ext denotes some possibly time-dependent external potential, and the interaction potential v β is defined as v β (x) : for some bounded, spherically symmetric and compactly supported function v : R d → R. In the following, we will make use of the abbreviation Note that the prefactor (N − 1) −1 in front of v β is chosen such that the interaction energy and the kinetic energy per particle are of the same order. The mean inter-particle distance is of order N − 1 d and therefore much smaller than the range of the interaction, which scales as N −β . Hence, on average, every particle interacts with many other particles, and the interactions are weak since (N − 1) −1 N dβ → 0 as N → ∞. This implies that we consider a mean-field regime. In particular, the case β = 0 is known as the Hartree scaling regime.
In this paper, we study the time evolution of the N -body system for large N when the bosons initially exhibit Bose-Einstein condensation. We impose suitable conditions on the external potential V ext (t) such that H β (t) is self-adjoint on D(H β (t)) = H 2 (R dN ) for each t ∈ R. Consequently, H β (t) generates a unique family of unitary time evolution operators {U (t, s)} t,s∈R via the Schrödinger equation The N -body wave function at time t ∈ R is determined by for some initial datum ψ(0) = ψ 0 ∈ L 2 sym (R dN ). Due to the interactions, the characterisation of the time evolution U (t, s) is a difficult problem. Even if the system was initially in a factorised state, where all particles are independent of each other, the interactions instantaneously correlate the particles such that an explicit formula for U (t, s) is quite inaccessible.
To describe U (t, s) approximatively, one observes that the dynamics of the many-body system can be decomposed into • the dynamics of the condensate wave function ϕ(t) ∈ L 2 (R d ), and • the dynamics of the fluctuations around the (time-evolved) condensate.

A first order approximation to the N -body dynamics
A first approximation to the N -body dynamics is provided by the time evolution of the condensate wave function. Its dynamics yield a macroscopic description of the Bose gas, which, in the limit N → ∞, coincides with the true dynamics in the sense of reduced density matrices. In order to formulate this mathematically, one assumes that the system is initially in a Bose-Einstein condensate with condensate wave function ϕ 0 , i.e., lim N →∞ where γ (1) (t) := Tr 2,...,N |ψ(t) ψ(t)| is the one-particle reduced density matrix of ψ(t) at time t. Then it has been shown, see e.g. [1,2,11,13,18,19,32,53], that lim N →∞ Tr γ (1) (t) − |ϕ(t) ϕ(t)| = 0 (7) for any t ∈ R, where ϕ(t) is the solution of the Hartree equation with initial datum ϕ(0) = ϕ 0 and with Note that for β = 0, the equation (8) is the N -independent Hartree (NLH) equation. For β > 0, the evolution is N -dependent and converges to the non-linear Schrödinger (NLS) dynamics with N -independent coupling parameter v in the limit N → ∞. The parameter µ ϕ(t) is a real-valued phase factor, which we choose as for later convenience. For the convergence with respect to reduced densities, this phase is irrelevant since it cancels in the projection |ϕ(t) ϕ(t)|. One way to prove the convergence (7), and consequently to derive the NLH/ NLS equation from a system of N bosons, is via the so-called BBGKY 1 hierarchy, which was prominently used in the works of Lanford for the study of classical mechanical systems in the infinite particle limit [36,37]. The first derivation of the NLH equation via the BBGKY hierarchy was given by Spohn in [54], and this was further pursued, e.g., in [1,2,20,21]. About a decade ago, Erdős, Schlein and Yau fully developed the BBGKY hierarchy approach to the derivation of the NLH/NLS equation in their seminal works including [18,19]. Subsequently, a crucial step of this method was revisited by Klainerman and Machedon in [33], based on reformulating combinatorial argument in [18,19] and a viewpoint inspired by methods of non-linear PDEs. This, in turn, motivated many recent works on the derivation of dispersive PDEs, including [11,12,13,14,15,32,53]. In [52], Rodnianski and Schlein introduced yet another method for proving (7), which uses coherent states on Fock space and was inspired by techniques of quantum field theory and the pioneering work of Hepp [29].
In the context of the current paper, the most relevant works on the derivation of the NLH/NLS equation are due to Pickl [50,51], who introduced an efficient method for deriving effective equations from the many-body dynamics, transforming the physical idea behind the mean-field description of an N -body system into a mathematical algorithm. Instead of describing the condensate as the vacuum of a Fock space of fluctuations, this approach remains in the N -body setting and uses projection operators to factor out the condensate. This strategy was successfully applied to prove effective dynamics for N -boson systems in various situations, e.g., [4,8,17,30,31,34,40,41].
A much stronger notion of distance than the one expressed in (7) is provided by the L 2 (R dN )-norm. Whereas closeness in the sense of reduced densities implies that the majority of the particles (up to a relative number that vanishes as N → ∞) is in the state ϕ(t), the norm approximation requires the control of all N particles. In particular, this implies that the fluctuations around the condensate can no longer be omitted from the description. In this sense, the norm approximation of ψ(t) can be understood as next-to-leading order correction to the mean-field description.
For the dynamics U (t, s), a norm approximation in d dimensions was proven in [38] for β = 0 and V ext = 0 under quite general assumptions on the interaction potential v. In [44], this result was extended to the range β ∈ [0, 1 3 ) for the three-dimensional defocusing case, and in [45], the focusing case in dimensions one and two was treated for β > 0 and β ∈ (0, 1), respectively. In these works, the authors consider initial data of the form for some appropriate initial fluctuation vector It is then shown that there exist constants C, C > 0 such that where δ = 1 for β = 0, δ = 1 − 3β for the three-dimensional defocusing case with β ∈ [0, 1 3 ), and δ = 1 2 and δ < 1 3 (1 − β) for the one-and two-dimensional focusing case, respectively. The fluctuations Here, H Bog (t) denotes the Bogoliubov Hamiltonian 2 , an effective Hamiltonian in Fock space which is quadratic in the number of creation and annihilation operators. For three dimensions and scaling parameter β = 0, a similar result was obtained in [42,43] via a first quantised approach. More precisely, denote p ϕ(t) j := |ϕ(t, x j ) ϕ(t, x j )| 2 Written in second quantized form, HBog(t) is defined as where a * x and ax denote the operator-valued distributions corresponding to the usual creation and annihilation operators on F(L 2 (R d )). Besides, K1(t) := Q(t) K1(t)Q(t) with Q(t) := 1 − |ϕ(t) ϕ(t)|, where K1 is the Hilbert-Schmidt operator on L 2 (R d ) with kernel K1(t, x, y) := ϕ(t, x)v β (x − y)ϕ(t, y). Further, K2(t) := (Q(t) ⊗ Q(t)) K2(t), where the kernel of the two-body function K2(t) is given by K2(t, x, y) := ϕ(t, x)v β (x − y)ϕ(t, y) (e.g. [44,Equation (31)]). and q ϕ(t) j The auxiliary N -particle Hamiltonian H ϕ(t) (t) is defined by subtracting from H β (t) in each coordinate the mean-field Hamiltonian h ϕ(t) (t) from (8), inserting identities (p on both sides of the difference, and discarding all terms which are cubic, C ϕ(t) , or quartic, Q ϕ(t) , in the number of projections q ϕ(t) (see Lemma 2.2). This yields which has a quadratic structure comparable to that of the Bogoliubov-Hamiltonian H Bog (t): (t), which form an effective two-body potential, contain exactly two projectors q ϕ(t) onto the complement of the condensate wave function, while H Bog (t) is quadratic in the creation and annihilation operators of the fluctuations. However, H ϕ(t) (t) is particle number conserving and acts on the N -body Hilbert space L 2 (R dN ), i.e., it determines the evolution of both condensate wave function and fluctuations. In contrast, H Bog (t) operates on the fluctuation Fock space F {ϕ(t)} ⊥ , does not conserve the particle number, and exclusively concerns the dynamics of the fluctuations with respect to the condensate wave function evolving according to (8).
Under appropriate assumptions on the initial datum ψ 0 , the time evolution U ϕ (t, s) generated by H ϕ(t) (t) approximates the actual time evolution U (t, s). More precisely, there exist constants C, C > 0 such that [42, Theorem 2.6]. Further, in the limit N → ∞, the fluctuations in U ϕ (t, 0)ψ 0 coincide with the solutions of the Bogoliubov evolution equation: let ξ ϕ 0 = ξ (k) ϕ 0 N k=0 denote the fluctuations around ϕ 0 ⊗N in the initial state ψ 0 under the decomposition (5), let ξ ϕ(t) = ξ (k) ϕ(t) N k=0 denote the fluctuations around ϕ(t) ⊗N in U ϕ (t, 0)ψ 0 , and let χ(t) = χ (k) (t) k≥0 denote the solutions of (13) with initial datum ξ ϕ 0 for 0 ≤ k ≤ N and ξ (k) [42, Lemma 2.8]. Hence, the combination of (15) and (16) yields (12), with a different timedependent constant but the same N -dependence. Beyond the mean field regime, a statement similar to (12) was shown in [46] for the range β ∈ [0, 1 2 ). For larger values of the scaling parameter, the evolutions of ϕ(t) and ξ ϕ(t) do not (approximately) decouple any more as a consequence of the short-scale structure related to the two-body scattering process. For β ∈ (0, 1), an accordingly adjusted variant of (12) for appropriately modified initial data was obtained in [9] in the three-dimensional defocusing case. Similar estimates for the many-body evolution of appropriate classes of Fock space initial data have been obtained in [6,16,22,23,25,26,27,28,35,52] for various ranges of the scaling parameter. A related result for Bose gases with large volume and large density was proved in [49].

Higher order approximations to the N -body dynamics
In this paper, we introduce a novel method for deriving a more precise characterisation of the dynamics, which approximates the N -body wave function to arbitrary order in powers of N −1 . This is achieved by constructing a sequence of N -body wave functions, which are defined via an iteration of Duhamel's formula with the time evolution U ϕ (t, s) generated by H ϕ(t) (t). We work in the first quantized framework as was the case, e.g., in [42].
• In our first result, we estimate the growth of the first A moments of the number of fluctuations when the system evolves under the dynamics U (t, s) or U ϕ (t, s). Estimates of this kind are often needed to derive effective descriptions of the dynamics of interacting bosons, e.g., in [5,6,10,42,49,52] In the remainder of the paper, we assume that for some A ∈ {1, ..., N }, the first A moments of the number of fluctuations in the initial state are sub-leading (see Assumption A3). More precisely, let γ ∈ (0, 1]. We assume that for all a ∈ {0, ..., A}, there exists some constant C(a) depending only on a such that Here, ξ ϕ 0 denotes the fluctuation vector corresponding to the initial state ψ 0 as in (5), and N ϕ 0 is the number operator on the Fock space F ≤N {ϕ 0 } ⊥ of fluctuations around ϕ 0 ⊗N . Note that γ = 0 corresponds to the trivial bound ξ ϕ 0 , N a ϕ 0 ξ ϕ 0 ≤ N a . In this sense, our assumption states that the expected number of fluctuations in ψ 0 is sub-leading. Clearly, the larger we choose γ, the stronger is the assumption.
• Under these conditions, we show in Corollary 2.5 that at any time t and for sufficiently large N , the first A moments of the number of fluctuations remain sub-leading, and the N -dependence N (1−γ)a in (17) is replaced by N c(β,γ)a for some (1 − γ) ≤ c(β, γ) < 1.
• In our second and main result (Theorem 1), we prove higher order corrections to the norm approximation (12) for the scaling regime β ∈ [0, 1 4d ). This is to be understood in the following sense: we construct a sequence of N -body wave functions {ψ for some time-dependent constant C(t). The exponent δ(β, γ) is positive, depends on β and γ and is determined in Theorem 1.
The first element of the approximating sequence {ψ For a = 1, the estimate (18) is thus well known since it coincides with the norm approximation (15) and consequently with (12). To obtain the next higher correction with respect to N , we add an appropriate correction term to ψ (1) ϕ (t). We expand the difference U (t, s) − U ϕ (t, s) ψ 0 using Duhamel's formula, identify the leading order contribution, and approximate it by replacing U (t, s) with U ϕ (t, s). This leads to the second element For the third element, we expand the difference U (t, s) − U ϕ (t, s) ψ 0 to the next order, using Duhamel's formula twice, and subsequently follow the same strategy as before. In this way, we construct all higher elements of the sequence as Duhamel expansions with finitely many terms, each of which exclusively contains ψ 0 , the auxiliary time evolution U ϕ (t, s), and the cubic and quartic interaction terms C ϕ(t) and Q ϕ(t) . The precise definition of ψ (a) ϕ (t) for any a, as well as a more detailed explanation of the construction, is provided in Definition 2.2 and the preceding discussion.
We note that higher order approximations of the reduced density matrices were obtained by Paul and Pulvirenti in [47] for β = 0 and factorized initial data, based on the method of kinetic errors from the paper by Paul, Pulvirenti and Simonella [48]. For j ∈ {1, ..., N }, the authors of [47] construct a sequence {F N,n j (t)} n∈N of trace class operators on L 2 (R jd ), which approximate the j-particle reduced density matrix γ (j) (t) of the system with increasing accuracy up to arbitrary precision. The approximating operators F N,n j (t) can be determined by a number of operations scaling with n. They depend on the initial data as well as the knowledge of the solution of the Hartree equation and its linearization around this solution.
Due to different methods used, it is not straightforward to compare the results of [47] with the results of this paper. However, we list some features of our paper that differ from the operator-based method of kinetic errors [47,48]. In contrast to the approach in [47], we derive approximations directly for the time-evolved N -body wave function. Our construction is in terms of the Bogoliubov time evolution U ϕ instead of the linearized Hartree flow, and it is implemented as a robust algorithm that requires an a-dependent, N -independent number of explicit calculations to compute the a'th order approximation. Moreover, the results obtained in this paper cover more generic initial data satisfying (17) and include positive values of β.
Notation. In the following, any expression C that is independent of both N and t will be referred to as a constant. Note that constants may depend on all fixed parameters of the model such as ϕ 0 , ψ 0 , v and V ext (0). Further, we denote A B and A B to indicate that there exists a constant C > 0 such that A ≤ CB, resp. A ≥ CB, and abbreviate Finally, we use the notation 2 Main results

Framework and assumptions
Let us first recall from [39,42,43] the explicit decomposition of an N -body wave function ψ in terms of a condensate ϕ ⊗N and k-particle fluctuations around this condensate. To this end, recall the following projections introduced in [50]: and the corresponding projection operators on L 2 (R dN ) For 0 ≤ k ≤ N , define the many-body projections and P ϕ k = 0 for k < 0 and k > N . Further, for any function f : We will in particular need the operators n ϕ and m ϕ corresponding to the weights The part of ψ in the condensate is given by P ϕ 0 ψ, and the part of ψ corresponding to k particles fluctuating around the condensate is precisely P ϕ k ψ for k ≥ 1. By construction, for some ξ To determine the explicit form of ξ (k) ϕ , observe that by Definition 2.1, where, by definition of the symmetric tensor product, ϕ is symmetric under permutations of all of its coordinates, and ξ (k) ϕ is orthogonal to ϕ in every coordinate, i.e., for every j ∈ {1, ..., k}. Hence, ξ (6). The relation between the N -body state ψ and the corresponding fluctuation vector ξ ϕ is given by the unitary map where ξ ϕ is defined by (20). The vacuum (1, 0, ..., 0) of F ≤N {ϕ} ⊥ corresponds to the condensate ϕ ⊗N , and the probability of k particles being outside the condensate equals by (20). The number operator N ϕ on F ≤N {ϕ} ⊥ , counting the number of fluctuations, is defined by its action (N ϕ ξ ϕ ) (k) := k ξ (k) ϕ . The expected number of fluctuations around the condensate ϕ ⊗N in the state ψ is thus given by with n ϕ from Definition 2.1.
Let us now state our assumptions on the model (1) and on the initial data.
A1 Interaction potential. Let v : R d → R be spherically symmetric and bounded uniformly in both be normalised. Let γ ∈ (0, 1] and A ∈ N. Assume that for any a ∈ {0, ..., A}, there exists a set of non-negative, a-dependent constants {C a } 0≤a≤A with C 0 = 1 such that, for sufficiently large N , Our analysis is valid as long as the solution ϕ(t) of the non-linear equation (8) and depends on the dimension d, the sign of v ϕ(t) , and the regularity of the external trap V ext (t).
Assumptions A1 and A2 are rather standard in the rigorous treatment of interacting manyboson systems. Note that we make no assumption on the sign of the potential or its scattering length and thus cover both repulsive and attractive interactions. Besides, we admit a large class of time-dependent external traps V ext , with basically the only restriction that V ext (t) must not obstruct the self-adjointness of H β (t) on H 2 (R dN ).
The third assumption provides a bound on the expected number of fluctuations around the condensate ϕ 0 ⊗N in the initial state ψ 0 . Note that while γ = 0 is the trivial bound, the condition becomes more restrictive as γ increases. We have chosen this particular formulation of A3 for later convenience 3 . However, its physical meaning is better understood from one of the following two equivalent versions of A3: Assume that for any a ∈ {0, ..., A}, there exists a set of nonnegative, a-dependent constants {C a } 0≤a≤A with C 0 = 1 such that, for sufficiently large N , Assume that for any a ∈ {0, ..., A}, there exists a set of non-negative, a-dependent constants {C a } 0≤a≤A with C 0 = 1 such that, for sufficiently large N , The equivalence A3 ⇔ A3 ⇔ A3 follows immediately from Lemma 2.1, whose proof is postponed to Section 3.1.
3 Note that the operators n ϕ and m ϕ are equivalent in the sense that they are related via (36), namely hence all results in terms of n ϕ can be translated to the corresponding statements in terms of n ϕ . We chose to work with m ϕ instead of n ϕ because this makes in particular Proposition 2.4 easier to write. For example, in terms of n ϕ , Proposition 2.4b reads which contains an additional term N −j+n . Since the proof of our main result requires an iteration of this proposition, the version with n ϕ is more practicable.
Hence, A3 can be understood as follows: Let A ∈ N and consider sufficiently large N such that A = O(1) with respect to N , i.e. A 1. Then we assume that for any a ≤ A, the part of the wave function with any a particles outside the condensate is at most of order N −γa .
Equivalently, A3 states that the first A 1 moments of the number of fluctuations must be sub-leading with respect to the particle number; for γ = 1, they must even be bounded uniformly in N . Here, "sub-leading" means that the moments of the relative number of fluctuations, i.e., the expectation values of (N ϕ(t) /N ) A , vanish as N → ∞. This, in turn, provides a bound on the high components of the fluctuation vector: for example, In other words, it must be very unlikely that significantly many particles are outside the condensate, whereas we impose no restriction on fluctuations involving only few particles (with respect to N ).
As soon as a becomes comparable to N , i.e., a N , the constants C ( , ) a are N -dependent and the assumption is trivially satisfied. However, note that we demand that N be large enough that A 1.
The simplest example of an N -body state satisfying A3 is the product state ψ = ϕ 0 ⊗N . Whereas the ground state of non-interacting bosons (v = 0) is of this form, the ground state as well as the lower excited states of interacting systems are not close to an exact product with respect to the L 2 (R dN )-norm due to the correlation structure related to the interactions.
Besides, it seems reasonable to expect that states exhibiting Bose-Einstein condensation satisfy A3 for some (possibly very small) γ, as it is well known that . Note, however, that we require a certain minimal size of γ, which is strictly greater than 2 3 (Theorem 1). To the best of our knowledge, there exists only one rigorous result in [43,Chapter 3] that identifies situations where a Bose gas satisfies assumption A3. This work concerns a homogeneous Bose gas on the d-dimensional torus and is restricted to the scaling β = 0. For this case, it is shown that the ground state as well as the lower excited states fulfil assumption A3 with γ = 1 (and consequently for all γ ∈ (0, 1]). More precisely, let ϕ 0 be the minimizer of the Hartree functional on the torus with ground state energy E 0 , and let ψ n denote the n'th excited state with energy E n . Then the author proves that there exist constants C, D > 0 such that P ϕ 0 a ψ n 2 ≤ Ce −Da for all (E n − E 0 ) ≤ a ≤ N . As a corollary of this statement, it is shown that there exists C a > 0 such that which implies that assumption A3 is satisfied. Let us conclude the discussion of our assumptions with a remark on the relation between A3 and the so-called Wick property of quasi-free states 4  For a quasi-free state χ on a Fock space F, it is known (e.g. [44,Lemma 5]) that for every a ≥ 1, there exists a constant C a > 0 such that Hence, A3 ( , ) holds with γ = 1 for quasi-free states. Since it is somewhat similar to the Wick property, it is referred to as quasi-free type property in [43].
Finally, let us recall from (14) the Hamiltonian H ϕ(t) (t) introduced in [42,43], which can be understood as first-quantised analogue of a Bogoliubov Hamiltonian. As pointed out in the introduction, H ϕ(t) (t) differs from H β (t) precisely by terms with three or four projectors q ϕ(t) , denoted by C ϕ(t) and Q ϕ(t) . In this sense, it is a quadratic Hamiltonian comparable to H Bog (t). Proof.
) before and after the expression in the brackets and uses the relations which concludes the proof.
The time evolution generated by H ϕ(t) (t) is denoted by U ϕ (t, s), and its well-posedness is recalled in the following lemma.
and generates a unique family of unitary time evolution operators U ϕ (t, s). U ϕ (t, s) is strongly continuous jointly in s, t and leaves H 2 (R dN ) invariant. For an initial datum ψ 0 ∈ L 2 sym (R dN ), the corresponding N -body wave function at time t ∈ R will be denoted by Proof. As a consequence of the Sobolev embedding theorem (e.g. [3, Theorem 4.12, Part IA]), . Hence, the statement of the lemma follows from [24].

Control of higher moments of the number of fluctuations
In our first result, we prove bounds on the growth of the expected number of fluctuations under the time evolution. We consider both the actual N -body dynamics U (t, s) and the dynamics U ϕ (t, s) generated by the Hamiltonian H ϕ(t) (t). The estimates are stated for ( m ϕ ) a ψ 2 as these expressions are required for the proof of our main theorem. By Lemma 2.1, they easily translate to bounds on the corresponding quantities q 1 · · · q a ψ 2 and ξ ϕ , N a ϕ ξ ϕ . The proofs of Proposition 2.4 and Corollary 2.5 are postponed to Section 3.2.
and let ϕ(t) be the solution of (8) with initial datum ϕ(s). Then it holds for t ∈ s, s + T ex d,v,V ext and j ∈ {1, ..., N } that where C t,s j := j! 3 j(j+1) e 9 j t s ϕ(s 1 ) 2 Proposition 2.4 provides an extension to positive β of [42, Lemma 2.1], where a comparable statement was shown for β = 0, γ = 1 and d = 3 with a similar method. Under the additional assumption A3 on the initial data, this implies the following estimates: Corollary 2.5. Assume A1 -A2 and A3 with γ ∈ (0, 1] and A ∈ {1, ..., N }. Let ψ(t), ψ ϕ (t) and ϕ(t) denote the solutions of (4), (26) and (8) with initial data ψ 0 and ϕ 0 from A3. Then for t ∈ 0, T ex d,v,V ext , sufficiently large N and a ∈ {0, ..., A}, it holds that At the threshold γ = 1−dβ, the leading order terms in the sums in Proposition 2.4 change, hence we obtain two different estimates. The additional restrictions on β and γ in part (a) stem from the second sum in Proposition 2.4a. Only if either β < 1 2d or γ > dβ, it is possible to choose b sufficiently large that the first sum dominates for large N .
By Lemma 2.1, Corollary 2.5 yields estimates on the growth of the first A moments of the fluctuation number, given A3 with parameters A and γ. Let Then, for sufficiently large N and for all a ∈ {0, ..., A}, we obtain for β ∈ [0, 1 2d ) where we estimated a, C a , C a 1 for the sake of readability. For β = 0, both time evolutions preserve the property A3 exactly, i.e., with the same power γ of N , up to a constant growing rapidly in t and a. For β > 0, the conservation is exact only for small γ, whereas one looses some power of N for larger γ. Further, note that for the range γ ∈ (0, dβ), we do not obtain a non-trivial estimate for the fluctuations ξ ϕ(t) in U (t, 0)ψ 0 .

Higher order corrections to the norm approximation
Based on the estimates obtained in Proposition 2.4, our main result establishes corrections of any order to the norm approximations (12) and (15): under assumption A3 on the initial data, we construct a sequence {ψ for some δ(β, γ) > 0, which may depend on β as well as on the parameter γ from assumption A3. For reasons given below, our analysis is restricted to the scaling regime β ∈ [0, 1 4d ). As explained in the introduction, it is well known that the actual time evolution ψ(t) is close to the evolution ψ ϕ (t) from (26) in norm. Hence, the first element of the approximating sequence {ψ (a) ϕ } a∈N is determined by Using Duhamel's formula, the difference between U (t, s)ψ and U ϕ (t, s)ψ can be expressed as for any ψ ∈ L 2 (R dN ). Consequently, by the triangle inequality and as a consequence of the unitarity of U (t, s). The leading order contribution in (28) is the term containing C ϕ(s) because the cubic interaction terms are larger than the quartic ones in the following sense: Lemma 2.6. Let ψ ∈ L 2 sym (R dN ) and denote by ϕ(t) the solution of (8) with initial datum The proof of this lemma is postponed to Section 3.3. For j = 0, it gives a bound on the cubic and quartic terms; the more general statement j ≥ 0 is included for later convenience.
When applying Lemma 2.6 to (28), we obtain expressions of the form ( m ϕ(s) ) j U ϕ (s, 0)ψ 0 2 . To be able to use assumption A3 on the initial data, we need to interchange, in a sense, the order of U ϕ (s, 0) and ( m ϕ(s) ) j . This is where Proposition 2.4 comes into play: from part 2.4b, it follows for sufficiently large N that As in Corollary 2.5, the size of γ determines the leading order term in the sum: for γ ≥ 1 − dβ, the dominant contribution issues from n = 3, whereas otherwise the addend corresponding to n = 0 is of leading order. Consequently, To ensure that (29) converges to zero as N → ∞, we have restricted the range of parameters γ admitted by assumption A3 to γ ∈ ( 2+dβ 3 , 1]. Besides, in the first case, the bound is only small for β < 1 4d , and the second case is anyway only possible for β < 1 4d . Hence, we can only cover the parameter regime β ∈ [0, 1 4d ). Analogously to (29), we also obtain Note that β < 1 4d implies that −2 + 6dβ < −1 + 4dβ, and besides, it follows from γ > 2+3d 3 and β < 1 4d that 2 + 2dβ − 4γ < 2 + dβ − 3γ. Consequently, the contribution with C ϕ(s) dominates in (28) for sufficiently large N , which leads to the estimate for some constant c(1) > 0 and with This yields (18) for n = 1.
To construct the second element ψ . As a consequence of Lemma 2.6, we define which equals the leading order contribution in (27) but with the true time evolution U (t, s) replaced by U ϕ (t, s). Put differently, the leading order contribution is cancelled but for the difference between U (t, s) and U ϕ (t, s). Since this difference is evaluated on C ϕ(s) U ϕ (s, 0)ψ 0 , which is small in norm, this is an improvement compared to the first order approximation ψ (1) ϕ (t). To verify this, let us compute the difference between ψ(t) and ψ (2) ϕ (t). Using twice Duhamel's formula, we obtain Due to the unitarity of U (t, s), we obtain with the triangle inequality The leading order term in (33) can be estimated as As before, considering the two ranges of γ separately yields for sufficiently large N (32), where we have used that C t a is increasing in a. Analogously, the second term can be estimated as and the third term was already treated in (30). Combining all bounds, we obtain for some c(2) > 0, which yields (18) for n = 2.
Iterating Duhamel's formula (27) (a − 1) times, we construct ψ (a) ϕ (t) as an expansion with a − 1 terms, where the last term contains the true time evolution U (t, s) and all others exclusively contain U ϕ (t, s). Consequently, to construct ψ (2) ϕ (t), we iterate (27) once more, which yields The leading order contributions issue from the first integral and from the expression with two cubic interaction terms. Analogously to above, they determine the next element ψ and similar calculations as before yield ψ(t) − ψ γ) . Continuing the iteration of (27), we obtain for any a ≥ 1 and s 0 = 0 the expansion All products are to be understood as ordered, i.e. L =0 P := P 0 P 1 · · · P L for L ∈ N and any expressions P . Extracting the leading contributions in each order, we construct the sequence {ψ  Theorem 1. Let β ∈ [0, 1 4d ) and assume A1 -A3 with A ∈ {1, ..., N } and γ ∈ ( 2+dβ 3 , 1]. Let ψ(t) and ϕ(t) denote the solutions of (4) and (8) with initial data ψ 0 and ϕ 0 from A3, respectively, and let ψ Hence, given any desired precision of the approximation, there exists some a ∈ N such that the corresponding function ψ   (t), an a-dependent number of steps is required, as well as the knowledge of the first quantised Bogoliubov time evolution. Put differently, all higher order corrections to the norm approximation follow from the (first order) norm approximation U ϕ (t, 0)ψ 0 after an N -independent number of operations. We cover initial states where the first A moments of the number of fluctuations are sub-leading, where A depends on a but is independent of N .

Preliminaries
where h ϕ(t) j (t) denotes the one-particle operator h ϕ(t) (t) from (8) acting on the j th coordinate.
Proof. For part (a), see, e.g., [51,Lemma 4.1] and note that ϕ(t) L ∞ (R d ) by the Sobolev embedding theorem. Part (b) can be shown as in the proof of [51, Lemma 6.2].

Proof of Proposition 2.4
Proof of Proposition 2.4. The proof of this proposition is essentially an adaptation of the proof of [49,Corollary 4.2]. We begin with part (a). Let ψ ∈ L 2 (R dN ) symmetric, s ∈ R and f : N 0 → R + 0 some weight function. Define and Let us for the moment abbreviate U (t, s)ψ =: ψ t . By Lemma 3.1b, where we have inserted 1 = (p ϕ(t) 2 ) on both sides of the commutator and used Lemma 3.3a. Since q = 0, we conclude that (40) equals zero. From now on, we will for simplicity drop the superscripts ϕ(t). Let Since, for example, f − f −2 By Lemmas 3.1 and 3.2 and since v β 2 N dβ , the first term in (44) leads to To obtain the estimate in the last line, note first that , which leads to by Lemma 3.1. The second term in (44) can be estimated as Now we choose for f the family of weight functions w j λ : k → (w λ (k)) j given by for some 0 < λ ≤ 1 − dβ and j ∈ {0, ..., N }. The set corresponding to L f from (43) is called L w j λ . To bound the operators in L w j λ , note that for any a, b ∈ N 0 , a > b, where we have used in the second line that for every 1 ≤ m ≤ j − 1, and that a j ≥ a − b for any 1 ≤ ≤ j and j ≥ 1 (the statement is trivial for j = 0). Since w λ (k) ≤ k+1 N λ for all k, especially also if k > N λ − 1, we conclude that Besides, one computes analogously to above that (k Finally, w λ (k) = 1 for k > N λ − 1, hence the above estimates imply For all other values of k, the differences yield zero. Thus, every element of L w j λ can be bounded, in the sense of operators, by the operator corresponding to the weight function Besides, since l j λ (k) = 0 for k > N λ + 1, one obtains Inserting (48) to (50) into (45) and (46) and using that λ Now we apply Grönwall's inequality, for now on using the abbreviations α ψ,ϕ,s (w j λ ; t) =: α j (t) and I t := ≤ e j3 j It α j (s) + j3 j e j(3 j +3 j−1 )It I t N dβ−λ α j−1 (s) where we have used that all integrands are non-negative and thus the upper boundary of all integrals could be replaced by t. Written explicitly, this gives where we have estimated I j t e 2j3 j It < e 9 j It . To relate this estimate to m j ψ 2 , observe that for for any b ∈ N. Inserting these estimates into (52) yields To minimise the second term, we choose the maximal λ = 1 − dβ, which concludes the proof of part (a).
The proof of part (b) is much simpler since we now consider the time evolution U ϕ (t, s). The term corresponding to (42) vanishes, which implies that we may directly consider the weights m 2j (k) instead of taking the detour via w j λ (k). Analogously to (38), we define α ψ,ϕ,s (f ; t) := U ϕ (t, s)ψ, f ϕ(t) U ϕ (t, s)ψ .
We will now abbreviate U ϕ (t, s)ψ =: ψ t . In this notation, We now evaluate this expression for the weight m 2j (k), i.e.
Proof of Corollary 2.5. From Proposition 2.4a and the assumptions on the initial data, we conclude that for every b ∈ N and sufficiently large N , If γ ≥ 1 − dβ, the leading order terms in both sums are the ones with maximal n, hence If one chooses b > a 1−dβ 1−2dβ for fixed β < 1 2d , the second term is for sufficiently large N dominated by the first one. For γ < 1 − dβ, the leading order terms are those with n = 0, hence which yields a non-trivial bound only for γ > dβ. Part (b) follows analogously from part (b) of Proposition 2.4 without the restrictions on β and γ that are due to the second sum.

Proof of Theorem 1
Proof of Lemma 2.6. We use the abbreviation (39), and drop all superscripts ϕ(t) in p ϕ(t) , q ϕ(t) and m ϕ(t) for simplicity. By Lemma 3.3a, by (35), and analogously i,j,k,l m a ψ, q i q j q k q l m a ψ N m 4+a ψ 2 , This implies part (a). For part (b), note that by Lemma 3.3a, Consequently, similarly to the estimate of m a Q ϕ(t) ψ . The last inequality follows because by Lemma 3.1a, to Young's inequality and since v β 2 in the sense of operators. As in the estimate of Q ϕ(t) , we thus obtain for ∈ {−1, 1} and analogously q 1 q 2 m a ψ < 4 a N m a+3 ψ 2 and q 1 q 2 q 3 m a ψ < 4 a m a+3 ψ 2 . Together, this implies part (b).