A Sinatra

Y Castin

E Witkowska

Non-diffusive phase spreading of a Bose-Einstein condensate at finite temperature

Keywords: numbers: 03, 75, Kk, 03, 75, Pp

We show that the phase of a condensate in a finite temperature gas spreads linearly in time at long times rather than in a diffusive way. This result is supported by classical field simulations, and analytical calculations which are generalized to the quantum case under the assumption of quantum ergodicity in the system. This super-diffusive behavior is intimately related to conservation of energy during the free evolution of the system and to fluctuations of energy in the prepared initial state.

I. INTRODUCTION

Phase coherence is one of the fundamental properties of Bose-Einstein condensates. It is also a key feature in the present developments of the research on condensates which, ten years after the first experimental realization, go in the direction of integrating this powerful tool into other branches of physics, of which metrology and quantum information are two promising examples [1].

The problem of the condensate phase dynamics due to atomic interactions at zero temperature has been analyzed by different authors in theory [2] and in experiment [3,4,5]. It is now well understood that an initially prepared relative phase between two condensates will spread in time due to the corresponding uncertainty in the relative particle number as the relative phase and the relative particle number are conjugate variables. The phase dynamics of a two component condensate in realistic situations including harmonic traps, non stationarity and fluctuations in the total number of particles was analyzed in [6], where a comparison to the experiments of [4] is also performed. An important conclusion was that the zero temperature theory could not account for the coherence times observed in experiment, which raises the question of the role of the non-condensed fraction.

In this paper we address the fundamental problem of phase spreading of a Bose-Einstein condensate in a finite temperature atomic gas. In order to obtain simple and general results, we consider the ideal case of a spatially uniform condensate at thermodynamic equilibrium, and we assume that one has access to the first order temporal correlation function a † 0 (t)a 0 of the component a 0 of the atomic field in the condensate mode. In real life, the situation is more complex: the atoms are trapped in harmonic potentials, and the measurement of phase coherence is a delicate procedure, usually relying on the interference between two condensates [4]. In the literature two well distinct predictions exist for the long time spreading of the condensate phase at finite temperature, either a diffusive behavior (variance growing linearly in time) [7,8,9,10] or a ballistic behavior (variance growing quadratically in time) [11]. We study this problem first with a classical field model [12,13,14,15,16,17], where exact numerical simulations can be performed. We then explain the numerics analytically, and extend the analytical approach to the quantum case.

The important result that we obtain is that the variance of the phase increases quadratically in time. This is at variance with the prediction of phase diffusion from the "quantum optics" open system approaches of [7,8,9,10] assuming the condensate to evolve under the influence of Langevin short memory fluctuating forces. Our prediction results from two ingredients, (i) the system is prepared in an initial state with an energy fluctuating from one experimental realization to the other, here sampling the canonical ensemble, and (ii) the system is isolated in its further evolution and therefore keeps a constant energy. As we shall see, the combination of these two ingredients prevents some temporal correlation functions to vanish at long times. Our prediction qualitatively agrees with the one of [11], but not quantitatively, as we obtain a different expression for the long time limit of the variance of the phase over the time squared. This difference is due to the fact that we take into account ergodicity in the system resulting from the interactions among Bogoliubov modes such as the Beliaev-Landau processes.

In section II we present the classical field model; numerical predictions for this model are presented in section III, and analytical results reproducing the numerics at short or long times are given in section IV. These analytical results are extended to the case of the quantum field in section V. We conclude in section VI.

II. THE CLASSICAL FIELD MODEL

In this section we develop a classical field model that has the advantage that it can be exactly simulated numerically. This will allow us to understand the physics governing the spreading of the condensate phase and to test the validity of various approximations, paving the way to the quantum treatment.

We consider a lattice model for a classical field ψ(r) in three dimensions. The lattice spacings are l 1 , l 2 , l 3 along the three directions of space and dV = l 1 l 2 l 3 is the volume of the unit cell in the lattice. We enclose the atomic field in a spatial box of sizes L 1 , L 2 , L 3 and volume V = L 1 L 2 L 3 , with periodic boundary conditions. The discretized field has the following Poisson brackets ih{ψ(r 1 ), ψ * (r 2 )} = δ r1,r2 dV

where the Poisson brackets are such that df /dt = {f, H} for a time-independent functional f of the field ψ. The field ψ may be expanded over the plane waves

