Functional Integral and Stochastic Representations for Ensembles of Identical Bosons on a Lattice

Regularized coherent-state functional integrals are derived for ensembles of identical bosons on a lattice, the regularization being a discretization of Euclidian time. Convergence of the time-continuum limit is shown for various discretized actions. The focus is on the integral representation for the partition function and expectation values in the canonical ensemble. The connection to the grand-canonical integral, and a number of differences, are discussed. Uniform bounds for covariances are proven, which simplify the analysis of the time-continuum limit and can also be used to analyze the thermodynamic limit. The relation to a stochastic representation by an ensemble of interacting random walks is made explicit, and its modifications in presence of a condensate are discussed.


Introduction
The operator-algebraic ('second-quantized') formulation of quantum statistical mechanics [12] provides a precise mathematical framework for the study of quantum many-body models. For the analysis of equilibrium states and dynamics at (small) positive temperatures, functional integral representations (FR) and stochastic representations, in particular interacting Brownian motions (IBM), have proven very useful alternative approaches. In theoretical physics, formal (i.e. nonrigorous) FIR have almost become the method of choice; stochastic representations also have a long history, going back to Feynman [21,22]. In mathematical physics, FIR have also been studied already for a long time. Clearly, there are subtleties involved when dealing mathematically with infinite-dimensional integrals. A traditional way of obtaining natural regularizations that also provide a strict relation to the algebraic formulation proceeds via Suzuki-Trotter product formulas and coherent-state integrals, to get a discrete-time functional integral, the time-continuum limit of which can then be analyzed. For fermionic systems, this has (in combination with multiscale techniques) led to mathematical results that have not yet been proven by any other method, see for instance [10,13,14,16,19,28,29,36]. Bosonic systems pose the additional technical problem of dealing with the unboundedness of operators, which in the FIR translates to the unboundedness of fields. This has been dealt with using decompositions into large and small fields. Using these has led to proofs of existence of massless scalar field theories with a fixed short-distance regulator [26,18]. More recently, the existence of the time-continuum limit of the Bose system (again with a fixed short-distance regularization) has been proven [7]. IBM representations start from the Feynman-Kac formula. Given this formula for a single particle, it is straightforward to generalize it to N particles. Unlike the coherent-state formalism, where (anti)symmetrization is built in, Bose (resp. Fermi) statistics then needs to be imposed by an explicit (anti)symmetrization operation, given by a (signed) sum over permutations of the labels of the N particles. The IBM allows to bring large-deviation techniques to bear on the problem [1,2,3]. It is naturally suited for the canonical ensemble (but can also be used in the grand canonical one).
For Bose systems, it is a major outstanding problem in mathematical physics to prove the occurrence of Bose-Einstein condensation (BEC) in the thermodynamic limit at fixed nonvanishing density. The case of the Gross-Pitaevskii limit for N bosons in a trapping potential of fixed extension has been treated in [33,31,32], where BEC was proven in the sense that the reduced density matrix has an eigenvalue of order N , associated to a condensate state. In the FIR, BEC shows up in the behaviour of the 'zero mode' of the field (the constant part of the field in a spatially homogeneous system). At small enough temperatures and high enough density, the most likely value for the zero mode is nonvanishing, and condensation takes place when it acquires a nonzero expectation value, which breaks the global U (1) symmetry of the action spontaneously. In the IBM representation and its subsequent large-deviation analysis, BEC is characterised as the occurrence of infinite cycles (as N → ∞) in the permutations dominating the above-mentioned symmetrization sum.
Since both FIR and IBM represent the same mathematical objects (the standard traces for partition functions and expectation values over Fock space), it is clear that they are related. Indeed, this relation is essentially Symanzik's random-walk representation of quantum field theory, the rigorous version of which has led to breakthrough results in mathematical quantum field theory [4,23]. (A detailed exposition is in [20].) It is interesting to make this relation explicit for specific models, to see how far techniques can be combined. (For the systems considered in this paper, this is done here in Section 5.) Another motivation for this is as follows.The IBM (and in the lattice setting, its discrete version, the interacting random walk (IRW)) corresponds to a one-particle kinetic energy given by the Laplacian, i.e. p 2 2m , where p is the momentum of the particle. In presence of a condensate, Bogoliubov [11] predicted a spectrum of low-lying states above the ground state which is linear in |p| at small |p|, namely E B (p) ≈ c|p|, where c is the velocity of sound in the condensate: the excitations are sound waves. It is certainly an interesting question how this spectrum arises in a collection of Brownian motions when an infinite cycle forms. In the FIR, there is a short answer: once a condensate has formed, the fluctuation fields describing deviations from the condensate have a matrixvalued covariance, the eigenvalues of which exhibit the Bogoliubov spectrum (see, e.g. [9]). It is thus natural, given that a relation between FIR and IRW representation can be derived without reference to a condensate, to look if one can obtain a stochastic representation also in presence of a condensate. First steps in this direction are described in Section 9. It is shown there that, at least in the time-continuum limit and for a local interaction, the kinetic term does generate a stochastic process, albeit one with long-range jumps (corresponding to c|p|). A main difference is in the interaction term, which differs from the case without condensate by terms causing branching and coalescence of random walks.
The question about the relation of FIR and IRW also motivates taking a closer look at the FIR for the canonical ensemble of bosons, and this is the main focus of this paper. The methods developed here, and some results, also apply to the grand canonical ensemble. The latter has been analysed mathematically in a series of papers [5,6,7,8], in which essential foundations for an analysis of Bose condensation using methods of constructive quantum field theory, in particular mathematically rigorous Wilsonian renormalization group methods, have been laid, and the flow has been studied into the symmetry-breaking regime. Further papers in that series can be obtained from [15]. More recently, the FIR for boson systems and its relation to IBM representations have been studied in [24,25], with a focus on the grand-canonical ensemble.
In this paper, I give a detailed discussion of the FIR for the canonical ensemble of identical bosons on a lattice, with a careful derivation from the operator formulation and a detailed discussion of the time continuum limit and the freedom of using different discrete-time actions that may be suitable for different aspects of the problem. The main points and results are (similarly to [24,25]) based on using auxiliary (Hubbard-Stratonovitch) fields, in terms of which the integral can be rewritten in terms of a positive Gaussian measure, with an interaction given by the inverse of a determinant (grand-canonical case) or a large permanent (canonical case) of a boson covariance in the background field given by the auxiliary field. The permanent implies the explicit symmetrization mentioned above. In this representation, there is no oscillatory part of the Gaussian integral. The auxiliary fields are ultralocal, i.e. a lattice-regularized version of white noise.
Some details of the grand-canonical and the canonical integrals are actually quite different, both concerning the covariance in a background field and the factorization properties of inverse determinant vs. permanent (see Section 8).
A main technical result shown here is that the fully regularized model has properties that allow for bounds that are uniform in the auxiliary fields, very much in the spirit of the loop vertex expansion [34]. This allows for completely elementary and technically very simple proofs of the existence ot the timecontinuum limit for a variety of different actions. They do not require any mutliscale analysis and in that sense, they are as simple as the corresponding ones in the fermionic case [17,35]. The uniformity of the bounds also allows to analyze the infinite-volume limit and the properties of correlation functions in that limit, in cases where renormalization is not required. This will be done in another paper [40]. Although technically comparatively simple, the hypotheses and proofs contain a few subtle points (to be discussed below, see in particular Sections 8 and 10), and it remains to be seen whether a similar uniformity also holds in genuinely multiscale situations (for a multiscale version of the loop vertex expansion, see [30]).
Since one aim of this paper is to provide a technically simple approach, I have added appendices to make it essentially self-contained, given a knowledge of the operator algebraic formulation of quantum statistical mechanics (as exposed, e.g., in Section 5.2 of [12]).

