Weak Edgeworth expansion for the mean-field interacting Bose gas

We consider the ground state and the low-energy excited states of a system of $N$ identical bosons with interactions in the mean-field scaling regime. For the ground state, we derive an Edgeworth expansion for the fluctuations of bounded one-body operators, which yields corrections to a central limit theorem to any order in $1/\sqrt{N}$ . For suitable excited states, we show that the limiting distribution is a polynomial times a normal distribution, and that higher order corrections are given by an Edgeworth-type expansion.


Introduction
A quantum mechanical system of N identical bosons is described by a wave function Ψ that is square integrable and symmetric under the exchange of any two particles, i.e., for d ≥ 1 the spatial dimension of the system and where ⊗ sym denotes the symmetric tensor product. We study the statistics of measurements described by self-adjoint operators on H N . In particular, we consider one-body operators on H N , i.e., operators of the form for bounded self-adjoint operators B on H. Since we consider indistinguishable bosons, we study symmetrized operators, i.e., operators of the form N j=1 B j . An example is the number of particles in a bounded volume V ⊂ R d , described by the operator N j=1 χ V (x j ) , (1.4) where χ V denotes the characteristic function on V . The goal of this article is to better understand the statistics of such operators.
Due to the permutation symmetry (1.1), the family of one-body operators {B j } N j=1 defines a family of identically distributed random variables, whose distribution is determined by the wave function Ψ via the spectral theorem. The probability that the corresponding random variable B j takes values in a set A ⊂ R is given by where χ A denotes the characteristic function of the set A and where ·, · denotes the inner product of H N . Functions of self-adjoint operators are defined via the functional calculus. Note that the operators j B j are formally the analogue of sample averages, which, in probability theory, are often interpreted as repeated measurements. This interpretation does not apply in our setting: the operator j B j does not describe N single-particle measurements on N copies of the system (these measurements would always be independent of each other).
If the N -body wave function is a product state, i.e., if Ψ = ϕ ⊗N for some ϕ ∈ H, the random variables B j are independent and identically distributed (i.i.d.). Consequently, N −1 j B j satisfies the law of large numbers (LLN), and the fluctuations around the expectation value are, in the limit N → ∞, described by the central limit theorem (CLT). Moreover, for large but finite N , the fluctuations can be expanded in an asymptotic Edgeworth series, providing higher order corrections to the central limit theorem to any order in 1/ √ N (see Section 3.5 for a more detailed discussion).
A factorized wave function Ψ = ϕ ⊗N describes the ground state of an ideal Bose gas, i.e., a system without interactions between the particles. In this work, we are interested in the situation where the bosons interact weakly with each other. We consider a system of N bosons in R d described by the many-body Hamiltonian acting on H N sym , under suitable assumptions on the interaction v and the external trapping potential V (see Section 1.1). This describes a Bose gas in the so-called mean-field (or Hartree) regime, where the interactions are weak and long-ranged. We consider the ground state Ψ gs N and suitable low-energy excited states Ψ ex N of the Hamiltonian H N , i.e., and where E gs N := inf spec(H N ) is the ground state energy and E ex N denotes a suitable excited eigenvalue of H N (see Definition 2.1). Due to the interactions between the particles, these states are no product states but correlated. Consequently, the family {B j } j of one-body operators defines a family of (weakly) dependent random variables. In fact, one can deduce from [4] that their covariance is where E Ψ N [·] := Ψ N , · Ψ N . Despite this dependence, the family {B j } j satisfies a LLN, which is comparable to the situation of i.i.d. random variables (see Section 3.2). Moreover, one can prove a CLT (see, e.g., [1,6,29,28]), which is a result of the formal analogy of quasi-free states and Gaussian random variables. Due to the dependence of the random variables {B j }, the variance of the limiting Gaussian in the CLT is not given by Var ϕ [B] but differs by O(1) (see Section 3.3). 1 In this work, we prove that the statistics of bounded one-body operators with respect to the N -body ground state Ψ gs N admit a weak Edgeworth expansion, which differs from the expansion for the i.i.d. case due to the interactions. Moreover, we prove an Edgeworth-type expansion for a class of low-energy excited states Ψ ex N .

