text:
# Untitled Document

THE PHASE OF THE SCATTERING MATRIX

José M. Gracia-Bondía1,2

1Department of Theoretical Physics I, University Complutense, Madrid 28040, SPAIN

2Department of Theoretical Physics, University of Zaragoza, Zaragoza 50009, SPAIN

Vacuum polarization in external fields is treated by way of calculating —exactly and then perturbatively— the phase of the quantum scattering matrix in the Shale–Stinespring approach to field theory. The link between the Shale–Stinespring method and the Epstein–Glaser renormalization procedure is highlighted.

1. Introduction

Renormalization theory is bound to suffer a reappraisal in the light of the reconstruction of Zimmermann’s forest formula in Hopf-algebraic terms, together with the interpretation of the dimensional regularization method given by Connes and Kreimer , and the birth of quantum field theory on noncommutative spaces , together with the evidence that Yang–Mills theories on noncommutative manifolds are ultraviolet divergent —see  and the other references in that paper.

Now, part of the advantages and new popularity (see, for instance, \[5–8\]) of the Epstein–Glaser renormalization method \[ 9,10\] stems from the fact that it is locally defined, and so in principle applies to models on nonflat manifolds. There is however some contention on whether, as claimed by some practitioners , the Epstein–Glaser method is a fundamental one.

In cases like these, it sometimes helps to look at a simpler problem, for which an absolutely reliable method is known, to see how the marketed procedures fare in its respect.

The chosen problem is that of vacuum polarization in external fields and the chosen reliable method is the Shale–Stinespring approach to linear quantum field theories. This is an entirely rigorous algebraic method; in  its prima facie applicability to “implementable” theories on noncommutative as well as commutative spaces was brought to the fore.

After reviewing the Shale-Stinespring theorem in Section 2, in the body of this article we show, by a refinement of its technique, that the phase of the scattering matrix is well defined and finite for implementable linear theories. We give explicit formulae for the phase.

We then exploit in QED a perturbative version of this approach, which leads in a direct way to the formulas reached by Scharf  in his account of vacuum polarization in QED by the Epstein–Glaser method. Some simplification of his calculations results from using gauge-invariant variables. The contention by Scharf and followers that the Epstein–Glaser renormalization procedure is a fundamental one is vindicated to some extent.

2. A reminder on Shale–Stinespring theory

Let $`=^+^{}=:P_+P_{}`$ be a Hilbert space, graded by the projections $`P_\pm `$ on the positive and negative spectral subspaces for a free Dirac operator. Operators on $``$ are presented in block form:

$$A=\left(\begin{array}{cc}A_{++}& A_+\\ A_+& A_{}\end{array}\right).$$

We have in mind particularly the classical (or “first-quantized”) scattering matrix

$$S=\left(\begin{array}{cc}S_{++}& S_+\\ S_+& S_{}\end{array}\right).$$

Introduce a nomenclature for its even and odd parts:

$$p_S:=S_{\mathrm{even}}=\left(\begin{array}{cc}S_{++}& 0\\ 0& S_{}\end{array}\right),q_S:=S_{\mathrm{odd}}=\left(\begin{array}{cc}0& S_+\\ S_+& 0\end{array}\right).$$

Similarly,

$$S^1=S^{}=\left(\begin{array}{cc}S_{++}^{}& S_+^{}\\ S_+^{}& S_{}^{}\end{array}\right),p_S^{}=\left(\begin{array}{cc}S_{++}^{}& 0\\ 0& S_{}^{}\end{array}\right),q_S^{}=\left(\begin{array}{cc}0& S_+^{}\\ S_+^{}& 0\end{array}\right).$$

Unitarity of $`S`$ gives the identities

$$\begin{array}{cc}\hfill S_{++}S_{++}^{}+S_+S_+^{}& =S_{++}^{}S_{++}+S_+^{}S_+=P_+,\hfill \\ \hfill S_{}S_{}^{}+S_+S_+^{}& =S_{}^{}S_{}+S_+^{}S_+=P_{},\hfill \\ \hfill S_{++}S_+^{}+S_+S_{}^{}& =S_{++}^{}S_++S_+^{}S_{}=0,\hfill \\ \hfill S_{}S_+^{}+S_+S_{++}^{}& =S_{}^{}S_++S_+^{}S_{++}=0.\hfill \end{array}$$

It is clear that $`p_S^1`$ exists if and only if $`S_{++}`$ and $`S_{}`$ are invertible as operators on $`^+`$ and on $`^{}`$, respectively; this is the generic case, that always holds when $`S`$ is close to the identity, and will be the only one considered in the sequel. We then define the skewadjoint operators

$$T_S:=q_Sp_S^1=\left(\begin{array}{cc}0& S_+S_{}^1\\ S_+S_{++}^1& 0\end{array}\right),\widehat{T}_S:=T_S^{}=\left(\begin{array}{cc}0& S_{++}^1S_+\\ S_{}^1S_+& 0\end{array}\right).$$

Consider the Fock space constructed on $``$ with the new scalar product

$$\eta |\phi :=\eta _+|\phi _++\phi _{}|\eta _{},$$
$`(1)`$
where $`\eta _\pm :=P_\pm \eta `$. Let $`\{\varphi _k\}`$ and $`\{\psi _k\}`$ denote arbitrary orthonormal bases for $`^+`$ and $`^{}`$, respectively; we shall abbreviate $`b_k:=b(\varphi _k)`$, $`d_k:=d(\psi _k)`$ in the notation of the “particle” and “antiparticle” annihilation operators, and similarly for the creation operators $`b_k^{}`$, $`d_k^{}`$. For any operator $`A`$ on $``$ we have the quantum (or “second-quantized”) counterpart, acting on Fock space:

$$d\mathrm{\Lambda }(A):=b^{}A_{++}b+b^{}A_+d^{}+dA_+b+:dA_{}d^{}:.$$

Here, for instance, $`:dA_{}d^{}:`$ denotes $`_{j,k}d_k^{}\psi _j|A_{}\psi _kd_j`$, the double colon meaning, as usual, a normally ordered product, and $`b^{}A_+d^{}`$ is $`_{j,k}b_k^{}\varphi _k|A_+\psi _jd_j^{}`$; the other cases should be clear.

This rule corresponds to the infinitesimal spin representation ; it is independent of the drafted orthonormal bases and makes sense only when $`A_+`$, $`A_+`$ are Hilbert–Schmidt.

The rule is mainly applied to selfadjoint operators, and yields (at least formally) selfadjoint operators in turn. For instance, in QED the free Dirac equation is written as

$$i\frac{}{t}\psi =\beta m\psi i\stackrel{}{\alpha }\frac{}{\stackrel{}{x}}\psi =:D_0\psi ,$$

where the Dirac matrices, say in the chiral representation, are given by

$$\stackrel{}{\alpha }:=\gamma ^0\stackrel{}{\gamma }=\left(\begin{array}{cc}\stackrel{}{\sigma }& 0\\ 0& \stackrel{}{\sigma }\end{array}\right),\beta :=\gamma ^0:=\left(\begin{array}{cc}0& 1_2\\ 1_2& 0\end{array}\right).$$

To this classical selfadjoint operator $`D_0`$ corresponds the quantum free Hamiltonian

$$𝐇_0:=d\mathrm{\Lambda }(D_0)=b^{}D_0b^{}+:dD_0d^{}:,$$

which is a positive operator. The classical interaction Hamiltonian is of the form

$$H(t)=e(A^0(t)\stackrel{}{\alpha }\stackrel{}{A}(t)),$$

where $`e`$ denotes the electromagnetic coupling constant and $`(A^0,\stackrel{}{A})=:A`$ the electromagnetic vector potential, a real $`c`$-number function. Note the covariant form $`H(t)=e\gamma ^0A/`$, with $`A/:=\gamma ^\mu A_\mu `$. In the interaction picture one considers

$$V(t):=e^{iD_0t}H(t)e^{iD_0t}.$$
$`(2)`$
The quantum interaction Hamiltonian is then

$$𝐕(t):=d\mathrm{\Lambda }(V(t))=b^{}V_{++}(t)b^{}+b^{}V_+(t)d^{}+d^{}V_+(t)b^{}+:dV_{}(t)d^{}:.$$
$`(3a)`$
This can be rewritten in terms of the formal fermion field $`\mathrm{\Psi }`$, as

$$𝐕(t)=d^3x:\overline{\mathrm{\Psi }}(x)A/(x)\mathrm{\Psi }(x):,$$
$`(3b)`$
with $`x=(t,\stackrel{}{x})`$ and the bar meaning the Dirac adjoint. For that, just write the fermion field in the form

$$\mathrm{\Psi }(x)=\underset{k}{}(b_k\varphi _k(x)+d_k^{}\psi _k(x)),$$

where $`\varphi _k(x)=e^{iD_0t}\varphi _k(\stackrel{}{x})`$, and similarly for the $`\psi `$’s.

For the quantum scattering matrix $`𝐒`$ we need instead the global spin representation , which we call $`\mathrm{\Lambda }`$. It is given by

$$𝐒:=e^{i\theta }\mathrm{\Lambda }(S)=e^{i\theta }|0_{\mathrm{in}}|0_{\mathrm{out}}|:\mathrm{exp}d\mathrm{\Lambda }(I):=0_{\mathrm{in}}|0_{\mathrm{out}}:\mathrm{exp}d\mathrm{\Lambda }(I):.$$

Here $`0_{\mathrm{in}}`$ denotes the incoming vacuum, $`0_{\mathrm{out}}:=𝐒\mathrm{\hspace{0.17em}0}_{\mathrm{in}}`$, and

$$I:=\left(\begin{array}{cc}(S_{++}^{})^11& S_+^{}S_{}^1\\ S_{}^1S_+^{}& 1S_{}^1\end{array}\right).$$

Again, this makes sense if and only if $`S_+`$, $`S_+`$ are Hilbert–Schmidt, and then we say that $`S`$ is implementable. The absolute value of the vacuum persistence amplitude $`0_{\mathrm{in}}|0_{\mathrm{out}}`$ is given by

$$\begin{array}{cc}\hfill |0_{\mathrm{in}}|0_{\mathrm{out}}|& =det^{1/4}(1T_S^2)=det^{1/2}(S_{}^{}S_{}^{})=det^{1/2}(S_{++}^{}S_{++}^{})\hfill \\ & =det^{1/2}(1S_+^{}S_+^{})=det^{1/2}(1S_+^{}S_+^{}).\hfill \end{array}$$

For our present purposes, this is the content of the Shale–Stinespring theorem . The phase $`\theta `$ is in principle undetermined and conventionally taken equal to zero; this is all that is needed to compute transition probabilities.

The phase of the vacuum persistence amplitude does matter physically, however: the current density is modified with respect to the free field situation by the vacuum polarization effect (that bears on the radiative correction to the photon propagator in the nonlinear theory), and the interacting current density is found by functional derivation of $`𝐒`$ with respect to the gauge potential, in which the phase intervenes .

The question is then to find an appropriate and computable definition for the phase of the quantum scattering matrix.

3. Computing the phase in the Shale–Stinespring framework

The difficulty comes from the fact that the global spin representation is projective. Let $`U(s,t)`$ be the classical unitary propagator in the interaction representation, which interpolates between the identity and $`S`$. Then $`U(s,t)`$ solves the equation

$$U(s,t)=1i_t^sV(u)U(u,t)𝑑u;$$
$`(4)`$
an explicit form being given by the Dyson expansion

$$U(s,t)=1+\underset{n=1}{\overset{\mathrm{}}{}}(i)^n_t^sV(t_1)_t^{t_1}V(t_2)\mathrm{}_t^{t_{n1}}V(t_n)𝑑t_n\mathrm{}𝑑t_2𝑑t_1,$$

from which follow the propagator properties:

$$U(t,t)=1,U(t,s)U(s,r)=U(t,r).$$

Here $`S=U(\mathrm{},\mathrm{})`$.

Now, for unitary implementable operators $`U_1`$, $`U_2`$, in general

$$\mathrm{\Lambda }(U_1)\mathrm{\Lambda }(U_2)=c(U_1,U_2)\mathrm{\Lambda }(U_1U_2),$$

where the cocycle $`c`$ (with the vanishing phase convention) is given  by

$$c(U_1,U_2)=\mathrm{exp}\left(i\mathrm{arg}det^{1/2}(1T_{U_2}\widehat{T}_{U_1})\right)=\mathrm{exp}\left(i\mathrm{arg}det^{1/2}(p_{U_1}^1p_{U_1U_2}p_{U_2}^1)\right).$$
$`(5)`$
This is not a trivial cocycle because the determinants of the $`p_U`$ operators are not individually defined in general. Assume that the interpolating family $`U(s,t)`$ is implementable and (strongly) differentiable with respect to its parameters —this will happen if the external field is sufficiently well behaved. Then