Identical Bosons on a Lattice
This section contains the detailed setup of the class of models considered here, namely lattice models of identical bosons with a short-range density-density interaction, in the second-quantized formulation.

Kinematics.
2.1.1. Space. To avoid technical issues in the setup, I start from a lattice system. Let η > 0, L be a very large integer multiple of η, and X is a lattice with spacing η > 0 and periodic boundary conditions, i.e. a discrete torus. The limit η → 0 at fixed L (i.e. L/η → ∞) is usually called the continuum or ultraviolet limit, and the limit L → ∞ is called the thermodynamic or infrared limit. In principle, one would like to take both limits, but this paper is mainly concerned with the thermodynamic limit at fixed η > 0. It is expected that in this limit (and with an appropriate choice of kinetic term 3 and interaction) Bose-Einstein condensation happens at low enough temperatures, even in presence of an ultraviolet regulator. For maps x → f x and x → g x from X to some C-algebra A, let For f and g arising as restrictions of continuous functions on R d /LZ d , this is a Riemann sum approximation to an integral that arises in the limit η → 0, which motivates the alternative notation (f | g) X = x∈X f x g x = X f g, as well as the notation x F x = η d x∈X F x . This notation will also be used at fixed η > 0, e.g. at η = 1.

2.1.2.
Hilbert spaces for bosons and the CCR algebra. The quantum mechanical Hilbert space for a spinless particle on X is H = L 2 (X, C). In our setting it is simply the finite-dimensional space C X of all maps f from X to C, with inner n! π∈Sn f (x π(1) , . . . , x π(n) ). Here S n denotes the set of all permutations on {1, . . . , n}. The bosonic Fock space is defined as Recall the standard C * algebra formulation of quantum many-particle systems given, e.g., in [12]. The algebra for boson systems is generated by bosonic annihilation operators a and creation operators a † satisfying the canonical commutation relations (CCR) (4) a x a y − a y a x = 0 and a x a † y − a † y a x = η −d δ x,y for all x, y ∈ X . Let is given by the vectors where the sequences ν = (ν x ) x∈X ∈ N X 0 satisfy the condition x∈X ν x = N .
Define the local density operators and the number operator N by Every n x maps the N -particle space F to itself, and on this space, it has operator norm n x = N . The restriction of N to F (N ) B is N times the identity operator. The n x all commute: (9) n x n y = n y n x for all x, y ∈ X . 4 2.2. Hamiltonians. For a single particle, the Hamiltonian is given by a selfadjoint matrix E, acting as (Ef ) x = y E x,y f y , i.e. a hopping matrix. Typical choices will be the discrete Laplacian −∆, or −∆ + m 2 , with m 2 > 0, or, for particles in an external potential W , −∆ + W (x). This operator will be assumed to be nonnegative. The particles interact by a two-body potential, i.e. v x,y ∈ R is the interaction energy contributed by a pair of particles at sites x and y. The interaction is symmetric: v x,y = v y,x . Thus v can be regarded as a self-adjoint operator on C X . v is called translation-invariant if v x+z,y+z = v x,y for all z, and in this case I use the notation v x,y = v(x − y).
The Hamiltonian on the bosonic N -particle Hilbert space F and where v xm,xn acts as a multiplication operator. In second-quantized form, it can be written as Inserting the notation (2), gives H = H 0 + V with (13) H 0 = (a † | E a) X and V = 1 2 (n | v n) X The case v = 0 describes free (i.e. noninteracting bosons).

2.3.
Ensembles. In the following, the Hamiltonian H is assumed to be nonnegative; in the above setup this is the case if both operators E ≥ 0 and v ≥ 0.
2.3.1. Canonical ensemble. For β > 0 define (14) [A] (N,β,X) = Tr F B e −βH A P N , Then Z (N,β,X) c is the canonical partition function and A (N,β,X) is the normalized canonical expectation value of A, for a system of N bosons at inverse temperature β on configuration space X. The range of the projection operator P N is finite-dimensional, so P N is trace class on F B , and the trace in (14) exists whenever e −βH A is a bounded operator. Because P 2 N = P N and the trace is cyclic, (16) [A] (N,β,X) = Tr F B P N e −βH A P N .
Because H conserves particle number,

Canonical Bose Ensemble
Thus it suffices to consider expectations of operators F that map F (N ) B to itself. It is convenient to introduce a generating function with insertion F as (19) Z (N,β,X) Normalized expectation values can then be obtained by differentiation as The dependence on H has been made explicit in the notation. In general, for Hamiltonians preserving particle number, normalized canonical expectation values are unchanged under shifts in the energy, as follows. For any α ∈ C and ψ ∈ F B , (21) e αN P N ψ = e αN P N ψ As a consequence (22) The prefactor drops out in the quotient defining the normalized canonical expectation value, so Grand-canonical ensemble. The grand canonical partition function at chemical potential µ is and the grand canonical expectation value is defined similarly. Convergence of the sum over N in general requires conditions on µ and H (e.g. µ < inf spec H 0 when V = 0) which will be specified when appropriate.
Returning to the shift invariance (23), it is very important to note that this invariance does not hold for the partition functions. Instead, the prefactor in (22) implies that in the sum for the grand canoncial partition function, the parameter µ gets shifted to µ − α. Depending on α, this shift may cause the sum over N to diverge. Thus, when one attempts to obtain the usual relation between canonical and grand canonical ensembles, arbitrary shifts are no longer allowed. The convergence condition will show up prominently in the coherent-state integral derived in the next section.

