Quantum Fluctuations of Many-Body Dynamics around the Gross-Pitaevskii Equation

We consider the evolution of a gas of $N$ bosons in the three-dimensional Gross-Pitaevskii regime (in which particles are initially trapped in a volume of order one and interact through a repulsive potential with scattering length of the order $1/N$). We construct a quasi-free approximation of the many-body dynamics, whose distance to the solution of the Schr\"odinger equation converges to zero, as $N \to \infty$, in the $L^2 (\mathbb{R}^{3N})$-norm. To achieve this goal, we let the Bose-Einstein condensate evolve according to a time-dependent Gross-Pitaevskii equation. After factoring out the microscopic correlation structure, the evolution of the orthogonal excitations of the condensate is governed instead by a Bogoliubov dynamics, with a time-dependent generator quadratic in creation and annihilation operators. As an application, we show a central limit theorem for fluctuations of bounded observables around their expectation with respect to the Gross-Pitaevskii dynamics.


Introduction
In the Gross-Pitaevskii regime, we consider systems of N bosons confined by an external field in a volume of order one (after appropriate choice of the length unit) and interacting through a repulsive potential with small effective range of the order 1/N . The corresponding Hamilton operator is given by with V ext (x) → ∞ as |x| → ∞ and V ≥ 0 compactly supported. According to bosonic statistics, H trap N acts as a self-adjoint operator on L 2 s (R 3N ), the subspace of L 2 (R 3N ) consisting of wave functions that are symmetric with respect to permutations of the N particles. In [45,49] it was proven that, to leading order in N , the ground state energy E trap N of (1.1) satisfies with the Gross-Pitaevskii energy functional Here a > 0 denotes the scattering length of the interaction potential V , which is defined through the solution of the zero-energy scattering equation with the boundary condition f (x) → 1, as |x| → ∞, by requiring that f (x) = 1 − a/|x| outside the support of V . Let ϕ GP ∈ L 2 (R 3 ) denote the (unique, up to a phase) normalized minimizer of (1.3). It turns out that the ground state vector of (1.1) and, in fact, every sequence of approximate ground states, exhibit complete Bose-Einstein condensation in the oneparticle state ϕ GP . In other words, let us consider a normalized sequence ψ N ∈ L 2 s (R 3N ) as N → ∞ (ie. ψ N is a sequence of approximate ground states). Let γ N denote the oneparticle reduced density matrix associated with ψ N , which is defined as the non-negative operator on L 2 (R 3 ) with the integral kernel γ N (x; y) = dx 2 . . . dx N ψ N (x, x 2 , . . . , x N )ψ N (y, x 2 , . . . , x N ) , normalized so that tr γ N = 1. Then, as first proven in [43,44,49], The convergence (1.4) implies that, in the states ψ N , the fraction of particles orthogonal to ϕ GP vanishes, in the limit N → ∞.
Recently, the estimates (1.2), (1.4) have been improved in [48,18,50,19] for integrable interaction potentials, through a rigorous version of Bogoliubov theory. In the translation invariant setting, these improvements have been previously achieved in [8,9]. This approach determines the ground state energy E trap N of (1.1), up to errors vanishing as N → ∞. Moreover, it gives precise information on the low-energy excitation spectrum of (1.1) and on the depletion of the Bose-Einstein condensate (in particular, it shows that the rate of convergence in (1.4) is proportional to 1/N ) and it provides a norm-approximation for the ground state vector. It is interesting to remark that, in the Gross-Pitaevskii regime, such precise estimates on the spectrum of the Hamiltonian cannot be obtained restricting the attention on quasi-free states. Instead, it is important to take into account corrections that can be described by the action of a unitary operator on the Fock space of excitations of the condensate, given by the exponential of a cubic expression in creation and annihilation operators. For related recent results concerning equilibrium properties of Bose gases in the (translation invariant) Gross-Pitaevskii regime, beyond the Gross-Pitaevskii regime and in the thermodynamic limit, see [38,1,15,3,4,37,12,21,22,31,32,5].
Gross-Pitaevskii theory is not only useful to predict the ground state energy of Bose gases described by the Hamilton operator (1.1). It can also be used to approximate their time-evolution. From the point of view of physics, it is relevant to study the dynamics of an equilibrium state of (1.1), after the external fields are switched off (so that the system is no longer at equilibrium). At (or very close to) zero temperature, it is therefore interesting to study the solution of the time-dependent Schrödinger equation for initial data ψ N,0 approximating the ground state of (1.1). In [26,27,28,29], it was first proven that the time-evolution ψ N,t of an initial data ψ N,0 exhibiting Bose-Einstein condensate in a one-particle state ϕ ∈ L 2 (R 3 ) still exhibits Bose-Einstein condensation, in a new one-particle state ϕ t , given by the solution of the nonlinear time-dependent Gross-Pitaevskii equation with the initial data ϕ t=0 = ϕ. More precisely, denoting by γ N,t the one-particle reduced density associated with the solution ψ N,t ∈ L 2 s (R 3N ) of the Schrödinger equation (1.5), it turns out that lim N →∞ ϕ t , γ N,t ϕ t = 1 (1.8) for any fixed t ∈ R, if (1.8) holds true at time t = 0. Analogous stability results have been later established in [51,7]. In [17], the convergence (1.8) is shown to hold with the optimal rate 1/N , for every fixed t ∈ R. It is easy to check that (1.8) implies the convergence γ N,t → |ϕ t ϕ t | in the trace-class topology (and also the convergence of the k-particle reduced density matrix associated with ψ N,t towards the product |ϕ t ϕ t | ⊗k , for every fixed k ∈ N). However, (1.8) does not provide an approximation for the many-body wave function ψ N,t in the strong L 2 (R 3N ) topology. To obtain a normapproximation, it is not enough to approximate the evolution of the condensate. It is instead crucial to take into account the evolution of its orthogonal excitations.
Norm-approximations for many-body dynamics have been derived in the mean-field setting, where particles are initially trapped in a volume of order one and interact weakly through a potential whose range is comparable with the size of the trap (so that every particle interacts effectively with all other particles in the system). In this case, as shown in [39,35,54,36,24,41], the solution of the many-body Schrödinger equation can be approximated, after removing the condensate wave function (whose evolution is described here by the nonlinear Hartree equation, see also [55,30,2,33,34]), by a unitary dynamics on the Fock space of excitations, with a generator quadratic in creation and annihilation operators, acting as a family of time-dependent Bogoliubov transformations (an approximation to arbitrary precision has been recently obtained in [14]). A similar norm-approximation has been derived in [10,16] for the many-body evolution generated by Hamilton operators having the form with β ∈ (0; 1), interpolating between the mean-field and the Gross-Pitaevskii scaling (analogous results have been also obtained in [40,47], for β < 1/2; higher order estimates have been derived, for sufficiently small β > 0, in [13]). To achieve this goal, it was important to combine the unitary dynamics with quadratic generator, describing the evolution of excitations on macroscopic scales, with an additional Bogoliubov transformation generating the correct microscopic correlation structure. In the present paper, we prove a norm-approximation for the many-body dynamics generated by the Hamilton operator (1.6), in the Gross-Pitaevskii regime (ie. for β = 1). Compared with the techniques developed in [10,16] for β < 1, there is an important difference in the construction of the approximating wave function. In fact, as already observed in [8,50,19] in the time-independent setting, for β = 1 microscopic correlations cannot be precisely modelled only through a Bogoliubov transformation; they require instead an additional unitary conjugation with a phase cubic in creation and annihilation operators. This makes our analysis significantly more involved. While the inclusion of the cubic phase is crucial to compare the generators of the full many-body evolution and of the quadratic dynamics (and thus to establish convergence for the corresponding evolutions), at the end it does not substantially change the L 2 (R 3N )-norm of the approximation and it can therefore be removed, providing a quasi-free norm-approximation to the many-body evolution, similar to those obtained in [10,16] for β ∈ (0; 1).