Assumptions
It is well known that the ground state Ψ gs N as well as the low-energy excited states Ψ ex N of H N exhibit Bose-Einstein condensation (BEC), i.e., lim N →∞ for any k ≥ 0. Here, |ϕ ϕ| denotes the projector onto ϕ ∈ H, i.e., the operator with integral kernel ϕ(x)ϕ(y), and γ (k) N denotes the k-particle reduced density matrix of Ψ N ∈ {Ψ gs N , Ψ ex N }, whose integral kernel is defined as (1.11) The condensate wave function ϕ is given by the minimizer of the Hartree energy functional, The corresponding Hartree energy is denoted by (1.14) We make the following assumptions on the interaction potential v and the trap V , which, in particular, ensure that ϕ is unique and can be chosen real-valued: be measurable, locally bounded and non-negative and let V (x) tend to infinity as |x| → ∞, i.e., and v ≡ 0, and assume that there exists a constant C > 0 such that, in the sense of operators on Q(−∆) = H 1 (R d ), Besides, assume that v is of positive type, i.e., that it has a non-negative Fourier transform.
Assumption 3. Assume that there exist constants C 1 ≥ 0 and 0 < C 2 ≤ 1, as well as a function ε : in the sense of operators on D(H N ).
Assumption 1 ensures that V is a confining potential; an example is the harmonic oscillator potential, V (x) = x 2 . Assumption 3 ensures that low-energy eigenstates of H N exhibit complete BEC in the Hartree minimizer, with a sufficiently strong rate. Assumptions 2 and 3 are, for example, satisfied by any bounded and positive definite interaction potential v, and by the repulsive three-dimensional Coulomb potential, v(x) = 1/|x|. Assumptions 1 to 3 are precisely the assumptions made in [4]. They ensure that we can expand the low-energy eigenstates of H N and the corresponding energies in an asymptotic series in 1/ √ N (see Section 2.3), which is crucial for deriving the Edgeworth expansions. Our main result holds for the ground state Ψ gs N of H N and for a class of excited eigenstates

