Direct simulation Monte Carlo method for cold-atom dynamics: Classical Boltzmann equation in the quantum collision regime

In this paper, we develop a direct simulation Monte Carlo method for simulating highly nonequilibrium dynamics of nondegenerate ultracold gases. We show that our method can simulate the high-energy collision of two thermal clouds in the regime observed in experiments [Thomas et al. Phys. Rev. Lett. 93, 173201 (2004)], which requires the inclusion of beyond $s$-wave scattering. We also consider the long-time dynamics of this system, demonstrating that this would be a practical experimental scenario for testing the Boltzmann equation and studying rethermalization.

Figure 1

Ultracold atom collider: (a) schematic of the precollision arrangement of two clouds at $\sim 200\mathrm{nK}$ approaching at a collision energy of $\sim 200\mu \mathrm{K}$; (b) schematic of a postcollision system. (c) and (d) experimental images of postscattering density for two collision energies spanning the $d$-wave shape resonance. (e) and (f) show the theoretical calculations matching the experimental results using the direct simulation Monte Carlo (DSMC) method developed in this paper.

Figure 2

A two-dimensional schematic of the cells used for a swarm of test particles. (a) The rectangular master cells are all of the same size and are chosen to ensure all particles lie within the boundaries of this grid. Cell boundaries are indicated by lines, and particles are indicated by dots. (b) An enlargement of two master cells showing their adaptive subdivisions into smaller LAC subcells. The number of subcells is determined by the number of particles within the master cell.

Figure 3

A two-dimensional schematic of how collisions are performed within the TASC scheme. A single cell (outer boundary line) and the distribution of test particles (black dots) are shown in (a) and (b) for two different random collisions. The finer grid of internal lines represents the boundaries of the TASC subcells. The first particle of the collision pair is selected at random from all the particles in the cell. In (a), the first particle occupies a TASC subcell that contains other particles, and the second participant in the collision is chosen at random from these other particles. In (b), the first particle (which occupies the central subcell) is the sole occupant of a TASC subcell. In this case, we check to see if there are any particles in layer 1, and if so, the collision partner is chosen at random from these other particles. If there were no particles in layer 1, we would then check layer 2, and so on.

Figure 4

An example of the sequence of steps in a DSMC evolution. (a) A simple DSMC scheme where the whole system evolves according to a single global time step $\mathit{\Delta }t$. (b) An example of a cell using LATS. In this example, the global time step ($\delta t$) is held constant, while the local time step ($\delta t_c$) is shown to vary. Collisionless evolution occurs at each global time step. A collision step is performed at the global step when, at least, $\delta t_c$ has passed since the last collision step. At global time $t_3$, we show a collision step, at which the local time counter ($t_c$) is updated and a new local time step ($\delta t_c^{\prime \prime }$) is established. Here, the local time step decreases, showing two further collision steps that follow shortly after the first.

Figure 5

The relative error of the total collision rate, Eq. (18), for the equilibrium distribution $f_{\mathrm{eq}}(\mathbf{p},\mathbf{r})$ against $\gamma$ with ${\mathcal{N}}_T={10}^7$ is shown for the cases without (solid line) and with cell adaption where $N_{\mathrm{th}}=2$ (dotted line), 150 (dashed-dotted line), and 500 (dashed line). The results shown here are averaged over 200 initial conditions, while the error bars give the standard deviation. Without adaption, the error increases with increasing $\gamma$, since $f_{\mathrm{eq}}(\mathbf{p},\mathbf{r})$ becomes more coarsely grained. However, with the inclusion of adaption, this behavior is combatted as the LAC subcells adapt accordingly. We obtain the initial conditions for the test particles from $f_{\mathrm{eq}}(\mathbf{p},\mathbf{r})$ using the Monte Carlo acceptance-rejection method. System parameters: The harmonic potential is chosen to be the same as that used in the ultracold collider experiment with ${\omega }_x={\omega }_y=2\pi \times 155\mathrm{Hz}$ and ${\omega }_z=2\pi \times 12\mathrm{Hz}$, ${\mathcal{N}}_P=2\times {10}^5$ and $T=600\mathrm{nK}$.

Figure 6

Schematic of the ultracold collider used in Sec. 4c. Two clouds initially separated by a distance of $2r_0$ collide at a relative momentum of $2p_0$. The number of atoms that have scattered out of the clouds, after they have passed through each other, is referred to as ${\mathcal{N}}_{sc}$.

Figure 7

The fraction of scattered atoms due to the collision of two equilibrium clouds as a function of $a_{sc}$. Equation (25) [solid green (gray) line] has poor agreement with the solution of Eq. (24) (solid black line) for ${\mathcal{N}}_{sc}/{\mathcal{N}}_P>0.05$, since approximation (a4) is no longer valid. Group (i) scattered atoms (dashed black line) and total scattered atoms [groups (i) and (ii)] (dashed-dotted black line) from the DSMC solution. The system parameters are given in Fig. 5, and for the DSMC simulation, $\gamma =0.2,{\mathcal{N}}_T={10}^7$, and $N_{\mathrm{th}}=2$. The standard deviation error is not shown as it is on the order of the linewidth.

