Simulation method for evaporative cooling of trapped Bose gases at finite temperatures

Experimentally, atomic Bose-Einstein condensates (BEC) are created by evaporatively cooling trapped Bose atoms. The evaporative cooling technique consists of the selective removal of high-energy atoms and collisional relaxation of the remaining gases. This is the finite-temperature dynamics in which the condensation and the thermal cloud coexist. The experiments at finite temperatures have observed many phenomena that can be attributed to the existence of the thermal cloud, such as vortex-lattice creation and Landau damping. In particular, an interesting experiment of vortex-lattice creation was reported in 2001 by the JILA group [1]. This experiment created vortex lattices by first spinning a normal Bose gas and then cooling the rotating gas through the BEC transition. This experiment indicated that vortex nucleation can arise from angular momentum transfer between the thermal cloud and the condensate. Ref.  [1] also showed that the above procedure can create vortex lattices more efficiently than creating BEC and then rotating it.

Motivated by the experiment of Ref. [1], in this paper, we study the dynamics of evaporative cooling of rotating Bose gases at finite temperatures. For this purpose, we first develop a simulation method for evaporative cooling. The dynamics of evaporative cooling has been theoretically studied by using the quantum kinetic theory within the so-called “ergodic approximation” [2, 3]. However, in order to clarify the dynamical effect of the thermal cloud, it is important to simulate the full dynamics of the thermal cloud. In this paper, we use the formalism of Zaremba, Nikuni, and Griffin (ZNG) to include the generalized GP equation and a semiclassical Boltzmann equation for the thermal cloud, which treats the excitations semiclassically within the Hartree Fock (HF) approximation. The generalized GP equation includes the mean field due to thermal could and source term associated with collisions between the condensate and the thermal cloud. We note that in Ref[4], we simulated the ZNG equations in the presence of a rotating trap potential, and showed that vortex creation occurs due to the transfer of angular momentum between the condensate and the thermal cloud. However, vortex creation by evaporative cooling of rotating gas has not been studied theoretically.

The paper is organized as follows. In Sec. 2, we summarize ZNG equations. In Sec. 3, we explain our approach to evaporative cooling. In Sec. 4, we discuss the evaporative cooling of a rotating Bose gas. Finally, we conclude in Sec. 5.

The ZNG formalism consists of a generalized GP equation for the condensate order parameter $\Phi({\bf r},t)$ and the semiclassical Boltzmann kinetic equation for the thermal cloud distribution function $f({\bf p},{\bf r},t)$. The derivation of the ZNG formalism can be found in Ref. [5] The coupled ZNG equations are given by

where $n_{\rm c}=|\Phi|^{2}$ is the condensate density, $\tilde{n}=\int\frac{d{\bf p}}{(2\pi\hbar)^{3}}f({\bf p,r},t)$ is the noncondensate density, and $V_{\rm ext}({\bf r},t)$ is a time-dependent trap potential. The $C_{22}$ term describes two-body collisions between noncondensate atoms, and the $C_{12}$ term describes collisions between condensate and noncondensate atoms. The detailed expressions for the $C_{22}$ and $C_{12}$ terms are given in Ref. [5]. In Eq. (2), $U=V_{\rm ext}+2gn_{\rm c}+2g\tilde{n}$ is the effective potential for the noncondensate, which includes the self-consistent Hartree–Fock mean field. In the ZNG model, noncondensate atoms have the the local energy $\tilde{\epsilon}_{p}({\bf r},t)=p^{2}/(2m)+U({\bf r},t)$, and the condensate has the local energy $\epsilon_{c}=\mu_{c}+\frac{1}{2}mv_{c}^{2}$, where ${\bf v}_{c}\equiv-\hbar(\Phi^{*}\nabla\Phi-\Phi\nabla\Phi^{*})/(2imn_{c})$ and $\mu_{c}\equiv\hbar^{2}\nabla^{2}\sqrt{n_{c}}/(2m\sqrt{n_{c}})+V_{\rm ext}+gn_{c}+2g\tilde{n}$ is the local condensate chemical potential. The source term $R$ appearing in Eq. (1) is directly related to the $C_{12}$ collision term through $R=(\hbar/2n_{c})\int d{\bf p}/(2\pi\hbar)^{3}C_{12}$.

