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 Nd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d$$\end{document}-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)-1Ndβv(Nβ·)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(N-1)^{-1}N^{d\beta }v(N^\beta \cdot )$$\end{document} for β∈[0,14d)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \in [0,\frac{1}{4d})$$\end{document}. We derive a sequence of N-body functions which approximate the true many-body dynamics in L2(RdN)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2({\mathbb {R}}^{dN})$$\end{document}-norm to arbitrary precision in powers of N-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^{-1}$$\end{document}. 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 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.
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 d N ) 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 d N ). 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 excitations from the (time-evolved) condensate. The evolution of ϕ(t) approximates the N -body dynamics ψ(t) in the sense of reduced densities (see Sect. 2.1). Moreover, if the dynamics of the excitations are suitably approximated and added to the description, one obtains an approximation of ψ(t) with respect to the L 2 (R d N ) norm, in the sense that ≤ C(t)N −δ for some power δ ∈ (0, 1] depending on the choice of β (see Sect. 2.2).
In this paper, we introduce a novel method for deriving a more precise characterisation of the dynamics. 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 an auxiliary Hamiltonian H ϕ(t) (t) (see (21) for a precise definition). Under the assumption that sufficiently high moments of the number of excitations in the initial state are subleading, we prove higher order corrections to the norm approximation 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 {ψ (a) ϕ (t)} a∈N ⊂ L 2 (R d N ) such that, for sufficiently large N , for some time-dependent constant C(t). The positive exponent δ(β, γ ) is determined in Theorem 1. It depends on β and on a parameter γ which is related to the initial number of excitations (see assumption A3). Let us remark that the approximating functions ψ (a) ϕ (t) are N -body wave functions. Since they are explicitly given in terms of the dynamics U ϕ (t, s) related to the (first order) norm approximation, the functions ψ (a) ϕ (t) are much more accessible than the true dynamics ψ(t).
In particular, all higher order corrections can be obtained from the norm approximation U ϕ (t, 0)ψ 0 by computing an N -independent number of integrals.
Moreover, if the initial excitation vector is quasi-free, one can extend the method introduced in this paper. We conjecture that it is possible to approximate all n-point correlation functions of the full dynamics ψ(t) to arbitrary precision by expressions depending only on approximations of the 2-point correlation functions, whose computation reduces to solving two coupled linear one-body equations ( [28,] and [47,Equation (34)]). This would mean a huge simplification of the full N -body problem (3) corresponding to the Hamiltonian (1) for certain initial states, in particular with regard to a numerical analysis, and we plan to show this in a separate paper.
Finally, we note that higher order approximations of the reduced density matrices were obtained by Paul and Pulvirenti [50] for β = 0 and factorised initial data, based on the method of kinetic errors from the paper by Paul et al. [51]. For j ∈ {1, . . . , N }, the authors of [50] 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 linearisation around this solution.
Due to different methods used, it is not straightforward to compare the results of [50] with the results of this paper. However, we list some features of our paper that differ from the operator-based method of kinetic errors [50,51]. In contrast to the approach in [50], we derive approximations directly for the time-evolved N -body wave function, and our construction 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 include positive values of β and cover more generic initial data than [50], where the initial state is assumed factorised, i.e., with zero initial excitations. 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 ≤ C B, resp. A ≥ C B, and abbreviate Finally, we use the notation for r ∈ R.

Leading Order Approximation: Reduced Densities
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 →∞ Tr γ (1) is the one-particle reduced density matrix of ψ(t) at time t. Then it has been shown, see e.g. [1,2,13,15,20,21,35,57], that lim N →∞ Tr γ (1) 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 Eq. (7) 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 (6), 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 [39,40]. The first derivation of the NLH equation via the BBGKY hierarchy was given by Spohn [58], and this was further pursued, e.g., in [1,2,22,23]. 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 [20,21]. Subsequently, a crucial step of this method was revisited by Klainerman and Machedon [36], based on reformulating combinatorial argument in [20,21] and a viewpoint inspired by methods of non-linear PDEs. This, in turn, motivated many recent works on the derivation of dispersive PDEs, including [13][14][15][16][17]35,57]. In [55], Rodnianski and Schlein introduced yet another method for proving (6), which uses coherent states on Fock space and was inspired by techniques of quantum field theory and the pioneering work of Hepp [32].
In the context of the current paper, the most relevant works on the derivation of the NLH/NLS equation are due to Pickl [53,54], 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 excitations, 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,10,19,33,34,37,43,44].

