Quick search Find article
Quick search
Find article
J. Phys. B: At. Mol. Opt. Phys. 42 No 8 (28 April 2009) 081001 (5pp)
doi:10.1088/0953-4075/42/8/081001

FAST TRACK COMMUNICATION

Stochastic field equation for the canonical ensemble of an ideal Bose gas

Sigmund Heller and Walter T Strunz

Institut für Theoretische Physik, Technische Universität Dresden, D-01062 Dresden, Germany

Received 6 February 2009, in final form 11 March 2009
Published 30 March 2009

Abstract. We present a novel norm preserving stochastic evolution equation for a Bose field at finite temperature. Ensemble averages are quantum expectation values in the canonical ensemble. This numerically very stable equation suppresses high-energy fluctuations exponentially, which prevent cutoff problems from occurring. We present 3D simulations for an ideal gas in various trapping potentials and focus on ground-state occupation numbers and spatial correlation functions for a wide range of temperatures above and below the critical temperature. Although rigorously valid for non-interacting bosons only, we argue that weakly interacting Bose gases may also be amenable to this approach, in the usual mean-field approximation.

Ultracold quantum gases in traps are currently being investigated to a hitherto unknown precision and in a variety of circumstances [1–3]. The fascinating experimental possibilities to manipulate relevant parameters such as trap geometry, temperature, particle number and even interaction strength show that these gases are ideal quantum systems to investigate and revisit fundamental concepts of quantum many-body and statistical physics.

In these experiments, after cooling, traps contain a gas with a certain finite number of particles; no particle reservoir is present. Thus, the choice of a canonical ensemble is to be preferred over a grand-canonical description. Furthermore, for these finite systems the two ensembles give different predictions; as an example, the fluctuations of thermodynamic quantities may behave very differently (which may persist even in the thermodynamic limit [4]). For a thorough analysis of condensate fluctuations of a finite trapped ideal gas in the canonical and the microcanonical ensemble, see [5–9].

We present in this work a stochastic evolution equation for a c-number field ψ(x, t) such that quantum-statistical expectation values in the canonical state can be evaluated through an ensemble mean over these stochastic fields. To be specific, we will determine canonical ground-state occupation numbers and correlation functions for an N-particle system, for instance G_1(x,x^{\prime})=\langle\skew2\hat{\psi}^{\dagger}(x)\skew2\hat{\psi}(x^{\prime})\rangle_N= \langle\!\langle \psi^*(x,t) \psi(x^{\prime},t)\rangle\!\rangle_{\rm {ensemble}} . The latter expression \langle\!\langle \cdots \rangle\!\rangle_{\rm {ensemble}} denotes the ensemble mean over the stochastic field ψ(x, t) which will be clarified in detail below. Such (spatial) correlation functions have recently been measured in impressive experiments [10–14].

Various stochastic field approaches to ultracold Bose gases at finite temperature can be found in the literature; their relations to more traditional approaches in many-body quantum theory were reviewed very recently in [15]. Such stochastic methods have a long history in quantum optics [16] and they are intimately connected to phase-space representations of quantum states and operators. Not only do they provide the equilibrium state asymptotically, but very often they provide important insight when used to describe the transition from a non-equilibrium initial state to final equilibrium in a phenomenological manner. In the field of ultracold Bose gases, the use of stochastic field equations (SFE) has so far concentrated on the grand-canonical ensemble [15, 17–21] (taking into account two-body interactions). We will comment on these equations and their relation to our result towards the end of this work.

We stress two crucial properties of our novel equation that overcome problems of earlier approaches: first, in the high-energy region, our noise is exponentially suppressed which prevent ultraviolet problems from occurring. Second, the equation is norm preserving, which reflects the canonical ensemble. Both these properties ensure a very stable numerical solution such that full 3D canonical ensembles can be tackled. Moreover, our equation may be implemented in position space which allows us to describe gases in arbitrary trapping potentials without any difficulty. Eventually, interactions may be taken into account.

For the main part of this work, however, we consider an ideal gas of N particles in a trap with single-particle Hamiltonian H=\frac{p^2}{2m}+V(x)= \sum\nolimits_k \epsilon_k |\epsilon_k \rangle\langle\epsilon_k | and eigenenergies epsilonk. For the determination of N-particle mean values langle · · · rangleN =  tr (\cdots \skew4\hat{\rho}_N ) we start with the canonical density operator

Equation (1)

