Projected Gross-Pitaevskii equation for harmonically confined Bose gases at finite temperature

We extend the projected Gross-Pitaevskii equation formalism of Davis et al. [Phys. Rev. Lett. 87, 160402 (2001)] to the experimentally relevant case of thermal Bose gases in harmonic potentials and outline a robust and accurate numerical scheme that can efficiently simulate this system. We apply this method to investigate the equilibrium properties of the harmonically trapped three-dimensional projected Gross-Pitaevskii equation at finite temperature and consider the dependence of condensate fraction, position, and momentum distributions and density fluctuations on temperature. We apply the scheme to simulate an evaporative cooling process in which the preferential removal of high-energy particles leads to the growth of a Bose-Einstein condensate. We show that a condensate fraction can be inferred during the dynamics even in this nonequilibrium situation.

Figure 1

Thermalized classical fields for two values of energy. Density slices taken in the $y=0$ plane for a system in an anisotropic trap with ${\lambda }_x=\sqrt{8}$ and ${\lambda }_y=1$. (a) Low-energy case with an average energy per particle of $E=10\hslash {\omega }_z$. (b) Higher-energy case with an average energy per particle of $E=24\hslash {\omega }_z$. For both simulations $C_{\mathrm{nl}}=2000$ and the ground-state energy is $E_g=8.54\hslash {\omega }_z$. The single-particle cutoff energy is $E_{\text{cut}}=31\hslash {\omega }_z$, below which 1739 single-particle states remain in the coherent region.

Figure 2

Equilibrium properties calculated using the PGPE. (a) The condensate fraction $\left(f_c\right)$ and (b) temperature, calculated as a function of energy for $C_{\mathrm{nl}}=1000$. Other simulation parameters are the same as in Fig. 1. The averaging was performed using $N_s=1400$ samples over a evolution period of $T=1200∕{\omega }_z$. The unit of temperature is $T_0=N_C\hslash {\omega }_z∕k_B$ [9], where $k_B$ is the Boltzmann constant.

Figure 3

Time-averaged column densities in momentum (a)–(c) and position (d)–(e) space of classical field simulations with various energies. Cases: (a) and (d) $E=24\hslash {\omega }_z$, (b) and (e) $E=18\hslash {\omega }_z$, and (c) and (f) $E=10\hslash {\omega }_z$. Other simulation parameters are the same as in Fig. 1. The averaging was performed using $N_s=1400$ samples over a evolution period of $T=1200∕{\omega }_z$.

Figure 4

Normalized coherence functions $g_2\left(r\right)$ and $g_3\left(r\right)$ (see text) as a function of radial position in the $x=0$ plane (averaged over angle about the symmetry axis) and for various condensate fractions. The simulation parameters are the same as in Fig. 1, and all results are for the case $C_{\mathrm{nl}}=2000$. For reference the coherent value $(g_n=1)$ and thermal value $(g_n=n\text{!})$ of the coherence functions are indicated as dotted and dashed lines in the plots, respectively.

Figure 5

(Color online) Evolution of an evaporatively cooled matter-wave field. (a)–(d) Momentum-space density in the $k_y=0$ plane, (e)–(h) corresponding position-space density in the $y=0$ plane. The + signs are used to indicate the zero coordinate in the plots for reference. The initial state for the simulation is shown in Fig. 1b, and the evaporation is applied until $t=64\pi ∕{\omega }_z$ by setting $\Psi \left(\mathbf{x}\right)=0$ for $\mid z\mid >9x_0$ during the simulation. Other parameters are as in Fig. 1.

Figure 6

Growth of condensate and loss of energy and normalization during evaporative cooling process. (a) The condensate fraction of the total remaining classical field. The open circles are calculated by diagonalizing the one-body density matrix estimated by time averaging over two trap periods, while the solid curve is calculated by fitting bimodal distributions to single shots as described in the text. (b) The classical field energy (dashed line) and normalization (solid line) as a function of time. The point of time at which evaporation is stopped is marked by a vertical dotted line. Simulation parameters are explained in Fig. 5.

Figure 7

Example of the bimodal fitting procedure for a single shot of a classical field with a condensate fraction of about 0.24. (a) The image on the left is the column density of the classical field, while the image on the right is the bimodal fit. (b) A slice through another column density of the same field. The open circles are the data points and the solid line is the fitted curve.