Main Results
Here I state the main integral formulas and convergence theorems.
By definition, E generates a stochastic process iff for all τ ≥ 0 and all x and y ∈ X, the x, y matrix element of e −τ E is a nonnegative real number: Theorem 1. Let H have kinetic term E and interaction v, as in (13). Let E and v be translation invariant, E > 0 and v ≥ 0 as operators, and E generate a stochastic process. Then the canonical partition function has the integral representation Here X includes a discrete time lattice of spacing ǫ = β ℓ , i.e. X = T × X, where T = ǫ{0, . . . , ℓ} 1 The integral runs over complex fields a : X → C, Da denotes the volume form on C X , and the action S X is The integral notation τ is shorthand for ǫ τ ∈T , ∂ τ denotes the discrete forward time derivative with the boundary condition a(β + ǫ, x) = 0, and The limit remains unchanged if (E (ǫ) a)(τ, x) is replaced by (Ea)(τ + ǫ, x) or by (Ea)(τ, x) in the action.
Two important features of this representation are characteristic for the canonical ensemble: (1) the N -dependent polynomial factor (a(β) | a(0)) X N appears in the integrand, and (2) the time derivative operator in the quadratic part of the action does not have a periodic boundary condition at τ = β. (See also Theorem 6, in particular (79)). The 'mass term'v(0) may be surprising at first sight. It arises in the proof in an integration by parts step which is necessary to show convergence. Its occurrence is also related to that of the prefactor e β 2 v(0)N , which is necessary to fulfil certain positivity conditions in intermediate steps of the proof. Another variant of the interaction, where the |a| 4 term is replaced by a term of at most quadratic growth for |a| → ∞, is discussed in Section 4. Theorem 1 is proven in Section 7, using the preparations made in Sections 6. Translation invariance is assumed here only to simplify the statement. In Section 7.2, a more general statement, Theorem 16, where translation invariance is not assumed, is proven. Theorem 1 applies in particular to a E of the form E = E 0 − µ1, with E 0 > 0 and µ ≤ 0, and it implies that the limit of the integral (27) equals (up to the explicit prefactor) the canonical partition function Z (N,β,X) c . By (22), Z (N,β,X) c is analytic in µ for all µ ∈ C. If v > 0, the quartic term in (29) is bounded below by δ τ,x |a(τ, x)| 4 for some δ > 0, so the integral (27) converges absolutely for all µ ∈ C as well, and it defines an analytic function of µ. Moreover, for µ 0 < 0, the convergence as ℓ → ∞ is uniform on compact subsets of {µ ∈ C : Re µ < µ 0 }, hence the limiting function is analytic in µ there. Thus it coincides with Z (N,β,X) c for all µ by the identity theorem, and hence the integral representation extends to all µ.
The action (29) is such that it admits a rewriting by an integral over auxiliary fields. Indeed, the theorem is proven by using such a representation, which is the basis of all the analysis done here, hence more fundamental. It is given as follows.
Theorem 2. Assume that E > 0 and v ≥ 0, and let F be any operator on F (N ) The integration over h runs over R X , and µ V is the normalized, centered The operator Q(h), given in (78) below, has a positive hermitian part, and the integral over a converges absolutely. The vectors v a ∈ F B are unnormalized coherent states, given in (38). For F = 0, the matrix element is The covariance C(h) = Q(h) −1 is entire analytic in h on C X . If E generates a stochastic process, then C(h) satisfies the uniform bound Theorem 2 follows from Theorem 6, Lemma 7, and (86). It will be used in further work to study properties of correlation functions.
Theorem 3. Let µ ∈ R and assume that E − µ > 0. Then the grand canonical partition function (24) has the representation Z where K(h) is given by Q(h), but with E replaced by E − µ and a periodic boundary condition in time (see (182) below). The grand-canonical covariance G(h) = K(h) −1 is analytic in h on a neighbourhood of R |X| . It E generates a stochastic process, then G(h) satisfies the uniform bound G(0) is the time-ordered Green function for free bosons with kinetic term E −µ.
Theorem 3 follows directly from Theorem 17.

Coherent states.
Here I summarize some well-known properties of coherent states, in the form that I shall use them in the following. For convenience of the reader, their (elementary) proofs are given in Appendix A.
Coherent states are eigenvectors of the annihilation operators: for any function a : X → C there is a nonzero vector v a ∈ F B such that for all x ∈ X (37) a x v a = a x v a . The vector v a is explicitly given by (38) v a = e (a|a † ) X Ω where Ω is the vacuum vector in the bosonic Fock space (see (5)). The right hand side of (38) is defined by the series expansion of the exponential. Theseries is norm-convergent for all a, so the coherent state v a is entire analytic in a.
Coherent states are grand canonical objects, since a maps F Their inner product on Fock space is Thus v a ∈ F has norm (40) v a = e 1 2 (ā|a) X .
Denote the unit vector the orthogonal projector on the subspace spanned by κ a by |κ a κ a |, , and C R = {z ∈ C : |z| ≤ R}. Then, in the sense of strong convergence, and for trace class operators A, (see Theorem 2.26 of [5]).

4.2.
The integral for the canonical ensemble. By (44), and with the notation (48), A priori, the application of (43) results in two limits R → ∞ and R ′ → ∞, but since both exist, the limit can be taken keeping R = R ′ . This will again be done later, when there are many insertions of type (43).
To streamline notation further, I will write By (41), and with the notation (48) 4.2.1. Matrix element of the N -particle projection.
Proof. Let N ∈ N 0 . Because the expansion (38) for the coherent state v a is norm convergent for any a and the projection P N is continuous, Because a † is the adjoint of a and because v a ′ is an eigenvector of a, The inner product v a ′ | Ω = 1.
and the limit R → ∞ is given by an absolutely convergent integral.
Proof. The integral representation follows by inserting (51) into (19). If the norm of e −βH is bounded by M > 0, then (40), the absolute value of the integrand is bounded by Note that if H commutes with N, then (57) tr e −βH P N = tr P N e −βH P N = tr e −βP N HP N P N so for H that is unbounded below, a formula similar to (54), with H replaced by H N = P N HP N , still holds.
By the residue formula (r > 0), the canonical partition function then becomes (59) It is this form which allows to connect conveniently to the other ensembles in the limit of large N by a stationary-phase argument.
A little calculation gives This is easy to interpret: when a y is applied to the N -particle state P N v a , the result is an N − 1-particle state, so the power decreases to N − 1 in that term.
It is straightforward to generalize this to operators of the form 2 Noninteracting bosons. It is instructive to see how this leads to the standard formulas for the canonical partition function for free bosons. In the absence of interaction, H = H 0 , so (45) can be used to rewrite the canonical partition function as Let r < 1, then the matrix in the quadratic form in the exponent has a positive definite hermitian part, and therefore the Gaussian integral is absolutely convergent, hence the limit R → ∞ can be taken. The Gaussian integral over a and a ′ then gives the inverse of its determinant, so that Evaluated in an eigenbasis of E, this determinant factorizes, and the (convergent) expansion in z then gives the sum which follows straightforwardly from the definition of the canonical ensemble as a trace. Here the eigenbasis is indexed by the set A, i.e. (E α ) α∈A are the eigenvalues of E.
Without the contour integral, the canonical partition function can be written as a permanent: Again, the limit R → ∞ can be taken since the real part of the quadratic form is strictly positive. Writing out each factor (ā ′ , a) X = xā ′ x a x and doing the Gaussian integral (noting that this time, the determinant is 1 and the inverse obviously 1 e −βE 0 1 ), this gives Again by (44) and (43), and with the new notation a 0 = a and a ℓ = a ′ , this becomes By (43), As discussed above, the limit may be taken with R ′ = R. By (41) this can be rewritten in terms of the unnormalized coherent states as 4.4.2. Auxiliary fields. As noted before, the interaction v is assumed to be a nonnegative operator, and the n x generate an abelian Banach algebra, so the exponential of the interaction can be written as a Gaussian integral where dµ v denotes the normalized Gaussian measure on R X with mean zero and covariance v. This identity is used for every  in the product (72). The Gaussian integral is absolutely convergent, hence can be taken out of the inner product and the integral over the a fields. This, and (45), then gives The collection of additional integration variables h ,x indexed by (, x) ∈ {1, . . . , ℓ} × X, is called the auxiliary field in the following. In (75) and similar expressions, e i √ ǫh denotes the multiplication operator, With the notations the following statement holds. Then If E > 0, then the quadratic form Re (ā | Q(h) a) X is strictly positive, so the Gaussian integral over a converges absolutely and the limit R → ∞ can be taken by dropping the restriction on the integral domain. In particular, for F = 0, Proof. The explicit formula (80) has already been derived; recall also (51) and (61). The positivity of the real part is equivalent to operator positivity of , and it holds because for E > 0, the diagonal part of Q S (h) strictly dominates the nondiagonal parts (note that e −ǫE e i √ ǫ h = e −ǫE < 1). For E = 0 and h = 0, Q S (0) = −∆, the Laplacian on {0, . . . , ℓ}, which is nonnegative. In general, set K = e −ǫE and U  = e i √ ǫh . Then K is a positive operator because E is hermitian, and U  is unitary for all  because h  is real-valued. Consequently, is a positive definite scalar product on C X , and  (78), is a matrix Therefore, det Q(h) = 1, and the inverse of Qh, the covariance exists and is also upper triangular, and given by a terminating Neumann series as . Because the inverse is given by a finite sum, C(h) is entire analytic in h. In fact, for every finite ℓ, it is a trigonometric polynomial in h.