Figure 8

The relative error of ${\mathcal{N}}_{sc}$ as calculated by the DSMC solution of Eq. (24) [see text] to that of our pseudospectral solution of Eq. (24). Here, we show the two cases ${\mathcal{N}}_T={10}^5$ (black) and ${10}^7$ [green (gray)] with $\gamma =0.2$, $N_{\mathrm{th}}=25$, and the system parameters given in Fig. 5. We use ${\eta }_{\mathrm{coll}}=0.01,{\eta }_{\max }=0.1$, and ${\eta }_{\mathrm{tr}}=0.01$. The results shown here are averaged over 200 simulations, while the error bars give the standard deviation.

Figure 9

(a) Numerically calculated $s$-wave (dotted line) and $d$-wave (dashed line) phase shifts of Ref. [32]. (b) $s$-wave (dotted line), $d$-wave (dashed line), and total (solid line) cross sections.

Figure 10

Column densities at time $\pi /2{\omega }_x$ after the clouds reach the center of the trap. $T_{\mathrm{coll}}=200\mu \mathrm{K}$ [(a) and (c)] is a regime of $s$- and $d$-wave interference, while $T_{\mathrm{coll}}=300\mu \mathrm{K}$ [(b) and (d)] is a $d$-wave regime. (c) and (d) only show the scattered atoms. The initial conditions for the clouds are chosen as in Sec. 4c, and the system parameters are given in Fig. 5, while the simulation parameters are $\gamma =0.2,{\mathcal{N}}_T={10}^5$, and $N_{\mathrm{th}}=2$. The results were averaged over 200 simulations. Note, we have compared these results to simulations with ${\mathcal{N}}_T={10}^7$ also averaged over 200 runs, and we find that the number of scattered particles and the angular scattering distributions agreed to within $1\%$.

Figure 11

Column densities illustrating the long-time dynamics of rethermalization. (a) At time $t{\omega }_z=7.04$ after the clouds have passed through each other twice. The colliding clouds are still visible (density peaks). When the colliding clouds are depleted, the system continues to evolve through collective oscillations that are illustrated by the images (b) and (c) at $t{\omega }_z=18.85$ and 19.60, respectively. The decay of these collective oscillations occurs on a slower time scale than the depletion of the colliding clouds, and the distribution does not take on the equilibrium distribution until much later times as seen in (d) at $t{\omega }_z=500.02$. The trap frequencies are ${\omega }_z=2\pi \times 50\mathrm{Hz}$ and ${\omega }_x={\omega }_y=2{\omega }_z$, and each of the initial clouds has ${\mathcal{N}}_P={10}^6$ and $T=600\mathrm{nK}$. We use an isotropic differential cross section with $a_{sc}=10\mathrm{nm}$. The initial separation is chosen such that there is insignificant overlap of the clouds. The momenta are chosen to give $T_{\mathrm{coll}}=32.4\mu \mathrm{K}$, giving a final equilibrium temperature of $T=6\mu \mathrm{K}$. Note, for an isotropic trap, the system does not completely thermalize without mean-field effects, since the breathing mode does not damp [12].

Figure 12

Envelope of the oscillations (as seen in the lower inset) of the root-mean square of $r$. The rapid decay of the envelope within the first ten trap cycles is attributed to the depletion of the colliding clouds, while the slower decay is the decay of the collective modes. The upper inset shows the mean number of collisions per atom.

Figure 13

Collective oscillation induced by a small contraction of the radial trap confinement for a system with ${\omega }_z=2\pi \times 50\mathrm{Hz},{\omega }_x={\omega }_y=10{\omega }_z,{\mathcal{N}}_P={10}^4$, and $T=600\mathrm{nK}$. The legend gives ${\mathcal{N}}_T$ used, and as this increases, the results for a single run become increasingly indistinguishable from the ensemble-averaged result for ${\mathcal{N}}_T={10}^7$.

Figure 14

Relative error of the numerical total collision rate in the case of the equilibrium distribution (15) for the choices of $M_c$, where the error bars indicate the standard deviations of 500 averages. Here, the original DSMC algorithm is implemented with ${\mathcal{N}}_T={10}^7$ and bin parameter (23) $\gamma =0.2$. The choice $M_c^a$ is seen to diverge for low ${\mathcal{N}}_T$ as $\overline{N_c}$ in Eq. (A3) becomes appreciable, which agrees well with the theoretically calculated error ($M_c^a$ theory) using Eq. (A3). Using $M_c^b$ removes this divergence, and the error is seen to agree well with the expected error for this discretization ($M_c^b$ theory). The final data set ($M_c^b$ no fix) demonstrates the error that arises when ${\tilde{M}}_c$, having noninteger values after rescaling, is not accounted for (see Sec. 3c3) [i.e., not including the ceiling function in Eq. (11)]. The system parameters are given in Fig. 5.