Setting and Main Results
We aim at approximating the solution ψ N,t of the many-body Schrödinger equation (1.5), for a class of initial data exhibiting complete Bose-Einstein condensation in a normalized one-particle wave function ϕ ∈ L 2 (R 3 ). In view of (1.4), from the point of view of physics it is interesting to choose ϕ as the minimizer of a Gross-Pitaevskii functional of the form (1.3). Here, we will keep the choice of ϕ open, requiring only sufficient regularity. First of all, we need to approximate the evolution of the condensate. While (1.7) provides a good approximation at the level of reduced density matrices, to derive a norm-approximation it is more convenient to consider a slightly modified, N -dependent, nonlinear equation. In the modified equation, the interaction potential appearing in the Hamilton operator (1.6) is corrected, to take into account correlations among particles. In order to describe correlations, we fix ℓ ∈ (0; 1/2) and we consider the ground state solution of the Neumann problem on the ball |x| ≤ N ℓ (we omit here the N -dependence in the notation for f ℓ and for λ ℓ ; notice that λ ℓ scales as N −3 ), with the normalization f ℓ (x) = 1 for |x| = N ℓ. We extend f ℓ to R 3 , setting f ℓ (x) = 1 for all |x| > ℓ and we also introduce the notation w ℓ (x) = 1 − f ℓ (x). To describe correlations created by the rescaled interaction appearing in (1.1) and in (1.6), we will use the functions f N, . By scaling, we observe that on the ball |x| ≤ ℓ. Some important properties of λ ℓ , f ℓ , w ℓ are collected in Lemma A.1 in Appendix A. With f N,ℓ , we can now define the condensate wave function at time t ∈ R as the solution ϕ t of the modified Gross-Pitaevskii equation with initial data ϕ t=0 = ϕ. As discussed in Lemma A.1, we have It is therefore easy to check that, as N → ∞, ϕ t converges to the solution of the limiting Gross-Pitaevskii equation (1.7) (with the same initial data). This convergence is part of the statement of Prop. A.2, in Appendix A, where we also collect some standard properties of the solutions of (1.7) and of (2.3) which will be used throughout the paper. As explained in the introduction, to obtain a norm-approximation for the timeevolution it is not enough to approximate the evolution of the condensate; we also need to take into account its excitations. To this end, it is convenient to factor out the (timedependent) condensate wave function introducing, for every t ∈ R, the unitary map N,t }, corresponding to the unique decomposition ⊗sn is symmetric with respect to permutation and orthogonal to ϕ t in each coordinate. Denoting by a(f ), a * (g) the usual creation and annihilation operators, the action of U N,t is characterized (see [42,17]) by the rules 4) where N = a * x a x dx is the number of particles operator, and b * (f ), b(g) are modified creation and annihilation operators satisfying the commutation relations x the corresponding operator valued distributions, we also find for all x, y, z ∈ R 3 . After factoring out the condensate with the unitary operator U N,t , we need to approximate the evolution of its orthogonal excitations in F ≤N ⊥ ϕt . Here, we need to distinguish between microscopic excitations, varying on small length scales between 1/N and ℓ (which is chosen of order one) and macroscopic excitations, varying on scales of order one. Let us first worry about the microscopic excitations, characterising all low-energy states. It is natural to include them on the initial data and to propagate them along the time-evolution. These excitations only depend on time through the time-dependence of the condensate wave function ϕ t . We are going to describe them through a (generalized) Bogoliubov transformation.
We define the integral kernel With Lemma A.1, we find k t ∈ L 2 (R 3 × R 3 ) (with bounded norm, uniform in N ). Hence, k t defines a Hilbert-Schmidt operator on L 2 (R 3 ), which we will again denote by k t . To obtain a Bogoliubov transformation acting on the Hilbert space F ≤N ⊥ ϕt , defined on the orthogonal complement of the condensate wave function ϕ t , we set q t = 1 − | ϕ t ϕ t | and We also denote µ t = η t − k t . With η t , we define the antisymmetric operator and we consider the unitary operator e Bt on F ≤N ⊥ ϕt . An important consequence of the bound η t 2 ≤ k t 2 ≤ C, uniformly in N ∈ N and t ∈ R, is the fact that e Bt does not substantially change the number of excitations. The proof of the following lemma can be found, for example, in [17, Lemma 3.1].
Lemma 2.1. Let B t be the anti-symmetric operator defined in (2.9). Then for every n ∈ Z there exists a constant C > 0 (that depends only on η t ) such that e −Bt (N + 1) n e Bt ≤ C(N + 1) n as an operator inequality on F ≤N .
On states with few excitations (on which the commutation relations (2.5) are almost canonical), e Bt approximately acts as a Bogoliubov transformation. In fact, we can write where we introduced the notation for ♯ = * , ·, for all f ∈ L 2 (R 3 ), x, y ∈ R 3 and n ∈ Z. Some bounds on the operators η t , µ t , γ t , σ t are collected in Lemma A.3, in Appendix A. We are interested in the time-evolution of initial data having the form under appropriate conditions on the excitation vector ξ N ∈ F ≤N ⊥ϕ (we will make assumptions on moments of number of particles and kinetic energy operators, in the state ξ N ). This allows us to consider initial data which are expected to describe the ground state vector of trapped Hamiltonian like (1.1) (see the Remark after Theorem 2.2 for more details).
We consider the many-body evolution ψ N,t = e −iH N t ψ N of (2.13), generated by the translation invariant Hamilton operator (1.6). At time t = 0, we expect ψ N,t to have again approximately the form (2.13), but now with ϕ replaced by the solution ϕ t of (2.3). For this reason, we introduce the excitation vector ξ N,t ∈ F ≤N ⊥ ϕt requiring that In other words, ξ N,t =Ū N (t; 0)ξ N , with the fluctuation dynamics (2.14) While e Bt takes care of the microscopic excitations of the condensate, to derive a norm approximation for e −iH N t ψ N we still need to take into account the evolution of the macroscopic excitations. To reach this goal, we introduce a unitary dynamics U 2,N (t; 0), whose generator is quadratic in creation and annihilation operators (a timedependent Bogoliubov transformation), approximating the fluctuation dynamicsŪ N (t; 0) and providing therefore an approximation of ξ N,t =Ū N (t; 0)ξ N . To this end, let us introduce the "projected" modified creation and annihilation operators where q x (y) = q t (y, x) = δ(x−y)− ϕ t (x) ϕ t (y) is the kernel of the projection orthogonal to the condensate wave function. Then, we define the time-dependent self-adjoint operator (2.17) and Here we used the kernels η t , µ t , γ t , σ t introduced in (2.8), (2.11) and, additionally, we defined p t = γ t − 1 and r t = σ t − η t . Moreover, for s ∈ [0; 1], γ  t are defined as γ t , σ t , but with η t replaced by sη t . Furthermore, we used the notation f x (y) = f t (y, x) for kernels of operators acting on L 2 (R 3 ) (here we drop the label t, to keep the notation as light as possible). Also, we set V N (x) = V (N x) and we used the notation for the kinetic energy operator. The self-adjoint operator J 2,N (t) generates a twoparameter family of unitary transformations U 2,N , satisfying the equation with U 2,N (s; s) = 1 for all s ∈ R (the well-posedness of (2.19) is part of the statement of the next theorem; it will be established in Prop. 4.3). In our first main theorem, we show that the Bogoliubov dynamics U 2,N can be used to describe the evolution of macroscopic excitations of the condensate, providing a normapproximation for the solution of the many-body Schrödinger equation.
be non-negative, spherically symmetric and compactly supported. Let ϕ ∈ H 6 (R 3 ). Let η t be defined as in (2.8), with parameter ℓ > 0 small enough. Then (2.19) defines a unique 2-parameter strongly continuous unitary group Moreover, let B t be defined as in (2.9) and (2.20) Consider a sequence of normalized initial data ψ N ∈ L 2 (R 3 ) ⊗sN , with excitation vectors 22) uniformly in N ∈ N. Then there exist constants C, c > 0 such that for any t ∈ R, and N ∈ N large enough.
Remark. From (2.21), we have ψ N = U * N,0 e B 0 ξ N . Thus, (2.23) is equivalent to with the fluctuation dynamics (2.14). In other words, the proof of Theorem 2.2 reduces to the comparison of the evolutionŪ N with its quadratic approximation U 2,N .
Remark. Observe that the choice ξ N = e B Ω, with and with a kernel τ ∈ (q 0 × q 0 )H 2 (R 3 × R 3 ) (with q 0 = 1 − |ϕ ϕ| projecting orthogonally to the condensate wave function) is compatible with the condition (2.22). This allows us to consider initial many-body wave functions of the form ψ N = U * N,0 e B 0 e B Ω which are expected to approximate (in the L 2 (R 3N ) norm which is preserved over time) the ground state vector of Hamilton operators of the form (1.1), describing trapped Bose gases in the Gross-Pitaevskii regime. This fact has been proven in [8] in the translation invariant setting, for a gas confined on the unit torus. Notice that, in [8, Eq. (6.7)], the norm-approximation of the ground state vector includes also the unitary operator e A , with A cubic in creation and annihilation operators; in fact, this cubic phase can be removed, at the expense of another error of order smaller than N −1/4 (arguing similarly as we do below, in the proof of Theorem 2.2, to show (4.13)). If the trapping potential is sufficiently smooth, it is also easy to verify the condition ϕ ∈ H 6 (R 3 ), with ϕ minimizing the functional (1.3).
In (2.23), after factoring out the evolving Bose-Einstein condensate and the microscopic correlation structure, we approximate the evolution of the macroscopic correlations by the Bogoliubov dynamics U 2,N (t), which still depends on N . It is thus natural to ask whether U 2,N (t) approaches a limiting, N -independent, quadratic evolution U 2,∞ (t), as N tends to infinity. To answer this question, we start by defining the pointwise limit of N w N,ℓ (x), setting 25) and the corresponding limiting integral kernel k t, where ϕ t is the solution to the limiting (N -independent) Gross-Pitaevskii equation (1.7). Similarly to (2.8), we define η t,∞ = (q t ⊗ q t )k t,∞ , projecting along q t = 1 − |ϕ t ϕ t |, orthogonally to ϕ t and µ t,∞ = η t,∞ − k t,∞ . We also introduce the notation With this notation, we can introduce the generator of the limiting quadratic evolution, setting whereâ,â * denote annihilation and creation operators, projected on F ⊥ϕt and where (2.28) Here γ t,∞ are defined like γ t,∞ , σ t,∞ , but with η t,∞ replaced by sη ∞,t . The twoparameter unitary evolution generated by J 2,∞ (t) is denoted by U 2,∞ (t; s). It satisfies the Schrödinger equation with initial condition U 2,∞ (s, s) = 1 for all s ∈ R (the well-posedness of (2.29) is part of the statement of next theorem). Notice that U 2,∞ (t; s) maps the Fock space F ⊥ϕs into F ⊥ϕt ; here, there is no truncation on the number of particles.
we find C, c > 0 (only depending on the expectation (2.30)) such that for any t ∈ R, and N ∈ N large enough. In (2.31), both U 2,N (t; 0)ξ and U 2,∞ (t; 0)ξ are thought of as vectors in the full Fock space F, and . denotes the norm in this space.
Because its generator is quadratic in creation and annihilation operators, the evolution U 2,∞ acts as time-dependent Bogoliubov transformations. Thus, its action on annihilation and creation operators can be calculated explicitly.
It can be written as Remark. Differentiating the action of U 2,∞ on A(f, g) yields the differential equation ]; using the notation introduced in the last two lines of (2.27), the bounds G t,∞ op , H t,∞ 2 < ∞, which are assumed in [52] follow from the analysis in the proof of Prop. 5.2.
Using the approximation in terms of the Bogoliubov dynamics U 2,∞ , we can establish a central limit theorem for the evolution of initial data approximating ground states of the trapped Hamiltonian (1.1).
Theorem 2.5. Under the same assumptions as in Theorem 2.2, consider initial wave functions having the form where B 0 is defined as in (2.9) and (3.2) respectively, and where O (j) = 1 ⊗ · · · ⊗ O ⊗ · · · ⊗ 1 is the operator O acting only on the j-th particle. We set, with U, V indicating the linear maps introduced in (2.32), Then there exists c > 0 and, for all −∞ < a < b < ∞, C > 0 such that, in the state for all N large enough. Here G t is a centered Gaussian with variance f t 2 .
Remark. The bound (2.35) follows through standard arguments (see [20, Proof of Corollary 1.2]) from an estimate of the form for the expectation of the random variable g(O N,t ), valid for any g ∈ L 1 (R), with (1 + s 4 )ĝ(s) ∈ L 1 (R). The proof of (2.36), which is based on the approximation (2.23) of the many body evolution, is given in Section 6. Similarly to [20,53], we could also extend (2.36) to a multivariate central limit theorem, proving that expectations of products of observables of the form (2.34) approach a Gaussian limit, as N → ∞.

