| J. Phys. B: At. Mol. Opt. Phys. 42 No 8 (28 April 2009) 081001 (5pp) |
| doi:10.1088/0953-4075/42/8/081001 |
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
. The latter expression
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
and eigenenergies
k. For the determination of N-particle mean values
· · ·
N = tr
we start with the canonical density operator
in second quantization with the corresponding energy
, the canonical partition function ZN and the projector
onto the N-particle subspace. The number states are
, where nk is the occupation number of the kth single-particle eigenstate |
k
. We also make use of the usual bosonic field operator
.
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),
with the product measure
. Coherent states |{z}
= |z0
|z1
· · · |zk
· ·· with
are used for all modes and
(see [16, 22]). For the projector on the N-particle subspace we find with [23]
With (2) and (3) we may express (normally ordered) quantum correlation functions as phase-space integrals, for instance
with the weight functions
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
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)
x|
k
which is the eigenvalue of the field operator when applied to the then stochastically evolving coherent state |{z(t)}
, i.e.
. With the notation ψ(x, t) =
x|ψ(t)
the desired SFE takes the form of a nonlinear, norm preserving stochastic Schrödinger equation, here using Stratonovich calculus [24]
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
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) =
x| dξ(t)
which represents white noise with correlations
dξ(t)|x![]()
x '| dξ(t)
= δ(x – x ') dt. Note that the white noise is always acted upon by the operator
from (6). For those field components with low energies
H
kT the operator K acts simply as the multiplication with the thermal energy kT. However, for high energy components with
H
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
—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
; the stochastic simulation requires a distribution of values for
. 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
is sufficient. As preliminary results suggest [25], this ceases to be true for the condensate fluctuations, for which the full
-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
. The liberty to choose other values for E0 (and other norms
, 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
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
. 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
. 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. 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-
and second-order
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. 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
and
(with μ the chemical potential and
the number operator), equation (5) reduces to
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
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
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
, where g = 4πa
2/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.
ReferencesSigmund Heller and Walter T Strunz 2009 J. Phys. B: At. Mol. Opt. Phys. 42 081001
T Moiseev et al 2009 J. Phys. D: Appl. Phys. 42 072003
Syud A Ahmed et al 2009 Environ. Res. Lett. 4 034004
Anh-Thu Le et al 2008 J. Phys. B: At. Mol. Opt. Phys. 41 081002
G I Fann et al 2005 J. Phys.: Conf. Ser. 16 461
P Bolognesi et al 2008 J. Phys. B: At. Mol. Opt. Phys. 41 051003
Ken B. Henisey et al. 2009 ApJ 706 705
S.-B. Qian et al 2009 ApJ 706 L96
J. I. González Hernández et al. 2009 ApJ 706 866
J A Souza and R F Jardim 2009 J. Phys. D: Appl. Phys. 42 032006