The numerical procedure for solving ZNG equations was developed in Ref. [6]. In this approach, the dynamics of the thermal cloud is calculated by using $N$-body simulations. The dynamics of the condensate is determined by numerically propagating the GP equation using a split-operator fast Fourier transform method. We use this approach to simulate evaporative cooling.

Before explaining our approach to evaporative cooling, we briefly mention the earlier studies in Refs.[2, 3]. Theoretical descriptions of the evaporative cooling in Refs.[2, 3] are based on the ergodic approximation, in which the noncondensate distribution function is given as a function of the single-particle energy. This approximation allows one to derive a simplified kinetic equation for the distribution function. The evaporative cooling process is modeled by using a Bose distribution function truncated at a cutoff energy as an initial state. Albeit approximate, this approach provides a physically reasonable description of evaporative cooling.

In the present study, we propose a simulation method for evaporative cooling based on the Jackson-Zaremba approach [6]. The main difficulty of applying this approach to evaporative cooling is that one does not directly work with the distribution function. Instead, one solves the classical $N$-body problem for test particles and constructs the phase-space distribution function. In this case, it is not obvious how to incorporate the evaporation process into the simulation. In the present study, we model evaporative cooling by eliminating the test particles with energy above a cutoff energy. One can then effectively prepare a truncated distribution.

More explicitly, we simulate the time evolution of test particles following the procedure described in Ref. [6]. After updating the phase-space coordinates $({\bf r}_{i},{\bf p}_{i})$ of test particles at each timestep of calculating the collisional term, we eliminate the test particle that has the single-particle energy greater than the cutoff energy: $\tilde{\epsilon}_{p_{i}}({\bf r}_{i},t)>\epsilon^{{\rm max}}_{\rm cut}(t)$. In order to describe the evaporative cooling, cutoff energy must be lowered in time. We adopt the following time dependence of the cutoff energy: $\epsilon_{\rm cut}(t)=\epsilon^{0}_{\rm cut}\exp(-\alpha t^{2}).$ In order to determine the initial cutoff energy $\epsilon^{{\rm max}}_{0}$, we first find the maximum single-particle energy $\tilde{\epsilon}^{\rm{max}_{0}}_{p}({\bf r},0)=\max\tilde{\epsilon}_{p_{i}}({\bf r}_{i},t=0)$ in the initial state $t=0$. We then set the initial cutoff energy as $\epsilon_{\rm cut}^{0}=\gamma\tilde{\epsilon}^{{\rm max}_{0}}_{p}$. The optimal values for $\alpha$ and $\gamma$ have been chosen such that evaporative cooling works the most efficiently. We varied the value of $\alpha$ in increments of 0.01 and the value of $\gamma$ from $0.5$ to $0.9$ in increments of 0.05 and observed the cooling behavior. We found that the values $\alpha=0.01$ and $\gamma=0.85$ resulted in the fastest attainment of a temperature of approximately 0.5 times the BEC transition temperature ($T\sim 0.5T_{c}$). That is, if the maximum cutoff energy is too large the evaporation does not occur. On the other hand, if the maximum cutoff energy is too small, the system quickly loses the particles and one cannot retain the phase-space density. As for the time dependence of $\epsilon_{\rm cut}(t)$, we explored the other monotonically decreasing functions, such as linear, quadratic, and logarithmic functions. However, we found that the exponential function allows for the fastest attainment of $T\sim 0.5T_{c}$ without any numerical uncertainty.