4.6.
Remarks on the structures arising from (81). Eq. (81) is useful in various respects. It represents the partition function as a Gaussian integral both in the original boson field a and in the auxiliary field h. Integrating over the auxiliary field will give a regularized version of the formal functional integral for boson systems; see Subsection 4.6.2. Instead integrating over the a-field gives a representation in terms of the auxiliary field only. This will be used to derive a random-walk representation which is similar to that of [3, 1, 2] (see Section 5) but, being an integral with respect to a real Gaussian measure it also provides a starting point for a field-theoretical analysis by convergent expansions.
4.6.1. Stochastic field equation. One can also derive a complex stochastic equation for a by integration by parts in the integral. Let F (a) denote the normalized canonical expectation value, then using and integrating by parts gives for k < ℓ the 'Schwinger-Dyson' equation Q(h)a = 0, which, in a formal continuum limit corresponds to the Gaussian white noise average over h of the equation The complex integral obtained after integrating out the auxiliary field. For simplicity I consider here the case of an on-site interaction, i.e. v x,y = vη −d δ x,y , with v ∈ R, v > 0. More general v's can be treated as well in terms of expansions, but the essential points are visible already in this case. With this choice of v, the h-integral factorizes both over the time slice index  and over the space index x, as follows. The h-dependent part of the quadratic form withã −1 = e −ǫE a −1 , so its exponential factorizes; the assumption that v is on-site implies that also the Gaussian integral factorizes. Calling z = a −1,x a ,x , the integral over h = h ,x is then simply The function V(z) = − ln Φ(z) is then the interaction term in the action for the a-integral. Obviously, |Φ(z)| ≤ e |z| for all |z|, so the large-field behaviour is bounded by a quadratic term. Moreover V(z) = z + 1 2 ǫvz 2 for small |z|, i.e. V contains a quartic interaction term for small fields.
Clearly,ã −1,x a ,x is complex, and therefore this representation loses the positivity of v. However, it is possible to obtain a representation with a positive interaction term, in the sense that one can prove that the action can be changed to one with a positive interaction term without changing the limit ǫ → 0. This will be done below in Section 7. 4.6.3. The probabilistic expectation obtained after integrating over the a-fields. The a-integral in (81) is Gaussian as well. Performing this integral, the polyomial of degree N in the integrand results in a permanent of order N , specifically The integrand is still complex because C(h) is a complex function, but this function will be shown to be bounded. The Gaussian measure itself is real and normalized, i.e. a probability measure). This motivates the notation Introduce the additional notation for the X-averaged permanent, A similar expression holds for the generating function with F insertion. This representation will be the basis for most of the bounds given in this paper.

The Random-Walk Representation
By (86), the k th factor in the product in (92) is given by .
In the basis of position eigenstates, e −i √ ǫh j k is diagonal, and thus writing out all matrix products gives for the product of transition amplitudes along the walk. Let X = {1, . . . , ℓ}×X. Define the density associated to y (k) by ν (k) : The Gaussian integral over h can now be done by the standard formula and the bilinear form in the exponent then worked out in terms of summations over k and k ′ ∈ {1, . . . , N } using the definition of ν. With the notations the canonical partition function then becomes in the translation-invariant case 3 Thus the partition function is a sum over collections of N random walks y (1) , . . . , y (N ) that have a local-in-time interaction v. This is a discrete analogue of the representation as symmetrized interacting Brownian motions of [3]. The convergence of this partition function for interacting random loops as ǫ → 0 holds as a direct consequence of the Lie product formula (68). While it is clear that the representations are related since they equal the same partition function, this provides an explicit illustration of the role of the factor (51) in the integral representation of the canonical partition function, and a natural finite-dimensional approximation of the Wiener integral over Brownian motions.

Analysis 1: Covariances
This section contains the basic estimates for the covariance. They are elementary and the proofs are easy, but they will be crucial because they are uniform in the auxiliary field h.
From now on, the spatial lattice constant will be set to unity, η = 1.
6.1. Stochasticity and bounds for the covariance.
Lemma 7. If E generates a stochastic process (see (25)), then for all  ′ ≥ , all x, x ′ ∈ X, and all h Proof. By (86), and writing out the matrix products in the random-walk notation of (97), All h-dependent factors have absolute values equal to 1, and, by hypothesis, the matrix elements of e −ǫE are all nonnegative, hence unaffected by taking absolute values. Thus 6.2. Modified covariances. In this section, a number of modifications of the operator Q(h) are introduced, and their properties, as well as those of their inverses, are discussed. These modifications will be related to different actions that all reproduce the same limit for the canonical partition function as ǫ → 0. 17 is invertible, and If E generates a stochastic process, then for all h and all , Proof. Eq. (112) holds because B 0 = 1. By definition, Every B  is diagonal in the spatial indices, so B(h) −1 is diagonal, and Thus (113) follows by matrix multiplication and Lemma 7.
6.2.2. Q 2 and Q 3 . Consider furthermore the following two operators motivated by a small-ǫ expansion. Let u : X → R and letẼ be a self-adjoint linear operator on C X , and regard, as with h, functions of u as multiplication Let h ∈ R X and u ∈ R X .
(a) Q 2 (h, u) is invertible and its inverse (HereK = e ǫ 2 u e −ǫẼ e ǫ 2 u and b  = e − ǫ 2 u a  .) (c) IfẼ generates a stochastic process, then For any δ > 0, there is ℓ 0 = ℓ(δ, β, Ẽ , u ) so that for all ℓ > ℓ 0 , all ,  ′ ∈ {0, . . . , ℓ} and all x, x,x ′ < δ Consequently, for ℓ > ℓ 0 , the right hand side of (121) is bounded uniformly in ,  ′ and ℓ: Proof. (a) Viewed as a matrix in ,  ′ , Q 2 is upper triangular, and the diagonal part D of Q 2 , given by is invertible for all u, h ∈ R X . For all h ∈ T ǫ , the inverse D(h, u) −1 is analytic in h because it only contains factors of (1−i √ ǫ h j,x ) e −ǫux ) −1 . The inverse of Q 2 is given by a Neumann series that terminates after at most ℓ+1 terms, hence has the same analyticity properties. The first equality in (b) holds because the hdependent term drops out of the hermitian part of Q; the inequality the follows by the hypothesis made in (b). (c) The proof of (121) is a straightforward adaptation of that of Lemma 7, writing out the finite Neumann sum for C 2 as a matrix product, and using that by (124), To see (122), note that an easy variation of the standard proof of the Lie product formula (as given, e.g., in [37]) results in and an estimate of the last factor gives O(ǫ 2 ), hence (see Appendix C) Given δ > 0, ℓ 0 is defined by the condition R ≤ δ. Then use |R x,y | ≤ R , and | −  ′ | ≤ ℓ, so that | −  ′ |ǫ ≤ β, hence e −(− ′ )ǫ(Ẽ−u) x,x ′ ≤ e βE min . Thus for ℓ > ℓ 0 , (123) holds.

Lemma 10.
Let h ∈ R X and u ∈ R X .
(c) IfẼ x,y ≤ 0 for all x = y and if ǫ is so small that 1 − ǫẼ x,x ≥ 0 for all x, then Moreover, the right hand side of (129) is bounded uniformly in ℓ.
Proof. Similar to the proof of Lemma 9.
6.2.3. Q 4 . To place the kinetic termẼ in the time-local ( =  ′ ) term, consider The bounds for this covariance are analogous to the previous ones; the argument is slightly different, and will be given in several steps in the following, without collecting it in a single Lemma.
IfẼ > 0, then Lemma 11. Let B be the algebra of bounded operators on a Hilbert space and Q, H ∈ B. Furthermore, let Q 0 be a positive operator and H = H † , and let Then Q −1 exists , and Proof. Because Q 0 > 0, it is invertible and has a positive square root Q 1 2