Fluctuation Dynamics
While the wave function e −i t 0 κ N (s)ds U * N,t e Bt U 2,N (t)ξ N appearing in (2.23) provides a good norm-approximation for the full evolution e −iH N t ψ N = e −iH N t U * N,0 e B 0 ξ N , the difference of their energy does not converge to zero, as N → ∞. For this reason, it seems difficult to show (2.23) directly. To circumvent this problem, we introduce an alternative approximation for the many-body evolution, this time having the correct energy. At the end, we will show that the two approximations are close in norm.
To define the new approximation of the many-body evolution, we will use a cubic phase. First, we define, as in [50], a cutoff Θ : N → R in the number of particles, setting Throughout this paper we always assume that N is sufficiently large so that this last condition holds true. Next we define the kernel with m = N −α and we introduce the antisymmetric operator where we recall the definition (2.15) of the projected operatorsb Notice that, in contrast with the quadratic kernel (2.7), we cut off (3.1) on a length scale m = N −α , vanishing as N → ∞. This makes sure that, as discussed in Lemma A.5, ν 2 ≤ C √ m ≤ CN −1/4 is small and therefore it allows us to compute the action of the phase e At , which will be used in combination with the generalized Bogoliubov transformation e Bt to generate the microscopic correlation structure, expanding the exponential to first and second order (all higher order contributions will be negligible).
In the following, we will choose α = ε = 1/2. An important observation is that conjugation with e At does not substantially change the number of excitations and their energy. The proof of the next lemma is deferred to Section 7.
Then, for every k ∈ N, there is C > 0 such that holds true as an operator inequality on F ≤N ⊥ ϕt . Remark. In (3.5) it is crucial that the exponent k is the same on both sides of the inequality. In fact, this is the reason for the introduction of the cutoff Θ(N ) in (3.2). For cubic phases without cutoff, as the one used in the time-independent setting in [19], an estimate similar to (3.5) holds true, but only with an additional term (N + 1) k+2 on the r.h.s. of the equation; see [19,Lemma 5.8].
With the cubic phase, we define a new fluctuation dynamics U N , setting To prove Theorem 2.2, we will first show that the difference betweenŪ N (t; 0)ξ N and U N (t; 0)ξ N is small in norm, in the limit N → ∞. Afterwards, we will prove (2.24), but withŪ N replaced by U N . To this end, we will need some properties of the cubically renormalized fluctuation dynamics U N , which we establish in the rest of this section. First of all, we need to control the growth of the number of particles and of the energy along the evolution U N . The proof of the next proposition is based on the estimates in [17, Prop. 6.1] for the dynamicsŪ N .
Under the same assumptions as in Theorem 2.2, let U N be defined as in (3.6). Then there exists C, c > 0 such that for all ξ ∈ F ≤N ⊥ϕ and all t ∈ R. Proof. From Lemma 3.1 (with k = 0), we find From [17, Prop. 6.1 and following remark] and applying again Lemma 3.1, we conclude While controlling the growth of the expectation of N , H N with respect to the fluctuation dynamics U N is enough to establish convergence of the reduced density, to obtain a norm approximation for the dynamics U N we need more precise information on its generator. To this end, we remark that (3.6) satisfies the Schrödinger type equation with the time-dependent generator J N (t) given by In the next proposition, we compute J N (t) up to errors vanishing in the limit N → ∞.
The proof of this proposition is deferred to Section 8.
where κ N (t) and J 2,N (t) are defined as in (2.20), (2.16) and where the error term E J N (t) satisfies