ψ(r) = k a k e i k•r √ V , (2)
where k is restricted to the first Brillouin zone, k α ∈ [-π/l α , π/l α [ where α labels the directions of space. We assume that, in the real physical system, the total number of atoms is fixed, equal to N . In the classical field model, this fixes the norm squared of the field:

dV r |ψ(r)| 2 = N. (3)
Equivalently the density of the system

ρ = N V (4)
is fixed for each realization of the field. The evolution of the field is governed by the Hamiltonian

H = k Ẽk a * k a k + g 2 r dV ψ * (r)ψ * (r)ψ(r)ψ(r), ( 5
)
where Ẽk is the dispersion relation of the non-interacting waves, and the binary interaction between particles in the real gas is reflected in the classical field model by a field self-interaction with a coupling constant g = 4πh 2 a/m, where a is the s-wave scattering length of two atoms.

In general, we expect the predictions of a classical field model to be cut-off dependent, i.e. the predictions of our model may depend on the lattice spacings l α . We use here a refinement to the usual classical field model, which makes it cut-off independent for some observables like the condensate fraction, a quantity expected to play an important role here. An obvious example of a quantity which will remain cut-off dependent is the mean value of the Hamiltonian H in thermal equilibrium.

Let us consider first the non-interacting case (g = 0) in presence of a condensate. For a thermalized classical field the occupation numbers of the excited plane wave modes are given by the equipartition formula

a * k a k = k B T Ẽk . (6)
We adjust the dispersion relation Ẽk in order to reproduce the Bose law for the occupation numbers of the quantum field in the Bose-condensed regime:

1 e βh 2 k 2 /2m -1 = k B T Ẽk . (7)
For all modes with large occupation number Ẽk ≃ h2 k 2 /2m, while the occupation of modes with h2 k 2 /2m ≫ k B T , whose quantum dynamics is not well approximated by the classical field model anyway, is exponentially suppressed as in the quantum theory.

In the interacting case, one could adapt the same trick of a modified dispersion relation, by including the fact that the relevant spectrum is not h2 k 2 /2m but the Bogoliubov spectrum [START_REF]To get the equation defining Ẽk , one would then replace[END_REF]. The resulting Ẽk would now start growing exponentially with k when the Bogoliubov energy

[(h 2 k 2 /2m)(2ρg + h2 k 2 /2m)] 1/2 reaches k B T .
In the classical field model we restrict our analysis to the regime k B T ≫ ρg so that at energies of the order of k B T , the Bogoliubov energy is dominated by the kinetic term h2 k 2 /2m. One can then simply use in the Hamiltonian the modified dispersion relation Ẽk as given by Eq. (7). This is what we did in the simulations of this paper, so that the classical field ψ evolves according to the non-linear equation [START_REF]∆ stands here for the operator acting on functions on the lattice[END_REF]:

ih ∂ t ψ = k B T exp -β h2 2m ∆ -1 + g|ψ(r, t)| 2 ψ . (8)
In practice this equation is integrated numerically with the FFT splitting technique.

We then introduce the density and the phase of the condensate mode

a 0 = e i θ N 0 . (9)
In what follows, we concentrate on three physical quantities: the condensate amplitude correlation function

a * 0 (t) a 0 (0) , (10)
the condensate atom number correlation function δN 0 (t) δN 0 (0) where δN 0 = N 0 -N 0 , (11) and the variance of the condensate phase change during t:

Var ϕ(t) = ϕ(t) 2 -ϕ(t) 2 where ϕ(t) = θ(t)-θ(0). (12)
The averages are taken over stochastic realizations of the classical field, as the initial field samples a thermal probability distribution.

III. CLASSICAL FIELD: NUMERICAL RESULTS

We consider a gas of N = 4 × 10 5 atoms with ρg = 700 h2 /mV 2/3 in a box of non commensurable square lengths to guarantee efficient ergodicity in the system, in the ratio

L 2 1 : L 2 2 : L 2 3 = √ 2 : (1 + √ 5)/2 : √ 3.
We choose the number of the lattice points in a temperature dependent way, such that the maximal Bogoliubov energy [(h 2 k 2 /2m)(2ρg + h2 k 2 /2m)] 1/2 on the lattice is equal to 3k B T .

To generate the stochastic initial values of the classical field we proceed as follows. (i) For each realization, we generate a non condensed field ψ ⊥ (r) at temperature T in the Bogoliubov approximation as explained in [START_REF] Sinatra | [END_REF]. In practice we generate complex numbers {b k } for each vector k on the grid according to the probability distribution

P (b k ) = 1 π ǫk k B T e -(|b k | 2 ǫk /kB T ) (13)
where ǫk = [ Ẽk ( Ẽk + 2ρg)] 1/2 . With a set of {b k } for a given realization we build the non condensed field

ψ ⊥ (r) = e iθ k =0 b k Ũk e ik•r √ V + b * k Ṽk e -ik•r √ V ( 14
)
where the initial value of the condensate phase θ is randomly chosen with the uniform law in [0, 2π[, and where the real amplitudes Ũk , Ṽk , normalized as Ũ 2 k -Ṽ 2 k = 1, are given by the usual Bogoliubov theory, here with the modified dispersion relation, so that

Ũk + Ṽk = Ẽk Ẽk + 2ρg 1/4 . (15)
(ii) We create the classical field with the constraint that the total number of atoms N is fixed:

ψ(r) = a 0 √ V + ψ ⊥ (r) ( 16
)
where a 0 = √ N -N ⊥ e iθ , N ⊥ is the number of non condensed atoms,

N ⊥ = r dV |ψ ⊥ (r)| 2 . ( 17
)
(iii) We let the field evolve for some time interval with the Eq.( 8) to eliminate transients due to the fact that the Bogoliubov approximation used in the sampling does not produce an exactly stationary distribution. After this 'thermalization' period we start calculating the relevant observables, as ψ evolves with the same Eq.( 8). First we investigate the mean condensate phase change ϕ (t). We find a linear dependence with time, with a slope slightly different from the value -ρg/h naively expected, e.g. from the zero temperature Gross-Pitaevskii equation. The slope difference is temperature dependent and is expected physically to correspond to the discrepancy between the zero temperature chemical potential ρg and the actual finite temperature one µ(T ). This we shall confirm using Bogoliubov theory in Sec. IV (see also [21]).

In figure 1, we show the real part of the amplitude correlation function of the condensate a * 0 (t)a 0 (0) as a function of time, for a temperature T = 0.17T c , where T c is the critical temperature of the ideal gas. The zerotemperature evolution e iρgt/h is removed so that the oscillations in the figure are due to the above mentioned effect µ(T ) = ρg. Due to the finite temperature in the system, the correlation function of the condensate amplitude is smeared out at long times. (10) normalized to its t = 0 value, and divided by the zero temperature evolution e iρgt/h , as a function of time: in the vertical axis label, ã0(t) stands for a0(t)e iρgt/h . (a) Short times behavior and (b) long times behavior. In solid line from an average over 500 solutions of Eq.( 8), in dashed line (red) the Bogoliubov approximation (40). Here the temperature is kBT = 3077.3 h2 /mV 2/3 = 0.1711Tc, where Tc is the critical temperature kBTc = (2πh 2 /m)(ρ/ζ(3/2)) 2/3 of the ideal gas, the number of particles is N = 4 × 10 5 and the coupling constant is such that the zero-temperature chemical potential is ρg = gN/V = 700h 2 /mV 2/3 .

Correspondingly the standard deviation of the condensate phase change increases with time, as we show in figure 2 for five different values of the temperature, up to T = 0.65T c . In all cases, at long times, we observe a quadratic growth of Var ϕ contrarily to the phase diffusion behavior ∝ t predicted in the literature [7,8,9,10].

To complete the physical picture, we show in figure 3 the correlation function of the condensate atom number (11). At very short times, see the beginning of the curves in Fig. 3a, the simulation (square symbols) confirms the Bogoliubov prediction (dashed oscillating line); at long times, see Fig. 3b, the correlation function drops to a value significantly smaller than the Bogoliubov prediction (fast oscillations are not shown in the figure); a key point is that this long time value of the correlation function of the condensate atom number is not zero.

One may fear at this stage that the classical field model is missing some source of damping in the dynamics of the system. However it is a well established fact that the classical field model is able to simulate damping processes, including the finite temperature Beliaev-Landau processes [22,23,24,25], since the interaction among the Bogoliubov modes is included in this model [12,13,15,16,[START_REF] Sinatra | [END_REF]21,26,27,28,29]. More quantitatively we now check that the damping times due to the Beliaev-Landau processes in the simulation are much shorter than the evolution times considered here. To this end, we extract from the simulations the temporal correlation

functions b * k (t)b k (0) and |b k | 2 (t)|b k | 2 (0) -|b k | 2 2
, obtained by projecting the classical field over the corresponding Bogoliubov mode and averaging over many realizations. We show these correlation functions for the lowest energy Bogoliubov mode and for an excited Bogoliubov mode in Fig. 4.

We come then into a paradox. On one side, the various Bogoliubov oscillators b k decorrelate at long times. On the other side, the variance of the phase change ϕ of the condensate varies quadratically at long times, which implies, as we shall see in Sec. IV, that the derivative of the phase φ does not decorrelate at long times, although it is a function of the b k 's; similarly, the fluctuations of the number of condensate atoms δN 0 , which are functions of the b k 's, do not decorrelate at long times. This paradox will be explained in Sec. IV, and quantitative predictions for long times behavior of the condensate atom number correlation function and of the variance of the condensate phase change will be derived. Anticipating these analytical results, we show in Fig. 5a the long time limit of (Var ϕ) 1/2 /t as a function of T /T c , from the results of the classical field simulations, but also from the predictions of the Bogoliubov approximation Eq.(38), and of the ergodic theory of Sec. IV. In figure 5b we show the same results and predictions for the asymptotic value of the condensate atom number correlation function.

X = X * (t)X(0) -| X | 2 .
The purple dashed-dotted line is an exponential function of t, given by Eq.( 24). It reproduces well the simulation results, and we have checked that the agreement is good for all vectors k lying in an arbitrarily chosen plane in k-space. The parameters N , T and ρg have the same values as in Fig. 1.

IV. CLASSICAL FIELD: ANALYTICAL RESULTS

The general procedure used here to obtain analytical results is the following. First one expresses the quantity of interest (the number of condensate atoms or the time derivative of the condensate phase) in terms on the amplitudes b k of the field ψ over the Bogoliubov modes,

b k (t) = dV r Ũk e -ik•r √ V e -iθ(t) ψ ⊥ (r, t) + Ṽk e ik•r √ V e iθ(t) ψ * ⊥ (r, t) (18)
where ψ ⊥ is the component of ψ orthogonal to the condensate mode. Second one evaluates the correlation functions of products of b k in various physical limits.

A. Correlation function of the condensate atom number

As the total number of particles is fixed, it is equivalent to calculate the correlation function of δN 0 in Eq.( 11) and of the number of non-condensed particles N ⊥ . Injecting the expansion Eq.( 14) for the time dependent non-condensed field ψ ⊥ over the Bogoliubov modes into Eq.( 17) we obtain

N ⊥ (t) = k =0 | Ũk b k (t) + Ṽk b * -k (t)| 2 . ( 19
)
Bogoliubov theory: In Bogoliubov theory interaction among the Bogoliubov modes is neglected so that at all times

b k (t) = b k (0) e -iω k t with ω k = ǫk /h . ( 20
)
As Wick's theorem applies for the initial thermal distribution we obtain for the correlation function of the condensate atom number:

δN 0 (t) δN 0 (0) Bog = k =0 ñ2 k | Ũ 2 k e iω k t + Ṽ 2 k e -iω k t | 2 +( Ũk Ṽk ) 2 |e iω k t + e -iω k t | 2 (21)
where ñk = k B T /ǫ k is the Bogoliubov mean occupation number of a mode for the classical field. At very short times, a good agreement of the Bogoliubov prediction with the simulation is observed in Fig. 3a. Smearing out the terms oscillating rapidly at Bohr frequencies 2ω k , we obtain a prediction directly comparable to the coarse grained numerical result of Fig. 3b:

δN 0 (t) δN 0 (0) Bog non osc = k =0 ( Ũ 2 k + Ṽ 2 k ) 2 ñ2 k . ( 22
)
This amounts to considering the correlation function of

N non osc ⊥ (t) = k =0 ( Ũ 2 k + Ṽ 2 k )b * k (t)b k (t), (23)
deduced from ( 19) by eliminating the oscillating terms such as b k b -k . As can be seen in Fig. 3b, Bogoliubov theory fails at long times. Note that in the thermodynamic limit, where the above sum is dominated by the low k terms, one may approximate Ṽk ∼ -Ũk , so that Eq.( 22) is roughly half of the t = 0 value of Eq.( 21); in other words, it is approximately half of the variance of the condensate number. In the numerical result of Fig. 3, the correlation function drops by much more than a factor 2.

Gaussian theory: A possible approach to improve Bogoliubov theory consists in assuming that the b k are Gaussian variables with a finite time correlation due to the Beliaev-Landau mechanism:

| b * k (t)b k (0) | 2 = ñ2 k e -2γ k |t| (24)
where γ k is calculated with time dependent perturbation theory including the discrete nature of the spectrum as in [START_REF] Sinatra | [END_REF]. This amounts to weighting each term of Eq.( 22) by exp(-2γ k |t|). This assumption is supported by numerical evidence for a single mode, see Fig. 4, and by an analytic derivation in the thermodynamic limit for one or two modes, see Appendix A. Nevertheless, the resulting prediction for the correlation function of N 0 , while looking promising at short times, see Fig. 3a, is in clear disagreement with the simulation at long times, see Fig. 3b. Since the assumption of a long time decorrelation of b * k (t) with b k (0) is physically reasonable, one may suspect that the Gaussian hypothesis is not accurate when a large number of modes are involved as for the correlation function of N 0 . This is indeed the case, as we now show. Ergodic theory: A systematic way to calculate the long time limit of the correlation function is to assume that the non-linear dynamics generated by Eq.( 8) is ergodic: at long times, the b k (t)'s for a given realization of the field explore uniformly a fixed energy surface in phase space [30]. In the Bogoliubov approximation for the energy, this means that the b k (t)'s sample the unnormalized probability distribution

P ∞ ({b k }) = δ   E - k =0 ǫk b * k b k   (25)
where the Bogoliubov energy E is fixed by the initial value of the field:

E = k =0 ǫk b * k (0)b k (0). ( 26
)
First, for a given initial condition of the field, we calculate the expectation value of N ⊥ (t) as given by Eq.( 19) over the ergodic distribution Eq.( 25), which is equivalent to the temporal average of N ⊥ (t) over an infinite time interval. The terms of the form b b or b * b * have a zero mean, since the phases of the b k 's are uniformly distributed over 2π, according to Eq.( 25). To calculate the expectation value of the b * b terms, it is convenient to introduce rescaled variables

B k = ǫk E 1/2 b k . (27)
According to Eq.( 25) the real parts and the imaginary parts of all the B k are uniformly distributed over the unit hypersphere in a space of dimension 2M, where M = V /dV -1 is the number of Bogoliubov modes so that we obtain |B k | 2 = 1/M where the overline stands for the average over the ergodic distribution (25). As a consequence the ergodic average of N ⊥ is

N ⊥ = 1 M k =0 ( Ũ 2 k + Ṽ 2 k ) E ǫk . ( 28
)
Note that this ergodic average depends on the t = 0 value of the b k 's via (26). Second, we average the product N ⊥ N ⊥ (0) over the thermal canonical distribution for the initial values b k (0). This gives the long time limit of the correlation function of the number of condensate atoms:

δN 0 (t → +∞)δN 0 (0) ergo = 1 M   k =0 ( Ũ 2 k + Ṽ 2 k )ñ k   2 .
(29) This prediction is in good agreement with the simulations at long times, see Fig. 3b for a fixed value of the temperature, and Fig. 5b as a function of temperature. Note that, according to Schwartz inequality, the ergodic value is lower than the coarse grained Bogoliubov prediction Eq.( 22), as was expected physically.

This clearly shows that the existence of infinite time correlations in the number of condensate atoms is a consequence of the conservation of energy during the free evolution of the system.

To understand the failure of the Gaussian model, we give the ergodic prediction of the long-time limit of the correlation function of the Bogoliubov mode occupation numbers

n k = |b k | 2 , n k (t → +∞)n k ′ (0) ergo -n k n k ′ = ñk ñk ′ M . (30)
This long-time value is non-zero, contrarily to the Gaussian model prediction. One may argue that the value Eq.(30) tends to zero in the thermodynamic limit, so that the error in the Gaussian model looks negligible for a large system. However, in calculating the correlation function of a macroscopic quantity such as N ⊥ , a double sum over the Bogoliubov modes appears, so that the small deviations Eq.( 30) from the Gaussian model prediction sum up to a macroscopic value. In other words, in the calculation of a given correlation function, one is not allowed to take the thermodynamic limit before the end of the calculation.

B. Variance of the condensate phase change

To reproduce the approach of the previous subsection for the phase, one should express the phase change ϕ(t) of the condensate amplitude a 0 as a function of the b k 's. It turns out that the quantity easily expressed in terms of the b k 's is the time derivative φ. The variance of ϕ is then related to the correlation function C of φ:

Var ϕ = t 0 dτ t 0 dτ ′ C(|τ -τ ′ |) (31)
where time translational invariance in steady state imposes for a classical field that C depends only on |τ -τ ′ |:

C(|τ -τ ′ |) = φ(τ ) φ(τ ′ ) -φ(τ ) φ(τ ′ ) . (32)
If C(τ ) → 0 fast enough when τ → ∞ then Var ϕ grows linearly in time. On the other hand, if C(τ ) has a nonzero limit at long times, then Var ϕ grows quadratically in time [31].

To express φ in terms of the b k 's, we write the equation of motion for a 0 :

ih ȧ0 = ih{a 0 , H} = ∂ a * 0 H = g √ V r ψ * (r)ψ 2 (r) ( 33
)
where we used ∂ a * 0 ψ * (r) = 1/ √ V obtained from Eq.( 2). We split ψ as in Eq.( 16); we eliminate the condensate amplitude in the resulting expression for ȧ0 /a 0 (i) by using |a 0 | 2 = N -N ⊥ , where N ⊥ is a function of the b k 's, see Eq. [START_REF]∆ stands here for the operator acting on functions on the lattice[END_REF], and (ii) by introducing the field [32] Λ(r) = e -i θ ψ ⊥ (r) [START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School, Quantum Gases in Low Dimensions[END_REF] which is a function of the b k 's only according to Eq.( 14). This leads to

ih ȧ0 a 0 = ρg + g V r dV Λ(r) 2 + |Λ(r)| 2 + g √ V r dV Λ * (r)Λ 2 (r) √ N -N ⊥ . ( 35
)
The real part of the above equation gives -h θ, which is also -h φ.

Restricting to a weak non-condensed fraction, we drop the cubic terms in Eq. [START_REF] Mora | [END_REF], to obtain [33] h φ ≃ -ρg -

1 2 g V r dV [Λ(r) + Λ * (r)] 2 = -ρg - 1 2 g V k =0 ( Ũk + Ṽk ) 2 |b k + b * -k | 2 . ( 36
)
It turns out that the products b k b -k generate oscillating terms which do not contribute to a coarse grained time average. It is thus useful to define

h φnon osc = -ρg - g V k =0 ( Ũk + Ṽk ) 2 |b k | 2 . ( 37
)
Bogoliubov theory: By using ( 20) and Wick's theorem we calculate the correlation function of Eq. (36). By temporal integration we obtain the variance of the condensate phase change

(Var ϕ) Bog = g hV 2 k =0 ( Ũk + Ṽk ) 4 ñ2 k t 2 + sin 2 ω k t (2ω k ) 2 .
(38) Qualitatively Bogoliubov theory correctly predicts a quadratic growth of the variance of ϕ at long times. As we show in Fig. 5a, however, it is not fully quantitative: it does not reproduce the value of the dephasing rate obtained from the simulations. This is not surprising as in the full non linear theory the b k 's interact and do not follow Eq. [START_REF] Sinatra | [END_REF].

To be complete, we also give the Bogoliubov approximation for the correlation function of the condensate amplitude a 0 . Neglecting the fluctuations of the modulus of a 0 , one can set

a * 0 (t) a 0 (0) ≃ N 0 e -i ϕ(t) . ( 39
)
Dropping the oscillating terms in b k b -k and b * k b * -k in φ(t), which give a small contribution, we get

a * 0 (t) â0 (0) Bog ≃ N 0 k =0 1 1 + i g hV ( Ũk + Ṽk ) 2 ñk t .
(40) The resulting expression is plotted as a dashed line in Fig. 1 against the result of the simulation. Gaussian theory: If we add by hand a decorrelation of the b k 's and assume Gaussian statistics, we get a diffusive spreading of the condensate phase change, with the variance of ϕ non osc growing linearly at long times:

(Var ϕ) Gauss = g hV 2 k =0 ( Ũk + Ṽk ) 4 ñ2 k e -2γ k t -1 2γ 2 k + t γ k , (41)
in clear contradiction with the numerical simulations. This prediction corresponds to a correlation function C vanishing at long times, whereas the correct correlation function has a finite limit, see Fig. 6. Ergodic theory: as in subsection IV A we calculate the long time value of the correlation function for φ using the ergodic assumption. The various steps of the calculation are rigorously the same as in Sec. IV A and lead to

C(τ → +∞) ergo = g hV 2 1 M   k =0 ( Ũk + Ṽk ) 2 ñk   2 .
(42) This prediction is in excellent agreement with the simulations: it gives the correct asymptotic value of C, see Fig. 6, and from the asymptotic expression Var ϕ ≃ C(+∞)t 2 it gives the correct values of the long time limit of (Var ϕ) 1/2 /t, see Fig. 5a, as a function of temperature.

V. QUANTUM TREATMENT: ANALYTICAL RESULTS

So far the classical field model was very useful in revealing the physical processes governing the long time behavior of the phase and atom number fluctuations in the condensate. However it is not a fully quantitative theory, as the long time limits of the correlation functions considered here depend on the precise choice of the energy cut-off, that is on the number of Bogoliubov modes M in the simulation, as is apparent on Eqs. (29,42). In this section, we therefore adapt the previous physical reasonings to the quantum field case.

A. The quantum model

We use a straightforward generalization of the classical field lattice model, taking here for simplicity a cubic lattice, as discussed in [START_REF] Sinatra | [END_REF][START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School, Quantum Gases in Low Dimensions[END_REF][START_REF] Mora | [END_REF]. The bosonic field ψ evolves according to the Hamiltonian

H = k h2 k 2 2m â † k âk + g 0 2 r dV ψ † ψ † ψ ψ, (43)
where âk annihilates a particle of wavevector k in the first Brillouin zone. The dispersion relation of the wave is now the usual one. The total number of atoms is fixed, equal to N . The coupling constant g 0 depends on the lattice spacing l in order to ensure a l-independent scattering length for the discrete delta interaction potential among the particles [START_REF] Castin | Lecture notes of the 2003 Les Houches Spring School, Quantum Gases in Low Dimensions[END_REF][START_REF] Mora | [END_REF]:

V (r 1 -r 2 ) = g 0 dV δ r1,r2 . (44)
Since we consider here the weakly interacting regime, we can restrict to a lattice spacing much larger than the scattering length a so that g 0 is actually very close to g = 4πh 2 a/m. To be able to use Bogoliubov theory as we did in the classical field reasoning, we restrict to the low temperature regime T ≪ T c with a macroscopic occupation of the condensate mode. We thus neglect the possibility that the condensate is empty, which allows us to use the modulus-phase representation of the condensate mode:

â0 ≃ e i θ N0 (45)
where N0 = â † 0 â0 and where θ is a Hermitian 'phase' operator obeying the commutation relation

[ N0 , θ] = i. (46)
This allows to consider the correlation of the condensate atom number fluctuation δ N0 ≡ N0 -N0 but also the variance of the condensate phase change φ(t) ≡ θ(t)θ(0), as we did for the classical field.

B. Correlation function of the condensate atom number

To predict the correlation function of δ N0 , we use Bogoliubov theory at short times and the quantum analog of the ergodic theory at long times.

In the number conserving Bogoliubov theory [32,36], written here for a spatially homogeneous system, one introduces the field conserving the total number of particles Λ(r) ≡ e -i θ ψ⊥ (r) (47)

where the non-condensed field ψ⊥ is obtained by projecting out the component of the field ψ on the condensate mode. The field Λ then admits the modal expansion on the Bogoliubov modes

Λ(r) = k =0 bk U k e ik•r √ V + b † k V k e -ik•r √ V ( 48
)
where the real amplitudes

U k , V k , normalized as U 2 k - V 2 k = 1
, are given by the usual Bogoliubov theory,

U k + V k = h2 k 2 /2m 2ρg 0 + h2 k 2 /2m 1/4 . ( 49
)
Since the total number of particles is fixed to N , it is equivalent to consider the fluctuations of N0 or of the number of non-condensed atoms N⊥ = r dV Λ † (r) Λ(r).

(50)

This, together with the expansion (48), expresses N⊥ as a function of the bk 's.

The equilibrium state of the system is approximated in the canonical ensemble by the Bogoliubov thermal density operator

ρBog (T ) = 1 Z Bog e - k =0 ǫ k b † k bk /kB T (51)
where the normalization factor Z Bog is the Bogoliubov approximation for the partition function, and where we have introduced the Bogoliubov spectrum

ǫ k = h2 k 2 2m h2 k 2 2m + 2ρg 0 1/2 . ( 52
)
Bogoliubov theory: In the Bogoliubov approximation for the time evolution, the bk merely accumulate a phase, at the frequency ω k = ǫ k /h, similarly to the classical field case. From Wick's theorem one then obtains

1 2 {δ N0 (t), δ N0 (0)} Bog = k =0 nk (n k + 1)(U 2 k + V 2 k ) 2 +2U 2 k V 2 k cos(2ω k t) [n 2 k + (n k + 1) 2 ] ( 53
)
where

nk (T ) = 1 exp(ǫ k /k B T ) -1 (54)
is the mean occupation number of the Bogoliubov mode k. Note that we have considered here the so-called symmetric correlation function (as {X, Y } stands for the anticommutator XY +Y X of two operators) which is a real quantity, equal to the real part of the non-symmetrized correlation function. The time coarse grained version of the prediction (53) is obtained by averaging out the oscillating terms, which amounts to considering the correlation function of the temporally smoothed operator number of non-condensed particles

N non osc ⊥ ≡ k =0 U 2 k + V 2 k b † k bk + V 2 k . (55)
Quantum Ergodic theory: Discarding from the start the oscillating terms in N⊥ , as in (55), we face here the problem of calculating the long time limit of A(t)A(0) , where A is a linear function of the Bogoliubov mode occupation numbers,

A = k =0 γ k b † k bk . ( 56
)
As the quantum state of the system is given by the Bogoliubov approximation Eq.( 51), we may inject a closure relation in the Bogoliubov Fock eigenbasis:

A(t)A(0) = 1 Z Bog {n k } e -β k =0 ǫ k n k   k =0 γ k n k   × {n k }|A(t)|{n k } , (57)
where the sum is taken over all possible integer values of the occupation numbers, not to be confused with the mean occupation numbers (54). The non-explicit piece of this expression is the matrix element of A(t), which may be reinterpreted as follows:

{n k }|A(t)|{n k } = Tr [Aσ(t)] ( 58
)
where the density operator σ, initially a pure state in the Bogoliubov Fock basis,

σ(0) = |{n k } {n k }| (59)
evolves during t with the full Hamiltonian H. We know that this evolution involves Beliaev-Landau processes that will spread σ over the various Fock states |{n ′ k } . This evolution is complex. But we need here the long time limit only, in which we may assume that an equilibrium statistical description is possible. Since the system is isolated during its evolution, we take for σ(t → +∞) the equilibrium density operator in the microcanonical ensemble [37], and we calculate the expectation value of A with σ(t → +∞) as we did for the classical field model. The calculation can be done in the thermodynamic limit. As shown in the Appendix B, one can calculate to leading order in this limit the difference between canonical and microcanonical averages.

Here the microcanonical ensemble has an energy

E = E Bog 0 + k =0 ǫ k n k , where E Bog 0
is the ground state Bogoliubov energy. We introduce the effective temperature T eff such that the mean energy in the canonical ensemble at temperature T eff is equal to E,

0 = H Bog (T eff ) -E = k =0 ǫ k [n k (T eff ) -n k ] (60)
where . . . stands for an average in the canonical ensemble and H Bog is the Bogoliubov Hamiltonian. Using the results of Appendix B one gets

Ā(E) -A (T eff ) = - 1 2 k B T 2 eff A ′ H Bog ′ ′ (T eff ) (61)
where Ā(E) is the microcanonical average of A at energy E and where the apex ′ stands for derivation with respect to temperature. We further use the fact that, in the thermodynamic limit, for typical values of the occupation numbers n k , T eff weakly deviates from the physical temperature T . We calculate T eff by expanding (60) up to second order in T eff -T [38]. Evaluating (61) with this value of T eff , keeping terms up to the relevant order [38], gives the desired result

{n k }|A(t → +∞)|{n k } = A +   k =0 ǫ k (n k -nk )   A ′ H Bog ′ + 1 2 A ′ H Bog ′ ′ ×      k =0 ǫ k (n k -nk ) 2 H Bog ′ -k B T 2      (62)
where all the canonical averages are now evaluated at the physical temperature T [39]. It remains to inject this expression into Eq.(57). The resulting average over n k leads to the long time value of the correlation function:

A(+∞)A(0) -A 2 = A ′ H Bog ′ 2 Var H Bog (63) = k =0 γ k ǫ k nk (n k + 1) 2 k =0 ǫ 2 k nk (n k + 1) , ( 64
)
where we used Wick theorem and the property

dn k /dT = ǫ k nk (n k + 1)/k B T 2 [40].
Using Schwartz inequality, one can show that this long time value of the correlation function is less than its zero time value γ 2 k nk (n k + 1). To be complete, we present an alternative derivation of our prediction (64) in the Appendix C, based on results obtained in [37]. We also note that the quantum ergodic calculation directly leads to a prediction of the long time limit for the correlation function of the Bogoliubov mode occupation numbers, see Eq.(A17).

Replacing in Eq.( 64) the coefficients γ k by their expression from Eq.(55), γ k = U 2 k + V 2 k , we obtain the long time value of the condensate atom number correlation function in the quantum ergodic theory. Note that, in the thermodynamic limit, this long time value scales as the volume V , whereas the t = 0 value scales as V 4/3 .

C. Correlation function of the time derivative of the condensate phase

As in the classical field case, we first look for an expression of the first order time derivative of the condensate phase operator θ in terms of the amplitudes of the field Λ on the Bogoliubov modes. Taking as a starting point in Heisenberg picture

ih d dt θ = [ θ, H], (65)
we split the quantum field in a condensate part and a non-condensed part,

ψ(r) = â0 √ V + ψ⊥ (r), ( 66
)
and we insert this splitting in the expression of H. Using the modulus-phase representation of â0 and the commutation relation Eq.( 46), we obtain, using â †

0 â0 = N -N⊥ , - h d dt θ = g 0 V N - 1 2 + r dV Λ † Λ + g 0 2V r dV   Λ N0 + 1/2 N0 ( N0 + 1) Λ + h.c.   + g 0 2 √ V r dV 1 N0 Λ † Λ2 + h.c. . (67)
The quantity ( N0 + 1/2)/ N0 ( N0 + 1) is actually

1 + O(1/ N 2 0
) so it can to a high accuracy be replaced by unity. Furthermore, as we did in the classical field model, we now keep the leading terms in Λ, under the assumption of a weak non-condensed fraction. We can also replace θ by φ under the temporal derivative, since θ(0) is time independent. We obtain [33]

-h d dt φ ≃ g 0 V N - 1 2 + r dV Λ † Λ + 1 2 Λ2 + 1 2 Λ †2 .
(68) Taking the expectation value of this expression over the thermal state in the Bogoliubov approximation leads to an expression coinciding with the value of the chemical potential predicted by Eq.( 103) of [START_REF] Mora | [END_REF], which includes in a systematic way the first correction to the pure condensate prediction ρg 0 [41]:

µ(T ) = g 0 V   N - 1 2 + k =0 (U k + V k ) 2 nk + V k (U k + V k )   .
(69) At this order of the expansion, this analytically shows that -h d θ/dt is the chemical potential of the system. We now turn to various predictions for the symmetrized

correlation function of d φ/dt, C S (τ ) = 1 2 d dt φ (τ ), d dt φ (0) - d dt φ 2 . ( 70
)
Bogoliubov theory: At a time short enough for the interactions between the Bogoliubov modes to remain negligible, one can apply Bogoliubov theory to get

C Bog S (τ ) = g 0 hV 2 k =0 (U k + V k ) 4 nk (n k + 1) + 1 2 cos(2ω k t) n2 k + (n k + 1) 2 . (71)
The temporal coarse grained version of this correlation function is obtained by averaging out the cosine terms, which amounts to considering a temporal derivative of φ freed from the oscillating terms bb and b †b † :

d dt φ nonosc = - g 0 hV   N - 1 2 + k =0 V k (U k + V k )   - g 0 hV k =0 (U k + V k ) 2 b † k bk . (72)
Quantum ergodic theory: We directly apply to the smoothed temporal derivative (72) the reasoning performed in the previous subsection. Up to an additive constant, Eq.( 72) is indeed of the form (56), with

γ k = -(g 0 /hV )(U k +V k ) 2 .
From (64) we therefore obtain the long time behavior of the phase derivative correlation function

C ergo S (+∞) = g 0 hV 2 k =0 (U k + V k ) 2 ǫ k nk (n k + 1) 2 k =0 ǫ 2 k nk (n k + 1)
.

(73) The long time limit of the variance of the phase difference is then [42] Var φ ∼ C ergo S (+∞)t 2 .

(74)

Although our conclusion of a ballistic behavior for the phase agrees qualitatively with [11], the explicit expression of the coefficient of t 2 differs from the one of [11] due the fact that we account for interactions among Bogoliubov modes such as the Beliaev-Landau processes leading to ergodicity in the system, while in [11] the many-body Hamiltonian is replaced by the Bogoliubov Hamiltonian in the last stage of the calculation. As can be seen from (73) using Schwartz inequality, ergodicity results in a reduction of phase fluctuations with respect to the Bogoliubov prediction.

In the thermodynamic limit, analytical expressions can be obtained for this ergodic prediction. In the low temperature limit k B T ≪ ρg,

C ergo S (+∞) ∼ 8π 4 15
a 2 ξ V k B T h 2 k B T ρg 3 , (75)
ϕ) 1/2 /t [ρg(a 2 ξ/V) 1/2
/h] FIG. 7: (Color online) In the regime T ≪ Tc for a weakly interacting Bose gas, quantum ergodic prediction (73) for the long time limit of (Var φ) 1/2 /t, in the thermodynamic limit.

When expressed in units of (a 2 ξ/V ) 1/2 ρg/h, (Var φ) 1/2 /t is a function of kBT /ρg only, that is readily calculated numerically (solid line) or that may be approximated by asymptotic equivalents (75,76) in the low temperature or high temperature limit (green dashed line). Note that the dimensionless quantity a 2 ξ/V may also be written as

ρa 3 /(N √ 4π).
where ξ is the healing length such that h2 /mξ 2 = ρg. This tends to zero at zero temperature [44]. In the high temperature limit k B T ≫ ρg,

C ergo S (+∞) ∼ 12ζ(3/2) 2 5ζ(5/2) a 2 λ V k B T h 2 (76)
where the thermal de Broglie wavelength obeys λ 2 = 2πh 2 /mk B T and where ζ is the Riemann Zeta function.

Here we have identified g 0 to g [43]. In Fig. 7 we give the quantum ergodic prediction for lim t→∞ (Var φ) 1/2 /t calculated numerically, which is a universal function of k B T /ρg when expressed in the right units.

VI. CONCLUSION

We have investigated theoretically the phase spreading of a finite temperature weakly interacting condensate. The gas is assumed to be prepared at thermal equilibrium in the canonical ensemble, and then to freely evolve as an isolated system. After average over many realizations of the system, we find in classical field simulations that the variance Var ϕ of the condensate phase change grows quadratically in time. This non-diffusive behavior is quantitatively explained by an ergodic theory for the Bogoliubov modes, the key point being that conservation of energy during the free evolution prevents some correlation functions of the field from vanishing at long times. We have extended the analytical treatment to the quantum field case and we have determined the coefficient of the t 2 term in the long time behavior of Var φ, see Eq.( 73). This analytical result holds at low temperature T ≪ T c and in the weakly interacting regime ρa 3 ≪ 1, for a large number of thermally populated Bogoliubov modes, and relies on the assumption that the (although weak) interaction among the Bogoliubov modes efficiently mixes them (quantum ergodic regime).

A physical insight in our result is obtained from the following rewriting Var φ(t) ∼ t 2 h2

∂ T µ ∂ T H 2 Var H (77)
where Var H is the variance of the energy of the gas, here in the Bogoliubov approximation and in the canonical ensemble, µ(T ) is the chemical potential of the system as given by Eq.( 69) and H (T ) is its mean energy in the Bogoliubov approximation. This formula may also be obtained from the following reasoning. For a given realization of the system, of energy E, the long time limit of the condensate phase change φ(t) can be shown to behave as

φ(t) ∼ -µ micro (E)t/h, ( 78
)
where µ micro is the chemical potential calculated in the microcanonical ensemble [45]. For a large system, canonical energy fluctuations around the mean energy H (T ) are weak in relative value so that one may expand µ micro (E) to first order in E -H . Taking the variance of φ (t) over the canonical fluctuations of E then leads to (77), since ∂ T µ/∂ T H ≃ ∂ E µ micro ( H ) for a large system. This reasoning shows that a necessary condition for the observation of an intrinsic diffusive spreading of the condensate phase change is a strong suppression of the energy fluctuations of the gas. To this end one may try to prepare the system in a clever way, starting with a pure condensate and giving to the system a well defined amount of energy, e.g. by a reproducible change of the trapping potential [46]. Alternatively one may try to follow a given experimental realization of the system, measuring the phase of the condensate in a non-destructive way and replacing ensemble average by time average.

We thank Fabrice Gerbier and Li Yun for useful comments on the manuscript, and Mikhail Kolobov for interest in the problem. We thank Francis Hulin-Hubard for giving us access to a multiprocessor machine. We are grateful to Anatoly Kuklov for pointing to us his work [11] and for useful discussions. One of us (E. W.) acknowledges financial support from QuFAR. Laboratoire Kastler Brossel is a research unit of University Paris 6 and Ecole normale supérieure, associated to CNRS. Our group is a member of IFRAF. Using the master equation approach developed in quantum optics [47,48], we calculate the temporal correlation function of the operator b † q bq giving the number of Bogoliubov excitations in the mode of wavevector q, in the thermodynamic limit and including the Beliaev-Landau coupling among the Bogoliubov modes. This is useful to motivate the Gaussian model introduced in section IV, and to estimate the time required for the correlation function A(t)A(0) -A 2 , where A is of the form Eq.( 56), to depart from its value predicted by the Bogoliubov theory. The idea of the master equation approach is to split the whole system in a small system S and a large reservoir R with a continuous energy spectrum. Treating the coupling W between S and R in the Born-Markov approximation one obtains a master equation for the density operator σ S of the small system. Here the small system is the considered Bogoliubov mode, with unperturbed Hamiltonian H S = ǫ q b † q bq , and the reservoir is the set of all other Bogoliubov modes, with unperturbed Hamiltonian H R = H Bog -H S . In the thermodynamic limit, the reservoir indeed has a continuous spectrum, whereas the small system has a discrete spectrum. The coupling W between S and R is obtained from the next order Bogoliubov expansion of the Hamiltonian, that is from the part of the Hamiltonian cubic in the field Λ,

H cub = g 0 ρ 1/2 r dV Λ † ( Λ + Λ † ) Λ. (A1)
Inserting the modal decomposition Eq.( 48) in H cub , we isolate the terms that are linear in b q [49]:

W = gρ 1/2 bq R † + b † q R (A2) R = V -1/2 k,k ′ =0,q δ -q,k+k ′ A k,k ′ b † k b † k ′ + δ q,k+k ′ B k,k ′ bk bk ′ + δ q,k ′ -k 2C k,k ′ b † k bk ′ (A3)
where the operator R acts on the reservoir only and the coefficients have the explicit expressions:

A k,k ′ = U q V k V k ′ + (U q + V q )(U k V k ′ + U k ′ V k ) + V q U k U k ′ B k,k ′ = U q U k U k ′ + (U q + V q )(V k U k ′ + U k V k ′ ) + V q V k V k ′ C k,k ′ = U q V k U k ′ + (U q + V q )(U k U k ′ + V k V k ′ ) + V q U k V k ′ .
As a consequence of momentum conservation for the whole system, the action of R (respectively R † ) changes the reservoir momentum by -hq (respectively hq).

Let us denote with a tilde the operators in the interaction picture with respect to the Bogoliubov Hamiltonian H S +H R . In the Born-Markov approximation [47, 48] the master equation, for the density operator of the small system in contact with the reservoir in an equilibrium state, reads [50]

d dt σS (t) = - +∞ 0 dτ h2 Tr R [ W (t), [ W (t -τ ), σS (t)σ eq R ]]
(A4) where Tr R denotes the trace over the modes of the reservoir and the equilibrium density operator of the reservoir σ eq R is supposed here to be the Bogoliubov thermal equilibrium at temperature T . We expand the double commutator; because of momentum conservation, the resulting terms that contain two factors bq or two factors b † q exactly vanish when one performs the corresponding traces over the reservoir. Coming back to Schrödinger's picture we finally obtain

d dt σ S = 1 ih [ǫ q b † q bq , σ S ] + Γ - q bq σ S b † q + Γ + q b † q σ S bq - 1 2 Γ - q b † q bq + Γ + q bq b † q , σ S , (A5)
where {, } is the anticommutator and the new mode frequency is ǫq = ǫ q +h∆ q . The effect of the reservoir on the small system is then characterized by a frequency shift ∆ q of the mode, whose explicit expression we shall not need here [51], and by two transition rates Γ + q and Γ - q given by the Fourier transform of reservoir correlation functions at the mode frequency:

Γ + q = g 2 0 ρ h2 +∞ -∞ dτ e -iǫqτ /h Tr R [ R † (τ )Rσ eq R ] (A6) Γ - q = g 2 0 ρ h2 +∞ -∞ dτ e iǫqτ /h Tr R [ R(τ )R † σ eq R ]. (A7)
Since the reservoir is here at thermal equilibrium, the two rates are not independent but Γ - q = e βǫq Γ + q . This results from the Bose law property 1 + nk = e βǫ k nk . The rates are then conveniently characterized by their difference Γ q ≡ Γ - q -Γ + q . One finds

Γ q = g 2 0 ρ (2π) 2 h d 3 k 4C 2 k,k ′ (n k -nk ′ )δ(ǫ q + ǫ k -ǫ k ′ ) +2B 2 k,k ′ (1 + nk + nk ′ )δ(ǫ k + ǫ k ′ -ǫ q ) (A8)
where k ′ stands for |k -q| in the integrand. From [24] one checks that Γ q is simply the standard Beliaev-Landau damping rate for the Bogoliubov mode q, the contribution in C 2 corresponding to the Landau mechanism and the one in B 2 to the Beliaev mechanism. We now proceed with the calculation of the temporal correlation function of two operators A S , B S of the small system, the whole system being at thermal equilibrium. The quantum regression theorem [53] states that

A S (t)B S = A S (t) ≡ Tr S A S σ eff S (t) (A9)
for t ≥ 0, where the effective density operator σ eff S is in general not hermitian nor of unit trace but evolves with the same master equation as σ S with the initial condition

σ eff S (0) = B S σ eq S (A10)
where σ eq S = e -βHS /Z S is the unit trace equilibrium solution of Eq.(A5). Using the invariance of the trace under a cyclic permutation we obtain

d dt A S = iǫ q h [ b † q bq , A S ] + Γ - q 2 [ b † q , A S ] bq + b † q [A S , bq ] + Γ + q 2 [ bq , A S ] b † q + bq [A S , b † q ] . (A11)
Specializing to A S = B † S = b † q or bq and A S = B S = nq ≡ b † q bq leads to linear first order differential equations for A S (t) that are readily solved: b † q (t)b q = nq e (iǫq-Γq/2)t (A12) b q (t)b † q = (n q + 1)e (-iǫq-Γq/2)t (A13) nq (t)n q -n2 q = nq (n q + 1)e -Γqt . (A14)

In the classical field limit, where nk + 1 is assimilated to nk , this justifies the Gaussian theory of section IV.

In both the classical and quantum cases, this shows that the occupation numbers decorrelate with the rates Γ q corresponding to the Beliaev-Landau processes. These rates have a non-zero value in the thermodynamic limit. The present calculation is readily extended to the inclusion of two Bogoliubov modes in the small system, of wavevectors q and q ′ = q. The coupling of the small system to the reservoir now takes the form

W 2 = gρ 1/2 bq R † q + bq ′ R † q ′ + h.c. , (A15)
where the operators R have the same structure as in the single mode case, except that the double sum over k, k ′ is restricted to values different from 0, q, q ′ . In the resulting master equation for the density operator of the two modes, the only issue is to see if there will be crossed terms between the two modes, involving e.g. the product of b † q with bq ′ . By calculating the trace over the reservoir of the corresponding product of operators R, e.g. Tr R [ Rq (τ ) R † q ′ σ eq R ], we find in general that all crossed terms vanish, because of momentum conservation [52]. The master equation therefore does not couple the two modes, and one obtains nq (t)n q ′ -nq nq ′ = 0, for q ′ = q, (A16) as is assumed in the Gaussian model for the classical field of section IV.

It is instructive to compare the long time limit of the predictions Eqs.(A14, A16) to the quantum ergodic prediction. Adapting the reasoning leading to Eq.(64), we obtain the quantum ergodic result nq (+∞)n q ′ -nq nq ′ = ǫ q ǫ q ′ nq (n q + 1)n q ′ (n q ′ + 1)

k =0 ǫ 2 k nk (n k + 1)
.

(A17) In the thermodynamic limit this tends to zero, as in the master equation approach.

APPENDIX B: DEVIATION OF MICROCANONICAL AND CANONICAL AVERAGES

We wish to calculate the thermal expectation value of an observable A in the microcanonical ensemble rather than in the canonical one. For convenience, we shall parametrize the problem by the temperature T of the canonical ensemble. Restricting to the thermodynamic limit, where k B T is much larger than the typical level spacing of the system, we calculate the first order deviation of the two ensembles.

We start with the usual integral representation of the canonical ensemble in terms of the microcanonical one:

A (T ) =
dE Ā(E)e S(E)/kB e -βE dE e S(E)/kB e -βE (B1)

where the density of states is written in terms of the exponential of the microcanonical entropy S(E), Ā(E) and A (T ) stand for the expectation value of A in the microcanonical ensemble of energy E and in the canonical ensemble of temperature T respectively, and β = 1/k B T .

In the thermodynamic limit we expect the integrand to be strongly peaked around the value E 0 (T ) such that

d dE S(E) k B -βE E=E0(T )
= S ′ [E 0 (T )]/k B -β = 0, (B2) where f ′ (x) stands for the derivative of a function f with respect to its argument x. We then expand u(E) ≡ S(E)/(k B ) -βE up to third order in E -E 0 and we approximate the integrand as e u(E) = e u(E0) e S ′′ (E0)(E-E0) 2 /2kB × × 1 + 1 6 (E -E 0 ) 3 S (3) (E 0 )/k B + . . . .(B3)

We also expand Ā(E) up to second order in E -E 0 .

Performing the resulting Gaussian integrals leads to

A (T ) -Ā[E 0 (T )] = k B 2|S ′′ (E 0 )| × × Ā′ (E 0 ) S (3) (E 0 ) |S ′′ (E 0 )| + Ā′′ (E 0 ) + . . . (B4)
This relation can be inverted to first order, to give the microcanonical average as a function of the canonical one; to this order, we can assume that Ā[E 0 (T )] = A (T ) in the right hand side of (B4). Furthermore, using the implicit equation (B2) one is able to express the derivatives with respect to E 0 in terms of derivatives with respect to T , e. We present in this section an alternative derivation of the ergodic result (63) for the correlation function of an hermitian operator A, here introduced in (56). The long time limit of the correlation function is rigorously defined in terms of the temporal average (C2) Here the p m = Z -1 exp(-βE m ) are the statistical weights defining the average in the canonical ensemble. Equation (C2), specialized for γ k = (g 0 /V )(U k +V k ) 2 , is equivalent to Eq.( 22) in [11] for the dephasing time, provided one replaces there H ′ by A. This makes the link between our approach and the one of [11].

The delicate point is now to relate the formal expression (C2) (involving the unknown exact eigenstates |m ) to an explicit expression treatable in the Bogoliubov approximation. If one directly approximates the exact eigenstates by eigenstates of the Bogoliubov Hamiltonian, |m ≃ |{n k } , as done in [11] (see Eq.(61) there), one obtains the Bogoliubov result

C Bog A (+∞) = γ 2 k nk (n k + 1), ( C3
)
which is a good approximation for the t = 0 value of the correlation function, but not for its long time limit. We argue that the exact eigenstates are in fact coherently spread over a large number of Bogoliubov eigenstates of very close energies, because of the Beliaev-Landau couplings among them. Following [37], we thus assume that

m|A|m ≃ Ā(E m ) (C4)
where Ā(E m ) is the microcanonical ensemble average at the energy E m , a thermodynamic quantity that is now treatable in the Bogoliubov approximation as we have already done in Eq.( 62) [56]. After average over the canon-ical distribution for the energy E m = E, we then obtain for the correlation function,

C A (+∞) ≃ Ā(E) -A 2 ≃ A ′ H Bog ′ 2
Var H Bog , (C5) where . . . stands for the canonical average at temperature T . We recover Eq.(63). ally on the lattice model considered here, the momentum along the direction α is conserved modulo 2πh/lα only by the interaction term. Since the energy cut-off in the simulations is of order of kBT , 2πh/lα is also of the order of the typical momentum of the populated modes of the field, so that we do not expect conservation of momentum. The simulations indeed confirm this expectation.

[31] This results e.g. from the identity: [32] Y. Castin, R. Dum, Phys. Rev. A 57, 3008 (1998).

FIG. 1 :

 1 FIG.1:(Color online) Real part of the condensate amplitude correlation function(10) normalized to its t = 0 value, and divided by the zero temperature evolution e iρgt/h , as a function of time: in the vertical axis label, ã0(t) stands for a0(t)e iρgt/h . (a) Short times behavior and (b) long times behavior. In solid line from an average over 500 solutions of Eq.(8), in dashed line (red) the Bogoliubov approximation (40). Here the temperature is kBT = 3077.3 h2 /mV 2/3 = 0.1711Tc, where Tc is the critical temperature kBTc = (2πh 2 /m)(ρ/ζ(3/2)) 2/3 of the ideal gas, the number of particles is N = 4 × 10 5 and the coupling constant is such that the zero-temperature chemical potential is ρg = gN/V = 700h 2 /mV 2/3 .

FIG. 2 :

 2 FIG. 2: (Color online) Standard deviation of the condensate phase change ϕ(t) (12) as a function of time for (a) T = 0.08245Tc (lower curve) and T = 0.1711Tc (upper curve), (b) T = 0.29467Tc, (c) T = 0.453Tc, (d) T = 0.6473Tc. Thick solid line (black): numerical solution from the classical field model Eq.(8) averaged over 500 realizations. Thin solid line (red): a linear fit. The parameters N and ρg have the same values as in Fig.1.

FIG. 3 :

 3 FIG. 3: (Color online) Correlation function of the condensate atom number (11). (a) Short times, (b) long times. The classical field results are obtained from an average over 500 solutions of Eq.(8); they are represented by symbols in (a) and a solid line in (b). The dashed lines (red) are the Bogoliubov approximations (21) (oscillating line) and (22) (horizontal line) in (a), and only (22) in (b). The dashed-dotted line (purple) is the Gaussian model. For clarity in (b) we washed out fast oscillations in the simulation result and in the Gaussian model, by averaging over consecutive points over a time width 0.45mV 2/3 /h. The horizontal dashed-dotteddotted line (blue) in (a) and (b) is the ergodic long time limit prediction, described in section IV. The parameters N , T and ρg have the same values as in Fig.1.

FIG. 4 :

 4 FIG. 4: (Color online) Correlation function of |b k | 2 (black solid line) and modulus squared of the correlation function of b k (green dashed line), as obtained from 2000 realisations of the classical field simulations, (a) for the mode k = (0, 0, 2π/L3) and (b) for the mode k = (2π/L1, 10π/L2, 4π/L3). We define the correlation function of a quantity X as Corr X = X * (t)X(0) -| X | 2 . The purple dashed-dotted line is an exponential function of t, given by Eq.(24). It reproduces well the simulation results, and we have checked that the agreement is good for all vectors k lying in an arbitrarily chosen plane in k-space. The parameters N , T and ρg have the same values as in Fig.1.

FIG. 5 :

 5 FIG.5: (Color online) (a) Long time limit of (Var ϕ) 1/2 /t, in units of h/mV 2/3 , as a function of the temperature over the ideal gas critical temperature. (b) Long time limit of the condensate atom number correlation function, as a function of T /Tc. In square symbols the numerical results of the classical field simulation (averaged over 500 realizations), in dashed line (red) the Bogoliubov prediction and in dashed-dotteddotted (blue) the ergodic prediction. The parameters N and ρg have the same values as in Fig.1.

FIG. 6 :

 6 FIG.6: (Color online) Correlation function Cnon osc of the quantity φnon osc defined in Eq.(37) calculated from the simulation (square symbols; the solid line is a guide to the eye), or using the Gaussian theory (purple dashed-dotted line going to zero at long times). In dashed line (red) the Bogoliubov prediction. In dashed-dotted-dotted line (blue) the long time prediction (42) of the ergodic theory. The parameters T , N and ρg have the same values as in Fig.1.

  g. S ′′ [E 0 (T )] = -1/[T 2 E ′ 0 (T )]. This leads to Ā[E 0 (T )]-A (T ) = -k B T A ′ (T )It is actually more convenient to parametrize the result in terms of the mean canonical energy H (T ) rather than in terms of E 0 (T ). Applying (B5) to A = H allows to calculate E 0 (T ) -H (T ) to first order. One then uses the first order expansionĀ[ H (T )] = Ā[E 0 (T )] + [ H (T ) -E 0 (T )] 0 (T )] + . . . (B6)In the first order term of this expression, we can replace Ā[E 0 (T )] by the canonical average A (T ), and we can identify E 0 (T ) with H (T ); we can do the same identification in the right hand side of (B5). We obtain[54] Ā[ H (T )] -A (T ) = -ALTERNATIVE DERIVATION OF THE LONG TIME LIMIT OF CORRELATION FUNCTIONS

m m|A|m 2 .

 2 τ )A(0) -A 2 . (C1) We then insert in (C1) a closure relation using the exact N -body eigenstates |m of the interacting system with eigenenergies E m . In the absence of degeneracies we obtain a single sum over m, C A (+∞) = m p m m|A|m 2m p

  (|τ -τ ′ |) = 2 t 0 (t -τ ) C(τ ) dτ.

HBog ′ (T ) .

As the n k are randomly distributed according to the canonical distribution of temperature T , this expression has a vanishing average, and a variance scaling as the inverse of the volume V . When V → ∞, A is of order V , the right hand side of (61) is of order V 0 , and T eff -T is of order V -1/2 , so that, to be consistent, one actually has to expand the left hand side of (60) up to second order in T eff -T to obtain the next order correction to T eff -T . In evaluating (61) we expand A (T eff ) up to second order in T eff -T , keeping only terms up to O(V 0 ); in the right hand side of (61) one can to this order replace T eff by T . One then obtains the microcanonical expectation value Ā up to the desired O(V 0 ) order, as given by Eq.( 62).

[39] As a consistency check one calculates the average of Ā, that is of Eq.( 62), over the canonical distribution for the occupation numbers n k : using dn k /dT = ǫ k nk (n k + 1)/kBT 2 , one indeed finds Ā (T ) = A (T ).

[40] Due to the fact that the third term of Eq.( 62) has a vanishing mean, its contribution to A(+∞)A(0) -A 2 is much smaller than the one of the second term in the thermodynamic limit and can be neglected. We then obtain Eq.(64).

[41] More precisely it exactly coincides if one replaces the value of the chemical potential µ by its lowest order value ρg0 in the Bogoliubov expectation values of products of Λ in the right hand side of Eq.(103) of [START_REF] Mora | [END_REF]. This is allowed at the order of accuracy of Eq.( 68).

[42] Taking the classical limit nk + 1

→ g, we recover the classical prediction Eq.( 42).

[43] From Eq.(42) in the thermodynamic limit, one finds that the ergodic value of the correlation function C for the classical field also behaves like Eq.( 76) at high temperatures, except for the numerical factor, which depends on the energy cut-off. For the particular choice of energy cut-off made in this paper, one finds a numerical factor (2π

12.920 . . . which, accidentally, differs from the quantum one 12.209 . . . by about 5% only.

[44] Note that there is no phase spreading of the condensate at zero temperature, since the correlation function â † 0 (t)â0(0) does not tend to zero at large times but oscillates at a frequency given by the chemical potential, see S.T. Beliaev, Zh. Eksp. Teor. Fiz. 34, 417 (1958) [Sov. Phys. JETP 34, 289 (1958)].

[45] On one hand, the microcanonical chemical potential is the derivative of the energy with respect to N for a fixed entropy, that is for a fixed number of quan-tum states within the energy shell defining the microcanical ensemble. In the Bogoliubov frame, we obtain, from the Hellmann-Feynman theorem, µmicro

k is the Bogoliubov ground state energy. We then calculate explicitly all the derivatives with respect to N . E.g., ∂N ǫ k = (g0/V )(U k +V k ) 2 . On the other hand, in the spirit of the ergodic method we calculate the expectation value of (d φ/dt)non osc in the microcanonical ensemble of energy E. The expressions of the chemical potential obtained in these two ways coincide. Finally we split (d φ/dt)non osc as the sum of its mean microcanonical value φ and of fluctuations. Assuming that the fluctuations have vanishing correlations over long times for the considered given realisation, we obtain Eq.(78) after temporal integration.

[46] J. Kinast, A. Turlapov, J.E. Thomas, Q. Chen, J. Stajic term in bq, and has a finite number of quadratic terms in bq which are thus negligible in the thermodynamic limit.

[50] The Born-Markov approximation implies that the temporal correlation functions of R and R † with R and R † have a temporal width much smaller than the damping time of σR. This is the case in the Bogoliubov limit, as defined in [32], where g0 → 0 for fixed values of ρg0, T , q, dV , since the correlation functions are then fixed and the damping rates Γ± tend to zero linearly in g0.

[51] The correct calculation of ∆q may involve the first order effect of the quartic perturbation to HBog, neglected here. See [25] for a complete calculation.

[52] The only exception is when q + q ′ = 0; the crossed terms in this case, however, involve in the interaction picture a product of bq(t) and b-q(t), or a product of their hermitian conjugates; the corresponding functions of t oscillate with the frequency ±2ǫq/h; in the Bogoliubov limit, as in [50], this frequency is much larger than the relaxation rates of σS(t), so that the crossed terms can be neglected in the so-called secular approximation.