. Set
IfẼ > 0, and because h is real, Q 4 (h) satisfies the hypotheses of Lemma 11 with Q 0 = Q 4 (0). In the estimate of products as above, then use When doing the terminating Neumann series for this inverse, use the spectral theorem forẼ and the elementary inequality that for all x ≥ 0 and  ≤ ℓ ∈ N 0 to see that the norm of Q 4 (0) −1 is bounded uniformly in ℓ.
Lemma 12. IfẼ ≥ 0 and the off-diagonal matrix elements ofẼ are negative, then C 4 (h) = Q 4 (h) −1 satisfies and the right hand side is bounded uniformly in ℓ.
Proof. The diagonal of Q 4 (h), as a matrix in the  indices, contains for all x and x ′ . By hypothesis, the geometric series for the inverse x,x ′ . Using the terminating Neumann series for Q −1 4 then concludes the proof.

Analysis 2: Time-Continuum Limit
The existence of the time-continuum limit ℓ → ∞ for the integrals (80), (91), and (92) for the canonical partition function holds simply by the Lie product formula, and it is already in this form very useful to do analysis. Sometimes, however, one is interested in simpler actions that resemble more closely those occurring in formal continuum-time functional integrals. Heuristically, they are obtained by dropping higher-order terms in ǫ, but proving that these terms can really be dropped is not straightforward. For many-fermion systems, it has been done using tree expansion techniques and determinant bounds, sometimes combined with multiscale analysis [10]. Later [17,35], it was done without multiscale techniques.All these proofs rely on determinant estimates which allow for convergent perturbation expansions under suitable conditions. The bosonic case is harder. In [7], existence of the time-continuum limit for a FIR of the many-boson system in the grand-canonical ensemble was proven by multiscale techniques (renormalization group using a decimation transformation in time). Ref. [7] contains much more than just the proof that the limit exists; it also gives a rigorous justification for an effective action with a fixed short-time cutoff, which can be used in the analysis of Bose-Einstein condensation in the thermodynamic limit [8]. The analysis required to do this involves decompositions in large and small fields, a conceptually clear and analytically powerful field theoretical method which, however, entails some technical overhead.
In this section, I prove that a variety of different actions give rise to the same limit ǫ → 0 of the partition function and the unnormalized expectation values of the canonical ensemble of bosons. The proof given here does not require any multiscale technique, only the uniform bounds on covariances given in the last section. This makes it possible to avoid a large-field analysis. As mentioned, it is not strictly necessary for the further analysis of the system to consider modified actions, but it serves to illustrate the application of the main observations made above, namely (i) the analyticity of the covariances in the auxiliary field h (ii) their uniform boundedness as functions of h, and (iii) the fact that the permanent in the canonical ensemble factorizes easily 21 over covariances and hence allows for a straightforward use of (i) and (ii) in remainder estimates. The factorization can be expressed in terms of the cycles of permutations, and then used to study the thermodynamic limit by convergent decoupling expansions [40].
Recall that ǫ → 0 always means ℓ → ∞ with ǫ = β/ℓ; this will be important in estimates since convergence or uniform boundedness as ǫ → 0 really means convergence or uniform boundedness as ℓ → ∞, and statements that constants are independent of ǫ will mean that they are independent of ℓ as well.
7.1. Convergence of modified partition functions. For simplicity, I consider only the case F = 1 here.
The following lemma implies that the h-dependence can be moved to the diagonal part of the Q-matrix, without changing the partition function.
It is interesting to note that formally, already this gives rise to an integral representation where the interaction term is local in time. Assuming for simplicity that the interaction is local, i.e. v x,y = v x δ x,y , where V = − ln Φ, with Φ given in (90). The interaction term then has the properties that it is local in the time index , hence positive, and moreover Φ can grow at most like e |a| 2 as |a| → ∞. Formally, one can also derive this representation directly from the trace using a modified coherent-state formula. That derivation is, however, formal because integrals that get exchanged do not converge absolutely. This is reflected here in that Re (ā | Q 1 (h) a) X is not positive for arbitrary h, hence Q 1 (h) cannot be used as a covariance for a Gaussian integral, so again, going from the right hand side of (143) to the right hand side of (144) involves non-convergent integrals. This problem may be bypassed by an additional regularization suppressing large a.
In a heuristic procedure, a further step is to expand the exponential to linear order, e −i √ ǫh,x = 1 − i √ ǫh ,x + O(ǫ) and drop the order ǫ remainder. The integral over h then becomes Gaussian, and results in a quartic interaction for the a fields. This does not really work, but, as shown in the following, a slight modification does, using integration by parts in the auxiliary field h, as given by Lemma 21 in Appendix D.

Lemma 14.
Assume that E generates a stochastic process. Let Q 2 (h, u) be defined as in (116), with u x = 1 2 v x,x andẼ = E, and 22

Canonical Bose Ensemble
Overview of the proof. -The permanent P N,X contains a product over n ∈ {1, . . . , N } of C 1 (h) (0,xn),(ℓ,x π(n) ) in Z c and a product of C 2 (h) (0,xn),(ℓ,x π(n) ) in Z 2 . By the discrete product rule and the resolvent formula, a difference will be made explicit. This term is of order ǫ for any u (e.g. for u = 0), but because a sum over ı ∈ {0, . . . , ℓ} arises as well and ǫℓ = β, a straightforward estimate only implies that the difference of partition functions is finite, but not that it vanishes. To get a better estimate, integrate by parts using Lemma 21, with a suitable (real) choice of u. Because E generates a stochastic process and u is real,Ẽ also generates a stochastic process. Thus the uniform bounds (113) and (121) hold, and they imply that the integral of the difference over h is bounded by O(ǫ 3 2 ), which, even after multiplication with ℓ, still vanishes as ǫ → 0.
Details of the proof. -The difference of partition functions is In the last equality, the discrete product rule was used. Now The difference of Q's is a diagonal matrix, so (151) (F depends on ǫ and on π, ı, x, y 1 , . . . , y N , etc.). Then The sum over ı is the one mentioned in the overview. All other summations run over index sets that are independent of ℓ. To show that this expression vanishes in the limit ℓ → ∞, it will suffice to use the bound
Thus, for all r ≥ 0 there is K r > 0 so that for all ℓ and all ı, x so that this contribution vanishes in the limit ℓ → ∞ (since ǫ = β ℓ ). Because it follows, again from (113)

7.2.
Recovering the a-integral. Under the stronger condition that the real part of the quadratic form defined by Q 2 (0, u) is strictly positive, one can now go backwards and rewrite the permanent as a Gaussian integral over a, up to an additional h-dependent factor that arises from the determinant of Q 2 . Because Q 2 is upper triangular in the time slice index  and because the diagonal part is diagonal in x as well, this determinant is AssumeẼ − u > 0 and recall the definition of Da in (76). By Lemma 9 (b), Q 2 (h, u) + Q 2 (h, u) † > 0, thus the Gaussian integral given by Q 2 is absolutely convergent and Da e −(ā|Q 2 (h,u)a) X = det Q 2 (h, u) −1 (see Lemma 23). By (271), Under the hypotheses of Lemma 14, Proof. The proof is an easy generalization of that of Lemma 14, the only difference being the presence of the factor This just gives additional contributions to the discrete product rule, namely is analytic in h ı,x and of order ǫ The prefactor involvingũ = |X| −1 x u x is included for later convenience. It makes only an inessential change because its limit as ℓ → ∞ is 1. To remove the u-dependence from the term coupling h to a, change variables in the integral to a ′ ,x = e − ǫ 2 ux a ,x . (This also cancels the u-dependence of Q 2 (0, u) in the quadratic form.) The Gaussian expectation over h can now be performed and gives (dropping the primes from the notation) (174) Expanding the quadratic form of the interaction term then gives the integral representation with a positive interaction term, as follows.
The integral representations involving the variants of the kinetic term in the action stated in Theorem 1 right after (31) are obtained similarly, by using Q 3 and Q 4 instead of Q 2 . By Lemmas 10, 11 and 12, the proofs that the limit ℓ → ∞ exists and is unchanged for these modified actions are straightforward adaptations of the ones given for Q 2 .