$$\mathrm{\Lambda }(U(s,t))\mathrm{\Lambda }(U(t,r))=c(s,t,r)\mathrm{\Lambda }(U(s,r)),$$
$`(6)`$
with an obvious notation. Still, $`c(t,t,r)=c(s,t,t)=1`$. On the other hand, it can be shown  that

$$i\frac{}{s}|_{s=t}\mathrm{\Lambda }(U(s,t))=𝐕(t).$$

We now seek to redefine $`\mathrm{\Lambda }(U(s,t))`$ by multiplying a phase factor $`e^{i\theta (s,t)}`$ so that the new quantum family $`𝐔(s,t):=e^{i\theta (s,t)}\mathrm{\Lambda }(U(s,t))`$ fulfils

$$𝐔(t,t)=1;𝐔(s,t)𝐔(t,r)=𝐔(s,r),$$
$`(7)`$
just like the classical propagator. If we manage that, then $`\theta (+\mathrm{},\mathrm{})`$ will have every right to be called the phase of the quantum scattering operator $`𝐒`$. Let $`c(s,t,r)=:\mathrm{exp}(i\xi (s,t,r))`$. Differentiating equation (6), one gets:

$$𝐕(t)\mathrm{\Lambda }(U(t,r))=i\frac{}{s}|_{s=t}c(s,t,r)\mathrm{\Lambda }(U(t,r))+i\frac{}{t}\mathrm{\Lambda }(U(t,r)).$$

Then, we redefine

$$𝐔(s,t):=\mathrm{exp}\left(i_t^s\frac{}{\lambda }|_{\lambda =\tau }\xi (\lambda ,\tau ,t)d\tau \right)\mathrm{\Lambda }(U(s,t)),$$

which clearly satisfies

$$𝐔(s,t)=1i_t^s𝐕(u)𝐔(u,s)𝑑u.$$

This equation is the quantized version of (4) and sports the same kind of iteration solution:

$$𝐔(s,t)=1+\underset{n=1}{\overset{\mathrm{}}{}}(i)^n_t^s𝐕(t_1)_t^{t_1}𝐕(t_2)\mathrm{}_t^{t_{n1}}𝐕(t_n)𝑑t_n\mathrm{}𝑑t_2𝑑t_1.$$
$`(8)`$
In the present case, the quantum Dyson expansion is rigorous: although $`𝐕`$ is an unbounded operator, it is a pretty tame one. Let $`E_m`$ denote the projector on states containing at most $`m`$ particles. Then, in view of (3), $`𝐕E_m`$ is a bounded operator from $`E_m`$ into $`E_{m+2}`$, and the norms of

$$𝐔(s,t)_{:n}E_m:=_s^t𝐕(t_1)_s^{t_1}𝐕(t_2)\mathrm{}_s^{t_{n1}}𝐕(t_n)𝑑t_n\mathrm{}𝑑t_2𝑑t_1E_m$$

can be easily estimated. To see that, one introduces the norm

$$|||V(t)|||:=V_{\mathrm{even}}(t)+V_{\mathrm{odd}}(t)_2,$$

where the latter is the Hilbert–Schmidt norm. By continuity and uniform boundedness, $`|||V(t)|||a(s,r)`$ holds for some finite function $`a(s,r)`$, when $`rts`$, and is not difficult to check that there are constants $`C_m`$ such that

$$𝐕(t)E_mC_ma(s,r).$$

Consult  for precise analyses of these bounds (encompassing also the boson case). On the other hand, from the integral equation (4),

$$|||U(s,r)|||1+(sr)a(s,r).$$

Putting both inequalities together, one gets the estimate

$$𝐔(s,t)_{:n}E_mC_mC_{m+2}\mathrm{}C_{m+2n2}\frac{(st)^na(s,t)^n}{n!},$$

with the result that the series in (8) indeed converges to a unitary operator, for $`st`$ small enough. Equation (7) then follows from (8) and allow us to extend the validity of the last conclusion—and, in turn, its own domain of validity. More detail on this is found in the important paper .

We can give an exact formula for the phase now. Use the standard identities

$$\frac{d}{dt}(\mathrm{arg}z(t))=\mathrm{}\frac{d}{dt}(\mathrm{log}z(t)),\frac{d}{dt}(\mathrm{log}detA(t))=Tr\left(A(t)^1\frac{dA(t)}{dt}\right),$$

which give, from (5),

$$\frac{}{\lambda }|_{\lambda =\tau }\xi (\lambda ,\tau ,t)=\frac{1}{2}\mathrm{}Tr\left(T_{U(\tau ,t)}\frac{}{\lambda }|_{\lambda =\tau }\widehat{T}_{U(\lambda ,\tau )}\right).$$

Therefore,

$$\theta (s,t)=\frac{1}{2}\mathrm{}_t^sTr\left(T_{U(\tau ,t)}\frac{}{\lambda }|_{\lambda =\tau }\widehat{T}_{U(\lambda ,\tau )}\right)d\tau =\frac{i}{4}_t^sTr[\frac{}{\lambda }|_{\lambda =\tau }\widehat{T}_{U(\lambda ,\tau )},T_{U(\tau ,t)}]d\tau .$$

(The last expression is more symmetrical; the trace of this commutator is not zero, because it is taken in Fock space, whereupon, in view of the form of the scalar product (1),

$$Tr\left(\begin{array}{cc}A_{++}& A_+\\ A_+& A_{}\end{array}\right)=TrA_{++}+TrA_{}^{}.)$$

The phase of the scattering matrix is then

$$\theta =\frac{1}{2}\mathrm{}_{\mathrm{}}^{\mathrm{}}Tr\left(T_{U(\tau ,\mathrm{})}\frac{}{\lambda }|_{\lambda =\tau }\widehat{T}_{U(\lambda ,\tau )}\right)d\tau .$$
$`(9)`$
The analogous formula, with the same notation, for the boson case was first given, to the best of our knowledge, by Várilly and the author ; it differs only by a sign. Then, equivalent formulae both for the boson and fermion cases were found by Langmann . The latter apply to charged fields, which are the ones considered in this paper. However, under the form (9) and with a suitable interpretation, the phase formula is applicable to Majorana fields, which are more general than charged fields .

Also note, before continuing, that

$$\frac{}{s}|_{s=t}\theta (s,t)=\frac{}{\lambda }|_{\lambda =s}\xi (\lambda ,s,s)=0.$$
$`(10)`$
In other words, there is no contribution from the coincidence points of $`\widehat{T}_{U(\lambda ,\tau )}`$ and $`T_{U(\tau ,t)}`$. This will prove to be the crucial remark.

One has simply $`/\lambda |_{\lambda =\tau }\widehat{T}_{U(\lambda ,\tau )}=iV_{\mathrm{odd}}(\tau )`$ in our present framework. Therefore, on calling $`T(\tau ,t):=T_{U(\tau ,t)}`$, finally:

$$\theta (s,t)=\frac{1}{2}_t^sTr(V_+(\tau )T_+(\tau ,t)T_+(\tau ,t)V_+(\tau ))d\tau ,$$

a rather elegant expression. Note that it differs from zero only at second order in perturbation theory. At that order,

$$\theta _{:2}(s,t)=\frac{1}{2}_t^sTr(V_+(\tau )U_{+:1}(\tau ,t)U_{+:1}(\tau ,t)V_+(\tau ))d\tau $$

with an obvious notation. Since $`U_{:1}(\tau ,t)=i_\tau ^tV(\tau )𝑑\tau `$, we set out to compute

$$\theta _{:2}(s,t)=\frac{i}{2}_s^tTr\left[V_+(\tau )\left(_\tau ^tV(\lambda )𝑑\lambda \right)\left(_\tau ^tV(\lambda )𝑑\lambda \right)V_+(\tau )\right]d\tau .$$

Of course, (10) still applies at this approximation. The total phase $`\theta _{:2}:=\theta _{:2}(\mathrm{},\mathrm{})`$ at this approximation is then

$$\theta _{:2}=\frac{i}{2}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\theta (t_1t_2)(V_+(t_1)V_+(t_2)V_+(t_2)V_+(t_1))𝑑t_1𝑑t_2.$$
$`(11a)`$
It should be clear that, at the same order of approximation, this is precisely

$$\frac{1}{2}\mathrm{}<0_{\mathrm{in}}|_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}T[𝐕(t_1)𝐕(t_2)]dt_1dt_20_{\mathrm{in}}>,$$
$`(11b)`$
where $`T`$ denotes the time-ordered product.

4. Perturbative calculation of the phase in quantum electrodynamics

In QED, from (2), or working like in the derivation of (3b), the integrals (11) are recast as

$$\begin{array}{ccc}\hfill \theta _{:2}=\mathrm{}\frac{e^2}{2}& \theta (t_1t_2)tr[A/(x_1)S^{}(x_1x_2)A/(x_2)S^+(x_2x_1)\hfill & \\ & A/(x_1)S^+(x_1x_2)A/(x_2)S^{}(x_2x_1)]d^4x_1d^4x_2.\hfill & (12a)\hfill \end{array}$$

Here $`S^\pm `$ denote the Wightman “functions”.

The first thing to remark is that, since

$$S^{}(x)\gamma ^\nu S^+(x)S^+(x)\gamma ^\nu S^{}(x)=S_{JP}(x)\gamma ^\nu S^+(x)S^+(x)\gamma ^\nu S_{JP}(x),$$

and $`S_{JP}`$ has support inside the lightcone, then the integrand has support inside the lightcone. That allows one to substitute for $`\theta (t_1t_2)`$ the covariant expression $`\theta ((v(x_1x_2)))=:\chi (x_1x_2)`$, where $`v`$ is an arbitrary timelike vector, which can thus be varied at will, and a parenthesis has been used to denote the Minkowski product $`g_{\mu \nu }y^\mu x^\nu =:(yx)`$ of two four-vectors $`y,x`$. Let then

$$\begin{array}{cc}\hfill \stackrel{~}{F}^{\mu \nu }(x)& :=tr[\gamma ^\mu S^{}(x)\gamma ^\nu S^+(x)\gamma ^\mu S^+(x)\gamma ^\nu S^{}(x)],\hfill \\ \hfill F^{\mu \nu }(x)& :=\chi (x)\stackrel{~}{F}^{\mu \nu }(x).\hfill \end{array}$$

The second thing to remark is that the last expression —and hence (12a)— is only formal: $`\chi \stackrel{~}{F}^{\mu \nu }`$ is actually undefined as a product of distributions in view of the singularities of $`\stackrel{~}{F}^{\mu \nu }`$ on the lightcone. Because the apparent trouble occurs at the coincidence points and since $`\chi (x)\stackrel{~}{F}^{\mu \nu }(x)`$ makes sense for $`x0`$, one can try to define $`\chi \stackrel{~}{F}^{\mu \nu }`$ as a distributional extension —or “regularization” in the terminology of — of the latter. The scaling degree  or singular order  of the integral of the product of the two Wightman functions is 2; therefore, distinct extensions of this quantity will differ by linear combinations of the delta function at the origin and its derivatives up to order two —i.e., by polynomials in $`k`$ of degree at most two in momentum space. The procedure is undoubtedly sound in the present case, as we know a priori the phase to be finite. Moreover, our framework will allow to select the “good” extension.

It is indeed convenient to work in momentum space. There,

$$\theta _{:2}=\frac{e^2}{2}(2\pi )^2\mathrm{}F^{\nu \mu }(k)A_\mu (k)\overline{A}_\nu (k)d^4k,$$
$`(12b)`$
where we take into account that $`A(k)=\overline{A}(k)`$, with the bar meaning here complex conjugation, because $`A(x)`$ is real; and, formally, $`F^{\mu \nu }(k)=(2\pi )^{1/2}\chi \stackrel{~}{F}^{\mu \nu }(k)`$ with $``$ denoting ordinary convolution. We recall that, for timelike $`k`$ and when choosing a frame in which $`k=(k^0,\stackrel{}{0})`$,

$$\chi (k)=\frac{i\delta (\stackrel{}{k})}{\sqrt{2\pi }(k^0i\epsilon )}=\frac{i\delta (\stackrel{}{k})}{\sqrt{2\pi }}\left(P\frac{1}{k^0}+i\pi \delta (k^0)\right).$$
$`(13)`$
To compute $`\stackrel{~}{F}^{\mu \nu }(k)`$, one looks at the Fourier transform of $`tr[\gamma ^\mu S^+(x)\gamma ^\nu S^{}(x)]`$. By using the well known expressions of the Wightman functions in momentum space, this is expressed as

$$\begin{array}{ccc}\hfill \frac{1}{(2\pi )^4}& tr[\gamma ^\mu (p/+m)\gamma ^\nu (q/m)]\theta (p^0)\theta (q^0)\delta (p^2m^2)\hfill & \\ & \times \delta (q^2m^2)\delta ^4(kpq)d^4qd^4p=:\frac{1}{(2\pi )^4}T^{\nu \mu }(k).\hfill & (14)\hfill \end{array}$$