Next-to-Leading Order: Norm Approximation
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 excitations from 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.
A norm approximation for initial coherent states on Fock space was first obtained in [30,31] by Grillakis, Machedon and Margetis. For initial states ψ 0 ∈ L 2 (R d N ) with fixed particle number, Lewin, Nam and Schlein proved in [41] a norm approximation for β = 0 and V ext = 0 under quite general assumptions on the interaction potential v. Nam and Napiórkowski extended this result in [47] to the range β ∈ [0, 1 3 ), in [49] to the range β ∈ [0, 1 2 ) for the three-dimensional defocusing case, and in [48] to the focusing case in dimensions one and two for β > 0 and β ∈ (0, 1), respectively. As proposed in [42], the authors decomposed the N -body wave function ψ(t) into condensate and excitations as is the truncated bosonic Fock space over the orthogonal complement in L 2 (R d ) of the span of ϕ ∈ L 2 (R d ). A definition of ξ (k) ϕ(t) will be given in (16). Further, ⊗ s denotes the symmetric tensor product, which is for ψ a ∈ L 2 (R da ), ψ b ∈ L 2 (R db ) defined as where S a+b denotes the set of all permutations of a + b elements. The addend k = 0 in (10) describes the condensate, while the terms k ∈ {1, . . . , N } correspond to the excitations. In the following, we will refer to ξ (k) ϕ (t) as k-particle excitation. In [41,[47][48][49], the authors consider initial data of the form for some appropriate initial excitation vector It is then shown that there exist constants C, C > 0 such that where δ = 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 excitations Here, H Bog (t) denotes the Bogoliubov Hamiltonian, 2 an effective Hamiltonian on 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 by Mitrouskas et al. [45,46] via a first quantised approach, where the splitting of ψ(t) into condensate and excitations is realised by means of projections as introduced in [53]. Since we will work in the first quantised setting, let us recall this approach and introduce some notation.
Define the orthogonal projections on L 2 (R d ) and the corresponding projection operators on L 2 (R d N ) 2 Written in second quantized form, H Bog (t) is defined as where a * x and a x denote the operator-valued distributions corresponding to the usual creation and annihilation operators on F (L 2 (R d )). Besides, and P ϕ k = 0 for k < 0 and k > N . Further, for any function f : N 0 → R + 0 and any j ∈ Z, define the operators We will in particular need the operators n ϕ and m ϕ corresponding to the weights Hence, for any ψ ∈ L 2 (R d N ), the part of ψ in the condensate ϕ ⊗N is given by P ϕ 0 ψ, and the part of ψ corresponding to k particles being excited from the condensate is precisely P for some ξ To determine the explicit form of ξ (k) ϕ , observe that by Definition 2.1, where, by definition of the symmetric tensor product, Obviously, ξ (k) ϕ is symmetric under permutations of all of its coordinates, and ξ (k) ϕ is orthogonal to ϕ in every coordinate, i.e.,  (11).
The relation between the N -body state ψ and the corresponding excitation vector ξ ϕ is given by the unitary map where ξ ϕ is defined by (16). 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 (16). The number operator N ϕ on F ≤N ⊥ϕ , counting the number of excitations, is defined by its action The expected number of excitations from the condensate ϕ ⊗N in the state ψ is thus given by with n ϕ from Definition 2.1.
In [46], the authors introduce an auxiliary N -particle Hamiltonian H ϕ(t) (t) by subtracting from H β (t) in each coordinate the mean-field Hamiltonian h ϕ(t) (t) from (7), inserting identities 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) : and with j ) before and after the expression in the brackets and uses the relations which concludes the proof.
The auxiliary Hamiltonian H ϕ(t) (t) has a quadratic structure comparable to that of the , 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 excitations. However, H ϕ(t) (t) is particle number conserving and acts on the N -body Hilbert space L 2 (R d N ), i.e., it determines the evolution of the N -body wave function consisting of excited particles and particles in the condensate ϕ(t), with ϕ(t) the solution of (7). In contrast, H Bog (t) operates on the excitation Fock space F ⊥ϕ(t) , does not conserve the particle number, and exclusively concerns the dynamics of the excitations with respect to the condensate wave function evolving according to (7).
The time evolution generated by H ϕ(t) (t) is denoted by U ϕ (t, s). For an initial datum Existence and uniqueness of of the time evolution U ϕ (t, s) are recalled in Lemma 3.1. Under appropriate assumptions on the initial datum ψ 0 , the time evolution U ϕ (t, s) approximates the actual time evolution U (t, s). More precisely, there exist constants C, C > 0 such that [46, Theorem 2.6]. Further, in the limit N → ∞, the excitations in U ϕ (t, 0)ψ 0 coincide with the solutions of the Bogoliubov evolution equation: let [46, Lemma 2.8]. Hence, the combination of (23) and (24) yields (13), with a different time-dependent constant but the same N -dependence. Finally, let us remark that for larger values of the scaling parameter, beyond the mean field regime, 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 (13) for appropriately modified initial data was proved by Brennecke et al. [11] 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 [8,18,24,25,[28][29][30][31]38,55] for various ranges of the scaling parameter. A related result for Bose gases with large volume and large density was proved in [52].