in second quantization with the corresponding energy \skew3\hat{H}= \sum\nolimits_k \epsilon_k\hat{a}^{\dagger}_k\hat{a}_k , the canonical partition function ZN and the projector \hat{P}_N=\sum\nolimits_{\sum n_k=N}|\{n_k\}\rangle\langle\{ n_k \}| onto the N-particle subspace. The number states are |\{n_k\}\rangle= \prod\nolimits_k\big(\hat{a}^{\dagger}_k\big)^{n_k}/\sqrt{n_k!}|0\rangle , where nk is the occupation number of the kth single-particle eigenstate |epsilonkrangle. We also make use of the usual bosonic field operator \skew4\hat\psi(x) = \sum_k \langle x|\epsilon_k\rangle \hat{a}_k .

In order to express quantum expectation values in terms of phase-space integrals, we start with the Glauber–Sudarshan P-representation [16, 22] of the Gaussian exponential in (1),

Equation (2)

with the product measure {\rm d}\mu\{z\}=\prod_k\big(\frac{{\rm d}\{{\rm Re} z_k\}\,{\rm d}\{{\rm Im} z_k\}}{\pi}\big) . Coherent states |{z}rangle  =  |z0rangle|z1rangle · · · |zkrangle · ·· with |z_k\rangle = \exp\big\{{-}|z_k|^2/2+z_k {\hat a}_k^\dagger\big\}|0\rangle are used for all modes and P(\{z\})=\prod\nolimits_k({\rm e}^{\epsilon_k/kT}-1)\cdot \exp\big({-}\sum\nolimits_k|z_k|^2({\rm e}^{\epsilon_k/kT}-1) \big) (see [16, 22]). For the projector on the N-particle subspace we find with [23]

Equation (3)

With (2) and (3) we may express (normally ordered) quantum correlation functions as phase-space integrals, for instance

Equation (4)

with the weight functions W_N(\{z\})=\frac{1}{N!} \big(\sum_k|z_k|^2\big)^{N} {\rm e}^{-\sum_k|z_k|^2} P(\{z\}) and the normalization constant CN  =  ∫dμ{z}WN({z}). Note that second (or higher) order correlations require the use of WN–2 (or lower index) in expression (4), while the CN remains.

Phase space integrals of the form (4) may be evaluated efficiently with the help of Monte Carlo methods based on Langevin equations for complex numbers zk(t). The stochastic equation is constructed in such a way that the weight functions Wn({z}) (n  =  N, N – 1,  ...) turn out to be stationary solutions of the corresponding Fokker–Planck equation and \big\langle \hat{a}_j^{\dagger}\hat{a}_k\big\rangle_N = \langle\!\langle z^*_j(t) z_k(t)\rangle\!\rangle_{\rm {ensemble}} holds (for large t). We succeed in deriving a novel norm preserving stochastic equation that satisfies this condition. The result is most conveniently written in terms of the c-number field ψ(x, t)  =  ∑k zk(t)langlex|epsilonkrangle which is the eigenvalue of the field operator when applied to the then stochastically evolving coherent state |{z(t)}rangle, i.e. \skew4\hat\psi(x)|\{z(t)\}\rangle = \psi(x,t)|\{z(t)\}\rangle . With the notation ψ(x, t)  =  langlex|ψ(t)rangle the desired SFE takes the form of a nonlinear, norm preserving stochastic Schrödinger equation, here using Stratonovich calculus [24]

Equation (5)

which is the central result of this paper. The given particle number N determines the norm of the field ψ (see below). The given temperature T enters through the important operator

Equation (6)

that obviously reflects the underlying Bose statistics.

Furthermore, in equation (5), Λ is a phenomenological damping parameter which sets the timescale for the transition to equilibrium. Using the relation to the SFEs derived for an interacting gas discussed at the end of this paper, it may be expressed in terms of microscopic parameters of the gas (e.g., scattering cross section). Reflecting a fluctuation–dissipation relation, the parameter Λ appears both in the damping term and as a square root, in the fluctuations. The latter are determined by the complex random field dξ(x, t)  =  langlex| dξ(t)rangle which represents white noise with correlations langledξ(t)|xranglelanglex '| dξ(t)rangle  =  δ(xx ') dt. Note that the white noise is always acted upon by the operator \sqrt{K} from (6). For those field components with low energies langleHrangle ll kT the operator K acts simply as the multiplication with the thermal energy kT. However, for high energy components with langleHrangle gg kT, the fluctuations are exponentially suppressed, which is a crucial feature of our novel equation to which we will come back towards the end of the paper.