In the numerical simulation shown below, we set the parameters of the system so as to make the ZNG equations in a rotating trap potential remarkably stable: a total of $N_{0}=6\times 10^{6}$ ⁸⁷Rb atoms, a radial trap frequency $\omega_{\rho}/2\pi=6.8$Hz and the axial tap frequency $\omega_{z}=13.6$. The simulation range is $10R_{\rm TF}$, where $R_{\rm TF}=a_{\rm ho}\left(8\pi Na/a_{\rm ho}\right)^{1/5}$with $a_{\rm ho}=\sqrt{\hbar/m\omega_{\rho}}$ is the Thomas-Fermi radius, and the spatial computational mesh for the condensate is $\sim 0.1\xi$ (where $\xi=\hbar/\sqrt{2ngn_{0}}$ is the healing length). Finally, we choose the total number of test particles as $N_{\rm test}=1.0\times 10^{8}$. We prepare an initial equilibrium state at $T\sim 0.95T_{\rm c}$ [4]. Figure 1 shows the simulation results for the time evolution of the condensate during the evaporative cooling. Figures 1(a)-(c) display the two-dimensional condensate densities sliced at $x=0$ from the complete three-dimensional data. The initial equilibrium density distribution at $t=0$ is shown in Fig. 1(a). The evaporate cooling is started at $t=0$. Fig. 1(b) and Fig. 1(c) show the growth of the condensate during evaporative cooling. In Fig. 1(d), we plot the time evolution of the condensate fraction $N_{\rm c}/N$. We can see that $N_{\rm c}/N$ increases as time evolves. From the condensate fraction, the final temperature at $t=t_{0}$ is estimated as $T\simeq 0.56T_{c}$

The pronounced oscillations in the behavior of the condensate fraction can be attributed to the fact that the ZNG simulation takes into account collisions between particles in the condensate and the noncondensed component. These oscillations arise from processes where the condensate transitions into the noncondensed component and vice versa. We confirmed that our procedure successfully simulates the growth of the condensate by evaporative cooling.

In this section, we discuss the evaporative cooling of a rotating Bose gas. In the JILA experiment, vortex lattices were created by first spinning a normal Bose gas and then cooling the rotating gas through the BEC transition. However, the ZNG formalism does not allow one to go through the BEC transition since this formalism is based on the mean-field description. Therefore, we start with the temperature below, but very close to the BEC transition temperature where the condensate fraction is very small.

The system parameters are the same as in the previous section, except that we now consider a weakly anisotropic potential rotating with time-dependent rotation frequency $\Omega(t)$:

where $x^{\prime}=x\cos[\Omega(t)t]+y\sin[\Omega(t)t],\ y^{\prime}=-x\sin[\Omega(t)t]+y\cos[\Omega(t)t].$ We first prepare the initial equilibrium state at $T=0.9T_{\rm c}$ in the static potential $V_{\rm ext}({\bf r},0)$. We then start rotating the anisotropic trap potential at $t=0$, gradually increasing the rotation frequency $\Omega(t)$. More explicitly, we change the rotation frequency according to

Throughout this paper, we set $\Omega_{\rm max}=0.8\omega_{\rm rad}$, $\epsilon=0.001$, $t_{\rm in}=t_{0}$, $t_{\rm out}=3t_{0}$, and $a=10/t_{0}^{2}$ with $t_{0}=10/\omega_{\rm rad}$. With the present choice of parameters, we effectively start with $\Omega(0)\approx 0$. We stop increasing the rotation frequency at $t=2t_{0}$. At $t=4t_{0}$, where the system reaches a stationary state, we stop rotating the trap potential and start evaporative cooling. Finally, we stop evaporative cooling at $t=8t_{0}$.

In Fig.2, we show the evolution of the condensate density during the evaporative cooling. The initial equilibrium density distribution at $t=0$ is shown in Fig. 2(a). Figure 2(b) shows the equilibrium density distribution in the presence of a rotating thermal cloud. We start the evaporate cooling at $t=4t_{0}$ and the time evolution of the condensate density profile during evaporate cooling is shown in Figs 2(c)-(f). We can see the growth of condensate as well as the excitation of surface mode. The vortices then begin to penetrate into the condensate. Figure 2(g) shows that the vortex lattice is formed as an equilibrium state. After stopping evaporative cooling at $t=8t_{0}$, the vortices decay and the condensate size is decreased as shown in Fig. 2(h).

We developed a simulation method for evaporative cooling of trapped BEC at finite temperatures using ZNG formalism. In our method, the evaporative cooling process is modeled by eliminating test particles describing the noncondensate atoms. The simulation results successfully reproduced growth of BEC by evaporative cooling. We applied our method to simulate vortex lattice formation by cooling the rotating Bose gas. Our simulation results showed nucleation and formation of vortex lattice during the evaporative cooling Since we used the ZNG formalism, we could not start with the normal gas above $T_{c}$ and go through the BEC transition. In order to improve the method so that it can go through the BEC transition, one will need a theory that can include the fluctuation of the condensate order parameter, such as the classical-field technique [7].