Main Result
We are interested in the statistics of the symmetrized operators j B j . After centering around the expectation value, we rescale by dividing by √ N . This scaling is chosen as it is the size of the standard deviation of j B j , which follows from (1.9) and (1.10) because This leads to the random variable for self-adjoint B ∈ L(H), where E Ψ N denotes the expectation value of a random variable with respect to the probability distribution determined by Ψ N . Moreover, we consider operators B such that the Hartree minimizer ϕ is not an eigenstate of B. This is equivalent to the statement that the standard deviation σ of the limiting Gaussian in the CLT (see our theorem below) is nonzero, see (3.12). Our main result is the following: Theorem 1. Let Assumptions 1 to 3 hold and let Ψ N ∈ C (η) N for some η ∈ N 0 , with C (η) N as in Definition 2.1. Let a ∈ N 0 and g ∈ L 1 (R) such that its Fourier transform g ∈ L 1 (R, (1 + |s| 3a+4 ). Then, for any self-adjoint bounded operator B ∈ L(H) such that the Hartree minimizer ϕ is not an eigenstate of B, for σ as in (3.12). Here, the functions p j (x) are polynomials of finite degree with real coefficients depending on B, V and v. The error can be estimated as for some C(a) > 0, where · op denotes the operator norm on L(H).
N , then p gs j is a polynomial of degree 3j which is even/odd for j even/odd. In particular, with α 3 as in (4.25) and where H 3 is the third Hermite polynomial (see (3.26)).
N for some η > 0, then p ex j is a polynomial of degree 3j + 2η which is even/odd for j even/odd. The leading order p ex 0 is computed in Proposition 4.7.
In particular, our result does not imply an asymptotic expansion of the probability If the N -body system is in its ground state Ψ gs N , Theorem 1 implies that B N admits a weak Edgeworth expansion although the random variables are not independent. However, the interactions affect the precise form of the Edgeworth series: the standard deviation σ of the Gaussian as well as the polynomials p gs j differ from the expansion for the non-interacting Bose gas (see Sections 3.5 and 3.6 for a detailed discussion). To prove Theorem 1, we expand the characteristic function φ gs N (s) := Ψ gs N , e isB N Ψ gs N in powers of N −1/2 . To leading order, φ gs N (s) is given by the expectation value of a Weyl operator with respect to a quasi-free state . Quasi-free states satisfy a Wick rule comparable to Wick's probability theorem for Gaussian random variables, and this formal analogy is the reason why we obtain a CLT for the ground state. Technically, we use an equivalent formulation of Wick's rule, namely the fact that a quasi-free state is a Bogoliubov transformation of the vacuum. This allows us to reduce the computation of φ gs N (s) to the computation of vacuum expectation values, which are non-zero only if they contain equal numbers of creation and annihilation operators.
For low-energy excited states, the leading order of the corresponding characteristic function φ ex N (s) is no longer given by an expectation value with respect to a quasi-free state, but rather a state with a finite number of creation/annihilation operators acting on a quasi-free state. Consequently, the limiting distribution is not a Gaussian but a Gaussian multiplied with a polynomial. One still obtains an Edgeworth-type expansion, but each order of the distribution is now the Gaussian times a (different) polynomial.
Theorem 1 is, to the best of our knowledge, the first derivation of an Edgeworth expansion for an interacting quantum many-body system. Asymptotic expansions for (weakly) dependent random variables have been derived in [23,24,18] for Markov processes, in [14] for stochastic processes which are approximated by a suitable Markov process, and in [7] in the context of dynamical systems. In [11], the authors prove the existence of Edgeworth expansions for weakly dependent random variables under fairly generic conditions, which includes random variables arising from dynamical systems and Markov chains but excludes our model 2 .
As discussed in Section 3.5 for the i.i.d. situation, Theorem 1 yields a very precise description in the center of the distribution. In contrast, it does not generally provide a good approximation of the tails of the distribution. For the dynamics generated by H N , large deviation estimates have been proven in [20,30].
We expect that Theorem 1 can be generalized to all situations where the N -body wave function admits an (explicitly known) asymptotic expansion in the spirit of Lemma 2.2. For example, it seems obvious that a dynamical Edgeworth expansion should exist, which provides corrections to [1]; moreover, generalizations to k-body operators as in [28] and to k one-body operators as in [6] seem feasible.
The remainder of the article is structured as follows: In Section 2, we summarize the quantum many-body framework and collect known results for the mean-field Bose gas which we require for the proof. Section 3 is a review of the probabilistic picture, including existing results on the CLT for the interacting Bose gas. In particular, we analyse the effect of the interactions on the Edgeworth series (Sections 3.5 and 3.6). Finally, Section 4 contains the proof of Theorem 1.

Many-body framework 2.1 Excitations from the condensate
We consider N -body states Ψ which exhibit complete BEC in the Hartree minimizer ϕ in the sense of (1.10). However, this does in general not imply that Ψ = ϕ ⊗N ; in fact, an O(1) fraction of the particles forms excitations from the condensate. To describe them mathematically, one recalls, e.g. from [22], that any Ψ ∈ H N sym can be decomposed as with the usual notation {ϕ} ⊥ := {φ ∈ H : φ, ϕ = 0}. The sequence of k-particle excitations forms a vector in the (truncated) excitation Fock space over {ϕ} ⊥ , and vectors in F ⊥ (resp. F ≤N ⊥ ) are denoted as φ (resp. φ ≤N ). The creation and annihilation operators on F ⊥ , a † (f ) and a(f ) for f ∈ {ϕ} ⊥ , are defined in the usual way as In [11], the authors consider a Banach space B and assume that the characteristic function is of the form φN (s) = (L N s v), where Ls : B → B is a family of bounded linear operators and where v ∈ B, ∈ B . Applied to our setting, we would identify v with the ground state ΨN , and with the projection on the ground state. However, e isB N is not of the form L N s for some N -independent Ls. Even if we would introduce Ls = e is 1 N B N , this operator would not satisfy the assumptions made in [11], which include that the spectrum of Ls is contained in the open disc of radius 1 for all s = 0, and that L N for k ≥ 1 and k ≥ 0, respectively, and φ ∈ F ⊥ . They can be expressed in terms of the operator-valued distributions a † x and a x , which satisfy the canonical commutation relations We denote the second quantization in F ⊥ (resp. F) of an operator A by dΓ ⊥ (A) (resp. dΓ(A)). The vacuum is denoted by |Ω and the number operator on F ⊥ is given by An N -body state Ψ is mapped onto its corresponding excitation vector χ by the unitary as identities on F ≤N ⊥ . We extend U N,ϕ trivially to a map to the full space F ⊥ . Analogously, elements of F ≤N ⊥ are naturally understood as elements of F ⊥ .

Bogoliubov theory
It was shown in [4] that the low-energy eigenstates of H N can be retrieved by perturbation theory around the eigenstates of the Bogoliubov Hamiltonian, which is given by (2.10) Here, and where we used the orthogonal projectors p := |ϕ ϕ| , onto the condensate and its complement.

Bogoliubov transformations
The Bogoliubov Hamiltonian H 0 can be diagonalized by Bogoliubov transformations (see, e.g., [32]), which are defined as follows: For F = f ⊕ g ∈ {ϕ} ⊥ ⊕ {ϕ} ⊥ , one defines the generalized creation and annihilation operators A(F ) and A † (F ) as where is called a (bosonic) Bogoliubov map. It can be written in the block form where U and V denote the operators with integral kernels U (x, y) and V (x, y), respectively.
The identity (2.14) leads to a transformation rule of creation/annihilation operators under a Bogoliubov transformation, In particular, powers of N ⊥ conjugated with U V can be bound as

Quasi-free states
Finally, we recall that a normalized state φ ∈ F ⊥ is called quasi-free if there exists a Bogoliubov for a ∈ {a † , a}, n ∈ N and f 1 , ..., f 2n ∈ {ϕ} ⊥ . Here, P 2n denotes the set of pairings where S 2n denotes the symmetric group on the set {1, 2, ..., 2n}.

Eigenstates of H 0
We denote by U V 0 : F ⊥ → F ⊥ the Bogoliubov transformation that diagonalizes H 0 , i.e., where D > 0 is a self-adjoint operator on {ϕ} ⊥ . It admits a complete set of normalized eigenfunctions, denoted as {ξ j } j≥0 . The ground state χ gs 0 of H 0 is unique and given by Any non-degenerate excited eigenstate χ ex 0 of H 0 can be expressed as

Low-energy eigenstates of H N
Assumptions 1 to 3 ensure that H N has a unique ground state and a discrete low-energy spectrum. We will consider the following class of eigenstates of H N : N iff all of the following are satisfied: (c) χ 0 is a state with η quasi-particles, i.e., it is given by (2.25) with η 0 + η 1 + · · · + η k = η.
In particular, i.e., the ground state is contained in the set C (η) N with zero quasi-particles. To keep the notation simple, we will indicate the quasi-particle number η only when it is inevitable to avoid ambiguities. If Ψ N ∈ C (η) N for some η ∈ N 0 , it was shown in [4, Theorem 3] that χ = U N,ϕ Ψ admits an asymptotic expansion in the parameter (N − 1) −1/2 , namely for some constant C(a) > 0 and for coefficients χ ∈ F ⊥ given in [4, Theorem 3, Eqn. (3.26)].
For the proof of Theorem 1, it is more convenient to have a full expansion of these states in powers of N −1/2 instead of (N − 1) −1/2 , which can be deduced from the results in [4] in a straightforward way.
N for some η ∈ N 0 and denote the corresponding excitation vector by χ = U N,ϕ Ψ.
(a) For any a ∈ N 0 , there exists a constant C(a) > 0 such that for some functions Θ (η) ,j ∈ L 2 sym (R dj ).
(b) For any , b ∈ N, there exists a constant C( , b) such that For any a ∈ N 0 , there exists some constant C(a) > 0 such that

32)
where the coefficients can be bounded as for some constants C( ) > 0. In particular, B (0) = ϕ, Bϕ , and The functions Θ (η) ,j can be computed using perturbation theory, and we refer to [4] for the explicit expressions. In a similar way, one obtains explicit expressions for B ( ) ; see [2].