Remarks about the Grand-Canonical Ensemble
Here I briefly outline the connection to the grand canonical ensemble. By definition (see (24)), the grand canonical partition function is the Laplace trans- . The setup of the canonical ensemble has a few technical advantages, in particular, on a lattice with nonzero spacing, the Hamiltonian, as well as its kinetic and interaction parts, become bounded operators when restricted to the N -particle space, so that one can regard fixing particle number as a natural regularization for defining bosonic models mathematically. One can then also try to study the grand canonical ensemble based on results for the canonical one. In the following, I discuss a few significant points where the ensembles differ, and how things fit together.
8.1. The periodic boundary condition in time. The integral over the a fields in the canonical ensemble does not have a periodic boundary condition in the Gaussian: there is no matrix element in Q that couples  = ℓ and  = 0. The reason for this is explicit from 51 -the projection P N removes the exponential, replacing it by the N 'th order term in its expansion. It is therefore plausible that the exponential, and hence the periodic boundary condition, will be reinstated by the sum over N in (24). But there are of course convergence issues when going beyond formal considerations. When (81) is inserted into (24) and the sum over N is formally exchanged with the limit and the integral, one gets This series converges for all values of βµ, a ℓ and a 0 . But the additional term created in the exponent spoils positivity of the real part of the quadratic form, hence leads to a divergent Gaussian integral, if µ > 0. This is explained in more detail in Appendix B.
The appearance of a factor e βµ only in the term coupling the time slices  and 0 is also unusual. The standard combination E − µ1 can be obtained simply by rearranging the product similarly to what is done in (24), which changes E to E µ = E − µ1 everywhere. That the sum over N in (24) converges for µ < 0 and E ≥ 0 is then a simple consequence of the nonnegativity of the interaction term v. If v > 0, the condition µ < 0 can be dropped.
8.2. The covariance. The covariance of the canonical ensemble is upper triangular because of the absence of a periodic boundary condition in time in the operator Q. With the just described absorption of µ into the kinetic term E µ , the summation over N gives

and (81) gets replaced by
(Here e −ǫEµ appears in Q.) As shown in Appendix B, the covariance is given for h = 0 by the standard time-ordered Green function for free bosons, eq. (228).
Theorem 17. Assume that E µ > 0. Then the integral ( (181)) is absolutely convergent. The grand-canonical covariance G(h) = K(h) −1 exists for all h, and has a norm bounded uniformly in h. It is analytic in h if |Im h ,x | < √ ǫe min , where e min is the smallest eigenvalue of E µ . If E µ generates a stochastic process, then for all  ′ ≥ , all x, x ′ ∈ X, and all h, the covariance G(h) = K(h) −1 for the grand canonical ensemble satisfies the bound Proof. This proof is an application of Lemma 20, and I will use the notations introduced in Appendix B. Let A  = e −ǫEµ e i √ ǫh . Because E µ > 0, A  = e −ǫEµ < 1 uniformly in h, and A 0 = 1, hence P = A 0 . . . A ℓ has norm P < 1. By Lemma 20, the inverse G exists, is given by (220), and has norm bounded uniformly in h ∈ R X because the same holds for all factors appearing in (220). In particular, the inverse (1 − P) −1 is given by a convergent geometric series. Expanding out this series gives a linear combination of products of A  's, with nonnegative coefficients. Thus the bound (183) follows by straightforward adaptation of the proof of Lemma 7. Analyticity in the ǫdependent neighbourhood of R X holds because the norm of e −ǫE e i √ ǫh remains strictly smaller than 1 on that set.
Thus, all arguments based on the uniform bounds for the covariance, integration by parts, etc., have their analogues for the grand canonical ensemble, and similar theorems about convergence of integrals with modified actions hold, in particular, there is a representation with an |a| 4 interaction term in the action. (I omit the statements and proofs here to avoid further repetition.) As mentioned in Section 3 for the the canonical ensemble, positivity of the interaction, i.e. v > 0, implies that the integral (27) converges absolutely also beyond µ = 0 and it defines an analytic function of µ. By the identity theorem, the continuation represents the partition function for all µ. A similar statement can be proven for the grand canonical ensemble.
8.3. The determinant. A further interesting and striking difference between the ensembles is the role played by the (inverse of the) determinant arising from the Gaussian a-integral. In the canonical case, the determinant is simply equal to one. In the grand canonical case, the determinant differs from one, and the expansion of ln det K(h) in powers of h creates the well-known expansion in terms of loops corresponding to virtual particle-antiparticle pairs with higher and higher order interactions. These virtual pair creation processes do not happen in the canonical ensemble because there is a restriction to fixed N . Thus the constraint of fixed particle number shows up explicitly as the determinant being one.
Moreover, there is no backward-in-time propagation in the canonical ensemble because the covariance is upper triangular in time. By contrast, the grand canonical covariance G(h) is not upper triangular. It is simply the time ordered boson two-point function in the external field h, hence allows forward and backward propagation. All these features make explicit that fluctuations are different in the canonical and grand canonical ensemble.
Finally, there is another, more technical point. The inverse of the determinant in the grand canonical ensemble does not factorize. The permanent in the canonical ensemble factorizes over cycles of permutations, and this makes a rather simple polymer expansion, in which the bounds on the covariances can be used directly, in analogy to the proofs given in Section 7 possible [40].