Moreover, $`tr[(p/+m)\gamma ^\mu (q/m)\gamma ^\nu ]=4(p^\mu q^\nu +q^\mu p^\nu ((pq)m^2)g^{\mu \nu })`$. Then one of the integrations in (14) is immediately disposed of, with the help of the $`\delta ^4`$-function. The other is easily performed, with the help of the remaining $`\delta `$-functions, again by choosing a frame in which $`k=(k^0,\stackrel{}{0})`$ so that $`\stackrel{~}{F}^{\mu \nu }`$ can be regarded as a function of only one variable; we obtain

$$T^{\mu \nu }(k)=\frac{2\pi }{3}\left(\frac{k^\mu k^\nu }{k^2}g^{\mu \nu }\right)\left[k^2(1+\gamma (k^2))(12\gamma (k^2))^{1/2}\theta (12\gamma (k^2))\theta (k^0)\right],$$

where $`\gamma (k^2):=2m^2/k^2`$. All this is found in many books .

Thus we have been led formally to compute $`\chi \stackrel{~}{F}^{\mu \nu }`$, where $`\stackrel{~}{F}^{\mu \nu }(k)`$ equals

$$\frac{1}{3(2\pi )^3}\left(\frac{k^\mu k^\nu }{k^2}g^{\mu \nu }\right)\left[k^2(1+\gamma (k^2))(12\gamma (k^2))^{1/2}\theta (12\gamma (k^2))sign(k^0)\right].$$

That indeed behaves as a polynomial of degree two at high momentum transfer, confirming that the singular degree of $`F^{\mu \nu }`$ is two.

The correct (unique) recipe to regularize the imaginary part of $`\chi \stackrel{~}{F}^{\mu \nu }`$ is selected by prescribing that the result $`F^{\mu \nu }`$ vanishes, together with derivatives up to order two, at zero momentum. This kills the delta function at the origin and its derivatives in configuration space, which otherwise would give a nonzero contribution to $`/s|_{s=t}\theta (s,t)`$, contradicting (10). It is clear now, from the $`\delta `$-function in (13), that

$$\mathrm{}F^{\mu \nu }(k)=\frac{1}{2}\stackrel{~}{F}^{\mu \nu }(k),$$

whereas the relation between the real and imaginary parts of $`F^{\mu \nu }`$ is then given by a subtracted (at the origin) dispersion relation

$$\mathrm{}F^{\mu \nu }(k^0)=\frac{(k^0)^3}{\pi }\mathrm{P}_{\mathrm{}}^{\mathrm{}}\frac{\mathrm{}F^{\mu \nu }(\zeta )}{\zeta ^3(\zeta k^0)}𝑑\zeta .$$
$`(15)`$
The prescription that $`F^{\mu \nu }`$ possess a zero of the third (indeed, fourth) order at $`k=0`$ can be independently justified by an argument that, although heuristic, we deem very strong in the present context. On invoking the Maxwell equations (in the Lorentz gauge) $`A^\mu (k)=ȷ^\mu (k)/k^2`$ to conjure up the source $`ȷ`$ of the classical field, and on introducing

$$G(k):=\frac{1}{k^2}(1+\gamma (k^2))(12\gamma (k^2))^{1/2}\theta (12\gamma (k^2))sign(k^0),$$

one gets

$$\frac{e^2}{2}(2\pi )^2\mathrm{}F^{\mu \nu }(k)A_\mu (k)\overline{A}_\nu (k)d^4k=\frac{e^2}{24\pi }(ȷ(k)\overline{ȷ}(k))G(k)d^4k.$$
$`(16)`$
The continuity equation $`(ȷ(k)k)=0`$ has been employed to simplify the result. This simple expression exhibits only gauge-invariant variables.

Now, we remark that (classical) gauge transformations are in general not implementable in $`1+3`$ dimensions. This is a very good indicator of the existence of ultraviolet divergences in the nonlinear theory, and indeed it was used by J. C. Várilly and the author in  to point out that QFT theories on noncommutative manifolds had to be ultraviolet divergent. On the other hand, in the linear theory selfinteraction is absent, so we would not expect ultraviolet divergences on physical grounds. That nonimplementability is the only source of spurious divergence difficulties. We therefore expect to be able to express the phase in terms of gauge-invariant variables, in a similar way to (16):

$$\theta _{:2}=\frac{e^2}{24\pi }(ȷ(k)\overline{ȷ}(k))H(k)d^4k,$$
$`(17a)`$
with $`H(k)`$ is regular at $`k=0`$; this is equivalent to $`F^{\mu \nu }`$ having the aforementioned behaviour at the origin in momentum space.

Taking into account that $`G`$ is odd, equation (15) leads immediately to a simple form for $`H`$:

$$H(k)=\frac{1}{\pi }\mathrm{P}_{4m^2}^{\mathrm{}}\frac{(1+\gamma (\lambda ))(12\gamma (\lambda ))^{1/2}}{\lambda (\lambda k^2)}𝑑\lambda .$$
$`(17b)`$
The restriction to timelike $`k`$ can be removed by analytic continuation. Making the change of variable $`\lambda =:4m^2/(1v^2)`$, we get

$$H(k)=\frac{3}{\pi k^2}_0^1\frac{v^2v^4/3}{v^21+4m^2/k^2}𝑑v,$$

which is essentially the expression one finds in textbooks \[24, pp. 249–252\] for the vacuum polarization functional, after renormalization. The last integral can be easily carried out analytically, and it is then an instructive exercise to check that the function $`H`$ is perfectly smooth at $`k=0`$ (at $`k^2=4m^2`$, the onset of the absorptive part, $`H`$ has a cusp). We shall not go into the details.

Before rushing to the conclusions, a comment is in order: we have more or less treated the $`A(x)`$ as test functions, guaranteeing implementability of the interpolating operators, for the sake of the argument. However, it is clear that the final formula for the phase (9) is acceptable with only the milder requirement of the implementability of the scattering operator; consult  for a very efficient removal of technical conditions on the potentials, for this last purpose.

5. Conclusions

In this paper we have performed what amounts to an ab initio finite calculation of the “bubble” diagrams in linear quantum field theory.

Now, at first significant order in QED this is essentially the same as the one-loop vacuum polarization or “photon self-energy” diagram (see in this respect \[26, pp. 195–196\]). The computation done here does not appear to have been pushed to that finish line before now, although the tools have been there since the seventies at least . (Besides Várilly and the author , Langmann and Mickelsson  came close in the nineties.) The main point is that the “local causality condition” (10) selects the correct prescription among all the (finite) regularizations, with recourse to neither heuristic arguments  nor the extremely long and complicated “nonperturbative proof” in .

On the other hand, it will not have escaped the reader’s attention that, in order to avoid pitfalls, we reorganize the calculation in the same way as . The whole procedure is thus in the spirit of the Epstein–Glaser renormalization procedure, where there are Feynman graphs, but the Feynman rules do not necessarily apply —we avoided rewriting (12a) in terms of Feynman propagators. It is remarkable that, in the boson case, the quantum scattering matrix was found long ago by the Epstein–Glaser method by Bellissard  —without the phase.

From formulae (16) and (17) one can easily verify a posteriori Bogoliubov’s causality condition

$$\frac{\delta }{\delta A(x_1)}<0_{\mathrm{in}}|𝐒^{}\frac{\delta 𝐒}{\delta A(x_2)}0_{\mathrm{in}}>=0\text{for}x_1^0>x_2^0,$$

which was the starting point of Epstein and Glaser.

Acknowledgments

I am happy to thank Ricardo Estrada, Héctor Figueroa, Steve Fulling, Carmelo P. Martín and Joe Várilly for discussions, and my students Esteban Araya and Pablo Blanco for help with drawing $`H`$. The warm hospitality of the departments of theoretical physics at the Complutense in Madrid and at Zaragoza, and the financial support of Sabáticos UCM 1999 are most gratefully acknowledged.

References

 A. Connes and D. Kreimer, “Renormalization in quantum field theory and the Riemann–Hilbert problem I: the Hopf algebra structure of graphs and the main theorem”, hep-th/9912092, IHES, 1999.

 T. Filk, “Divergences in a field theory on quantum space”, Phys. Lett. B 376 (1996), 53–58.

 J. C. Várilly and J. M. Gracia-Bondía, “On the ultraviolet behaviour of quantum fields over noncommutative manifolds”, Int. J. Mod. Phys. A 14 (1999), 1305–1323.

 S. Minwalla, M. Van Raamsdonk and N. Seiberg, “Noncommutative perturbative dynamics”, hep-th/9912072, Princeton, 1999.

 R. Brunetti and K. Fredenhagen, “Microlocal analysis and interacting quantum field theories: renormalization on physical backgrounds”, math-ph/9903028, DESY, Hamburg, 1999.

 D. R. Grigore, “On the uniqueness of the non-abelian gauge theories in the Epstein–Glaser approach to renormalization theory”, hep-th/9806244, Bucharest, 1998.

 T. Hurth and K. Skenderis, “The quantum Noether condition in terms of interacting fields”, hep-th/9811231, 1998.

 G. Pinter, “The action principle in Epstein–Glaser renormalization and renormalization of the $`𝐒`$-matrix of $`\mathrm{\Phi }^4`$-theory”, hep-th/9911063, DESY, Hamburg, 1999.

 H. Epstein and V. Glaser, “The role of locality in perturbation theory”, Ann. Inst. Henri Poincaré A XIX (1973), 211–295.

 P. Blanchard and R. Sénéor, “Green’s functions for theories with massless particles (in perturbation theory)”, Ann. Inst. Henri Poincaré A XXIII (1975), 147–209.

 G. Scharf, Finite Quantum Electrodynamics: The Causal Approach, 2nd edition, Springer, Berlin, 1995.

 J. M. Gracia-Bondía, J. C. Várilly and H. Figueroa, Elements of Noncommutative Geometry, Birkhäuser, Boston, 2000.

 J. M. Gracia-Bondía and J. C. Várilly, “QED in external fields from the spin representation”, J. Math. Phys. 35 (1994), 3340–3367.

 D. Shale and W. F. Stinespring, “Spinor representations of infinite orthogonal groups”, J. Math. Mech. 14 (1965), 315–322.

 N. N. Bogoliubov and D. V. Shirkov, Introduction to the Theory of Quantized Fields, Wiley, New York, 1959.

 H. Grosse and E. Langmann, “A super-version of quasi-free second quantization: I. Charged particles”, J. Math. Phys. 33 (1992), 1032–1046.

 J. T. Ottesen, Infinite Dimensional Groups and Algebras in Quantum Physics, Lecture Notes in Physics: Monographs m27, Springer, Berlin, 1995.

 E. Langmann, “Cocycles for boson and fermion Bogoliubov transformations”, J. Math. Phys. 35 (1994), 96–112.

 J. C. Várilly and J. M. Gracia–Bondía, “$`S`$-matrix from the metaplectic representation”, Mod. Phys. Lett. A 7 (1992), 659–667.

 F. Scheck, Electroweak and Strong Interactions, Springer, Heidelberg, 1996.

 R. Estrada, “Regularization of distributions”, Int. J. Math. & Math. Sci. 21 (1998), 625–636.

 V. B. Berestetskii, E. M. Lifshitz and I. P. Pitaevskii, Quantum Electrodynamics (volume 4 of the Landau and Lifshitz Course of Theoretical Physics), Pergamon Press, Oxford, 1980.

 C. Itzykson and J.-B. Zuber, Quantum Field Theory, McGraw-Hill, New York, 1980.

 W. Greiner and J. Reinhardt, Quantum Electrodynamics, Springer, Berlin, 1992.

 E. Langmann and J. Mickelsson, “Scattering matrix in external field problems”, J. Math. Phys. 37 (1996), 3933–3953.

 J. M. Jauch and F. Rohrlich, The Theory of Photons and Electrons, Springer, Berlin, 1976.

 S. N. M. Ruijsenaars, “Charged particles in external fields. II. The quantized Dirac and Klein–Gordon theories”, Commun. Math. Phys. 52 (1977), 267–294.

 G. Scharf and W. F. Wreszinski, “The causal phase in Quantum Electrodynamics”, Nuovo Cim. A 93 (1986), 1–27.

 J. Bellissard, “Quantized fields in interaction with external fields I. Exact solutions and perturbative expansions”, Commun. Math. Phys. 41 (1975), 235–266.


format: text

text:
# The 𝑁 boson time dependent problem: an exact approach with stochastic wave functions

## 1 Introduction

Since the experimental realization of the first atomic gaseous Bose-Einstein condensates a few years ago , the physics of dilute Bose gases has been considered with a renewed interest. One fascinating aspect of these new systems is the possibility to accumulate in a single quantum state a large fraction of the atoms confined in a trap.

At very low temperature, a simple theoretical description of the dynamics of these systems is obtained by neglecting the uncondensed atoms, and by considering the wave function of the condensate, which obeys a Schrödinger equation with a non linear term originating from the mean-field interactions between the atoms. Such an approach neglects two- and more-particle correlations and is valid under a weak-interaction condition which is usually stated in terms of the density $`n`$ and the scattering length $`a`$ of the gas as $`(na^3)^{1/2}1`$. Current gaseous condensates satisfy such a condition. Nevertheless effects beyond the Gross-Pitaevskii equation may be considered at zero temperature; also finite temperature phenomena are not accounted for by the pure state mean field approach.