The derivation of our SFE (5) is based on the proof that the weight functions Wn in equation (4) are stationary solutions of the corresponding Fokker–Planck equation. Now it is necessary to establish carefully the relation between the norm \langle\psi|\psi\rangle = \sum_k|z_k|^2 = {\cal N} —which remains constant for a single realization for all times—and the particle number N. As is apparent from equation (4), the phase-space integral has to be extended over fields that are not restricted to hypersurfaces with a single norm {\cal N} ; the stochastic simulation requires a distribution of values for {\cal N} . Due to the strongly localizing character of the Poisson distribution arising from the projector (equation (3)), however, it turns out that for particle numbers much larger than one and for all the temperatures and observables of interest in this work a simulation with a single norm {\cal N} is sufficient. As preliminary results suggest [25], this ceases to be true for the condensate fluctuations, for which the full {\cal N} -distribution appears to be required. Also, there is a surprising subtlety: the centre of the distribution depends on the absolute value of the ground-state energy E0. It is for E0  =  0 only that we have \langle {\cal N} \rangle = N . The liberty to choose other values for E0 (and other norms {\cal N} , accordingly) can be used with benefit to achieve faster convergence in the numerical implementation (see [25]).

Two further remarks are called for. First, the numerical results presented below are ample evidence that it is possible to replace the ensemble mean by a time average \langle\skew2\hat{\psi}^{\dagger}(x)\skew2\hat{\psi}(x^{\prime})\rangle_N= \lim\limits_{t\rightarrow\infty}\frac{1}{t}\int\nolimits_0^{t}{\rm d}s \{\psi^{\ast}(x,s)\psi(x^{\prime},s)\} over a single realization ψ(x, s). As the propagation is based on the highly nonlinear equation (5), this is not too surprising. Second, we chose to propagate with the term (Λ  +  i) in our equation (5) such that for Λ  =  0 the remaining complex unit `i' describes `real' dynamics. In this way, we can simulate the transition from a non-equilibrium to an equilibrium state in a phenomenological manner—see recent experiments [26].

We now turn to applications. While a gas in a box is best treated in momentum space, a general implementation of equation (5) in position space is advantageous, since it can be adjusted easily to any trapping potential (and in a next step mean-field atomic interactions may be included—see later). The use of position space requires the generation of spatially correlated noise \langle x|\sqrt{K}|\,{\rm d}\xi\rangle . For the simulations presented here we use a Wigner–Weyl representation [16, 22] of the operator K and consider only terms of lowest order in hslash. This is an excellent approximation as the examples presented here and further examinations [25] show.

The true quality of our equation is to be verified by calculating various characteristic quantities. The implementation in position representation allows us to treat the Bose gas in any trapping potential easily without the need to determine any eigenenergies or eigenfunctions. In order to make contact to previous results for the canonical ensemble, however, we restrict ourselves in the following to a 3D box and a 3D harmonic oscillator potential.

First we show the ground-state occupation in figure 1. The results for the harmonic oscillator (with N  =  200 particles, plus signs) are obtained by propagating equation (5) on a position grid. No use is made of the known spectrum and eigenfunctions. We compare with an analytical approximation (full line) [27], which is known to be in good agreement with exact results, and the thermodynamic limit (dashed-dotted line). Next we show results for a 3D box (N  =  100 particles, crosses) obtained from our SFE (here computed in momentum space) compared with a result based on a path integral approach [28] (dashed line) and find very good agreement. Temperature is scaled to the critical temperature of the thermodynamic limit [1, 2]. Note that finite size effects are very significant for the box as seen when comparing with the thermodynamic limit (dotted line).

Figure 1

Figure 1. Ground-state occupation as a function of temperature for an ideal 3D Bose gas of 200 particles in a harmonic isotropic trap (plus signs) and for 100 particles in a box potential (crosses). We compare with an analytical calculation for the harmonic (in a low temperature approximation) [27] (full line) and the box potential [28] (dashed line). We also display results for the thermodynamic limit (harmonic: dashed-dotted line, box: dotted line).