Quadratic evolution and proof of Theorem 2.2
In this section, we study the quadratic evolution defined in (2.19). First of all, we establish important properties of the time-dependent generator J 2,N .
To show Prop. 4.1 (and some of the other bounds discussed in this section), the following lemma will be useful. Its proof is a straighforward adaptation of [10, Lemma 3.4, 3.6].
Lemma 4.2. Let F be a bounded operator and J 1 , J 2 two Hilbert-Schmidt operators on L 2 (R 3 ). We also denote by F, J 1 , J 2 the integral kernels of the three operators ( where # = * , ·. Then, we have for all p ∈ Z and ξ 1 , ξ 2 ∈ F ≤N . Moreover, (4.8) Remark. The bounds in Lemma 4.2 continue to hold true if we replace the operator b, b * (or the corresponding operator-valued distributions) with the projected operators b,b * , introduced in (2.15) and used in (2.16) to define the generator J 2,N (t). In fact, it is easy to see that switching from b, b * tob,b * corresponds to multiplying the operators J 1 , J 2 , F with the orthogonal projection q t = 1 − | ϕ t ϕ t | on the right and/or on the left; this does not increase the norms J 1 2 , J 2 2 , F op appearing on the r.h.s. of (4.7), (4.8). Furthermore, (4.7), (4.8) (and their proof) hold true also for operators A 1 , A 2 , A 3 on the full Fock space F (without truncation to N ≤ N ), defined like A 1 , A 2 , A 3 , but with b, b * replaced by the standard creation and annihilation operators a, a * or by their projected versionâ,â * , used in the definition of the limiting generator (2.27).
Proof of Prop. 4.1. With the notation introduced in the last two lines on the r.h.s. of (2.16), we claim that G t op , H t 2 ≤ Ce c|t| . For most terms this follows easily from Lemma A.1, Prop. A.2 and Lemma A.3. In fact, all contributions (to G t and to H t ) arising from the kinetic part (2.17) that involve derivatives of p, r or µ have L 2 -norm (and therefore also operator norm) less than Ce c|t| . Also the term To handle contributions arising from the potential part (2.18), we observe that (0)), that γ t , p t and σ t are also bounded operators on L 2 (R 3 ) and that all off-diagonal terms (contributing to H) involve at least one factor of p t or σ t (with bounded Hilbert-Schmidt norm, uniformly in N Thus, the L 2 -norm of the l.h.s. is bounded, uniformly in N . This concludes the proof of the bounds G t op , H t 2 ≤ Ce c|t| . From Lemma 4.2 (and from the remark after the lemma) we arrive at (4.1) and, using (4.8), to (4.2) (the second term on the r.h.s. of (2.16), the one proportional to a * x a x , can be handled in the same way). Since we also conclude (4.5).
As for (4.3) and (4.4), we observe that contributions to the time-derivativeJ 2,N (t) have the same form as contributions to J 2,N (t), either with a factor ϕ t replaced by˙ ϕ t , or with one of the kernel η t , γ t , σ t , p t , r t replaced byη t ,γ t ,σ t ,ṗ t ,ṙ t , or with one operator b,b * replaced by its time-derivative (the projection depends on time). Using the similar formula for ∂ tbx and the bounds in Lemma A.1, Prop. A.2 and Lemma A.3, we conclude thatJ 2,N (t) can be written as the sum of terms of the form (4.6) 1 (with two projected operatorsb ♯ or with oneb ♯ and one b ♯ ). For this reason, Lemma 4.2 also implies (4.3) and (4.4) (again, the term proportional to a * x a x on the r.h.s. of (2.16) can be handled similarly).
With the help of Prop. 4.1, we obtain well-posedness of the equation (2.19) (existence and uniqueness of the unitary quadratic evolution U 2,N ) and control on the growth of kinetic energy and number of particles. Compared with the bounds obtained in Prop. 3.2 for the full fluctuation dynamics, here we can derive stronger estimates, controlling arbitrary moments of the number of particles operator N and the second moment of the kinetic energy operator K. These improvements (which we can only show for U 2,N and not for the full fluctuation dynamics U N ) will play a crucial role in the proof of Theorem 2.2. (4.10) Remark. From [10, Lemma 3.10] (extended trivially to the case β = 1), we have which also implies H N ≤ C(K 2 + N 2 ) and Proof. To prove the well-posedness, we proceed as in [41,Theorem 7]. Note that (4.1), (4.3) and (4.5) are precisely the inequalities shown in [41,Lemma 9] and needed to apply the abstract result in [41,Theorem 8].
As for the second bound in (4.10), we first apply (4.2) to estimate where, in the last inequality, we applied the first bound in (4.10). To control the first term on the r.h.s. of the last equation, we observe that (4.4) and with the first bound in (4.10), we conclude that Proof of Theorem 2.2. As observed in the remark after Theorem 2.2, we have Next, we introduce cubic phases to pass fromŪ N to the new fluctuation dynamics U N .
To this end, we estimate Writing we obtain, recalling the definition (3.2) and the estimate ν 2 ≤ C √ m ≤ CN −1/4 from Lemma A.5 and applying Cauchy-Schwarz, With Lemma 3.1, we arrive at Similarly, using also Prop. 4.3, we find Hence, we conclude that (4.14) We now compute 3.3 and with the bound (4.11), we obtain Using V N ≤ CKN (which follows from Sobolev inequality, since Integrating over t, using the assumption (2.22) and combining with (4.14), we find (2.23).

Limiting quadratic evolution and proof of Theorem 2.3
In this section we show the well-posedness of the limiting Schrödinger equation (2.29), we control the growth of the number of particles w.r.t. the limiting quadratic evolution U 2,∞ and we show the convergence of U 2,N to U 2,∞ in the limit N → ∞, as stated in Theorem 2.3.  To prove the convergence of U 2,N towards U 2,∞ , we bound the difference of the two generators. Since J 2,N (t) is only defined on the truncated Fock space F ≤N , our estimate is restricted to this space.
Proposition 5.2. Under the assumptions of Theorem 2.2, we have, for every ξ 1 , ξ 2 ∈ F ≤N and for every t ∈ R, Proof. From (2.16) and (2.27), we write Observe that, with respect to (2.16), we extracted the operator on the second line from theb * (γ x )b(γ y )-contribution in the last summand in (2.18), denoting by G ′ t the difference between G t (as appearing in the last line of (2.16)) and this term (we isolate this term because it will require some additional care). In the expansion for J 2,∞ , we extracted the corresponding term, proportional toâ * xâx , from G t,∞ . Thus, we have where we recall that q t = 1 − | ϕ t ϕ t | and q t = 1 − |ϕ t ϕ t | and where we introduced the notation P t , P t,∞ for the operators with the integral kernels P t (x, For any x ∈ R 3 , we have Estimating Let us now consider the term II. Again with Lemma 4.2, we can estimate Going through the several contributions to G ′ , G ′ ∞ in (2.16) and (2.27) (all diagonal terms) and applying the bounds in Prop. A.2 (part iv)), Lemma A.3 and Lemma A.4, we find that G ′ op ≤ Ce c|t| (as already discussed in the proof of Prop. 4.3) and that In fact, to compare contributions from (2.18) with the corresponding contributions in (2.28), we often need to control the convergence of Here, it is important to observe that, in all terms contributing to G ′ t (which are compared with terms in G ′ t,∞ ), the factor N 3 V (N (x − y))f ℓ (N (x − y)) appears in convolution with a kernel σ x , σ y , p x , p y (this is not the case for III; that's why this term has to be handled separately). To further illustrate this point, consider for example the (this difference arises from the term proportional tob * (γ y )b(σ x ) in the second summand in (2.18)). With (A.2), we can bound With (5.2), we obtain Similarly, we can also estimate To bound the last term, we proceed as in [10,Lemma 5.2]. We find This is the only contribution where the kinetic energy is needed (exactly because, in contrast with contributions in G ′ t , here the difference N 3 V (N (x − y))f ℓ (N (x − y)) − 8πaδ(x − y) acts directly on the operators a * x a y , without convolution; therefore, some regularity of ξ 1 , ξ 2 is needed).
Finally, to control the term IV, we bound Using that q t − q t op ≤ Ce ce c|t| /N , that and that, going through the several contributions to H t , H t,∞ (the off-diagonal terms) in (2.16) and (2.27) and applying the bounds in Prop. A.2 (part iv)), Lemma A.3 and which concludes the proof of the proposition.
We can now proceed with the proof of Theorem 2.3.
Proof of Theorem 2.3. First of all, we observe that While we cannot move J 2,N (t) to the left of the projection ½(N ≤ N ), we can move J 2,∞ (t) to its right, generating a commutator. Thus With Prop. 5.2 and recalling the expression in the last two lines of (2.27) for the limiting generator J 2,∞ (t), we find d dt Integrating over t and with the assumption (2.30), we arrive at Inserting on the r.h.s. of (5.4) proves the desired estimate.
6 Central Limit Theorem: Proof of Theorem 2.5 Following the remark after Theorem 2.5, in this section we aim at proving that for every g ∈ L 1 (R) withĝ ∈ L 1 (R, (1 + s 4 )ds). For the initial wave function ψ N = U * N,0 e B 0 e B Ω with B 0 defined as in (2.9) and with B given by (2.33), with τ ∈ (q 0 ⊗ q 0 )H 2 (R 3 × R 3 ), we find that (2.22) is satisfied, with ξ N = e B Ω. Thus, Theorem 2.2 provides the norm approximation  Next, we conjugate the observable e is O N,t with the unitary operators defining the norm approximation. With the rules (2.4), we find denotes the second quantization of the one-particle operator R and q t = 1 − | ϕ t ϕ t |.
When inserting in (6.2), the contribution of the first term on the r.h.s. of (6.3) is small. Proceeding as in [53, Step 1 in Proof of Theorem 1.1], we arrive at for the action of the modified Weyl operator e iφ(f ) . From (2.10) and (2.12), the action of e Bt is given by and with an error D satisfying Proceeding as in [53, Step 2 in Proof of Theorem 1.1], from (6.4) we therefore arrive at Before proceeding with the last two unitary conjugations, we replace now the field To this end, we observe that for all ξ 1 , ξ 2 ∈ F ≤N (the operators b, b * are only defined on the truncated Fock space).
Thus, for ξ ∈ F ≤N , we have, with the notation ½ ≤N = ½(N ≤ N ), Inserting in (6.6) we find, with Lemma 2.1 and Prop. 4.3, (6.8) Next, we apply Theorem 2.3 to replace the quadratic evolution U 2,N with its limit U 2,∞ . Moreover, we replace h t with where γ t,∞ , σ t,∞ are defined as in (2.26) and q t = 1 − |ϕ t ϕ t |. From Prop. A.2 and Lemma A.4, we find From (6.8), we therefore obtain The action of U 2,∞ (t; 0) on the operators a, a * appearing in φ a is explicit and described by Prop. 2.4. Setting n t = U (t; 0)h ∞,t + V (t; 0)h ∞,t , we find Finally, we need to compute the action of e B . To this end, we first replace the operator B in (2.33) with proceeding as we did above to replace φ(h t ) with φ a (h t ) to show that e B Ω − e Ba Ω ≤ C/ √ N . Then we use the explicit formula for the action of the Bogoliubov transformation e Ba , which implies, setting f t = cosh(τ )n t + sinh(τ )n t , that e Ba Ω, e isφa(gt) e Ba Ω = Ω, e isφa(ft) Ω = Ω, e −s 2 ft 2 /2 e isa * (ft) e isa(ft) Ω = e −s 2 ft 2 /2 . From (6.9), we obtain which immediately implies (6.1). The statement of Theorem 2.5 now follows by standard arguments (see, for example, [20, Corollary 1.2]).

Control of action of A t and proof of Lemma 3.1
In this section, we consider the action of the cubic phase e At on number and energy of excitations. To this end, we compute commutators of A t with the Hamilton operator Lemma 7.1. Recall the definition of A t in (3.2), with parameter M = m −1 = N 1/2 , and recall the notation H N = K + V N , with K, V N the kinetic and potential energy operators Furthermore, for all n ∈ Z.
Remark. To apply Lemma 7.1 in the proof of Lemma 3.1, it is important that the total exponent of H N +N on the r.h.s. of (7.2) is one (so that (3.5) follows by Grönwall's Lemma). To reach this goal, we inserted the cutoff Θ(N ) in the definition (3.2) of A t (in [19,Lemma 5.7], where the cubic phase does not have a cutoff, the exponent of N increases).
Proof. We proceed similarly as in the proof of [19,Lemma 5.7]. We define A 1 t as A t in (3.2), but with b, b * replaced by b, b * . Using [K, a * x ] = −∆ x a * x and that K commutes with N , we have With Lemma A.3 and Lemma A.5, we find Writing ∇ x ν t (x, y) = −N ∇ y w N,m (x − y) ϕ t (y) and integrating by parts, we can bound M 4 can be controlled analogously. Also M 5 satisfies the same bound, since ∇ 1 σ ≤ Ce c|t| √ N by Lemma A.3. As for M 1 , the scattering equation (2.2) yields where we used Hardy's inequality in the first integral and the choice m = N −1/2 . As for M 11 , we write Finally, we have to control the term arising from the difference A t − A 1 t . We observe that, on Here we used the fact that, because of the projection in the kernel η t (x, y) and the definition of p t and σ t in (2.11), we have ϕ t , σ x = 0 = ϕ t , p x . Note that as quadratic forms on F ≤N ⊥ ϕt × F ≤N ⊥ ϕt , we have, for any operator D, with the shorthand notation Applying this with D = K, contributions arising from the commutator of this difference with K can be controlled as before, using the bounds of Prop. A.2 for ϕ t and its derivatives. We conclude that Proceeding as in (7.4) to remove the cutoff Θ(N ), we find Noticing that [a * s a s , b * (σ x )] = σ t (x; s) b * s , and applying Lemma A.3 to prove that |σ t (x; s)| ≤ CN , we find The term N 3 , N 4 can be bounded similarly. As for N 5 , N 6 , they can be estimated by Cauchy-Schwarz, using Lemma A.3 to show that sup z p z , sup z σ z ≤ Ce c|t| . We find we argue as we did above to handle [K, A t − A 1 t ], using the identity (7.5). Combining all the estimates above, we obtain (7.1) (in particular, the large term on the r.h.s. of (7.1) emerges summing the r.h.s. of (7.4) and the r.h.s. of (7.6).
The second claim in the Lemma follows analogously, noticing that powers of (N + 1) can be moved freely from one norm to the other, that in the bounds (7.3), (7.7) we can use the cutoff to remove the operator (N + 1) 1/2 applied to ξ 2 , at the expense of an additional factor M 1/2 ≤ CN 1/4 and that also the main terms on the r.h.s. of (7.4) and (7.6) can be controlled (with Cauchy-Schwarz) by the r.h.s. of (7.2). Now we are ready to show Lemma 3.1.
where, as above, A 1 t is defined as A t , but withb,b * replaced by b, b * and where Let us now consider (3.5). For ξ ∈ F ≤N ⊥ ϕt and s ∈ [0, 1] we define f ξ (s) = ξ, e −sAt H N (N + 1) k e sAt ξ .
We have For k = 0, Lemma 7.1 implies, together with (3.3), that The desired bound follows therefore from Grönwall's Lemma. For k ≥ 1, the contribution of the term proportional to [H N , A] is bounded similarly; using (7.2), with ξ 1 = e sAt ξ, ξ 2 = (N + 1) k e sAt ξ and n = k, we find To handle the second contribution on the r.h.s. of (7.9), we use (7.8). For k = 1, we write Contributions to the second term on the r.h.s. of (7.11) have the form (7.12) The two commutators contributing to A 1,γ t K, A 1,γ t V N can be bounded similarly to the terms M 2 , M 3 and, respectively, N 4 , N 5 in the proof of Lemma 7.1. All other terms can be bounded directly with Cauchy-Schwarz; for instance e sAt ξ, The hermitian conjugates of the terms in (7.12) satisfy similar estimates. We handle the last term on the r.h.s. of (7.11) analogously as we handled the terms proportional to A t − A 1 t in the proof of Lemma 7.1, using the identity (7.5). Observe here that on allowing to recover a commutator which has been estimated before. Thus, we find that (7.10) also holds true for k = 1. For k ≥ 2, we write and we argue similarly as we did for k = 1, after appropriately pulling factors of (N +1) 1/2 through the commutator [N , A t ]. We obtain again the estimate (7.10). The desired bounds follows by Grönwall's Lemma. In this section, we study the properties of the generator J N (t) of the full fluctuation dynamics (3.6). We start from the expression (3.7). A first step in the proof of Prop. 3.3 consists in applying the rules (2.4) to compute the generator where we recall that K and V N are the kinetic and the potential energy operators, as defined on F ≤N ⊥ ϕt in (3.4). Next, we have to consider the effect of the quadratic conjugation with the generalized Bogoliubov transformation e Bt . We define the renormalized generator In order to describe the operator G N (t), we define and the quadratic operator and The following proposition establishes, up to negligible errors, the form of G N (t), in terms of (8.3), (8.5), (8.6).
non-negative, spherically symmetric, and compactly supported. Let G N (t) be defined as in Eq. (8.2). Assume ℓ in (2.7) is small enough but of order one in N . Let where the phase κ G (t), G 2,N (t) and V N are given in Eq. (8.3),(8.4), (8.1) respectively, and the error term E G N (t) satisfies for any ξ 1 , ξ 2 ∈ F ≤N ⊥ ϕt .  [19] for all details (we give some more details for the term (i∂ t e −Bt )e Bt , which is absent in [19]).
As for the cubic term, we find with The proof of (8.14) is very similar to the proof of [19,Lemma 4.9]. Furthermore, following the strategy used in [19,Lemma 4.10], we obtain Finally, we claim that for any ξ 1 , ξ 2 ∈ F ≤N ⊥ ϕt , t ∈ R, and N ∈ N large enough. To prove (8.16), we use (2.10) to expand where the superscript s denotes that we consider cosh and sinh of sη t instead of η t . Here, E ∂t is given by where d x,s , d * x,s are the operator valued distributions associated with operators d (s) t , which are defined as d t in (2.10), but with η t replaced by sη t . With (2.12) and with the bounds from Lemma A.3, we have which concludes the proof of (8.17).
Furthermore, in the terms on the first line on the r.h.s. of (8.13), we split η t = µ t + k t and we combine the contribution associated with k t with the contributions extracted from the third summands on the r.h.s. of (8.12) and on the r.h.s. of (8.15) (expanding γ t = ½ + p t ). Observing that and that, from the scattering equation (2.2), we conclude that To conclude the proof of Prop. 3.3, we need therefore to control the action of A t on the terms on the r.h.s. of (8.8). We already have some information about the action of A t on the Hamilton operator H N = K + V N , thanks to Lemma 7.1. The action of A t on the quadratic terms in G 2,N (t) (excluding the kinetic energy operator K) is determined by the next lemma.

(8.22)
Moreover, assuming additionally that F (r, s) = F (s, r), we also have Finally, we need to control the action of A t on the cubic operator C N,t on the r.h.s. of (8.8). This is the aim of the next lemma. Lemma 8.3. Let A t be defined as in (3.2), with parameters M = m −1 = √ N and C N,t as in (8.7). Furthermore, let We will also use the short-hand notation B A lengthy but straightforward computation leads to Here we used that a(γ x )b * (σ r ) = (1 − N /N ) 1/2 a * (σ x )a y + (1 − N /N )σ t (x, r) and the fact that σ t 2 ≤ C. The other terms with 6 creation and annihilation operators, namely R 8 , R 10 and R 11 , can be bounded in the same way. Next, we consider terms with a contraction, quartic in creation and annihilation operators. Let us start with where we used Lemma A.3 to bound sup r p r 2 ≤ Ce c|t| . The term R 3 is estimated in exactly the same way. Furthermore, and R 7 can be bounded analogously. Also R 9 satisfies the same estimate, since | p r , σ x | ≤ p r σ x . Finally, we control and we observe that R 6 is essentially the same as R 2 , after renaming variables.
Let us now consider the second term on the r.h.s. of (8.26). Here, we have to compute the commutator . The computations are more involved than in (8.27), because now there can be multiple contractions, leading to contributions that are quadratic or even constant in creation and annihilation operators. The main contributions are those where b s and b r are contracted with b * x , b * y . There are two such contributions. Assuming that b, b * satisfy canonical commutation relations (it is easy to check that the corrections are negligible), they are given by with the error | ξ 1 , Eξ 2 | ≤ Ce c|t| N −1/4 (N + 1) 1/2 ξ 1 (N + 1)ξ 2 , needed to remove the cutoff Θ(N ) (arguing similarly as in (7.4), with the choice M = N 1/2 ). Terms involving contractions between b r , b s and B x (or between B * r and b * x , b * y ) are smaller, because they produce factors of σ t (r, x) or σ t (s, x), which are in L 2 . As an example, consider the term δ(s − x)σ t (r, x)B * r b * y (where we contracted b s with b * x and b r with B x (ignoring again corrections to the canonical commutation relations). It produces a contribution S to (8.26), which can be bounded by where we used |σ t (r, x)| ≤ CN | ϕ t (r) ϕ t (x)| from Lemma A.3 and ν t ≤ C √ m from Lemma A.5 (and the choice m = N −1/2 ). Terms arising from [B * r b s b r , b * x b * y B x ] containing 6 or 4 creation and annihilation operators can be bounded as we did above with R 1 , . . . , R 11 . The next term in (8.26) has the form The contribution T 2 is already normally ordered. It can simply be bounded by Cauchy-Schwarz. We find where we used sup x σ x 2 ≤ Ce c|t| from Lemma A.3 and sup x ν x 2 ≤ C √ me c|t| from Lemma A.5 (with m = N −1/2 ). All other normally ordered term emerging from (8.28) can be treated analogously. On the other hand, terms involving commutators (produced through normal ordering) are closely related with the contributions discussed above from the first two terms on the r.h.s. of (8.26). Due to the presence of the differences Θ(N + 1) − Θ(N ) (or similar), also the contributions where b s , b r are contracted with b * x b * y (arising from T 3 and T 4 ) are negligible, here (since (Θ(N + 1) − Θ(N ))ξ ≤ C/M ½(M/2 ≤ N ≤ M )ξ , we can gain a factor M −1 , arguing similarly as in (7.4)).
Finally, we deal with the commutator [C N,t , A t − A 1 t ] using the identity (7.5). The resulting terms can be treated analogously as we did with the contributions to [C N,t , A 1 t ] (but these terms are less singular and thus simpler to handle). They all satisfy the estimate (8.25). We skip the details.
We are now ready to proceed with the proof of Prop. 3.3. where, by Prop. 8.1,

Proof of
with the eror E N (t) satisfying the bound (8.9). From Lemma 3.1, we find With Lemma 7.1, Lemma 8.3 and Lemma 8.2, we claim that where Ξ 1 , Ξ 2 are defined in Eq. (8.24) and To prove (8.30), we start by observing that Going through the terms in On the other hand, with Duhamel's formula, we can write Similarly, Applying Lemma 7.1, Lemma 8.3 and Lemma 8.2 (noticing that the commutator of the quadratic operators Ξ 1 , Ξ 2 with A t is a sum of terms that can be bounded with (8.22), (8.23)) and propagating the estimates through the cubic phase with Lemma 3.1, we arrive at (8.30).
As for the first term on the r.h.s. of (8.29), we observe that, since sup x ν t,x ≤ Ce c|t| N −1/4 from Lemma A.5, We conclude that where | ξ 1 , E ′ ξ 2 | ≤ Ce c|t| N −1/4 (H N + N + 1) 1/2 ξ 1 (H N + N 3 + 1) 1/2 (N + 1)ξ 2 . (8.31) Next, we observe that in the terms Ξ 1 , Ξ 2 , as defined in (8.24), we can replace the parameter m = N −1/2 with the fixed, N -independent, parameter ℓ ∈ (0; 1), at the expense of small error. In fact, setting Ce c|t| √ N (N + 1) 1/2 ξ 1 (K + N + 1) 1/2 ξ 2 (8.33) for j = 1, 2. To show (8.33), we argue as we do in the proof of Prop. 5.2 to control the terms I, II, III (in that Proposition, we control convergence of N 3 V N f N,ℓ towards a δ-distribution, but the same argument implies convergence of N 3 V N f N,m towards a δ-distribution and therefore allows us to control the difference N 3 V N (f N,ℓ − f N,m ) = N 3 V N (w N,m −w N,ℓ ); note that, to handle N 3 V N f N,m we need to use (A.2) with ℓ replaced by m = N −1/2 , which makes some of the estimates, like the one for the second term on the r.h.s. of (5.3), worse). Combining the operator G 2,N (t), as defined in (8.4), with the terms Ξ ′ 1 , Ξ ′ 2 from (8.32) we obtain, as quadratic form on F ≤N ⊥ ϕt (so that the "projected" operatorsb,b * are the same as b, b * ), the operator J 2,N (t) in (2.16), up to a small error due to the term on the seventh line of (2.16), whose matrix elements can be bounded by A Properties of f ℓ , ϕ t , ϕ t , η t , η ∞,t , ν t In this appendix, we collect some analytic properties of functions and kernels that are used throughout the paper to construct the approximation of the many-body dynamics.
In the first lemma, whose proof can be found in [17,9,25], we consider the ground state solution of the Neumann problem (2.1) on the ball |x| ≤ N ℓ, with the normalization f ℓ (x) = 1 for |x| = N ℓ. ii) We have 0 ≤ f ℓ , w ℓ ≤ 1. Moreover there exists a constant C > 0 such that for all ℓ ∈ (0; 1/2) and N ∈ N.
iii) There exists a constant C > 0 such that for all x ∈ R 3 , ℓ ∈ (0; 1/2) and all N ∈ N large enough.
i) For ϕ ∈ H 1 (R 3 ), there exist unique global solutions t → ϕ t and t → ϕ t in C(R, H 1 (R 3 )) of the Gross-Pitaevskii equation (1.7) and, respectively, of the modified Gross-Pitaevskii equation (2.3) with initial datum ϕ. We have ϕ t = ϕ t = 1 for all t ∈ R. Furthermore, there exists a constant C > 0 such that ii) If ϕ ∈ H m (R 3 ) for some m ≥ 2, then ϕ t ,φ t ∈ H m (R 3 ) for every t ∈ R. Moreover there exist constants C depending on m and on ϕ H m , and c > 0 depending on m and on ϕ H 1 such that for all t ∈ R ϕ t H m , φ t H m ≤ Ce c|t| .
iii) Suppose ϕ ∈ H 4 (R 3 ). Then there exist constants C > 0 depending on ϕ H 4 , and c > 0 depending on ϕ H 1 such that for all t ∈ R φ t H 2 , φ t H 2 ≤ Ce c|t| .
Furthermore, if ϕ ∈ H 6 (R 3 ) there exist constants C > 0 depending on ϕ H 6 , and c > 0 depending on ϕ H 1 such that for all t ∈ R φ t H 4 , ˙ ϕ t H 4 ≤ Ce c|t| .
For ϕ ∈ H 6 (R 3 ) there are constants C, c > 0 such that Recall now the definition (2.8), depending on the parameters N, ℓ, of the kernel η t appearing in the generalized Bogoliubov transformation e Bt and the notation γ t = cosh η t , σ t = sinh η t . Furthermore, we set p t = γ t − ½, r t = σ t − η t and µ t = η t − k t (recall (2.7)). Several bounds for the operators η t , γ t , σ t , p t , r t (for their integral kernels) and for their time-derivatives are established in the next lemma, whose proof is a straightforward adaptation of [7,  Lemma A.3. Let ϕ t be the solution of (2.3) with initial datum ϕ ∈ H 4 (R). Let w ℓ = 1−f ℓ with f ℓ the ground state solution of the Neumann problem (2.1) and let ℓ ∈ (0; 1/2). Let k t , η t , µ t be defined as in (2.7), (2.8). Then there exist constants C, c > 0 depending only on ϕ H 4 (or lower Sobolev norms of ϕ) and on V such that the following bounds hold uniformly in ℓ, for all t ∈ R. i) We have η t ≤ Cℓ 1/2 and also ∇ j η t ≤ C √ N , ∇ j µ t ≤ C for j = 1, 2. With ∇ 1 η t and ∇ 2 η t we indicate the kernels ∇ x η t (x; y) and ∇ y η t (x; y), similar definitions hold for ∆ j η t , for j = 1, 2. Let σ t , p t , r t be defined as in (2.11) and after (2.18), we obtain σ t , p t , r t , ∇ j p t , ∇ j r t ≤ C , ∆ j p t , ∆ j r t , ∆ j µ t ≤ Ce c|t| , ∇ j σ t ≤ Ce c|t| √ N .
iii) Moreover we have where we indicate with sup x η x 2 = sup x |η t (x; y)| 2 dy and σ t,x , p t,x , r t,x ≤ Ce c|t| .
iv) For j = 1, 2 we have the following bounds for the time derivatives ∂ t η t , ∂ 2 t η t ≤ Ce c|t| , and also ∂ t ∇ j η t ≤ C √ N e c|t| , ∂ t ∇ j µ t ≤ Ce c|t| .
Furthermore ∂ t σ t , ∂ t r t , ∂ t p t , ∇ j ∂ t p t , ∇ j ∂ t r t , ∆ j ∂ t p t , ∆ j ∂ t r t , ∆ j ∂ t µ t ≤ Ce c|t| .
vi) Finally, we have ∂ t η x , ∂ t k t,x , ∂ t µ t,x ≤ Ce c|t| and ∂ t σ t,x , ∂ t p t,x , ∂ t r t,x ≤ Ce c|t| .
While the kernels η t , γ t , σ t , p t , r t considered in the last lemma are used in the definition of the fluctuation dynamics U N and of its quadratic approximation U 2,N , the limiting quadratic evolution U 2,∞ is defined in (2.29) in terms of limiting kernels η ∞,t , γ ∞,t , σ ∞,t , p ∞,t , r ∞,t . To show the well-posedness of U 2,∞ and to compare it with U 2,N , we need some bound on these limiting objects.