More complex theories have been developed in order to cope with effects beyond the Gross-Pitaevskii equation: Bogoliubov’s approach takes into account the next term in the $`(na^3)^{1/2}`$ expansion . Also quantum kinetic theories have been developed to study the formation of the condensate and to include the effect of the non-condensed particles . Unfortunately the corresponding calculations are quite heavy for 3D non-homogeneous systems such as trapped gases, and this constitutes a first limitation to the use of these methods. Also approximations used in some of these mean field theories are not under rigorous control, making it difficult to assess their domain of validity (for a review see e.g. ). Therefore a computational scheme capable to provide exact results can have a great importance both from a purely theoretical point of view and for a quantitative analysis of experimental data.

When the Bose gas is at thermal equilibrium such an exact numerical calculation of the properties of the gas is available, using the Quantum Monte-Carlo techniques, based on Feynman path integral formulation of quantum mechanics . The aim of this paper is to present an alternative exact and numerically tractable solution to the problem of the interacting Bose gas, a method not restricted to the case of thermal equilibrium but which allows for the study of the dynamics of the gas. The method is based on a stochastic evolution of Hartree states, in which all atoms have the same wave function, these Hartree states being either Fock states (fixed number of atoms) or coherent states. As a particular case of this solution with coherent states, we recover the stochastic scheme corresponding to the evolution of the density operator of the system in the positive $`P`$-representation . This evolution is known to lead to strong unstability problems, which fortunately do not show up for other implementations of the present method.

The outline of the paper is the following: In section 2, we present the stochastic formulation of the evolution of these Hartree states which, after average over the stochastic component, leads to the exact evolution. Section 3 is devoted to the presentation of two particular schemes implementing this stochastic formulation. We first present a ‘simple’ scheme, which minimizes the statistical spread of the calculated $`N`$-atom density matrix. We also investigate a more elaborate scheme in which the trace of the calculated density matrix remains strictly constant in the evolution. With this constraint, we recover for coherent states the known stochastic simulation associated with the positive $`P`$-representation . Finally we investigate in sections 4 and 5 two examples, a two-mode model system and a one-dimensional Bose gas respectively. These examples illustrate the accuracy and the limitations of the method. Generally speaking we find that the ‘simple’ scheme simulations are only limited by the computation power: the number of realizations needed for a good statistical accuracy increases exponentially with time for the simulation with Fock states. On the contrary the simulations with constant trace are subject to divergences of the norms of the stochastic wave functions in finite time, a phenomenon already known for coherent states in the context of the positive $`P`$-representation .

## 2 Stochastic formulation of the $`N`$-boson problem using Hartree functions

### 2.1 Model considered in this paper

The Hamiltonian of the trapped interacting Bose gas under exam can be written in terms of the Bose field operator $`\widehat{\mathrm{\Psi }}(x)`$ as:

$$=𝑑x\widehat{\mathrm{\Psi }}^{}(x)h_0\widehat{\mathrm{\Psi }}(x)+\frac{1}{2}𝑑x𝑑x^{}\widehat{\mathrm{\Psi }}^{}(x)\widehat{\mathrm{\Psi }}^{}(x^{})V(xx^{})\widehat{\mathrm{\Psi }}(x^{})\widehat{\mathrm{\Psi }}(x)$$
(1)
where $`x`$ is the set of spatial coordinates of a particle, $`h_0=\frac{\mathrm{}^2}{2m}^2+V_{\mathrm{ext}}(x)`$ is the single particle Hamiltonian in the external confining potential $`V_{\mathrm{ext}}`$ and where interactions are assumed to occur via a two-body potential $`V(xx^{})`$.

In practice we consider the dilute gas and the low temperature regimes, which correspond respectively to $`n|a|^31`$ and $`|a|\lambda `$ for a three-dimensional problem ($`\lambda =h/(2\pi mk_BT)^{1/2}`$ is the thermal de Broglie wavelength). The true interaction potential can then be replaced by a simpler model potential leading to the same scattering length $`a`$ provided that the range $`b`$ of this model potential is much smaller than the healing length $`\xi =(8\pi na)^{1/2}`$ and than $`\lambda `$ . This ensures that the physical results do not depend on $`b`$. For simplicity we will use here repulsive Gaussian potentials corresponding to a positive scattering length $`a>0`$.

### 2.2 A stochastic Hartree Ansatz with Fock states

From a mathematical point of view, the exact evolution of the $`N`$-body density matrix $`\rho `$ can be obtained from the Hamiltonian (1) using the quantum-mechanical equation of motion

$$\dot{\rho }(t)=\frac{1}{i\mathrm{}}[,\rho (t)]$$
(2)
but any concrete calculation is impracticable even for moderate particle numbers $`N`$, due to the multi-mode nature of the problem leading to a huge dimensionality of the $`N`$body Hilbert space.

For this reason approximate theories have been developed in order to get useful results at least in some specific ranges of parameters; the simplest one is the so-called mean-field theory, in which the $`N`$-particle density matrix is approximated by a Fock state Hartree ansatz

$$\rho (t)=|N:\varphi (t)N:\varphi (t)|.$$
(3)
The evolution of the normalized condensate wave function $`\varphi `$ is determined using either a factorization approximation in the evolution equation for the field operator  or a variational procedure . The result is the well-known mean-field equation

$$i\mathrm{}\frac{\varphi (x)}{t}=\left(\frac{\mathrm{}^2^2}{2m}+V_{\mathrm{ext}}(x)\right)\varphi (x)+(N1)\left(𝑑x^{}V(xx^{})\left|\varphi (x^{})\right|^2\right)\varphi (x).$$
(4)
For an interaction potential $`V(xx^{})`$ modeled by a contact term $`g\delta (xx^{})`$ (where $`g=4\pi \mathrm{}^2a/m`$ in a three-dimensional problem) it reduces to the Gross-Pitaevskii equation commonly used to analyze the dynamics of pure Bose-Einstein condensed gases.

A first attempt to improve the accuracy of the Hartree ansatz (3) is to allow for a stochastic contribution $`dB`$ in the evolution of the macroscopic wave function $`\varphi `$:

$$\varphi (t+dt)=\varphi (t)+Fdt+dB.$$
(5)
In all this paper the noise $`dB`$ is treated in the standard Ito formalism : it is assumed to have a zero mean $`\overline{dB}=0`$ and to have a variance $`\overline{dB^2}dt`$; a deterministic contribution is given by the “force” term $`Fdt`$. In this framework, the $`N`$-body density matrix would result from the stochastic mean over noise or, in other terms, from a mean over the probability distribution $`𝒫(\varphi )`$ in the functional space of the wave functions $`\varphi `$:

$$\rho (t)\stackrel{\mathrm{?}}{=}|N:\varphi (t)N:\varphi (t)|_{\mathrm{stoch}}=𝒟\varphi 𝒫(\varphi )|N:\varphi (t)N:\varphi (t)|.$$
(6)
An immediate advantage of this prescription over the pure state ansatz Eq.(3) is that it could deal with finite temperature problems . However as shown in §2.5, the simple generalization Eq.(5) of the Gross-Pitaevskii equation cannot lead to an exact solution of the $`N`$-body problem . Therefore we have to enlarge the family of dyadics over which we expand the density operator; more precisely we use Hartree dyadics in which the wave functions in the bra and in the ket are different:

$$\sigma (t)=|N:\varphi _1(t)N:\varphi _2(t)|.$$
(7)
The two wave functions $`\varphi _1(x)`$ and $`\varphi _2(x)`$ are assumed to evolve according to Ito stochastic differential equations:

$$\varphi _\alpha (t+dt)=\varphi _\alpha (t)+F_\alpha dt+dB_\alpha (\alpha =1,2).$$
(8)
The expansion Eq.(6) is then replaced by

$$\rho (t)=|N:\varphi _1(t)N:\varphi _2(t)|_{\mathrm{stoch}}=𝒟\varphi _1𝒟\varphi _2𝒫(\varphi _1,\varphi _2)|N:\varphi _1(t)N:\varphi _2(t)|.$$
(9)
We will see in the following that within this extended Hartree ansatz one can find a stochastic evolution for $`\varphi _{1,2}`$ reproducing the exact time evolution.

Actual calculations (see §4 and §5) will be performed with a Monte-Carlo technique, in which the evolution of the probability distribution $`𝒫`$ is simulated by a large but finite number $`𝒩`$ of independent realizations $`\varphi _{1,2}^{(i)}(t)`$, $`i=1,\mathrm{},𝒩`$. At any time the (approximate) density matrix $`\rho `$ is given by the mean over such an ensemble of wave functions:

$$\rho (t)\frac{1}{𝒩}\underset{i=1}{\overset{𝒩}{}}|N:\varphi _1^{(i)}(t)N:\varphi _2^{(i)}(t)|.$$
(10)
The expectation value of any operator $`\widehat{O}`$ is thus expressed by:

$$\widehat{O}\frac{1}{𝒩}\underset{i=1}{\overset{𝒩}{}}N:\varphi _2^{(i)}(t)|\widehat{O}|N:\varphi _1^{(i)}(t).$$
(11)
For an Hermitian operator one can equivalently consider only the real part of this expression since the imaginary part is vanishingly small in the large $`𝒩`$ limit.

Consider as an example the one-particle density matrix of the gas, usually defined as:

$$\rho ^{(1)}(x,x^{})=\widehat{\mathrm{\Psi }}^{}(x^{})\widehat{\mathrm{\Psi }}(x).$$
(12)
Inserting in this expression our form of the complete density matrix (10), we obtain the simple result

$$\rho ^{(1)}(x,x^{})=N\varphi _1(x)\varphi _2^{}(x^{})\varphi _2|\varphi _1^{N1}_{\mathrm{stoch}}$$
(13)
from which it is easy to obtain the spatial density $`n(x)=\rho ^{(1)}(x,x)`$ and the correlation function $`g^{(1)}(x,x^{})=\rho ^{(1)}(x,x^{})/(n(x)n(x^{}))^{1/2}`$. Also, the condensate fraction can be obtained from the largest eigenvalue of $`\rho ^{(1)}(x,x^{})`$.

Remarks:

1. The desired stochastic evolution, which has to satisfy Tr$`[\rho ]=1`$, cannot preserve the normalization of $`\varphi _{1,2}`$ to unity; we can write indeed

$$\text{Tr}[\rho (t)]=\varphi _2(t)|\varphi _1(t)^N_{\mathrm{stoch}}=1$$
(14)
which for $`|\varphi _1|\varphi _2`$ imposes $`\varphi _1\varphi _2>1`$.
2. The expansion Eq.(9) is always possible. Using the identity

$$\text{Id}_N=\underset{+\mathrm{}}{lim}\frac{1}{}\underset{j=1}{\overset{}{}}|N:\psi ^{(j)}N:\psi ^{(j)}|$$
(15)
where the functions $`\psi ^{(j)}`$ have a uniform distribution over the unit sphere in the functional space, we obtain

$$\rho =\underset{+\mathrm{}}{lim}\frac{1}{^2}\underset{j_1,j_2=1}{\overset{}{}}|N:\psi ^{(j_1)}N:\psi ^{(j_2)}|N:\psi ^{(j_1)}|\rho |N:\psi ^{(j_2)}.$$
(16)
We write the matrix elements $`N:\psi ^{(j_1)}|\rho |N:\psi ^{(j_2)}`$ as $`\xi _{(j_1,j_2)}^{2N}`$ and we set $`\varphi _1^{(j_1,j_2)}=\psi ^{(j_1)}\xi _{(j_1,j_2)}`$ and $`\varphi _2^{(j_1,j_2)}=\psi ^{(j_2)}\xi _{(j_1,j_2)}^{}`$. Putting $`𝒩=^2`$ and reindexing $`(j_1,j_2)`$ as a single index $`i`$ we recover the expansion Eq.(10). Note that this expansion is not unique and does not have the pretension to be the most efficient one. For instance if the system is initially in a Hartree state $`|N:\varphi _0`$, such a procedure is clearly not needed since one has just to set $`\varphi _1^{(i)}(t=0)=\varphi _2^{(i)}(t=0)=\varphi _0`$. This will be the case of the numerical examples in sections 4 and 5.

### 2.3 Stochastic evolution of a Fock state Hartree dyadic

In this subsection we calculate the stochastic time evolution during an infinitesimal time interval $`dt`$ of the dyadic $`\sigma (t)`$ given in Eq.(7). This will be used later in a comparison with the exact master equation.

After $`dt`$, the dyadic $`\sigma `$ has evolved into:

$$\sigma (t+dt)=|N:\varphi _1+d\varphi _1N:\varphi _2+d\varphi _2|,$$
(17)
where $`d\varphi _1`$ and $`d\varphi _2`$, defined according to (8), contain both the deterministic contribution $`F_\alpha dt`$ and the stochastic one $`dB_\alpha `$. Splitting each contribution into a longitudinal and an orthogonal component with respect to $`\varphi _\alpha `$ and isolating a Gross-Pitaevskii term in the deterministic contribution, we can write:

$`dB_\alpha (x)=\varphi _\alpha (x)d\gamma _\alpha +dB_\alpha ^{}(x)`$ (18)
$`F_\alpha (x)=F_\alpha ^{GP}(x)+\lambda _\alpha \varphi _\alpha (x)+F_\alpha ^{}(x).`$ (19)
Our choice of the Gross-Pitaevskii term is the following one:

$`F_\alpha ^{GP}(x)=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}\left[h_0+{\displaystyle \frac{(N1)}{\varphi _\alpha ^2}}{\displaystyle 𝑑x^{}V(xx^{})\left|\varphi _\alpha (x^{})\right|^2}\right]\varphi _\alpha (x)`$ (20)
$``$ $`{\displaystyle \frac{1}{i\mathrm{}}}\left[{\displaystyle \frac{(N1)}{2}}{\displaystyle \frac{\varphi _\alpha \varphi _\alpha |V|\varphi _\alpha \varphi _\alpha }{\varphi _\alpha ^4}}\right]\varphi _\alpha (x).`$
The first term gives the standard Gross-Pitaevskii evolution, including the kinetic term, the potential energy of the trap and the mean-field interaction energy; the second term, which arises naturally because we are considering Fock states (rather than coherent states as commonly done) takes into account the difference between the total mean-field energy per particle of the condensate and its chemical potential $`\mu `$ .

We split the field operator in its longitudinal and transverse components, keeping in mind that the wave functions $`\varphi _\alpha `$ are not of unit norm:

$$\widehat{\mathrm{\Psi }}^{}(x)=\frac{\varphi _\alpha ^{}(x)}{\varphi _\alpha ^2}\widehat{a}_{\varphi _\alpha }^{}+\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x)$$
(21)
with

$$\widehat{a}_{\varphi _\alpha }^{}=𝑑x\varphi _\alpha (x)\widehat{\mathrm{\Psi }}^{}(x).$$
(22)
The relevant bosonic commutation relations then read:

$$[\widehat{a}_{\varphi _\alpha },\widehat{a}_{\varphi _\alpha }^{}]=\varphi _\alpha ^2\text{and}[\widehat{a}_{\varphi _\alpha },\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x)]=0.$$
(23)
We will also need the projector $`𝒬_\alpha `$ onto the subspace orthogonal to $`\varphi _\alpha `$:

$$𝒬_\alpha ^{(x)}[\psi (x,x^{},\mathrm{})]=\psi (x,x^{},\mathrm{})\frac{\varphi _\alpha (x)}{\varphi _\alpha ^2}𝑑y\varphi _\alpha ^{}(y)\psi (y,x^{},\mathrm{}).$$
(24)
This projector arises in the calculation as we have introduced a component of the field operator orthogonal to $`\varphi _\alpha `$. Using $`𝑑x\varphi _\alpha (x)\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x)=0`$ we shall transform integrals involving $`\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x)`$ as follows:

$$𝑑x\psi (x,x^{},\mathrm{})\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x)=𝑑x𝒬_\alpha ^{(x)}[\psi (x,x^{},\mathrm{})]\delta \widehat{\mathrm{\Psi }}_\alpha ^{}(x).$$
(25)
Inserting these definitions in (17) the expression for $`\sigma `$ at time $`t+dt`$ can be written as

$`\sigma (t+dt)\sigma (t)=`$ $`S_1^{(0)}|N:\varphi _1N:\varphi _2|+\text{e.c.}`$ (26)
$`+`$ $`{\displaystyle }dxS_1^{(1)}(x)\delta \widehat{\mathrm{\Psi }}_1^{}(x)|N1:\varphi _1N:\varphi _2|+\text{e.c.}`$
$`+`$ $`{\displaystyle }{\displaystyle }dxdx^{}S_1^{(2)}(x,x^{})\delta \widehat{\mathrm{\Psi }}_1^{}(x)\delta \widehat{\mathrm{\Psi }}_1^{}(x^{})|N2:\varphi _1N:\varphi _2|+\text{e.c.}`$
$`+`$ $`{\displaystyle }{\displaystyle }dxdx^{}S^{(1,1)}(x,x^{})\delta \widehat{\mathrm{\Psi }}_1^{}(x)|N1:\varphi _1N1:\varphi _2|\delta \widehat{\mathrm{\Psi }}_2(x^{})`$
where the notation e.c. stands for the exchanged and conjugate of a quantity, that is the complex conjugate of the same quantity after having exchanged the indices 1 and 2. The explicit expressions for the $`S_\alpha ^{(i)}`$ are:

$`S_1^{(0)}=N{\displaystyle \frac{\varphi _1|F_1^{GP}}{\varphi _1^2}}dt+N\lambda _1dt+Nd\gamma _1+{\displaystyle \frac{N(N1)}{2}}d\gamma _1^2+{\displaystyle \frac{N^2}{2}}d\gamma _1d\gamma _2^{}`$ (27)
$`S_1^{(1)}(x)=\sqrt{N}\left\{𝒬_1^{(x)}\left[F_1^{GP}(x)\right]dt+F_1^{}(x)dt+dB_1^{}(x)+(N1)d\gamma _1dB_1^{}(x)+NdB_1^{}(x)d\gamma _2^{}\right\}`$ (28)
$`S_1^{(2)}(x,x^{})={\displaystyle \frac{\sqrt{N(N1)}}{2}}dB_1^{}(x)dB_1^{}(x^{})`$ (29)
$`S^{(1,1)}(x,x^{})=NdB_1^{}(x)dB_2^{}(x^{}).`$ (30)
Analogous expressions for $`S_2^{(0)}`$, $`S_2^{(1)}`$, $`S_2^{(2)}`$ are obtained by exchanging the indices 1 and 2. In the next subsection, we evaluate the exact evolution of the same dyadic during a time interval $`dt`$, so that we can determine the constraints on the force and noise terms entering into these equations.

### 2.4 Exact evolution of a Fock state Hartree dyadic

To make the stochastic scheme described in the previous sections equivalent to the exact dynamics as it is given by (1), the final result of the previous subsection (26)-(30) has to be compared with the exact evolution of the density matrix $`\sigma (t)`$. Consider a dyadic $`\sigma =|N:\varphi _1N:\varphi _2|`$ at time $`t`$; according to the equation of motion (2), after an infinitesimal time step $`dt`$ it has evolved into:

$`\sigma (t+dt)`$ $`=`$ $`\sigma (t)+{\displaystyle \frac{dt}{i\mathrm{}}}\left(\sigma (t)\sigma (t)\right)`$ (31)
$`=`$ $`\sigma (t)+E_1^{(0)}|N:\varphi _1N:\varphi _2|+\text{e.c.}`$
$`+`$ $`{\displaystyle }dxE_1^{(1)}(x)\delta \widehat{\mathrm{\Psi }}_1^{}(x)|N1:\varphi _1N:\varphi _2|+\text{e.c.}`$
$`+`$ $`{\displaystyle }{\displaystyle }dxdx^{}E_1^{(2)}(x,x^{})\delta \widehat{\mathrm{\Psi }}_1^{}(x)\delta \widehat{\mathrm{\Psi }}_1^{}(x^{})|N2:\varphi _1N:\varphi _2|+\text{e.c.}`$
where the $`E_\alpha ^{(i)}`$ are given by

$`E_1^{(0)}={\displaystyle \frac{Ndt}{i\mathrm{}}}\left[{\displaystyle \frac{\varphi _1|h_0|\varphi _1}{\varphi _1^2}}+{\displaystyle \frac{(N1)}{2}}{\displaystyle \frac{\varphi _1\varphi _1|V|\varphi _1\varphi _1}{\varphi _1^4}}\right]`$ (32)
$`E_1^{(1)}(x)={\displaystyle \frac{dt\sqrt{N}}{i\mathrm{}}}𝒬_1^{(x)}\left[\left(h_0+{\displaystyle \frac{(N1)}{\varphi _1^2}}{\displaystyle 𝑑x^{}V(xx^{})\left|\varphi _1(x^{})\right|^2}\right)\varphi _1(x)\right]`$ (33)
$`E_1^{(2)}(x,x^{})={\displaystyle \frac{dt\sqrt{N(N1)}}{2i\mathrm{}}}𝒬_1^{(x)}𝒬_1^{(x^{})}[V(xx^{})\varphi _1(x)\varphi _1(x^{})].`$ (34)
Analogous expressions for $`E_2^{(0)}`$, $`E_2^{(1)}`$, $`E_2^{(2)}`$ are obtained by exchanging the indices 1 and 2.

### 2.5 Validity conditions for the stochastic Fock state Hartree ansatz

The similarity of the structures of (26) and (31) suggests the possibility of a stochastic scheme equivalent to the exact evolution: to achieve this, it is necessary to find out specific forms of deterministic (19) and stochastic (18) terms for which the mean values of the $`S_\alpha ^{(i)}`$ equal the $`E_\alpha ^{(i)}`$:

$`\overline{S_1^{(0)}+S_2^{(0)}}=E_1^{(0)}+E_2^{(0)}`$ (35)
$`\overline{S_\alpha ^{(1)}(x)}=E_\alpha ^{(1)}(x)`$ (36)
$`\overline{S_\alpha ^{(2)}(x,x^{})}=E_\alpha ^{(2)}(x,x^{})`$ (37)
$`\overline{S^{(1,1)}(x,x^{})}=0.`$ (38)
From the last equation (38), it follows immediately why independent bras and kets are needed in the ansatz (7): in the case $`\varphi _1=\varphi _2=\varphi `$ such a condition would in fact lead to a vanishing orthogonal noise and finally to the impossibility of satisfying (37).

In terms of the different components, these conditions can be rewritten as:

$`\left(\lambda _1+\lambda _2^{}\right)dt+{\displaystyle \frac{(N1)}{2}}\left[\overline{d\gamma _1^2}+\overline{d\gamma _2^2}^{}\right]+N\overline{d\gamma _1d\gamma _2^{}}=0`$ (39)
$`F_1^{}(x)dt+(N1)\overline{dB_1^{}(x)d\gamma _1}+N\overline{dB_1^{}(x)d\gamma _2^{}}=0`$ (40)
$`F_2^{}(x)dt+(N1)\overline{dB_2^{}(x)d\gamma _2}+N\overline{dB_2^{}(x)d\gamma _1^{}}=0`$ (41)
$`\overline{dB_\alpha ^{}(x)dB_\alpha ^{}(x^{})}={\displaystyle \frac{dt}{i\mathrm{}}}𝒬_\alpha ^x𝒬_\alpha ^x^{}\left[V(xx^{})\varphi _\alpha (x)\varphi _\alpha (x^{})\right]`$ (42)
$`\overline{dB_1^{}(x)dB_{2}^{}{}_{}{}^{}(x^{})}=0.`$ (43)
As we shall discuss in detail in §3, several different stochastic schemes can be found satisfying (39)-(43); each of them gives an evolution identical in average to the exact one, but the statistical properties can be very different.

### 2.6 A stochastic Hartree ansatz with coherent states

Up to now we have worked out the case of a Fock state ansatz $`|N:\varphi _1N:\varphi _2|`$. Actually coherent states rather than Fock states are generally used, both in quantum optics and in condensed matter physics. We now show that our stochastic procedure also applies with a coherent state ansatz of the form:

$$\sigma (t)=\mathrm{\Pi }(t)|\text{ coh}:\varphi _1\text{ coh}:\varphi _2|,$$
(44)
with

$$|\text{ coh}:\varphi _\alpha =\mathrm{exp}\left(\overline{N}^{1/2}𝑑x\varphi _\alpha (x)\widehat{\mathrm{\Psi }}^{}(x)\right)|\mathrm{vac}$$
(45)
where $`\overline{N}`$ is the mean number of particles. We have included here a prefactor $`\mathrm{\Pi }(t)`$ which was absent in the case of the Fock state ansatz Eq.(7); in the Fock state case indeed such a prefactor could be reincluded into the definition of $`\varphi _1`$ and $`\varphi _2`$. The wave functions $`\varphi _\alpha (x)`$ and the prefactor factor $`\mathrm{\Pi }`$ evolve according to Ito stochastic differential equations

$`d\varphi _\alpha =F_\alpha dt+dB_\alpha `$
$`d\mathrm{\Pi }=fdt+db.`$ (46)
Splitting the field operator as

$$\widehat{\mathrm{\Psi }}(x)=\overline{N}^{1/2}\varphi _\alpha (x)+\delta \widehat{\mathrm{\Psi }}_\alpha (x)$$
(47)
and using

