Bose-Einstein condensation of an ideal gas is a famous example of a phase transition which relies purely on quantum statistics. The recent experimental achievement of Bose-Einstein condensation in a dilute gas of trapped alkali atoms [1, 2, 3, 4] triggered a revival of interest in various properties of the condensate. Some old but fundamental problems of statistical description of the Bose-Einstein condensate have been extensively studied. In particular, the issue of microcanonical and canonical fluctuations of a noninteracting condensate was solved successfully [5, 6, 7, 8, 9, 10, 11]. On the contrary, a closed system description of fluctuations of the weakly interacting condensate seems to be still a controversial problem at least at temperatures close to the critical one [12, 13, 14, 15]. The standard mean field description based on the Bogolubov approach fails [16] at the critical temperature, as then the order parameter vanishes. On the other hand many excited states become highly occupied particularly in the case of realistic systems consisting of a finite number of atoms. The complexity of the problem grows further when one is interested in an estimation of the characteristic time scale of the condensate fluctuations[17, 18]. The calculation of the two-time correlation function requires obviously a dynamical approach.

In the following, we present an approach that allows the study of the time evolution of an interacting Bose gas. The problem is very difficult and existing approaches lead to numerical algorithms which are hard to implement. The first one [19, 20] invokes quantum kinetic theory and leads to a quantum Boltzmann master equation. The second one [21, 22] is based on the separation of the “classical” mean field describing the condensate from its fluctuations. The fluctuations are quantal in nature and the dynamics couples averages of both their normal and anomalous parts to the condensate mean field. None of the two mentioned methods seems to be capable of handling a realistic case of a large number of particles at a relatively high temperature.

The first problem one encounters in describing the interacting condensate is the real N-body nature of the system. The link between this problem and the traditional “giant matter wave” description is not straightforward. The only exception is the case of a gas trapped in a box with periodic boundary conditions. Here the symmetry of the potential uniquely enforces the form of the eigenstates of the one-particle density matrix regardless the interaction. Due to the symmetry, the condensate can be associated with the ground state of the trap. Therefore, to avoid all ambiguities related to the identification of the condensate we are going to study here the system in a box with periodic boundary conditions [23].

We want to present our approach rooted in the theory of multi-mode lasers. Let us start with a general formulation of the N-body problem. The second-quantized Hamiltonian for the atomic system confined in a box with periodic boundary conditions and interacting via pairwise contact potential may be written in the following form:

where Φ is an atomic field operator, V=L ³ is a volume of the system (L being a size of the box) and $g=\frac{4\pi \hslash \alpha }{mV}$ characterizes the atom-atom interactions in the low-energy, s-wave approximation (a being the scattering length and m – the mass of the atom). The field Φ is expanded in natural modes of the system – the plane waves:

where a _k are bosonic annihilation operators, and $\mathbf{k}=\frac{2\pi }{L}\mathbf{n}$ with n_i =0,±1,±2,… (i=x,y,z). With this substitution the Hamiltonian assumes its final form:

where $\xi =\frac{\hslash }{2m}{\left(\frac{2\pi }{L}\right)}^2$. After elimination of a fast time dependence with the substitution a _k=exp(-iξn²t)α _k, the Heisenberg equations of motion for the operators α _k acquire the following form:

Solving the nonlinear operator equations (4) is difficult. The complexity of the problem obviously requires some approximation. A semiclassical approximation is particularly well suited for the description of a finite system at temperatures below the critical one, except the region close to the absolute zero.

The semiclassical approximation consists in replacing all operators α _k by c-number complex amplitudes (we are not going to introduce a separate notation for corresponding complex fields). At very low temperatures only the lowest lying states are macroscopically occupied and quantum fluctuations in excited states become important (see [21] for comparison). Therefore, the relevant range of temperatures for the applicability of our model excludes temperatures close to zero. From the viewpoint of the Bogolubov method [24], such an approach is legitimate as indeed many modes are macroscopically populated, i.e. their occupation is greater than quantum fluctuations.

The semiclassical approximation leads to nonlinear differential equations which must be solved numerically. The first observation is that Eqs.(4) can be viewed as a set of Hamilton equations for the complex degrees of freedom. Our approximate dynamics preserves the number of particles: N=∑_k ${\alpha }_{\mathrm{k}}^{*}$α_k, as well as the total energy of the system. It therefore corresponds to a genuine microcanonical description. Moreover, the resulting equations resemble the famous Fermi-Pasta-Ulam [25] problem of a system of harmonic oscillators coupled by a nonlinear interaction. A one-dimensional version of this dynamics has been studied recently [26] in the context of the pure Bose-Einstein condensate (T=0) reloaded from a harmonic into a rectangular trap. Equations studied here, however, in spite of a formal analogy, describe quite a different physical situation. Our complex amplitudes are not expansion coefficients of the condensate wave function in some convenient basis. They represent a number of coupled “mean fields” – a natural extension of the condensate mean field of the Bogolubov approach. Let us notice that if we start with 100% occupation of the k=0 mode and keep only this mode in the model we simply recover the standard Gross-Pitaevskii equation for the interacting condensate.

 

Fig. 1. Condensate occupation as a function of time for total energy per particle E/ħ=539Hz.

Download Full Size | PDF