The Effect of a Condensate in the Canonical Ensemble
In his famous work [11], Bogoliubov treated condensation of the weakly interating Bose gas by an ansatz in which the creation and annihilation operators of the kinetic-energy zero mode are replaced by commuting operators, and diagonalizing the resulting effective Hamiltonian by a canonical transformation, the Bogoliubov transformation. He was thus able to derive the linear spectrum of low-energy excitations of the condensate, which is known to be experimentally accurate. Even the elements of the Bogoliubov transformation could be identified in experiments [43].
At sufficiently low temperatures β −1 , weakly interacting bosons are thus expected to condense into a state resembling the lowest-energy state of the free boson system. The latter is determined by the zero mode of the one-particle kinetic term which, in a translation-invariant system, is spatially constant. In the FIR, this corresponds to splitting the integration variable into a 'condensate field' which is integrated over a low-dimensional subspace of field space, and all modes orthogonal to it. Condensation as a phase transition then corresponds to the condensate field acquiring a nonvanishing expectation value, which spontaneously breaks the U (1) symmetry of particle number conservation. In the FIR, this is the symmetry of action and integration measure under a τ,x → e iθ a τ,x (where θ is independent of τ and x). In the following, I do this field decomposition as a change of variables in the a integral, without, for the moment, making the assumption that the zero mode field is also constant in time (which is the usual zeroth order approximation for expansions). This decomposition can be applied irrespective of what ensemble is studied, but again, I will focus on the canonical ensemble here. It also makes only a technical difference which of the discretizations discussed above one uses; for convenience of presentation, I take the action with the time-local |a| 4 term and the quadratic form Q 4 for the kinetic term.
Assume that E is given by the discrete Laplacian. Then E has a nontrivial kernel, namely the constant functions on X. By (22), E ≥ 0 can be changed to E + ε1 with ε > 0, at the expense of having an explicit factor e βN ε in front of the partition function. Assume, moreover, that v is strictly positive, v > 0. Then the integral (27) converges absolutely for any ε, even ε < 0, and consequently, the limit ε → 0 can be taken. It is thus possible to proceed by introducing ε. Because E ε > 0, Theorem 1 applies, and one can also postpone 28 taking the limit ε → 0 to a convenient point. Thus the integral is and m =v(0). For simplicity of presentation, v x,y = vδ x,y (at η = 1) is taken in the following. Thenv(0) = v as well.
9.1. Orthogonal decomposition. The fundamental equation of it all is This is an orthogonal decomposition: (c(τ ) | b(τ )) X = 0 for all τ . With this, The contribution to the action S that is linear in b vanishes because of the orthogonality of b and c. For the same reason, and a binomial expansion gives The sum over K in (192) suggests that the condensate now plays the role of a particle reservoir -at fixed c, the ensemble described by the b integral is no longer at fixed N , but rather, all particle numbers K ∈ {0, . . . , N } contribute to the partition function. To be a true reservoir in a thermodynamical sense, the condensate fluctuations would have to be so small that it becomes approximately constant also in time τ , and that the likeliest values of K are small compared to N . This is expected in the thermodynamic limit at sufficiently 29 large β and with N growing together with |X| so that ρ = N |X| is fixed and large enough. However, in the present setup, c may still receive important contributions from the b-integral, especially because of the slow decay of the b-covariance (see Section 9.3 below). 9.2. Quadratic Form for the b Fields. By definition of the time derivative and by its boundary condition (see Theorem 1 and (130)), the operator Q(c) defining the action A 2 as a quadratic form in b is still upper triangular in the time slice indices  and  ′ . Explicitly, in the notation where times τ = ǫ are labelled by  ∈ {0, . . . , ℓ}, The c-dependent matrix appearing in (195) is obviously nonnegative (but does have an eigenvalue zero). The determinant of L(c  ) is no longer equal to one. However, and this is again peculiar to the canonical ensemble, det L(c) = 1 + O(ǫ), and the matrix Q corresponding to the quadratic form A 2 in b is therefore invertible if ǫ is small enough.
9.3. Time-independent condensate field. In the canonical partition function, there is no chemical potential that can be tuned to give a mexican hat potential. Instead it is the power (51) that introduces a tendency towards a nonzero value of c. A first understanding of this can be gained by assuming that c is constant in time, c(τ ) = c(0) = c 0 for all τ , and taking the contribution from K = 0 only. Then c(β)c(0) = |c 0 | 2 . Fixing ρ = N |X| means that |X| = σN with σ = 1 ρ . Then for large N , the factor in the integrand drives the maximum of the integrand in c 0 , hence also the most probable value of c 0 , away from zero. The above replacements amount to considering the simplified function The 1 term in the square brackets is due to the boundary condition for the discretized time derivative which implies that its integral does not evaluate to zero on a constant field. Explicitly, when inserting a -independent field into (79) (with E = 0 and h = 0): the difference of the two sums does not vanish, but gives |X||c 0 | 2 = N σγ. As discussed, the logarithm in F comes from the power in the integral for the canonical partition function, and it implies that F takes its minimum at a γ 0 = 0. As mentioned, the true ensemble average is more complicated than this. In particular, ρ 0 = N −K |X| is the condensate density, and whether condensation really takes place at a given temperature is a question of the distribution of K's. 30 9.4. Positivity in the continuum limit. If c 0 = 0, the eigenvalues of the matrix Q exhibit the typical Bogoliubov spectrum, leading to a kinetic term E B (p) = |p| w 2 + p 2 where w 2 = 2v|c 0 | 2 . (Here I have taken a spatial continuum limit as well.) The Fourier transform G t of e −t(|p|+ √ w 2 +p 2 ) is the convolution of that of e −t|p| and that of e −t √ w 2 +p 2 , which are both strictly positive functions, hence strictly positive, so that for all x and y. By iterated convolution, e −τ F (x, y) is positive for all x and y, and the same holds, again by convolution, for e −τ E B (x, y).
The decay is, however, slow since the kernels that get convolved only have power law decays. Specifically [41], This slow decay suggests that a stochastic representation as a random walk in a fluctuating condensate background has long jumps, hence differs from the one corresponding to Brownian motion. The cubic term A 3 also leads to branching and coalescence, corresponding to processes where one of two scattering particles emerges from, or gets absorbed in, the condensate.