$$\delta \widehat{\mathrm{\Psi }}_\alpha (x)|\text{ coh}:\varphi _\alpha =0$$
(48)
we find that the equivalence of the stochastic scheme and the exact evolution translates into the set of conditions:

$`f=0`$ (49)
$`F_1(x)dt+{\displaystyle \frac{1}{\mathrm{\Pi }}}\overline{dbdB_1(x)}={\displaystyle \frac{dt}{i\mathrm{}}}h_0\varphi _1(x)`$ (50)
$`F_2(x)dt+{\displaystyle \frac{1}{\mathrm{\Pi }^{}}}\overline{db^{}dB_2(x)}={\displaystyle \frac{dt}{i\mathrm{}}}h_0\varphi _2(x)`$ (51)
$`\overline{dB_\alpha (x)dB_\alpha (x^{})}={\displaystyle \frac{dt}{i\mathrm{}}}V(xx^{})\varphi _\alpha (x)\varphi _\alpha (x^{})`$ (52)
$`\overline{dB_1(x)dB_2^{}(x^{})}=0.`$ (53)
As we shall see in §3, such conditions are satisfied by several stochastic schemes. Very remarkably, the stochastic evolution deduced from the positive $`P`$-representation  arises naturally as one of them.

Within this coherent state ansatz the one-particle density matrix is evaluated using

$$\rho ^{(1)}(x,x^{})=\overline{N}\varphi _1(x)\varphi _2^{}(x^{})\mathrm{\Pi }(t)\mathrm{exp}(\overline{N}\varphi _2|\varphi _1)_{\mathrm{stoch}}.$$
(54)
In a practical implementation of the simulation it turns out to be numerically more efficient to represent $`\mathrm{\Pi }(t)`$ as the exponential of some quantity

$$\mathrm{\Pi }(t)=e^{\overline{N}S(t)}$$
(55)
and to evolve $`S(t)`$ according to the stochastic equation

$`dS`$ $`=`$ $`{\displaystyle 𝑑x[dB_1(x)\varphi _1^{}(x)+dB_2^{}(x)\varphi _2(x)]}`$ (56)
$`{\displaystyle \frac{\overline{N}dt}{2i\mathrm{}}}{\displaystyle 𝑑x𝑑x^{}V(xx^{})\left[|\varphi _1(x)|^2|\varphi _1(x^{})|^2|\varphi _2(x)|^2|\varphi _2(x^{})|^2\right]}.`$
## 3 Particular implementations of the stochastic approach

In the previous section we have derived the conditions that a stochastic scheme has to satisfy in order to recover the exact evolution given by the Hamiltonian (1); in the case of the Fock state ansatz (7), we get to the system (39)-(43), while in the case of the coherent state ansatz (44) we get to the conditions (49)-(53). As the number of these equations is actually smaller than the number of unknown functions there is by no mean uniqueness of the solutions, that is of the simulation schemes. We need a strategy to identify interesting solutions.

We therefore start this section by considering various indicators of the statistical error of the simulation (§3.1) which can be used as guidelines in the search for simulation schemes. These indicators are defined as variances of relevant quantities which are conserved in the exact evolution but which may fluctuate in the simulation. We show that these indicators are non decreasing functions of time; attempts to minimize the time derivative of a specific indicator will lead to particular implementations of the general stochastic method, such as the simple scheme (§3.2) and the constant trace scheme (§3.3).

### 3.1 Growth of the statistical errors

The first indicator that we consider measures the squared deviation of the stochastic operator $`\sigma (t)`$ from the exact density operator $`\rho (t)`$:

$`\mathrm{\Delta }(t)`$ $`=`$ $`\text{Tr}[(\sigma ^{}(t)\rho (t))(\sigma (t)\rho (t))]_{\mathrm{stoch}}`$ (57)
$`=`$ $`\text{Tr}[\sigma ^{}(t)\sigma (t)]_{\mathrm{stoch}}\text{Tr}[\rho (t)^2].`$
We now show that $`\mathrm{\Delta }(t)`$ is a non-decreasing function of time. When the stochastic scheme satisfies the validity conditions derived in the previous section, we can write the stochastic equation for $`\sigma `$ as:

$$d\sigma =\frac{dt}{i\mathrm{}}[,\sigma ]+d\sigma _s$$
(58)
where $`d\sigma _s`$ is a zero-mean noise term linear in $`dB_\alpha `$ (and $`db`$ for the coherent state simulation). In the case of simulation with Fock states it is given by

$$d\sigma _s=N^{1/2}\left\{𝑑x𝑑B_1(x)\widehat{\mathrm{\Psi }}^{}(x)|N1:\varphi _1N:\varphi _2|+𝑑x𝑑B_2^{}(x)|N:\varphi _1N1:\varphi _2|\widehat{\mathrm{\Psi }}(x)\right\}.$$
(59)
In the case of simulation with coherent states it is given by

$`d\sigma _s`$ $`=`$ $`db|\mathrm{coh}:\varphi _1\mathrm{coh}:\varphi _2|`$ (60)
$`+`$ $`\overline{N}^{1/2}\mathrm{\Pi }\left\{{\displaystyle 𝑑x𝑑B_1(x)\widehat{\mathrm{\Psi }}^{}(x)|\mathrm{coh}:\varphi _1\mathrm{coh}:\varphi _2|+𝑑x𝑑B_2^{}(x)|\mathrm{coh}:\varphi _1\mathrm{coh}:\varphi _2|\widehat{\mathrm{\Psi }}(x)}\right\}.`$
We calculate the variation of $`\mathrm{\Delta }`$ during $`dt`$, replacing $`\sigma `$ by $`\sigma +d\sigma `$ in Eq.(57) and keeping terms up to order $`dt`$. Using the invariance of the trace in a cyclic permutation and averaging over the noise between $`t`$ and $`t+dt`$ we finally obtain

$$d\mathrm{\Delta }=\text{Tr}[d\sigma _s^{}d\sigma _s]_{\mathrm{stoch}},$$
(61)
which is a positive quantity. Minimization of this quantity is the subject of §3.2. Physically $`d\mathrm{\Delta }0`$ means that the impurity of the stochastic density operator $`\sigma `$ always increases in average, while the exact density operator has a constant purity $`\text{Tr}[\rho ^2]`$.

The second kind of indicator that we consider measures the statistical error on constants of motion of the exact evolution. Consider a time independent operator $`X`$ commuting with the Hamiltonian. The stochastic evolution leads to an error on the expectation value of $`X`$ with a variance given by the ensemble average of

$`\mathrm{\Delta }_X(t)`$ $`=`$ $`\left|\text{Tr}[X\sigma (t)]\text{Tr}[X\rho (t)]\right|^2_{\mathrm{stoch}}`$ (62)
$`=`$ $`\left|\text{Tr}[X\sigma (t)]\right|^2_{\mathrm{stoch}}\left|\text{Tr}[X\rho (t)]\right|^2.`$
From Eq.(58) we obtain the variation after a time step $`dt`$ of the expectation value of $`X`$ along a stochastic trajectory:

$$d(\text{Tr}[X\sigma ])=\frac{dt}{i\mathrm{}}\text{Tr}\left(X[,\sigma ]\right)+\text{Tr}[Xd\sigma _s].$$
(63)
Using the invariance of the trace under cyclic permutation and the commutation relation $`[,X]=0`$ we find that the first term in the right hand side of Eq.(63) vanishes so that

$$d\mathrm{\Delta }_X=\left|\text{Tr}[Xd\sigma _s]\right|^2_{\mathrm{stoch}},$$
(64)
a quantity which is always non-negative.

Using expression (64) one can ‘design’ simulations preserving exactly the conserved quantity, the constraint to meet being $`\text{Tr}[Xd\sigma _s]=0`$: for instance in the Fock state simulation, one first chooses the transverse noises $`dB_\alpha ^{}`$ satisfying Eqs.(42,43); then one simply has to take for the longitudinal noise of $`\varphi _1`$:

$$d\gamma _1=\frac{1}{\sqrt{N}}(N:\varphi _2|X|N:\varphi _1)^1𝑑x𝑑B_1^{}(x)N:\varphi _2|X\delta \widehat{\mathrm{\Psi }}_1^{}(x)|N1:\varphi _1$$
(65)
and a similar expression for $`d\gamma _2`$; finally the force terms $`F_\alpha `$ are adjusted in order to satisfy Eq.(39-41). As natural examples of conserved quantities one can choose $`X=1`$ or $`X=`$; the former case is discussed in detail in §3.3.

### 3.2 The simple schemes

These schemes are characterized by the minimization of the incremental variation of the statistic spread of the $`N`$-particle density matrix $`\sigma (t)`$, a spread that we have already quantified in Eq.(57) by $`\mathrm{\Delta }(t)`$. To be more specific we assume that we have evolved a dyadic up to time $`t`$, and we look for the noise terms that minimize the increase of $`\text{Tr}[\sigma ^{}\sigma ]`$ between $`t`$ and $`t+dt`$.

#### 3.2.1 Simulation with Fock states

In the case of the Fock state ansatz, we calculate explicitly the variation of $`\text{Tr}[\sigma ^{}\sigma ]`$ from Eq.(59) and we get:

$$\frac{d\text{Tr}[\sigma ^{}\sigma ]}{N\text{Tr}[\sigma ^{}\sigma ]}=N\overline{|d\gamma _1+d\gamma _2^{}|^2}+\underset{\alpha =1,2}{}\varphi _\alpha ^2𝑑x\overline{|dB_\alpha ^{}(x)|^2}+\left[d\gamma _1+d\gamma _2^{}+\text{c.c.}\right].$$
(66)
We now look for the noise terms $`d\gamma _\alpha `$ and $`dB_\alpha ^{}`$ minimizing this quantity under the constraints Eqs.(39-43).

We first note that we can choose $`d\gamma _1=d\gamma _2=0`$ without affecting the transverse noises, as shown by Eqs.(39-43): the correlation function of the transverse noises do not involve the $`d\gamma _\alpha `$, and we can accommodate for any choice of $`d\gamma _\alpha `$ by defining appropriately the force terms $`F_\alpha ^{},\lambda _\alpha `$. In the particular case defining our simple scheme we take all these force terms equal to zero. Note that the choice of vanishing $`d\gamma `$’s immediately leads to a vanishing noise term in Eq.(66).

Secondly the terms involving the transverse noise in Eq.(66) are bounded from below: As the modulus of a mean is less than the mean of the modulus, we have

$$\left|\overline{dB_\alpha ^{}(x)dB_\alpha ^{}(x)}\right|\overline{|dB_\alpha ^{}(x)|^2},$$
(67)
with the left hand side of this inequality fully determined by condition Eq.(42).

We have found for the transverse noise a choice which fulfills Eqs.(42,43) and which saturates the inequality Eq.(67):

$$dB_\alpha ^{}(x)=\left(\frac{dt}{i\mathrm{}}\right)^{1/2}𝒬_\alpha ^{(x)}\left[\varphi _\alpha (x)\frac{dk}{(2\pi )^{d/2}}\left(\stackrel{~}{V}(k)\right)^{1/2}e^{ikx}e^{i\theta _\alpha (k)}\right]$$
(68)
where $`d`$ is the dimension of position space, $`\stackrel{~}{V}(k)`$ is the Fourier transform of the model interaction potential, supposed here to be positive:

$$\stackrel{~}{V}(k)=𝑑xV(x)e^{ikx}.$$
(69)
The phases $`\theta _\alpha `$ have the following statistical property:

$$\overline{e^{i\theta _\alpha (k)}e^{i\theta _\alpha (k^{})}}=\delta (k+k^{})$$
(70)
and $`\theta _1,\theta _2`$ are uncorrelated. In practice for half of the $`k`$-space (e.g. $`k_z>0`$) $`\theta _\alpha (k)`$ is randomly chosen between $`0`$ and $`2\pi `$; for the remaining $`k`$’s we take $`\theta _\alpha (k)=\theta _\alpha (k)`$. One can then check that this choice for the transverse noise reproduces the correlation function Eqs.(42,43).

We show now that the implementation (68) saturates the inequality Eq.(67), so that it leads to the minimal possible value for $`d\text{Tr}[\sigma ^{}\sigma ]`$ within the validity constraints of the stochastic approach. We calculate explicitly the right hand side of Eq.(67):

$`\overline{|dB_\alpha ^{}(x)|^2}`$ $`=`$ $`{\displaystyle \frac{dt}{\mathrm{}}}\left(𝒬_\alpha ^{(x)}𝒬_\alpha ^{(x^{})}[V(xx^{})\varphi _\alpha ^{}(x)\varphi _\alpha (x^{})]\right)_{x=x^{}}`$ (71)
$`=`$ $`{\displaystyle \frac{dt}{\mathrm{}}}|\varphi _\alpha (x)|^2\left[V(0)2{\displaystyle 𝑑y\frac{|\varphi _\alpha (y)|^2}{\varphi _\alpha ^2}V(xy)}+{\displaystyle \frac{\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha }{\varphi _\alpha ^4}}\right]`$
where $`𝒬_\alpha ^{}`$ projects onto the subspace orthogonal to $`\varphi _\alpha ^{}`$ and where we have used the positivity of the Fourier transform $`\stackrel{~}{V}`$ of the model interaction potential. The left hand side of Eq.(67) is calculated using Eq.(42):