Probabilistic picture
To illustrate the effect of the interactions, we compare in this section the random variables with probability distribution determined by Ψ N ∈ C To underline differencences between the ground state Ψ gs N ∈ C N and excited states Ψ ex N ∈ C (η) N for η > 0, we will indicate this in the notation by using the superscripts gs and ex when appropriate.

Random variables
A self-adjoint one-body operator B ∈ L(H) defines a family {B j } N j=1 of random variables with common probability distribution determined by the N -body wave function Ψ N . For Ψ iid N , the random variables are i.i.d., and the expectation value E ϕ [B], the variance Var ϕ [B] and the standard deviation σ iid are given by For an eigenstate Ψ N ∈ C due to the bosonic symmetry (1.1) of Ψ N . Note that by (1.10),

Law of large numbers
For the product state Ψ iid N , the weak LLN states that the empiric mean converges to its expectation value, i.e., lim N →∞ for any ε > 0. Abbreviating B := B − ϕ, Bϕ , Markov's inequality yields for the interacting gas (see, e.g., [1, Sec. 1]) hence (1.10) yields The LLN for Ψ N looks formally like the LLN for independent random variables. Let us stress that Ψ iid N is not the ground state of the ideal gas because ϕ is the minimizer of the Hartree energy functional, which depends on the interactions. In this sense, the interactions have an effect already on the level of the LLN.