Next we determine spatial correlation functions of first- G_1(x,x^{\prime})=\langle\skew2\hat{\psi}^{\dagger}(x)\skew2\hat{\psi}(x^{\prime})\rangle_N and second-order G_2(x,x^{\prime})=\langle\skew2\hat{\psi}^{\dagger}(x)\skew2\hat{\psi}^{\dagger}(x^{\prime}) \skew2\hat{\psi}(x)\skew2\hat{\psi}(x^{\prime})\rangle_N above and below the critical temperature. G1 can be measured through interference experiments (see [10]). G2 and related quantities have been measured more recently in impressive experiments [11–14]. For first-order correlations, the canonical results are close to the grand-canonical values. For second-order correlations and temperatures below the critical temperature, however, large deviations appear and the results of the canonical ensemble may be obtained approximately with the help of condensate corrections to the grand-canonical calculation [29]. In figure 2, we show G1(x, 0) (top) and G2(x, 0) (bottom) obtained from the SFE (5); here for a 3D Bose gas of 1000 particles in a box with periodic boundary conditions (left-hand side) and for a gas of 200 particles in an isotropic harmonic oscillator (right-hand side). Our values (plus signs, crosses, stars) are compared with a direct grand-canonical calculation (`gc' in figure 2) and the corrected description (`gc + c' in figure 2) based on the theory of [29] and find good agreement.

Figure 2

Figure 2. First- and second-order spatial correlations G1(x, 0) (top) and G2(x, 0) (bottom) for a Bose gas in a 3D box potential with periodic boundary conditions (left) and a 3D harmonic oscillator potential (right) for different temperatures. Results from the SFE (plus signs, crosses, stars) are compared to a direct calculation based on the theory of [29]: a grand-canonical description (`gc', solid, dashed and dashed-dotted lines) fails for G2, while a corrected theory (`gc + c' dashed-dotted and dotted lines) shows good agreement.

We see the main achievement of this paper in the fact that we established a numerically very robust (exact) SFE for the canonical state of an ideal Bose gas, amenable to an efficient solution in full 3D for arbitrary trapping potentials. Still, in the light of the wealth of activities involving ultracold quantum gases, it is certainly of great importance to also investigate the interacting case. Let us therefore relate equation (5) to previous stochastic equations constructed for the grand-canonical ensemble of an interacting Bose gas [17–21]. A good survey over these different findings, their relations and their limitations is given in [20]. Most notably, in several of these approaches [17–20], termed `classical field methods' in [20] (and see also [30]), due to ultraviolet problems, lowly occupied states must be cut off [18] or treated in a different formalism [19, 20].

Let us now turn to our SFE (5): if one omits the nonlinear terms and substitutes both H\rightarrow {\tildeH} = H-\mu{\hat N} and K \rightarrow \tilde{K} = {\tildeH}/({\rm e}^{{\scriptsize\tildeH}/kT} -1) (with μ the chemical potential and {\hat N} the number operator), equation (5) reduces to

Equation (7)

One can easily show that this equation indeed provides proper grand-canonical ensemble averages [31]. We hasten to stress that the canonical equation (5) is not merely a normalized version of equation (7). As discussed before, in (7) the operator \tilde{K} incorporates a natural high energy noise cutoff induced by proper quantum statistics.

Equation (7) for a grand-canonical ensemble establishes the link to those earlier stochastic approaches for an interacting gas: there, the operator \tilde{K} appears as the temperature kT, requiring a cutoff. On the positive side, in these approaches, the interaction is taken into account by a mean-field contribution {\tildeH}\rightarrow {\tildeH}+g|\psi(x)|^2 , where g  =  4πahslash2/m and a is the s-wave scattering length.

After these considerations it appears more than tempting to use equation (5) for a canonical ensemble even in the case of a weakly interacting gas, with V(x) → V(x)  +  g|ψ(x)|2. As argued, the resulting equation is free from ultraviolet problems and coincides (in the grand-canonical case) with previous (`classical') equations (including interaction and a high-energy cutoff). Moreover, it reduces to the (imaginary-time) Gross–Pitaevskii equation for kT → 0, describing a pure condensate (with the given particle number N).

Let us summarize our result: we present a norm-preserving stochastic field equation for the canonical state of a Bose gas, describing experiments with a finite number of atoms in an arbitrary trap. The appearance of the operator K (equation (6)) ensures that cutoff issues do not appear; the equation is numerically very stable. We are able to determine important quantities like spatial correlation functions and occupation numbers as a function of temperature in arbitrary traps. Finally, we relate our equation to stochastic Gross–Pitaevskii equations that exist for interacting gases in the grand-canonical ensemble, raising expectations that the new equation should be applicable to weakly interacting Bose gases as well.

Acknowledgments

We are grateful for inspiring discussions with Markus Oberthaler and Thimo Grotz. SH acknowledges support by the International Max Planck Research School for Dynamical Processes in Atoms, Molecules and Solids, Dresden. Computational resources were provided by the ZIH, TU, Dresden.

References
[1] 
Pethick C J and Smith H 2002 Bose–Einstein Condensation In Dilute Gases (Cambridge: Cambridge University Press)  
[2] 
Pitaevskii L and Stringari S 2003 Bose–Einstein Condensation (Oxford: Oxford University Press)  
[3] 
Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885 
CrossRef
[4] 
Ziff R M, Uhlenbeck G E and Kac M 1977 Phys. Rep. 32 169 
CrossRef
[5] 
Grossmann S and Holthaus M 1996 Phys. Rev. E 54 3495 
CrossRef
[6] 
Politzer H 1996 Phys. Rev. A 54 5048 
CrossRef
[7] 
Wilkens M and Weiss C 1997 J. Mod. Opt. 44 1801 
CrossRef
[8] 
Grossmann S and Holthaus M 1997 Phys. Rev. Lett. 79 3557 
CrossRef
[9] 
Holthaus M, Kalinowski E and Kirsten K 1998 Ann. Phys., NY 270 198 
CrossRef
[10] 
Bloch I, Hänsch T W and Esslinger T 2000 Nature 403 166 
CrossRefPubMed
[11] 
Hellweg D, Cacciapuoti L, Kottke M, Schulte T, Sengstock K, Ertmer W and Arlt J J 2003 Phys. Rev. Lett. 91 010406 
CrossRefPubMed
[12] 
Fölling S, Gerbier F, Widera A, Mandel O, Gericke T and Bloch I 2005 Nature 434 481 
CrossRefPubMed
[13] 
Öttl A, Ritter S, Köhl M and Esslinger T 2005 Phys. Rev. Lett. 95 090404 
CrossRefPubMed
[14] 
Schellekens M, Hoppeler R, Perrin A, Gomes J V, Boiron D, Aspect A and Westbrook C I 2005 Science 310 648 
CrossRefPubMed
[15] 
Proukasis N P and Jackson B 2008 J. Phys. B: At. Mol. Opt. Phys. 41 203002 
IOPscience
[16] 
Mandel L and Wolf E 1995 Optical Coherence and Quantum Optics (Cambridge: Cambridge University Press)  
[17] 
Stoof H T C 1997 Phys. Rev. Lett. 78 768 
CrossRef
Stoof H T C and Bijlsma M J 2001 J. Low Temp. Phys. 124 431 
CrossRef
Duine R A and Stoof H T C 2001 Phys. Rev. A 65 13603 
CrossRef
[18] 
Davis M J, Morgan S A and Burnett K 2001 Phys. Rev. Lett. 87 0160402 
CrossRefPubMed
Davis M J, Ballagh R J and Burnett K 2001 J. Phys. B: At. Mol. Opt. Phys. 34 4487 
IOPscience
[19] 
Gardiner C W, Anglin J R and Fudge T I A 2002 J. Phys. B: At. Mol. Opt. Phys. 35 1555 
IOPscience
[20] 
Bradley A S, Blakies P B and Gardiner C W 2005 J. Phys. B: At. Mol. Opt. Phys. 38 4259 
IOPscience
[21] 
Drummond P D, Deuar P and Kheruntsyan K V 2004 Phys. Rev. Lett. 92 040405 
CrossRefPubMed
Deuar P and Drummond P D 2006 J. Phys. A: Math. Gen. 39 1163 
IOPscience
Deuar P and Drummond P D 2006 J. Phys. A: Math. Gen. 39 2723 
IOPscience
Deuar P and Drummond P D 2007 Phys. Rev. Lett. 98 120402 
CrossRefPubMed
[22] 
Schleich W P 2001 Quantum Optics in Phase Space (Berlin: Wiley)  
CrossRef
[23] 
Mollow B R 1968 Phys. Rev. 168 1896 
CrossRef
[24] 
Gardiner C W 1983 Handbook of Stochastic Methods (Berlin: Springer)  
[25] 
Heller S and Strunz W T 2009  in preparation
[26] 
Ritter S, Öttl A, Donner T, Bourdel T, Köhl M and Esslinger T 2007 Phys. Rev. Lett. 98 090402 
CrossRefPubMed
[27] 
Kocharovsky V V, Scully M O, Zhu S-Y and Zubairy M S 2000 Phys. Rev. A 61 023609 
CrossRef
[28] 
Glaum K, Kleinert H and Pelster A 2007 Phys. Rev. A 76 063604 
CrossRef
[29] 
Naraschewski M and Glauber R J 1999 Phys. Rev. A 59 4595 
CrossRef
[30] 
Hohenberg P and Halperin B I 1977 Rev. Mod. Phys. 49 435 
CrossRef
[31] 
Heller S and Strunz W T 2008 Path Integrals—New Trends and Perspectives ed W Janke and A Pelster (Singapore: World Scientific)  




Please login to access our web services, or create an account if you don't yet have one.

You must have cookies enabled in your web browser to be able to login.

Username
Password

Forgotten your password? Get a new one here.