Assumptions
Let us state our assumptions on the model (1) and on the initial data.
Our analysis is valid as long as the solution ϕ(t) of the non-linear equation (7) exists in and depends on the dimension d, the sign of v ϕ(t) , and the regularity of the external trap V ext (t). Under assumptions A1 and A2 and for times s, t

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
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 [27].
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 d N ).
The third assumption provides a bound on the expected number of excitations from 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 non-negative, 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 3.2, whose proof is postponed to Sect. 4.1.
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 excitations 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 excitations, i.e., the expectation values of (N ϕ(t) /N ) A , vanish as N → ∞. This, 3 Note that the operators n ϕ and m ϕ are equivalent in the sense that they are related via (35) hence all results in terms of m ϕ 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 3.3 easier to write. For example, in terms of n ϕ , Proposition 3.3b 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 m ϕ is more practicable.
in turn, provides a bound on the high components of the excitation 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 excitations 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 d N )-norm due to the correlation structure related to the interactions.
Regarding interacting bosons, A3 is fulfilled for quasi-free states with subleading expected number of excitations, since for any quasi-free state ξ ∈ F and any a ≥ 1 there exists a constant C a > 0 such that [47,Lemma 5]). In [42,Theorem A.1], it was shown that the ground state of H Bog is a quasi-free state, which, via the map U ϕ N , defines an N -body state ψ Bog that converges to the actual ground state ψ 0 in norm as N → ∞ [42, Theorem 2.2]. Note that we require a certain minimal size of γ , which is strictly greater than 2 3 . Further, let us remark that assumption A3 for A = 1 means complete BEC, which was shown to be a sufficient condition for the validity of the Bogoliubov approximation [42]. It was shown in [56, Lemma 1] and [26, Lemma 1] that A3 with A = 1 and γ = 1 is satisfied for low-energy states of a d-dimensional Bose gas for β = 0 in the homogeneous and inhomogeneous setting, respectively. For the Gross-Pitaevskii scaling β = 1 in three dimensions, a comparable bound was proved in [6,7].
Besides, A3 with parameter γ < 1 is comparable to part of the assumption made in [49]. In this work, the authors assume that the initial excitation vector ξ ϕ 0 be a quasi-free state in F ⊥ϕ(t) such that for all ε > 0 and with κ ε > 0 independent of N , and prove a norm approximation for β ∈ [0, 1 2 ) with parameter δ = (1 − ε − 2β)/2. Since ξ ϕ 0 is quasi-free, the first part of this assumption is comparable to A3, whereas we do not impose any condition on the initial energy of the excitations.
Finally, Mitrouskas showed in [45,Chapter 3] that assumption A3 with γ = 1 (and consequently for all γ ∈ (0, 1]) is fulfilled by the ground state and lower excited states of a homogeneous Bose gas on the d-dimensional torus for β = 0. 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 As a corollary of this statement, it is shown that there exists C a > 0 such that which implies that assumption A3 is satisfied.

Control of Higher Moments of the Number of Excitations
In our first result, we estimate the growth of the first A moments of the number of excitations 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,8,12,46,52,55].
and let ϕ(t) be the solution of (7) with initial datum ϕ(s). Then it holds for t ∈ s, s Under the additional assumption A3 on the initial data, this implies that at any time t and for sufficiently large N , the first A moments of the number of excitations remain sub-leading: Corollary 3.4 Assume A1-A2 and A3 with γ ∈ (0, 1] and A ∈ {1, . . . , N }. Let ψ(t), ψ ϕ (t) and ϕ(t) denote the solutions of (4), (22) and (7) with initial data ψ 0 and ϕ 0 from A3. Let Then for t ∈ 0, T ex d,v,V ext , sufficiently large N and a ∈ {0, . . . , A}, it holds that (a) for the time evolution U (t, s) that or, equivalently, that or, equivalently, that with C t a := C t,0 a and where we estimated a, C a , C a 1 for the sake of readability.
At the threshold γ = 1 − dβ, the leading order terms in the sums in Proposition 3.3 change, hence we obtain two different estimates. The additional restrictions on β and γ in part (a) stem from the second sum in Proposition 3.3a. Only if either β < 1 2d or γ > dβ, it is possible to choose b sufficiently large that the first sum dominates for large N . 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 excitations ξ ϕ(t) in U (t, 0)ψ 0 .