Central limit theorem for the ground state
Let us first compare the ground state Ψ gs N of the interacting gas with the product state Ψ iid N . The fluctuations around the respective expectation values are described by the rescaled and centred random variables for U 0 and V 0 from (2.26) and q as in (2.13). In general, σ and σ N differ by an error of order O(1). Hence, the interactions have a visible effect on the level of the CLT: they change the variance of the limiting Gaussian random variable. The simplest way to understand this effect is via the characteristic functions of the random variables B iid N and B N , which are given by for the ideal gas, and by φ gs N (s) := Ψ gs N , e iB N s Ψ gs N (3.14) for the interacting gas. To compute the inner products in (3.13) and (3.14), one applies the map U N,ϕ from (2.8) to the N -body states ϕ ⊗N and Ψ N . Since ϕ ⊗N is the pure condensate, U N,ϕ maps ϕ ⊗N onto the vacuum |Ω of the excitation Fock space, whereas U N,ϕ Ψ gs Lemma 2.2). Conjugating B N with U N,ϕ and, for the interacting gas case, with U V 0 , leads to the identities (see Section 4.2 for the details). Since the inverse Fourier transform leads to the Gaussian probability densities as in (3.10) and (3.11). The mathematical derivation of quantum central limit theorems has first been studied in the 1970s in [9,17] and was followed by many works in different settings, e.g., [13,33,21,16,19,8]. For the ground state of an interacting N -body system, (3.11) was proven in [29] for interactions in the Gross-Pitaevskii regime. For the mean-field Bose gas, the corresponding dynamical problem was first studied in [1], where the authors consider the time evolution generated by H N of an initial product state. This was generalized in [6] to k one-body operators (corresponding to a multivariate setting), in [27] to singular interactions, and in [28] to k-body operators (corresponding to m-dependent random variables).

No Gaussian central limit theorem for low-energy eigenstates
So far, we have considered the situation where the interacting Bose gas is in its ground state. If, instead, it is in a low-energy eigenstate Ψ ex N ∈ C (see also [29, Appendix A]). The general case with n excitations is treated in Proposition 4.7.