Conclusion and Outlook
The coherent-state and auxiliary-field integrals for the boson system derived here provide a simple regularized version of the standard formal time-continuum functional integrals used in condensed-matter physics. The integral for the canonical partition function and expectation values exhibits a few significant differences to the usually considered grand-canonical case: there is no periodic boundary condition in time in the quadratic form defining the Gaussian integral, propagation is only in forward time direction, and the inverse determinant from Gaussian integration is equal to one. The integrals also allow to derive a straightforward random-walk expansion, which can be viewed as a regularized version of the Brownian motion representation used in [3,1,2].
It is often stated that the Wick rule fails for the canonical expectation values of products of field operators. The integral representations for the canonical expectation values provide an easy way to understand this -in Gaussian integrals, the Wick rule is derived by integration by parts, and the factor (ā ℓ | a 0 ) X that arose from the projection to the N -particle space F (N ) B , see Lemma 4, gives rise to additional terms in that procedure. In presence of an operator insertion F, this factor gets modified, see, e.g. (62), but it remains a polynomial in the fields. It is easy to derive a modified Wick rule using the 31 standard formalism of field Laplacians for representing Gaussian expectations (see Appendix E and [37]).
At the heart of the technical results of this paper are the pointwise bounds, uniformly in h, for the regularized covariances, which are similar to those used in [34] and subsequent works on the loop vertex expansion. Although it turns out that they are very simple to prove (see Lemma 7 for an example), it was not completely obvious a priori that they would be available in the discrete-time regularization. Above, I have used them to show that there is a considerable freedom in chosing the discretized action without changing the time-continuum limit of the partition function. Although this is expected, the point here is that the proofs require only very elementary bounds, such as the combination of resolvent equations like (150), Taylor expansion like in (156) and integration by parts as in (255. Moreover, these tehniques can be used to derive convergent expansions for expectation values, which allow to take the infinite-volume limit [40].
One should not jump to the conclusion that this already makes more involved machinery, such as a decomposition in large and small fields, unnecessary in the analysis of bosonic many-body problems, however. If renormalization becomes necessary, e.g. when regularizing determinants or performing any operation in which expansions in the fields are required, it is not yet clear whether uniform bounds will apply.
The bulk of this paper has dealt with the canonical ensemble, but (as discussed in some detail in Sections 3 and 8), many essential bounds carry over to the grand-canonical case. Moreover, one can also take the approach to use results for the canonical ensemble to make statements about the grand canonical one, by summation over N , provided one has good enough control over the N -dependence.
Finally, I have given a brief discussion coming back to the original motivation, namely the connection to stochastic representations. The energy of excitations in a Bose condensate is no longer quadratic, but linear in momentum, E(p) = c|p|, with c the velocity of sound. Nevertheless, in the formal time-continuum limit, this kinetic term still generates a stochastic process (albeit with longrange jumps). It will be interesting to see if the uniform bounds derived in this paper, for which this positivity was crucial, continues to hold in a regularized version of the functional integral.

Appendix A. Coherent-State Formulas
In terms of the orthonormal basis (6), the expansion of the coherent state is which is norm convergent for any a ∈ C X by inspection. Eq. (39) then follows by inserting the expansion.
By the commutation relations, a x a † Since the exponential series for v a is norm convergent, v For all x ∈ X, the operators a x and a † x have the common dense domain D = D √ nx , and v a ∈ D. Let w ∈ D. Then To see (43), let ψ ∈ F B . Then the basis coefficients ψ ν = e ν | ψ are square summable over ν ∈ N X 0 , and the sum converges absolutely and uniformly on C |X| R , hence can be exchanged with the a-integral. Writing a x = ρ x e iθx , The integral on the right hand side of (44) factorizes into a product over x of The θ-integral vanishes unless ν x = ν ′ x . The ρ x -integral is an incomplete Gamma function Γ(ν x + 1, R 2 ), which increases to ν x ! as R → ∞. Thus Thus S R = ν |ψ ν (1 − f ν (R))| 2 converges, and and f ν (R) → 1 as R → ∞. By the dominated convergence theorem, applied to the sum over ν, S R → 0 as R → ∞. Since the e ν are an ONB, S R = ψ − P R (ψ) 2 , so (43) holds.
To see (44), use that the series expansion converges absolutely and uniformly for a ∈ C X R , because A, as a trace class operator, is bounded hence | e ν ′ | A e ν | ≤ A . So the integral over a can be exchanged with the summation. By the same arguments as above, the sum reduces to a single summation over ν, and (213) Because A is trace class, the sum ν | e ν | A e ν | converges. Thus, again by the dominated convergence theorem, (44) follows.
Eq. (45) is an immediate consequence of the following Lemma.

Appendix B. Inversion Formulas
Let B be a Banach algebra with identity 1 and A 0 , . . . , A ℓ ∈ B. Define Q and where ,  ′ ∈ {0, . . . , ℓ}. The following lemma states the formulas for the inverses (218) for the general case where the A  need not commute.
Proof. Considered as a matrix in the indices  and  ′ , Q is upper triangular, more precisely of the form identity plus a nilpotent term R. So it is invertible, and the inverse given is by a terminating Neumann series in powers of R. This gives (219). Let D , ′ = δ ,ℓ δ 0, ′ A 0 so that K = Q − D. Due to the special form of D, the resolvent equation reads for the matrix elements Because A 0 C 0,ℓ = P, the equation for  ′ = ℓ reads (223) G ,ℓ = C ,ℓ + G ,ℓ P By hypothesis, 1 − P is invertible, so Inserting (219), and noting that the indicator functions in that expression evaluate to 1 in C ,ℓ for all , and in C 0, ′ , for all  ′ , one gets (220).
A sufficient condition for 1 − P to be invertible is that P < 1, which in turn holds if A ı ≤ 1 for all ı ∈ {0, . . . , ℓ} and A ı 0 < 1 for some ı 0 . For A  = e −ǫE e i √ ǫh with real h  , a sufficient condition is that E > 0. These conditions are not necessary for invertibility of 1 − P. But they are necessary for the convergence of the coherent-state functional integral of the grandcanonical ensemble. The series converges for all values of βµ, a ℓ and a 0 . Upon formal exchange of the summation and integration in (24), the resulting term in the exponent is just the difference D of Q and K, with A 0 = e βµ 1 of norm A 0 = e βµ . If µ > 0, this is bigger than one, and invertibility of P can fail for particular X, and even if it holds, it will typically not hold with a norm of the inverse that is uniform in X. Moreover, the quadratic form does not have a positive real part, and therefore the Gaussian integral diverges, which in turn invalidates the exchange of summation and integration. Similarly, if one wants to use the complex contour integral for enforcing that the particle number is N , one should take integration radius |z| < 1 (see (63)). 35 Appendix C. Baker-Campbell-Hausdorff Formulas The Baker-Campbell-Hausdorff expansion is well-known and has been considered in many works, often with the aim to simplify the recursion for generating more and more terms in the expansion [44,27,42]. The point here is, on the other hand, to avoid combinatorics by deriving and bounding remainder terms in integral form. I assume in the following that A and B are elements of some Banach algebra, and show that with H 2 (t) = O(t 2 ) and H 3 = O(t 3 ) analytic in A and B, and both having explicit bounds in terms of the norms of A and B. The application will be with t = ǫ ∼ ℓ −1 very small. It is then of course irrelevant whether the error term is in the exponent or 'downstairs', as it is small in norm compared to 1, so one can put it in the exponent simply by taking the logarithm. Indeed it is in general an advantage if the error term is not in the exponent. The way the expansion is generated also allows to go to any higher order, with an according remainder term. Variants of the Lie product formula including higher order terms in the exponent, and with accordingly improved convergence rate as ℓ → ∞, can be proven straightforwardly using these results. They are useful for generalizations of the systems discussed here, e.g. for certain Fermi-Bose mixtures.
Let Obviously this iteration can be continued to arbitrary orders and it gets more and more tedious as one goes on. But at each step one has an integral formula for the remainder which is useful for bounds.
By Taylor expansion of the exponential, ∆ ǫ (h) = − ǫ 2 h 2 + O(ǫ 3 2 h 3 ). The following elementary statements are used in the sequel: Proof. Call the left hand side of (258) I q (a, b). Because Re q > 0 this integral exists, and it is absolutely convergent, hence can be taken as the limit R → ∞ of the integral over C R = {z ∈ C : |z| ≤ R}, which will be denoted by I R q (a, b). At fixed R, expand inā and b. The compactness of the integration region justifies the exchange of summations and integral. This gives where γ 1 is the straight line from 0 to qL ∈ C. The integrand is an entire function of s, so the contour can be deformed to γ 2 + γ 3 , where γ 2 is the straight line from 0 to LRe q and γ 3 the straight line parallel to the imaginary axis from LRe q to Lq. The latter can be parametrized by t = Lu + iLvr, u = Re q, v = Im q and r ∈ [0, 1]. Then  Here a | b =ā ⊤ b denotes the standard inner product on C N .
Proof. The proof will by reduction to Lemma 22, using a diagonalization argument. Let Q 0 = 1 2 (Q+Q † ) and H = 1 2i (Q−Q † ) then Q = Q 0 +iH, H = H † , Q 0 = Q † 0 > 0. Thus Q satisfies the hypotheses of Lemma 11, and the same decomposition as in the proof of that lemma will be used here: let B = Q The Jacobian from the transformation with U is one because the determinants U and U † appear and multiply to 1. The w-integral factorizes because D is diagonal. The real part of every eigenvalue of D equals 1, so Lemma 22 applies. The product of eigenvalues is det D −1 , which combines with det B −2 to det Q −1 . The product of exponentials gives a quadratic form in the exponent as so (263) follows.
Permanent formulas, such as (67) and (92) As in Section 2.3 of [37], the generating function identity can be used to rewrite the Gaussian integral of polynomials in terms of the field space Laplacian. The result is stated in the following Lemma. (The operator Q in that lemma may be h-dependent.) Lemma 24. Let Q be a linear operator on C X and Q + Q † > 0. Then where C = Q −1 and With the notation δ δa,x = ǫ −1 η −d ∂ ∂a,x , the Laplacian can be rewritten as (270) ∆ C = ( δ δa | C δ δā ) X = (C δ δā | δ δa ) X Techniques using these Laplacians will be useful in doing convergent expansions, in analogy to [38].