Higher Order Corrections to the Norm Approximation
Based on the estimates obtained in Proposition 3.3, our main result establishes corrections of any order to the norm approximations (13) and (23): 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 (22) 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 d N ). Consequently, (27) by the triangle inequality and as a consequence of the unitarity of U (t, s). The leading order contribution in (27) is the term containing C ϕ(s) because the cubic interaction terms are larger than the quartic ones in the following sense: N ) and denote by ϕ(t) the solution of (7) with initial datum Then for any j ∈ N 0 and t ∈ 0, The proof of this lemma is postponed to Sect. 4.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 3.5 to (27), we obtain expressions of the form 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 3.3 comes into play: from part 3.3b, it follows for sufficiently large N that To enhance readability, we do not keep track of the different constants C(t) for now, but we specify it in more detail in Theorem 1. As in Corollary 3.4, 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 (28) converges to zero as N → ∞, we 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 (28), 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 (27) for sufficiently large N , which leads to the estimate This yields (5) for n = 1.
To construct the second element ψ (2) ϕ (t) of the approximating sequence, we need to extract from (26) the relevant contributions such that ψ(t) − ψ (2) . As a consequence of Lemma 3.5, we define which equals the leading order contribution in (26) 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 Duhamel's formula twice, we obtain Due to the unitarity of U (t, s), we obtain with the triangle inequality The leading order term in (32) can be estimated as As before, considering the two ranges of γ separately yields for sufficiently large N with δ(β, γ ) from (31). Analogously, the second term can be estimated as and the third term was already treated in (29). Combining all bounds, we obtain which yields (5) for a = 2. Iterating Duhamel's formula (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 ψ (3) ϕ (t), we iterate (26) 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) − ψ (3) γ ) . Continuing the iteration of (26), 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 {ψ where s 0 := 0. As above, the products are ordered. For n = k = 0, let T  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 (7) 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. We cover initial states where the first A moments of the number of excitations are sub-leading, where A depends on a but is independent of N .

Preliminaries
almost everywhere for some F : where h ϕ(t) j (t) denotes the one-particle operator h ϕ(t) (t) from (7) acting on the j th coordinate.

Lemma 4.2 Let
(c) In particular, this implies Proof For simplicity, let us drop all superscripts ϕ. Part (a) is shown e.g. in [54,Lemma 4.1]. For part (b), observe that for any 1 ≤ j ≤ N , q l ψ = q 1 · · · q j−1 nψ 2 by part (a). Since nψ is again symmetric, the statement follows by iteration.

Lemma 4.3 Denote by T i j an operator acting non-trivially only on coordinates i and j.
(a) Let ϕ ∈ L 2 (R d ), let f , g : N 0 → R + 0 be any weights and i, j ∈ {1, . . . ,

Proof of Lemma 3.2.
Let us for simplicity drop all superscripts ϕ. First, observe that hence in the sense of operators. The first part of (a) follows from Lemma 4.2b and the first line in (34). For the second part, Lemma 4.2a implies Further, note that max j={1,...,a−1} Part (b) follows from (20) and (35).

Proof of Proposition 3.3
Proof of Proposition 3.3. The proof of this proposition is essentially an adaptation of the proof of [52,Corollary 4.2]. We begin with part (a). Let ψ ∈ L 2 (R d N ) 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 4.1a, d dt α ψ,ϕ,s ( f ; t) where we have inserted 1 = ( p 2 ) on both sides of the commutator and used Lemma 4.3a. Since q = 0, we conclude that (39) equals zero. From now on, we will for simplicity drop the superscripts ϕ(t). Let Since, for example, f − f −2 By Lemmas 4.1 and 4.2 and since v β 2 N dβ , the first term in (43) leads to To obtain the estimate in the last line, note first that by Lemma 4.1. The second term in (43) 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 (42) is 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 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 (47) to (49) into (44) and (45) 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 : with where we have estimated I j t e 2 j3 j I t < e 9 j I t . To relate this estimate to m j ψ 2 , observe that for 0 ≤ k ≤ N , for any b ∈ N. Inserting these estimates into (51) 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 (41) vanishes, which implies that we may directly consider the weights m 2 j (k) instead of taking the detour via w j λ (k). Analogously to (37), we define We will now abbreviate U ϕ (t, s)ψ =: ψ t . In this notation, d dt α ψ,ϕ,s ( f ; t) We now evaluate this expression for the weight m 2 j (k), i.e.
This corresponds to w j λ (k) with the choice λ = 1 in (46). Consequently, we define l j (k) := j3 j N −1 m 2( j−1) (k) analogously to (47) and conclude that m 2 j (k) − m 2 j (k − 2) ≤ l j (k) and m 2 j (k + 2) − m 2 j (k) ≤ l j (k). Analogously to the estimate of the first term in (43) and using the relation (48) for λ = 1, we obtain d dt The same Grönwall argument which led to (51) concludes the proof.

Proof of Corollary 3.4.
From Proposition 3.3a 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 3.3 without the restrictions on β and γ that are due to the second sum.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.