$$\overline{dB_\alpha ^{}(x)^2}=\frac{dt}{i\mathrm{}}\varphi _\alpha ^2(x)\left[V(0)2𝑑y\frac{|\varphi _\alpha (y)|^2}{\varphi _\alpha ^2}V(xy)+\frac{\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha }{\varphi _\alpha ^4}\right].$$
(72)
As the expressions between square brackets in Eqs.(72,71) are real positive we deduce the equality in Eq.(67).

We can now calculate explicitly the variation of $`\text{Tr}[\sigma ^{}\sigma ]`$ by integrating Eq.(71) over $`x`$:

$$\frac{d\text{Tr}[\sigma ^{}\sigma ]}{N\text{Tr}[\sigma ^{}\sigma ]}=\frac{dt}{\mathrm{}}\left[2V(0)\underset{\alpha =1,2}{}\frac{\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha }{\varphi _\alpha ^4}\right].$$
(73)
This expression is particularly useful since it allows one to derive an upper bound on the increase of $`\text{Tr}[\sigma ^{}\sigma ]`$: as we assume here a positive Gaussian model potential $`V(xx^{})`$ the matrix element $`\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha `$ is positive so that the right hand side of Eq.(73) is smaller than $`2V(0)dt/\mathrm{}`$. After time integration, using Eq.(57) and the fact that the trace of the squared density operator $`\rho ^2`$ is a constant under Hamiltonian evolution we can deduce an upper bound on the squared statistical error $`\mathrm{\Delta }(t)`$:

$$\mathrm{\Delta }(t)+\text{Tr}[\rho ^2]\left[\mathrm{\Delta }(0)+\text{Tr}[\rho ^2]\right]e^{2NV(0)t/\mathrm{}}.$$
(74)
Note that it involves the model dependent quantity $`V(0)`$ and not only the physical parameters of the problem such as the chemical potential or the scattering length. It may be therefore important to adjust the model interaction potential $`V(xx^{})`$ in order to minimize the growth of the statistical error for given physical parameters.

To summarize the proposed simple scheme has several noticeable properties. The deterministic force acting on the $`\varphi _\alpha `$’s is simply the mean field contribution, so that the whole correction to the mean field evolution is provided by the transverse noises $`dB_\alpha ^{}`$. Also the evolutions of the two states $`\varphi _\alpha `$ are totally independent from each other.

#### 3.2.2 Simulation with coherent states

In the case of the coherent state ansatz, an explicit calculation of $`d\text{Tr}[\sigma ^{}\sigma ]`$ from Eq.(60) gives:

$`{\displaystyle \frac{d\text{Tr}[\sigma ^{}\sigma ]}{\text{Tr}[\sigma ^{}\sigma ]}}`$ $`=`$ $`\overline{\left|{\displaystyle \frac{db}{\mathrm{\Pi }}}+\overline{N}{\displaystyle 𝑑x𝑑B_1(x)\varphi _1^{}(x)}+dB_2^{}(x)\varphi _2(x)\right|^2}+\overline{N}{\displaystyle \underset{\alpha =1,2}{}}{\displaystyle 𝑑x\overline{|dB_\alpha (x)|^2}}`$ (75)
$`+\left[{\displaystyle \frac{db}{\mathrm{\Pi }}}+\overline{N}{\displaystyle 𝑑x𝑑B_1(x)\varphi _1^{}(x)}+dB_2^{}(x)\varphi _2(x)+\text{c.c}\right].`$
We now proceed to the minimization of the increment of $`\text{Tr}[\sigma ^{}\sigma ]`$ within the coherent state ansatz along the same lines as the previous subsection. First we optimize the noise $`db`$ on the normalization factor $`\mathrm{\Pi }`$:

$$db=\overline{N}\mathrm{\Pi }\left(𝑑x𝑑B_1(x)\varphi _1^{}(x)+dB_2^{}(x)\varphi _2(x)\right).$$
(76)
This choice leads to a vanishing noise term in Eq.(75). We insert this expression for $`db`$ in the validity conditions Eq.(50) and Eq.(51) and we obtain:

$$F_\alpha (x)=\frac{1}{i\mathrm{}}\left[h_0+\overline{N}𝑑x^{}V(xx^{})|\varphi _\alpha (x^{})|^2\right]\varphi _\alpha (x).$$
(77)
Finally minimization of the contribution of the noise terms $`dB_\alpha `$ with the constraint Eq.(52) is achieved with the choice

$$dB_\alpha (x)=\left(\frac{dt}{i\mathrm{}}\right)^{1/2}\varphi _\alpha (x)\frac{dk}{(2\pi )^{d/2}}\left(\stackrel{~}{V}(k)\right)^{1/2}e^{ikx}e^{i\theta _\alpha (k)}$$
(78)
where the phases $`\theta _\alpha (k)`$ are randomly generated as in Eq.(70).

The first equation Eq.(77) fixes the deterministic evolution to the usual mean-field equation (4). We note here that the mean-field term in Eq.(77) does not contain the normalization of the spatial density $`\overline{N}|\varphi _\alpha (x^{})|^2`$ by $`\varphi _\alpha ^2`$, a feature present in the Fock state simulation (see Eq.(20)). This is a disadvantage of the coherent state simulation since this normalization factor appearing in the Fock state simulation has a regularizing effect: the norms $`\varphi _\alpha `$ may indeed deviate significantly from unity in the stochastic evolution. The second equation Eq.(78) determines the stochastic noise on the wave functions in a way very similar to the Fock state case Eq.(68). In particular the evolutions of $`\varphi _1`$ and $`\varphi _2`$ are still uncorrelated. The only difference is the disappearance of the projector $`𝒬_\alpha `$ in the expression of the noise, which leads to an increased noise with respect to the simulation with Fock states.

As in the previous subsection we now estimate the squared error $`\mathrm{\Delta }`$. We calculate the variation of $`\text{Tr}[\sigma ^{}\sigma ]`$ for the choice of noise Eq.(78):

$$\frac{1}{\text{Tr}[\sigma ^{}\sigma ]}\frac{d\text{Tr}[\sigma ^{}\sigma ]}{dt}=\frac{\overline{N}V(0)}{\mathrm{}}\underset{\alpha =1,2}{}\varphi _\alpha ^2.$$
(79)
The average over all stochastic realizations of the norm squared of the wave functions can be calculated exactly:

$$\varphi _\alpha ^2_{\mathrm{stoch}}(t)=e^{tV(0)/\mathrm{}}\varphi _\alpha ^2_{\mathrm{stoch}}(0).$$
(80)
This leads to a remarkable identity on the trace of $`\sigma ^{}\sigma `$:

$$\mathrm{ln}\text{Tr}[\sigma ^{}\sigma ]_{\mathrm{stoch}}(t)=\mathrm{ln}\text{Tr}[\sigma ^{}\sigma ]_{\mathrm{stoch}}(0)+\overline{N}\left(e^{tV(0)/\mathrm{}}1\right)\underset{\alpha =1,2}{}\varphi _\alpha ^2_{\mathrm{stoch}}(0).$$
(81)
Using finally the concavity of the logarithmic function, leading to the logarithm of a mean being larger than the mean of the logarithm we obtain a lower bound on the squared error $`\mathrm{\Delta }`$ on the $`N`$body density matrix:

$$\mathrm{\Delta }(t)+\text{Tr}[\rho ^2]A\mathrm{exp}\left[2B\overline{N}\left(e^{tV(0)/\mathrm{}}1\right)\right]$$
(82)
where we have introduced the constant quantities