Edgeworth expansion for the product state
For the case of i.i.d. random variables, one can go beyond the order N −1/2 of the CLT and approximate the probability distribution of B iid N in an Edgeworth series, i.e., in a power series in powers of N −1/2 which is determined by the cumulants of the distribution. We follow the discussion from [12, Chapter 2]. The 'th cumulant of the distribution of B iid N is defined as for ∈ N, and one easily verifies that where we abbreviated B := B − ϕ, Bϕ . The first three cumulants coincide with the first three central moments; in particular, The basic idea of the Edgeworth series is to expand φ iid N around the characteristic function exp(−s 2 σ 2 iid /2) of the corresponding Gaussian random variable. Since the 'th cumulant is the 'th coefficient in the Taylor expansion of ln φ iid N (s) around zero, one (formally) computes with (3.23) φ iid N (s) = e ln φ iid where we abbreviated κ := κ [ B]. Applying the inverse Fourier transform leads to a series expansion for the probability density b iid N of the random variable B iid N , The functions H j are polynomials of degree j which are even/odd for j even/odd. The complete (formal) Edgeworth expansion is given by the formula The 'th Edgeworth polynomial p iid is a polynomial of degree 3 which is even/odd for even/odd and whose coefficients depend on the cumulants of B of order up to + 2. If the expansion is truncated after finitely many terms, the right hand side of (3.28) is in general no probability density since it may become negative for large values of |x| and is not necessarily normalized. The Edgeworth expansion is thus a local approximation, which is good in the center of the distribution but can be inaccurate in the tails. The expansion (3.28) was first formally derived by Chebyshev and Edgeworth in the end of the 20'th century, and the first proof is due to Cramér. Under the assumption that all relevant moments of the distribution exist, the rigorous statement is usually formulated as an asymptotic expansion of the cumulative distribution function or the probability density, with an error that is uniform in x, i.e., (see, e.g., [34,10,26,12,15] and the references therein). In general, one cannot take the limit a → ∞ since the series does usually not converge. Generalizations of Edgeworth expansions for i.i.d. random variables, for example to different statistics, the multivariate case or the situation when the leading order is not Gaussian, can be found in the literature mentioned above.

Edgeworth expansion for the interacting gas
Let us consider the ground state Ψ gs N of the interacting gas. Due to the dependence of the random variables, this situation is much more intricate than for the product state. In Theorem 1, we prove that the probability density b N of the random variable B N with probability distribution determined by Ψ N is given by in the weak sense of (1.21). Let us provide a formal derivation of this result. As a consequence of the interactions, the cumulants for σ as in (3.12) and α 3 as in (4.25). Note that the leading order of κ gs [B N ] for = 2, 3, 4 is N − /2+1 , which is the scaling behaviour of the corresponding cumulant in the i.i.d. case. Moreover, only even/odd powers of N −1/2 contribute for even/odd. Proving (3.33) in full generality for each ≥ 2 would be extremely tedious, which is why we refrain from following that route for a proof of Theorem 1. Assuming one could prove the (formal) identity for each ≥ 2, a computation along the lines of (3.24) (formally) yields

Weyl operators
As a preparation, we recall in this section the concept of Weyl operators (see, e.g., [31]) and collect some of their well-known properties. For any f ∈ H, the Weyl operator is defined as It is unitary with W * (f ) = W (−f ) and satisfies the shift property for all f, g ∈ H. Conjugation with a Bogoliubov transformation U V , V = U V V U , transforms a Weyl operator into another Weyl operator as The Baker-Campbell-Haussdorff formula yields The number operator transforms under a Weyl operator as for any ξ ∈ F.