For the initial conditions we assume “Bose-Einstein-like” occupation of different trap states. While calculating the energy in order to determine this distribution, however, we neglect the energy of interparticle interactions. This is evidently not the equilibrium distribution for the interacting system. Moreover, a population of individual states does not specify initial conditions uniquely. It defines the modulus of the corresponding amplitude but says nothing about its phase. In our approach each mode is assigned an initial, randomly chosen, phase. Any subsequent dynamics depends on the initial phases and, in a sense, a single simulation describes a single experimental realization. Microcanonical expectation values would require, therefore, an average over these initial phases. As we have checked, instead of doing this, it is enough to start with some random phases and trace the system evolution for a sufficiently long time. The observed self-averaging can be attributed to the ergodicity of the studied dynamics. Contrary to the 1-D analogue which is completely integrable [27], the 3-D version of the dynamics may be chaotic. We make use of this fact in our simulations and avoid averaging over phases of the complex amplitudes.

Values of the parameters in the model are ξ=71.373 Hz, N=10⁵ and g=0.018 Hz (the atomic mass and the scattering length are those of ⁸⁷Rb and the size of the box is equal to the Thomas-Fermi radius of a condensate of N atoms in a trap with frequency of ω ₀=2π 80 Hz). We performed our calculations for the model with 729 modes (n_i =-4,…,4, i=x,y,z). Further increasing of the number of modes does not lead to a substantial change in the results for the case studied in this paper. Our calculations show that after a time of the order of a few miliseconds the system reaches a dynamical equilibrium. The mean occupation of the condensate (k=0 mode) stabilizes at some value and on larger time scales (of the order of a second) it only fluctuates around this mean value – see Figure 1.

The similarity to the Fermi-Pasta-Ulam problem may cast some doubts on the genuine ergodicity of the dynamics investigated in this paper. Originally, Fermi, Pasta, and Ulam intended to test the ergodic hypothesis in the chain of coupled harmonic oscillators. However, they observed a quasiperiodic behavior identified by the returns of energy to the initial (lowest) mode. This kind of behavior has been discovered also by J.H. Eberly and co-workers [28] in a resonantly coupled system composed of a two-level atom and a single mode of a monochromatic electromagnetic field being initially in a coherent state. Occurrence of the so called revivals in the system proves the quantum nature of the electromagnetic field. In our calculations the largest time scale for which we have studied the dynamics was of the order of one second. On this time scale we did not observe any revivals but this, obviously, does not exclude the possibility of revivals on much larger time scales. In fact our numerical simulations involve a finite number of modes and so the numerical implementation inevitably leads to a quasiperiodic evolution.

Since the evolution (see Eqs.(4)) is Hamiltonian, the most natural choice of independent variables are the energy and the particle number. The number of particles was fixed through all our calculations (N=10⁵) and the energy of the system was the control parameter. Traditionally, however, the temperature, not the energy, is used as an independent thermodynamic variable. Calculation of the microcanonical temperature requires monitoring of the entropy of the system for different energies. Although, in principle, this can be done [29], we think that the energy is a much better characteristic of the system. The reason is twofold. First of all, it is not obvious whether the ergodic hypothesis can be applied in the studied case and therefore whether the notion of the temperature can be unambiguously introduced. Secondly, in realistic experiments it is rather the final energy of thermal atoms which is detected in destructive time-off-light measurements. The temperature is a parameter fitted to the observed velocity distribution.

In Figure 2 we present probability distributions of the condensate population. The different colors signify different total energies. Let us stress two main features which can be seen in Figure 2: (i) positions of the maxima of these distributions correspond to a mean (time-averaged) occupation of the Bose-Einstein condensate, (ii) the spread of these distributions around the most probable value is a measure of condensate fluctuations. The figure clearly indicates that the condensate population decreases monotonically with energy while fluctuations become larger.

 

Fig. 2. Probability distribution of the condensate (k=0 mode) population. Different colors signify different total energies per particle: blue E/ħ=2036 Hz, green E/ħ=1680 Hz, orange E/ħ=1424 Hz, yellow E/ħ=1164 Hz, red E/ħ=539 Hz.

Download Full Size | PDF

The occupation and the fluctuations of the interacting condensate are shown in detail in Figure 3, where the mean occupation of the condensate is depicted (in blue). We see that the condensate disappears at the energy per particle close to E/ħ=2000 Hz.

 

Fig. 3. Condensate occupation (in blue) and fluctuations (in red) versus total energy per particle.

Download Full Size | PDF

Figure 3 also presents fluctuations of the condensate (in red). Both the mean occupation and the fluctuations of the condensate are smooth functions of the energy. They do not show any discontinuity signifying a phase transition because they correspond to a finite system of N=10⁵ atoms. Let us notice, however, that fluctuations reach the maximum value at the energy per particle close to E/ħ=1450Hz. Moreover, this value of energy corresponds to the inflection point of the mean occupation of the condensate. Both curves distinguish the same characteristic value of the energy. Close to this energy the system undergoes rapid changes. This characteristic energy corresponds to the critical energy for the Bose-Einstein condensation.

In conclusions: first of all, we have proposed a very efficient approach to the time evolution of an interacting Bose-Einstein condensate. This approach is valid for a wide range of temperatures except for very small ones. We have shown that the method works very well in a realistic and relevant range of parameters: energy, number of particles, and interaction strength. We calculated the mean occupation as well as the fluctuations of the interacting condensate by averaging the corresponding time-dependent quantities. We did not explore all possible applications of our method. Our aim was rather to show its potential. We believe that the semi-classical approach presented here is perfectly suited for studying many properties of an interacting condensate which so far, due to their complexity, were beyond direct theoretical investigations. In particular, it may shed new light at the unresolved problem of nucleation of a condensate and other dynamical phenomena. The method is easily extendible to more than one field which allows to study coupled (via photoassociation or Feshbach resonances) atomic and molecular systems (see [30]) or other mixtures of Bose gases at finite temperatures.