$`A`$ $`=`$ $`\mathrm{exp}\left[\mathrm{ln}\text{Tr}[\sigma ^{}\sigma ]_{\mathrm{stoch}}(0)\right]`$ (83)
$`B`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1,2}{}}\varphi _\alpha ^2_{\mathrm{stoch}}(0)`$ (84)
It is quite remarkable that in the limit of times shorter than $`\mathrm{}/V(0)`$ this lower bound scales exponentially with time as the upper bound Eq.(74) obtained for Fock states. Consequently the simulation scheme with Fock states is likely to be more efficient that the simulation with coherent states. This will indeed be the case in the numerical examples given in the next sections.

### 3.3 The constant trace schemes

We have given the expression of the one-body density matrix $`\rho ^{(1)}`$ in terms of $`\varphi _\alpha (x)`$ for the simulation with Fock states Eq.(13). This expressions shows that $`\rho ^{(1)}`$ is very sensitive in the large $`N`$ limit to fluctuations of $`\varphi _2|\varphi _1`$. The same remark applies to two-body observables. In order to improve the statistical properties of the simulation one can consider the possibility of a simulation scheme with $`\varphi _2|\varphi _1=1`$ at any time. This actually corresponds to a conserved trace of each single dyadic $`\sigma (t)`$. This possibility is analyzed in §3.3.1; it is extended to the coherent state simulation in §3.3.2, leading to the well-known positive $`P`$-representation formalism.

#### 3.3.1 Simulation with Fock states

Within the Fock state ansatz, the conservation of the trace of the dyadic $`\mathrm{Tr}[\sigma ]`$ can be achieved by (i) choosing the transverse noises $`dB_\alpha ^{}`$ according to the formula Eq.(68) and (ii) using the expression Eq.(65) for the longitudinal noise with $`X=1`$. Point (ii) gives:

$$d\gamma _1=\varphi _2|\varphi _1^1𝑑x\varphi _2^{}(x)𝑑B_1^{}(x).$$
(85)
The forces terms $`\lambda _\alpha `$ and $`F_\alpha `$ are fixed by the conditions (39)-(41):

$$\lambda _1dt=\frac{N1}{2}\varphi _2|\varphi _1^2𝑑x𝑑x^{}\varphi _2^{}(x)\varphi _2^{}(x^{})\overline{dB_1^{}(x)dB_1^{}(x^{})}$$
(86)
and

$$F_1^{}(x)dt=(N1)\varphi _2|\varphi _1^1𝑑x^{}\varphi _2^{}(x^{})\overline{dB_1^{}(x^{})dB_1^{}(x)}.$$
(87)
The expressions for $`d\gamma _2,\lambda _2`$ and $`F_2^{}(x)`$ are obtained by exchanging the indices 1 and 2 in these results.

#### 3.3.2 Simulation with coherent states

In the case of the coherent state ansatz, the value of $`d\sigma _s`$, which is the zero-mean noise term entering the variation of the dyadic $`\sigma `$ during a time step $`dt`$, is given in Eq.(60). The requirement of a constant trace $`\mathrm{Tr}[\sigma ]=\mathrm{\Pi }e^{\overline{N}\varphi _2|\varphi _1}`$ leads to the following condition on the noise terms

$$db+\overline{N}\mathrm{\Pi }𝑑x\left(\varphi _2^{}(x)dB_1(x)+dB_2^{}(x)\varphi _1(x)\right)=0.$$
(88)
We choose the noise terms $`dB_\alpha `$ as in Eq.(78). The remaining parameters $`F_\alpha `$ are now unambiguously determined by (50)-(51):

$`F_1(x)={\displaystyle \frac{1}{i\mathrm{}}}\left[h_0+\overline{N}{\displaystyle 𝑑x^{}\varphi _2^{}(x^{})V(xx^{})\varphi _1(x^{})}\right]\varphi _1(x)`$ (89)
$`F_2(x)={\displaystyle \frac{1}{i\mathrm{}}}\left[h_0+\overline{N}{\displaystyle 𝑑x^{}\varphi _1^{}(x^{})V(xx^{})\varphi _2(x^{})}\right]\varphi _2(x).`$ (90)
This scheme exactly recovers the stochastic evolution in the positive $`P`$-representation, which was originally obtained with a different mathematical procedure .

## 4 Stochastic vs. exact approach for a two-mode model

In order to test the convergence of the stochastic schemes developed in the previous section we now apply this method to a simple two-mode system for which the exact solution of the $`N`$-body Schrödinger equation can also be obtained by a direct numerical integration. This allows (i) to check that the stochastic methods when averaged over many realizations give the correct result indeed, and (ii) to determine the statistical error for each of the four implementations of the stochastic approach (‘constant trace’ vs. ‘simple’, Fock vs. coherent states).

The toy-model that we consider describes the dynamics of two self-interacting condensates coherently coupled one to the other. It can be applied to the case of two condensates separated by a barrier  (Josephson-type coupling) or condensates in two different internal states coupled by an electromagnetic field  (Rabi-type coupling). In this model we restrict the expansion of the atomic field operator to two orthogonal modes,

$$\widehat{\psi }(x)=\widehat{a}u_a(x)+\widehat{b}u_b(x).$$
(91)
The Hamiltonian Eq.(1) takes the simple form:

$$=\frac{\mathrm{}\mathrm{\Omega }}{2}\left(\widehat{a}^{}\widehat{b}+\widehat{b}^{}\widehat{a}\right)+\mathrm{}\kappa \left(\widehat{a}^2\widehat{a}^2+\widehat{b}^2\widehat{b}^2\right)$$
(92)
where $`\widehat{a},\widehat{b}`$ annihilate a particle in modes $`u_a`$ and $`u_b`$, $`\kappa `$ characterizes the strength of the atomic interactions inside each condensate and $`\mathrm{\Omega }`$ is the Rabi coupling amplitude between the two condensates. Here we have restricted for simplicity to the case where (i) the condensates have identical interaction properties, (ii) the interactions between atoms in different wells are negligible, and (iii) the Rabi coupling is resonant. The most general two-mode case could be treated along the same lines.

The direct numerical solution of the Schrödinger equation is performed in a basis of Fock states $`|n_a,n_b`$ with $`n_{a,b}`$ particles in modes $`u_a,u_b`$. The numerical integration is simplified by the fact that $`n_a+n_b`$ is a quantity conserved by the Hamiltonian evolution. We start with a state in which all atoms are in mode $`u_b`$, either in a Fock state $`|n_a=0,n_b=N`$ (for the Fock state simulations) or in a coherent state $`\mathrm{exp}(N^{1/2}\widehat{b}^{})|0,0`$ (for the coherent state simulations). We watch the time evolution of the mean fraction of particles in mode $`u_a`$, $`p_a\widehat{a}^{}\widehat{a}/N`$.

Mean-field theory (the Gross-Pitaevskii equation), valid in the limit $`N1`$ with a fixed $`\kappa N/\mathrm{\Omega }`$ , predicts periodic oscillations of $`\widehat{a}^{}\widehat{a}/N`$; the peak-to-peak amplitude of the oscillations is equal to unity if $`\kappa N/\mathrm{\Omega }<1`$, and is smaller than one otherwise . In the exact solution the oscillations are no longer periodic due to emergence of incommensurable frequencies in the spectrum of $`H`$.

In the simulation method we evolve sets of two complex numbers representing the amplitudes of the functions $`\varphi _1(x)`$ and $`\varphi _2(x)`$ on the modes $`u_{a,b}(x)`$ (plus the $`\mathrm{\Pi }`$ coefficient in the coherent state case). The results are presented in figure 1 for $`N=17`$ particles and $`\kappa /\mathrm{\Omega }=0.1`$, together with the result of the direct integration of the Schrödinger equation.

The first row in the figure concerns the constant trace simulations. Figure 1a shows results of the simulation based on the positive $`P`$-representation, that is the constant trace simulation with coherent states. As well known  this scheme leads to divergences of the norm $`\varphi _1\varphi _2`$ for some realizations of the simulation. We have cut the plot in figure 1 at the first divergence. The same type of divergences occurs in the constant trace simulation with Fock states (figure 1b). Note however that the characteristic time for the first divergence to occur is somewhat longer. We have checked for these constant trace simulations that the probability distribution of $`\varphi _1\varphi _2`$ broadens with time, eventually getting a power law tail. The corresponding exponent $`\alpha `$ decreases in time below the critical value $`\alpha _{\mathrm{crit}}=3`$ for which the variance of $`\varphi _1\varphi _2`$ becomes infinite. This scenario is identical to the one found with the positive $`P`$-representation .

The simple simulation schemes plotted on the second row of figure 1 provide results which are at all time in agreement with the direct integration within the error bars. Contrarily to the constant trace schemes we do not observe finite time divergences in the simple schemes. For a given evolution time we have checked that the error bars scale as $`1/\sqrt{𝒩}`$ where $`𝒩`$ is the number of stochastic realizations. For a given $`𝒩`$ we found that the error bars increase quasi-exponentially with time.

The noise in the simple simulation schemes is investigated in more details in figure 2 which shows the error estimator $`\text{Tr}[\sigma ^{}\sigma ]_{\mathrm{stoch}}`$ as function of time, for coherent states in figure 2a and for Fock states in figure 2b. The coherent state result confirms the prediction Eq.(82). The Fock state result is found to be notably smaller than the upper bound Eq.(74). This is due to the fact that the terms proportional to $`\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha `$ in Eq.(73) are not negligible as compared to the term $`V(0)`$. We have checked these conclusions for various values of $`N`$ and $`\kappa /\mathrm{\Omega }`$.

For a large number of particles it is known  that the oscillations of $`\widehat{a}^{}\widehat{a}`$ experience a collapse followed by revivals. These revivals are purely quantum phenomena for the field dynamics and they cannot be obtained in classical field approximation such as the Gross-Pitaevskii equation. We expect to see a precursor of this phenomenon even for the small number of particles $`N=17`$. As the simple scheme simulation with Fock states is the most efficient of the four schemes for the investigation of the long time limit, we have pushed it to the time at which a “revival” can be seen, as shown in figure 3. This figure is obtained with $`𝒩=10^8`$ simulations.

## 5 Stochastic approach for a one-dimensional Bose gas

The interacting Bose gas is in general a multi-mode problem, and the simulation schemes may have in this case a behavior different from the one in a few-mode model such as in §4. We have therefore investigated a model for a one-dimensional Bose gas. The atoms are confined in a harmonic trap with an oscillation frequency $`\omega `$. They experience binary interactions with a Gaussian interaction potential of strength $`g`$ and range $`b`$:

$$V(xy)=\frac{g}{(2\pi )^{1/2}b}\mathrm{exp}\left[(xy)^2/(2b^2)\right].$$
(93)
At time $`t=0`$ all the atoms are in the same normalized state $`\varphi `$ solution of the time independent Gross-Pitaevskii equation

$$\mu \varphi (x)=\frac{\mathrm{}^2}{2m}\frac{d^2\varphi }{dx^2}+\frac{1}{2}m\omega ^2x^2\varphi (x)+(N1)𝑑x^{}V(xx^{})|\varphi (x^{})|^2\varphi (x).$$
(94)
At time $`t=0^+`$ the trap frequency is suddenly increased by a factor two, which induces a breathing of the cloud .

This expected breathing is well reproduced by the numerical simulations. The mean squared spatial width $`R^2`$ of the cloud as function of time is obtained by taking $`\widehat{O}=_{k=1}^N\widehat{x}_k^2/N`$ in Eq.(11) where $`\widehat{x}_k`$ is the position operator of the $`k`$th particle. The quantity $`R^2`$ is shown in figure 4 for the simulation schemes with Fock states. One recovers the key feature of the constant trace simulation, that is a divergence of the norm $`\varphi _1\varphi _2`$ in finite time for some realizations. Before the occurrence of the first divergence the stochastic variance of the size squared of the cloud, defined as

$$\text{Var}(R^2)_{\mathrm{stoch}}=\frac{1}{𝒩}\underset{i=1}{\overset{𝒩}{}}\left[R_i^2(t)R^2(t)\right]^2\text{with}R_i^2(t)=e\left[N:\varphi _2^{(i)}(t)|\widehat{O}|N:\varphi _1^{(i)}(t)\right]$$
(95)
is notably smaller in the constant trace scheme than in the simple scheme, as shown in figure 5a. This contrast between the two schemes for the statistical error on one-body observables was absent in the two-mode model of §4.

We have also investigated the noise on the $`N`$body density matrix characterized by $`\text{Tr}[\sigma ^{}\sigma ]_{\mathrm{stoch}}`$ (see figure 5b). As expected this error indicator is smaller with the simple scheme. For this simple scheme it varies quasi exponentially with time with an exponent $`\gamma 4\omega `$, which is smaller by a factor roughly 2 than the one of the upper bound Eq.(74). This difference is due to the fact that the range $`b`$ of the interaction potential is chosen here of the same order as the size $`R`$ of the cloud so that the terms $`\varphi _\alpha ,\varphi _\alpha |V|\varphi _\alpha ,\varphi _\alpha `$ neglected in the derivation of the upper bound are actually significant. We have checked for various ranges $`b`$ much smaller than $`R`$ that $`\gamma `$ then approaches the upper bound $`2NV(0)/\mathrm{}`$.

## 6 Conclusion and perspectives

In this paper we have investigated a general method to solve exactly the $`N`$body problem in the bosonic case. The principle of the approach is to add to the usual mean-field Gross-Pitaevskii equation a fluctuating term. We have determined the general conditions ensuring that the average over all possible realizations of this stochastic equation reproduces the exact $`N`$body Schrödinger equation.

This idea already received a particular implementation in quantum optics, in the frame of the positive $`P`$-representation. We recover here the scheme based on the positive $`P`$-representation as a particular case of a simulation evolving coherent states of the bosonic field with the constraint that the trace of the density operator should remain exactly equal to unity for each single realization. This provides a simple derivation of the stochastic evolution within this representation alternative to the usual one  based on analyticity properties.

Among the many possible implementations of the general stochastic approach we have also investigated schemes evolving Fock states (that is number states) rather than coherent states. This is well suited to situations where the total number of particles is conserved. In particular we have identified a scheme preserving exactly the trace of the density operator which is for number states the counterpart of the one based on the positive $`P`$-representation.

Schemes with constant trace are subject to divergence of the norm of some realizations in finite time. This effect, already known in the context of the positive $`P`$-representation , makes these schemes difficult to use.

In order to overcome this divergence problem we have investigated schemes in which the condition on the trace is relaxed. We have chosen instead to minimize the statistical spread on the $`N`$body density matrix, which gave rise to the ‘simple’ schemes, either with coherent states or Fock states. In this case the $`N`$body density operator is obtained as a stochastic average of dyadics such as $`|\text{coh}:\overline{N}^{1/2}\varphi _1\text{coh}:\overline{N}^{1/2}\varphi _2|`$ or $`|N:\varphi _1N:\varphi _2|`$, where the evolutions of $`\varphi _1`$ and $`\varphi _2`$ are fully decoupled. The deterministic parts are given by Gross-Pitaevskii equations, which preserves the norm of $`\varphi _{1,2}`$, contrarily to the case of constant trace schemes. The decoupling between the evolutions of $`\varphi _1`$ and $`\varphi _2`$ allows a reinterpretation of our representation of the $`N`$body density operator. If the initial density operator is given by $`\rho (t=0)=|N:\varphi _0N:\varphi _0|`$ it will be given at time $`t`$ by

$$\rho (t)=|\mathrm{\Psi }(t)\mathrm{\Psi }(t)|$$
(96)
with the $`N`$particle state vector

$$|\mathrm{\Psi }(t)=\text{lim}_𝒩\mathrm{}\frac{1}{𝒩}\underset{j=1}{\overset{𝒩}{}}|N:\varphi ^{(j)}(t).$$
(97)
In this expression $`\varphi ^{(j)}`$ are stochastic realizations with the initial condition $`\varphi ^{(j)}(t=0)=\varphi _0`$.

The ‘simple’ schemes have much better stability properties than the constant trace schemes: differently from the case of constant trace schemes, the deterministic evolution of the ‘simple’ schemes has a Gross-Pitaevski form and thus conserves the norms $`\varphi _{1,2}`$. This condition, together with the upper bound $`\zeta _{1,2}V(0)dt\varphi _{1,2}^2/\mathrm{}`$ on the eigenvalues $`\zeta _{1,2}`$ of the noise covariance operator $`\overline{dB_\alpha (x)dB_\alpha ^{}(x^{})}`$, can be used to prove that the stochastic equations possess a finite, non-exploding solution valid for all times (see  in page 94,  in §4.5).

We have numerically applied the simulation schemes to a two-mode model and to a one-dimensional Bose gas. In both cases we found that the constant trace schemes lead to some diverging realizations, while the simple schemes lead to a statistical spread on the $`N`$body density operator increasing exponentially with time with an exponent $`\gamma NV(0)`$. The simple schemes are therefore not well suited to determine small deviations from the mean-field approximation in the large $`N`$ limit but can be more efficiently applied to systems with a small number of particles, such as small atomic clouds tightly trapped at the nodes or antinodes of an optical lattice.

In the one-dimensional numerical example of this paper we have presented results for a simple one-body observable, the size of the atomic cloud. We have actually extended the calculations to more elaborate observables such as the first order and the second order correlation functions of the field. We have not presented the results here as the initial state of the gas was taken to be a (not very physical) Hartree-Fock state. We are presently working on the possibility to generate a more realistic initial state such as a thermal equilibrium for the gas, by extending our stochastic approach to evolution in imaginary time.

This work has several possible perspectives of extension. One can first use as building block a more sophisticated ansatz than the Hartree-Fock state $`|N:\varphi `$, such as Bogoliubov vacua (that is squeezed states of the atomic field) or a multimode Hartree-Fock ansatz (that is an arbitrary coherent superposition of number states in several adjustable modes of the field). One can also look for approximate rather than exact stochastic solutions to the $`N`$body problem but that would be better than mean-field approaches in some given situations. We hope to address some of these perspectives in the near future.