Strategy of proof
In this section, we give an overview of the proof of our main result, Theorem 1. We will in the following always assume that Assumptions 1 to 3 are satisfied and that Ψ N ∈ C (η) N for some η ∈ N 0 (see Definition 2.1). Moreover, we will use the notation χ = U N,ϕ Ψ and denote by χ n the coefficients of the asymptotic expansion (2.29). As above, we will only indicate the dependence on η in the notation where it is inevitable. Our goal is to compute the quantity for B N as in (3.9) and g : R → C some integrable and sufficiently regular function. As a first step, we use the excitation map U N,ϕ from (2.8) to re-write the characteristic function as where B denotes the operator on F defined by Applying the substitution rules (2.9) and expanding the square roots 1 − N ⊥ /N in N −1 leads to the following asymptotic expansion (see Section 4.3.1 for the proof): where for B ( ) as in (2.33), with c 0 = 1, for some constant C(a) > 0 and any ξ ∈ F.
Note that the estimate (4.15) is by far not optimal in the powers of N ⊥ except for a = 0, which determines the largest power of s in Proposition 4.4. In combination with Duhamel's formula, Lemma 4.2 leads to an expansion of e isB . Together with the asymptotic series (2.29) for χ, this yields the following expansion of (4.11), which is proven in Section 4.3.2: where T 0 (s) := e isB 0 , (4.18a) with To control the remainders of the expansion, it is crucial that B 0 = a † (qBϕ) + a(qBϕ), hence e iτ B 0 = W (iτ qBϕ) (4.21) is a Weyl operator. Moreover, the operators R and B can be bounded by powers of the number operator. Hence, applying Lemma 4.1 repeatedly and making use of the fact that moments of the number operator with respect to χ are bounded uniformly in N (Lemma 2.2b) yields an estimate of the error S a (s) (see Section 4.3.3 for a proof): where C B (a) ≤ C(a)(1 + B 3a+4 op ) for some constant C(a). The next step is to compute the coefficients in the expansion (4.17), which is done in Section 4.3.4. Since an explicit evaluation to any order is too complex to obtain in full generality, we focus on the dependence of the coefficients on s: where 25) and where Θ N with η > 0, we explicitly compute only the leading order polynomial. The proof of the following proposition is given in Section 4.5.

Proof of Lemma 4.2
We decompose B N as for any ξ ∈ F and b ∈ N 0 and where [ R 2b , N ⊥ ] = 0. Besides, by Lemma 2.2c, there exists some r B (a) ∈ R with |r B (a)| ≤ C(a) B op such that Consequently, for B as in (4.13) and where R a satisfies (4.15).

Proof of Proposition 4.3
From (4.13), it follows that B − B 0 = N −1/2 R 0 with Hence, (4.16) implies that a; (s) as in (4.18) and (4.19). This implies (4.17) by (2.29) because By (4.15), we find for any ξ, ξ ∈ F that where we used the notation f = qBϕ and abbreviated By definition (4.14) of the operators B and using Lemma 2.2c we find that Using this estimate repeatedly yields by definition (4.40) of γ . Moreover, k ≤ | | + 1, hence hence the operators B from (4.14) transform as We summarize these expressions in the following way, keeping only track on the τ -dependence and on the total number of creation/annihilation operators a : for some ξ µ ∈ L 2 (R dµ ). Here, we used the notation a −1 := a and a 1 := a † .
(a) Two polynomials (4.49) are equivalent with respect to the relation ∼ iff they have the same degree j and the number of operator-valued distributions a x in each summand is even/odd for j even/odd. We denote the representatives of the equivalence classes with respect to the relation ∼ by F j , i.e., (b) Two polynomials (4.49) are equivalent with respect to the relation ∼ j iff they have a degree ≤ j and the number of operator-valued distributions a x in each summand is even/odd for j even/odd. We denote the representatives of the equivalence classes with respect to the relation ∼ j by F ≤j , i.e., F ≤j ∼ j F j (4.51) for any j ≤ j. When using the notation F ≤j , we will drop the index j from ∼ j .
Computation of (4.65a). As above, we abbreviate ν = U 0 qOϕ + V 0 qOϕ. With (4.3) and (4.4), we find For any one-body operator A, any ONB (ϕ i ) of H ⊥ , and g ∈ H ⊥ , we have A ij a(V 0 ϕ i ) + V 0 ϕ i , g + a † (U 0 ϕ i ) + g, U 0 ϕ i × a(U 0 ϕ j ) + U 0 ϕ j , g + a † (V 0 ϕ j ) + g, V 0 ϕ j , (4.67) where we denoted A ij := ϕ i , Aϕ j . Consequently, expanding the exponential yields where c 1 , c 3 ∈ R are given by Computation of (4.65b). Using that Computation of (4.65c). We find, using first (4.5) and then Lemma 2.2c, that where H k (x) is the k-th Hermite polynomial as defined in (3.26). This yields (1.21) with polynomials p j (x) of degree 3j + 2η in x ∈ R which are even/odd for j even/odd. Note that the coefficients of the p j must be real-valued because E[g(B N )] ∈ R for real-valued g.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.