Quantitative simulation of a magneto-optical trap operating near the photon recoil limit

Abstract

We present a quantitative model for magneto-optical traps operating on narrow transitions, where the transition linewidth and the recoil shift are comparable. We combine a quantum treatment of the light scattering process with a Monte-Carlo simulation of the atomic motion. By comparing our model to an experiment operating on the $5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^3{\mathrm{P}}_1$ transition in strontium, we show that it quantitatively reproduces the cloud size, position, temperature and dynamics over a wide range of operating conditions, without any adjustable parameters. We also present an extension of the model that quantitatively reproduces the transfer of atoms into a far off-resonance dipole trap, highlighting its use as a tool for optimizing complex cold atom experiments.

Keywords:

• Magneto-optical trap
• Monte Carlo

1. Introduction

The advent of laser cooling and trapping [1–4] was a revolutionary advance leading to a plethora of ultra-cold atomic experiments. The magneto-optical trap [5] (MOT) is the workhorse for all experiments with cold neutral atoms and has been extended to molecules in recent years [6–9]. Whilst quantitative theories of laser cooling in the Doppler and sub-Doppler regimes have existed for many years [10], a quantitative model for the MOT is more challenging. The difficulty arises due to the complex three-dimensional polarized light field in the presence of a magnetic quadrupole field, and the effects of optical pumping. In MOTs operating on strong transitions, the re-scattering and absorption of light must also be taken into consideration [11]. Nonetheless, as the technique is extended to more complex systems such as molecules, there is a considerable interest in models with the ability to quantitatively predict MOT properties.

One approach typically used to simulate MOT dynamics is to make simplifying assumptions about the atomic system and perform a Monte-Carlo [12] integration of the classical equations of motion. This assumes that the atoms experience an average force from the laser beams. Wohlleben [13] and Chaudhuri [14] have used this method to accurately simulate the atomic trajectories of rubidium atoms in a 2D${}^{+}$ MOT. This method has also been used to simulate loading into optical traps [15,16]. For more complex systems where optical pumping must be included, such as molecular MOTs, this model breaks down. A more accurate model is produced using the optical Bloch equations [17]. By performing an adiabatic elimination of the density matrix coherences, one is left with a series of rate equations. Atutov [18] has shown this model to be accurate at modelling a sodium MOT involving optical pumping, whilst both Comparat [19] and Tarbutt [20] have utilized this method to study the formation of molecular MOTs.

Divalent atoms exhibit inter-combination lines that are spin-forbidden by the usual electric dipole selection rules, but which are weakly allowed through state mixing, leading to very narrow transitions. These narrow transitions enable the production of ‘narrow-line MOTs’ (nMOTs) where the dynamics are limited by photon recoil, leading to sub-$\mathrm{\mu }$K temperatures [21]. The ability to trap atoms at low temperatures has aided the field of precision spectroscopy as the Doppler broadening due to the motion of atoms is greatly reduced, leading to the development of atomic clocks operating in the optical domain [22]. The narrow transition also facilitates the creation of high-density atom clouds [23,24], since the radiation trapping that limits the density of conventional MOTs is suppressed. This has led to a variety of studies in the high-density regime, such as the study of quantum degenerate gases [23,25–28], multiple scattering [29–31] and Rydberg blockade [32] to name but a few.

As the transition linewidth in a nMOT is so small, the transition is often power-broadened, which prevents the conventional adiabatic elimination of the density matrix coherences used to obtain a rate equation model [19,20]. To address this, we develop a Monte-Carlo model which is based upon the steady-state solution of the optical Bloch equations.

The paper is structured as follows. We initially introduce the concept of a nMOT in Section 2 and subsequently detail the model in Section 3. The experimental configuration used to test the model is detailed in Section 4. In Section 5, we compare the model to experimental data. In Section 6, we explore future applications of this model by simulating the loading of atoms into a far off-resonance dipole trap (FORT). We conclude our findings in Section 7.

2. Narrow-line MOTs (nMOTs)

The experimental configuration for a nMOT is the same as that for a conventional MOT [5], with atoms of mass m cooled and confined by a combination of a quadrupole magnetic field and laser beams (wavelength $\mathit{\lambda }$) with the appropriate circular polarization. What makes nMOTs distinctive is the ratio $\mathit{\eta }=\mathrm{\Gamma }/{\mathit{\omega }}_{\mathrm{R}}$, where $\mathrm{\Gamma }$ is the natural linewidth of the cooling transition, and ${\mathit{\omega }}_{\mathrm{R}}=(4{\mathit{\pi }}^2ħ)/(2m{\mathit{\lambda }}^2)$ is the frequency shift due to the atomic recoil following the absorption or emission of a photon. In conventional magneto-optical traps operating on strong dipole allowed transitions $\mathit{\eta }\gg 1000$. In this regime a single scattering event does not significantly alter the subsequent probability to scatter a photon, and the effects of individual scattering events can be averaged out, leading to the conventional semi-classical theory of Doppler cooling [10].

Conversely, in a nMOT $\mathit{\eta }\approx 1$. This condition typically only occurs when cooling on narrow dipole-forbidden transitions. For example, the ${}^{88}\mathrm{Sr}$$5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^3{\mathrm{P}}_1$ transition that we consider in this paper has $\mathit{\eta }=1.6$. In this limit individual scattering events significantly alter the subsequent probability of absorption, and the ultimate limit of laser cooling is set by the recoil temperature rather than the Doppler temperature [21].

Loftus et al. [24] showed that the behaviour of atoms in a nMOT is governed by the scaled detuning $\mathit{\delta }=|\mathrm{\Delta }|/{\mathrm{\Gamma }}^{\prime }\left(S\right)$, where ${\mathrm{\Gamma }}^{\prime }\left(S\right)=\mathrm{\Gamma }\sqrt{1+S}$ is the power-broadened linewidth, and $\mathrm{\Delta }=\mathit{\omega }-{\mathit{\omega }}_0$ the laser detuning with ${\mathit{\omega }}_0$ and $\mathit{\omega }$ the angular frequencies of the atomic transition and the cooling laser, respectively. Here, the parameter $S=I/I_{\mathrm{S}}at$ is the intensity of the cooling light I normalized by the saturation intensity $I_{\mathrm{S}}at$. Three regimes can be identified.

Figure 1. Absorption images of the nMOT along with their associated force curves. (a) and (b) correspond to $S=250$ and $\mathrm{\Delta }=-2\mathit{\pi }\times 110\text{kHz}$, whilst (c) and (d) correspond to $S=20$ and $\mathrm{\Delta }=-2\mathit{\pi }\times 200\text{kHz}$. The force curve associated with $S=1$ is also shown in (d) for comparison. The dashed ellipse in (a) and (c) shows where $\mathrm{\Delta }{\mathit{\omega }}_z=\mathrm{\Delta }$ in the quadrupole field with a vertical gradient of $8\text{G/cm}$.

The regime that is closest to a conventional MOT occurs when $\mathit{\delta }\approx 1$ and $S>>1$. This is illustrated in Figure 1(a) and (b). Here, the power-broadened linewidth dominates and the cloud forms close to the quadrupole centre as atoms are resonant with all three pairs of laser beams. In this ‘Doppler’ regime (I) the power-broadened linewidth determines the temperature, and the atoms are comparatively hot.

If $\mathrm{\Delta }$ is increased such that $\mathit{\delta }\gg 1$ but $S>1$, then the trap no longer forms at the quadrupole centre, but is displaced to where the Zeeman shift $\mathrm{\Delta }{\mathit{\omega }}_z=\mathrm{\Delta }$. The resulting resonance condition forms an elliptical ‘shell’ around the quadrupole centre. Since gravity is strong compared to the light-induced forces, the atoms fall under gravity until the resonance condition is met, forming an elliptically shaped nMOT (shown in Figure 1(c)) where the atoms predominantly interact with the beam that directly opposes gravity. This is seen by the clearly separated force peaks displaced from the quadrupole zero in Figure 1(d). We refer to this as the ‘power-broadened regime’ (II).

Finally, the recoil dominated ‘quantum’ regime (III) occurs when $\mathit{\delta }\gg 1$, and $S\le 1$. As in the power-broadened regime, the MOT position is determined by $\mathrm{\Delta }$ and the magnetic field gradient. However, since a photon recoil is sufficient to tune an atom out of resonance with the nMOT beams, recoil effects dominate the scattering behaviour of the atoms. This regime enables the lowest temperatures, ultimately reaching half the photon recoil limit, which for ${}^{88}\mathrm{Sr}$ is $460\text{nK}$ [24].

3. Modelling the cloud

In this paper, we explicitly consider nMOTs operating on the lowest-lying intercombination lines in divalent atoms such as the alkaline earths and Yb. These $J=0\to J=1$ transitions are completely closed,¹ and there is no optical pumping or dark states, as shown in Figure 2(a). Despite this apparent simplicity, it is still very challenging to fully model the interaction of this four-level system with the spatially varying quadrupole field and laser polarization, since one must keep track of complex spatially varying phases between the laser beams that appear in the off-diagonal terms in the atomic density matrix. Therefore, we make a significant simplification and treat each Zeeman transition as an independent two-level system, as shown in Figure 2(b). This amounts to non-conservation of population in the limit of $S\gg 1$ and also neglects Raman-like transitions related to coherences between Zeeman sub-levels. We expect that this is a good approximation in regimes (II) and (III). In these regimes, the atoms fall under gravity until the resonance condition is met and predominantly interact with the laser beam which opposes gravity. The Zeeman splitting between the $m_j$ sublevels is much greater than the transition linewidth, effectively isolating the three Zeeman sublevels, of which the $m_j=-1$ state is most strongly driven due to the helicity of the laser beams. However, we expect the model to fail in the ‘Doppler’ regime (I), and we show that this is indeed the case.

Figure 2. (a) Energy level structure of the nMOT. $\mathrm{\Delta }{\mathit{\omega }}_z$ is the Zeeman splitting due to the quadrupole field. (b) Simplified energy level structure used in the simulation. (c) nMOT experimental schematic. The straight red and circular black arrows represent the laser beam propagation direction and helicity, respectively.

3.1. Mathematical formalism

We simulate the ${}^{88}\mathrm{Sr}$$5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^3{\mathrm{P}}_1$ nMOT using the conventions and definitions shown in Figure 2(c). The nMOT is constructed from three retro-reflected orthogonal laser beams in the laboratory co-ordinate system $\mathit{r}=\left(x,y,z\right)$ where the unit vector directions $\widehat{\mathit{x}}$, $\widehat{\mathit{y}}$ and $\widehat{\mathit{z}}$ are shown in Figure 2(c), and the origin of the coordinate system is taken to be the zero of the quadrupole magnetic field. Each circularly polarized laser beam i has angular frequency ${\mathit{\omega }}_i$ and wave-vector ${\mathit{k}}_{\mathit{i}}$ and the helicity of each beam is illustrated in Figure 2(c). The nMOT quadrupole field $\mathit{B}$ is defined as $\mathit{B}=\mathit{\gamma }\left(x\widehat{\mathit{x}}+y\widehat{\mathit{y}}-2z\widehat{\mathit{z}}\right)$ where $\mathit{\gamma }$ is the gradient of the magnetic field along the radial $\left(\widehat{\mathit{x}},\widehat{\mathit{y}}\right)$ directions. The magnetic field splits the ${}^3{\mathrm{P}}_1$ state into three Zeeman levels j with splittings $\mathrm{\Delta }{\mathit{\omega }}_z=g{\mathit{\mu }}_{\mathrm{B}}\left|\mathit{B}\right|/ħ$ where $g=3/2$ is the g-factor and ${\mathit{\mu }}_{\mathrm{B}}$ is the Bohr magneton.

To reproduce the macroscopic dynamics of the nMOT, the simulation is performed for $\sim 5000$ atoms. Initially, the atoms are uniformly, randomly placed into a user-defined ellipsoid in the laboratory frame with position $\mathit{r}$. The atoms are also assigned random velocity vectors $\mathit{v}$ taken from a 3D-Boltzmann distribution with a user-defined initial temperature. These initial conditions are chosen to be similar to the final nMOT to reduce the processing time.² Typically, an initial temperature of $T=1\mathrm{\mu }\text{K}$ is used.

The total simulation time $t_{\mathrm{tot}}$ is broken into time-steps of length $\mathit{\delta }t$. At each time-step, the probability that each atom scatters a photon from laser i via a transition j is given by:(1) $$\begin{array}{cc}\hfill P_{\mathit{ij}}& ={\mathrm{\Gamma }}_e{\mathit{\rho }}_{\mathit{ee}}^{\mathit{ij}}\mathit{\delta }t\hfill \end{array}$$(1) (2) $$\begin{array}{cc}& =\frac{{\mathrm{\Gamma }}_e}{2}\frac{W_{j}S\mathit{\delta }t}{1+W_{j}S+4{\left({\mathrm{\Delta }}_i-{\mathit{k}}_{\mathit{i}}·\mathit{v}-\mathrm{\Delta }{\mathit{\omega }}_z^j\right)}^2/{\mathrm{\Gamma }}_e^2},\hfill \end{array}$$(2)

where ${\mathit{\rho }}_{\mathit{ee}}^{\mathit{ij}}$ is the steady-state excited state population derived from standard two-level optical Bloch equations [17] and $\mathit{\delta }t=0.1/{\mathrm{\Gamma }}_e$ such that $P_{\mathit{ij}}\ll 1$.

The coupling strength $W_j$ is dependent on the local magnetic field and the polarization of the nMOT laser beam. Due to the spatial inhomogeneity of the magnetic field, $W_j$ must be calculated as a function of position for each laser beam. This is most easily performed by entering a local atomic frame for the calculation. This frame is defined such that the local z-axis ${\widehat{\mathit{z}}}^{\mathbf{\prime }}$ is directed along the local magnetic field vector. Firstly, the total electric field for each laser beam $\mathcal{E}$ is defined in the laboratory frame in the spherical basis ${\widehat{\mathit{ϵ}}}_{\mathit{q}}$ [33,34] as:(3) $$\begin{array}{c}\hfill \mathcal{E}=\sum _q{\mathcal{E}}^q{\widehat{\mathit{ϵ}}}_{\mathit{q}},\end{array}$$(3)

where $\mathcal{E}=\left({\mathcal{E}}^1,{\mathcal{E}}^0,{\mathcal{E}}^{-1}\right)$ and(4) $$\begin{array}{cc}\hfill {\mathcal{E}}^1& =-\frac{1}{\sqrt{2}}\left({\mathcal{E}}_x+i{\mathcal{E}}_y\right),\hfill \end{array}$$(4) (5) $$\begin{array}{cc}\hfill {\mathcal{E}}^0& ={\mathcal{E}}_z,\hfill \end{array}$$(5) (6) $$\begin{array}{cc}\hfill {\mathcal{E}}^{-1}& =\frac{1}{\sqrt{2}}\left({\mathcal{E}}_x-i{\mathcal{E}}_y\right),\hfill \end{array}$$(6)

where ${\mathcal{E}}_{x,y,z}$ is the electric field defined in cartesian coordinates in the laboratory frame. A rotation is then performed to enter the local co-ordinate system of each atom to determine which transitions can be driven along with the associated transition coupling strengths. The rotation matrix is given by:(7) $$\begin{array}{cc}\hfill M_q\left(\mathit{\theta }\right)& =UR\left(\mathit{\theta }\right)U^{†},\hfill \end{array}$$(7) (8) $$\begin{array}{cc}& =\frac{1}{2}\left(\begin{array}{ccc}1+cos\left(\mathit{\theta }\right)& -\sqrt{2}sin\left(\mathit{\theta }\right)& 1-cos\left(\mathit{\theta }\right)\\ \sqrt{2}sin\left(\mathit{\theta }\right)& 2cos\left(\mathit{\theta }\right)& -\sqrt{2}sin\left(\mathit{\theta }\right)\\ 1-cos\left(\mathit{\theta }\right)& \sqrt{2}sin\left(\mathit{\theta }\right)& 1+cos\left(\mathit{\theta }\right)\end{array}\right),\hfill \end{array}$$(8)

where U is the transformation from the cartesian basis to the spherical basis and $R\left(\mathit{\theta }\right)$ is the rotation matrix which maps ${\widehat{\mathit{k}}}_{\mathit{i}}$ onto $\mathit{B}$ by an angle $\mathit{\theta }$. This leads to a new polarization vector ${\mathcal{E}}^{\prime }=M_q\mathcal{E}$ where $W_j={\left|{\mathcal{E}}^{j\prime }\right|}^2$.

Once $P_{\mathit{ij}}$ is known, random sampling from an uniform distribution is used to determine whether a scattering event occurs or not. If no scattering event occurs, the atom follows its initial trajectory defined in the laboratory frame as(9) $$\begin{array}{cc}\hfill {\mathit{v}}^{\mathbf{\prime }}& =\mathit{v}+\mathit{g}\mathit{\delta }t\hfill \end{array}$$(9) (10) $$\begin{array}{cc}\hfill {\mathit{r}}^{\mathbf{\prime }}& =\mathit{r}+{\mathit{v}}^{\mathbf{\prime }}\mathit{\delta }t+\frac{1}{2}\mathit{g}\mathit{\delta }{t}^2,\hfill \end{array}$$(10)

where $\mathit{g}=g\left(0,0,-1\right)$ is the acceleration due to gravity and the prime notation represents the final atom position or velocity after a time step $\mathit{\delta }t$. If a scattering event does occur, the atom evolves as(11) $$\begin{array}{c}\hfill {\mathit{v}}^{\mathbf{\prime }}=\mathit{v}+\frac{ħ\left|{\mathit{k}}_{\mathit{i}}\right|}{M}\left({\widehat{\mathit{k}}}_{\mathit{i}}+{\widehat{\mathit{k}}}_{\mathit{s}}\right)+\mathit{g}\mathit{\delta }t,\end{array}$$(11)

where ${\mathit{k}}_{\mathit{i}}$ is the wavevector of the laser from which the atom initially absorbed a photon with ${\widehat{\mathit{k}}}_{\mathit{i}}$ its associated unit vector, and ${\widehat{\mathit{k}}}_{\mathit{s}}$ is a random unit vector representing the direction of spontaneous emission.³

During each time step, the atomic positions and velocities are recorded, yielding a complete trajectory of each atom. A simulated absorption image of the nMOT is constructed by histogramming the atomic positions in the $\widehat{\mathit{x}}-\widehat{\mathit{z}}$ plane and calculating the column density along $\widehat{\mathit{y}}$. This is then normalized such that comparisons between theory and experiment can be made. A vertical and horizontal temperature ($T_v$ and $T_h$) is associated with the motion in the x and z directions by fitting a Maxwell–Boltzmann distribution to the vertical and horizontal components of $\mathit{v}$. This allows us to obtain the spatial, thermal and temporal dynamics of the atom cloud.

Figure 3. Experimental (top row) and theoretical (bottom row) absorption image as a function of nMOT beam detuning with $S=9$. The dashed purple line shows the resonance condition where $\mathrm{\Delta }=\mathrm{\Delta }{\mathit{\omega }}_z$.

4. Experimental configuration

The experiment is described in detail elsewhere [35–38] and so it is briefly summarized here. Initially, atoms from a strontium oven were slowed using a Zeeman slower before a ‘blue-MOT’ was formed on the $5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^1{\mathrm{P}}_1$ transition. Atoms in the blue-MOT were cooled to a temperature of several mK. After initial cooling, the blue-MOT light was removed and cooling light at $689\text{nm}$, which addresses the $5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^3{\mathrm{P}}_1$ transition, was applied. This light was artificially broadened to match the Doppler-broadened profile of the atoms in the blue-MOT. After sufficient cooling, the broadening of the $689\text{nm}$ light was removed, leaving single-frequency light, and a cold, dense nMOT. The nMOT was imaged using absorption imaging on the $5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^1{\mathrm{P}}_1$ transition. Images are captured on a Pixelfly QE camera with a post-magnification pixel-size of $8\mathrm{\mu }\text{m}$. In order to measure the temperature, the cooling light and quadrupole field were switched off, and the atomic cloud was imaged after a variable time-of-flight. By fitting the variation of the cloud width with expansion time we determined the temperature in the $\widehat{\mathit{x}}$ and $\widehat{\mathit{z}}$ directions in the conventional way. The saturation parameter S is determined independently by measuring the power and beam profile of the nMOT beams after the vacuum chamber. The uncertainty in this measurement is estimated to be $\pm 3\%$.

We are able to achieve a range of nMOT sizes, densities and temperatures by varying the nMOT beam detuning, power and the initial loading rate. We are able to achieve nMOTs ranging in size from $1/{\mathrm{e}}^2$ radii of $30\mathrm{to}300\mathrm{\mu }\text{m}$, densities up to $1\times {10}^{12}{\text{cm}}^{-3}$ and 3D temperatures as low as $460\text{nK}$. The vertical magnetic field gradient is held at $8\text{G/cm}$ during the final stage of the nMOT.

5. Testing the model

As discussed in Section 2, the properties of a nMOT are significantly dependent on $\mathrm{\Delta }$ and S. This strong parameter dependence allows us to test the accuracy of the model in a wide variety of nMOT regimes. Firstly, we test the model operating in regime (II) where the width and position of the nMOT are strongly dependent on $\mathrm{\Delta }$. The top row of Figure 3 shows experimental absorption images of the nMOT at four different values of $\mathrm{\Delta }$. It is clear that the MOT ‘sags’ under gravity and forms at lower positions as $\mathrm{\Delta }$ is decreased. The lower row of Figure 3 shows the theoretical absorption images obtained from the simulation. We qualitatively observe excellent agreement in position and shape of the nMOT in the absence of fitting parameters.

In order to quantitatively compare the model to the experimental data, the mean vertical position $\overline{z}$ and full width at half maximum (FWHM) of the nMOT are extracted numerically (without fitting) from the experimental and theoretical data. The results are plotted as a function of $\mathrm{\Delta }$ in Figures 4 and 5, respectively. The residuals normalized to their estimated uncertainties $R_{\mathit{\nu }}$ [39], are shown below each figure. The vertical FWHM saturates as a function of $\mathrm{\Delta }$ as the width is determined by the temperature of the atoms. The horizontal FWHM on the other hand continually increases as the radius of the resonant ellipse is proportional to $\mathrm{\Delta }$. We observe excellent agreement between experiment and theory with no adjustable parameters. Although the nMOT position is largely determined by the resonance condition, the agreement that we observe indicates that our model also successfully predicts the offset that results from the interplay between the cooling and trapping forces and gravity.

Figure 4. Experimental (blue triangles) and theoretical (purple circles) vertical position of the nMOT as a function of nMOT beam detuning. The lower plots show the normalized residuals between theory and experiment, normalized to the error on the experimental data. The dashed black lines show $R_{\mathit{\nu }}=\pm 2$.

Figure 5. Experimental (blue triangles) and theoretical (purple circles) vertical (a) and horizontal (b) FWHM of the nMOT as a function of nMOT beam detuning. The lower plots show the normalized residuals between theory and experiment, normalized to the error on the experimental data. The dashed black lines show $R_{\mathit{\nu }}=\pm 2$.

A more stringent test of the model is provided by the temperature. Unlike the position, which is largely determined by the resonance condition, the nMOT temperature is strongly dependent on the intensity of the cooling beams. As S and hence the normalized detuning $\mathit{\delta }$ varies, the nMOT crosses between the different regimes identified in Section 2.

Figure 6. (a) The blue squares and purple circles represent the measured and simulated nMOT temperatures for $S=1.9$. The blue triangles and purple diamonds represent the measured and simulated nMOT temperatures for $S=60$. The dashed purple line is a guide to the eye. (b) and (c) are the theoretical absorption images for nMOT beam detunings of $-2\mathit{\pi }\times 220\text{kHz}$ and $-2\mathit{\pi }\times 40\text{kHz}$, respectively.

Figure 7. Experimental (blue triangles) and theoretical (purple circles) nMOT temperatures as a function of time after a decrease (a) or increase (b) in nMOT laser beam power. The two nMOT beam powers used here were $S=14$ and 31. The normalized residuals $R_{\mathit{\nu }}$, normalized to the error bars, are shown below each figure. The dashed lines show $R_{\mathit{\nu }}=\pm 2$.

The dependence of the nMOT temperature on $\mathrm{\Delta }$ for two different values of S is shown in Figure 6(a). Firstly, we consider a nMOT operating close to the quantum regime with $S=1.9$. The temperature is essentially independent of $\mathrm{\Delta }$, since as shown in Figure 4 the position of the nMOT just tracks the resonance condition, and the number of scattering events each atoms experiences remains largely unchanged. In this regime, our model is again in excellent agreement with the measurements with no adjustable parameters.

Figure 8. Theoretical nMOT dynamics of the simulated crossed FORT with a trap depth of $5\mathrm{\mu }\text{K}$. The top and bottom row shows a theoretical absorption image in the $\widehat{\mathit{x}}-\widehat{\mathit{z}}$ and $\widehat{\mathit{x}}-\widehat{\mathit{y}}$ plane, respectively. All images in each row have the same colour-scale in order to show particle dynamics.

At higher intensity ($S=60$) the nMOT operates in the power-broadened regime (II). As expected the cloud is hotter, and the temperature is also observed to be largely independent of detuning again. The model is in excellent agreement for $|\mathrm{\Delta }|>2\mathit{\pi }\times 140\text{kHz}$, but begins to deviate significantly from experiment close to resonance. Here, the power-broadened linewidth begins to approach the Zeeman splitting in the excited state. Thus, the nMOT crosses over into the conventional ‘Doppler’ regime (I) where the linewidth is dominant, forming near the quadrupole zero, as shown in Figure 6(c). As a result, our key assumption that the atoms scatter independently on each of the three Zeeman transitions no longer holds, and the model breaks down.

As well as the equilibrium properties, we have also considered whether our model can reproduce the out-of-equilibrium dynamics of the nMOT. To do this, we looked at the response of the temperature to a sudden increase or decrease in the power of the laser beams. Initially, the nMOT was allowed to reach equilibrium at $S=S_0$. At $t=0$, S is suddenly decreased (increased) to a new value $S^{\prime }$. Experimental measurements of the subsequent cooling (heating) are shown in Figure 7, along with the results of the simulation. The agreement is excellent in both cases. More quantitatively, the reduced chi-squared statistics [39] were ${\mathit{\chi }}_{\mathit{\nu }}^2=0.7$ and 1.8, respectively, illustrating that our technique quantitatively reproduces both the steady-state and dynamic properties of the nMOT.

In summary, in the regimes of interest for experiment, where the nMOT is cold and dense, the results in Figures 3–7 show that our approach yields highly accurate quantitative predictions for the position, size, temperature and dynamics of the nMOT, requiring knowledge only of the experimental control parameters (intensity, detuning and magnetic field gradient). The model breaks down gradually at high intensity and close to resonance, exhibiting significant deviations only when ${\mathrm{\Gamma }}^{\prime }(S)\approx \mathrm{\Delta }$, as expected.

6. Dipole trapping

Motivated by the quantitative agreement between theory and experiment reported in Section 5, we have extended our model to investigate the loading of atoms into a far off-resonance dipole trap (FORT). Optimizing the transfer of atoms into such conservative traps is useful in applications such as optical lattice clocks and Bose–Einstein condensation. Typically, the optimum parameters are found using an experimental exploration of a large parameter space. The ability to quantitatively model the transfer process would, therefore, be a useful tool.

Figure 9. Experimental (blue triangles) and theoretical (purple circles) temperatures of the atoms in the crossed FORT as a function of trap depth. The dashed line is a linear fit to the simulated atom temperature. Experimental data taken from [40].

To compare the model to experimental data, we simulated the experiment performed by Ido et al. [40]. Their FORT consisted of two crossed laser beams with $1/{\mathrm{e}}^2$ radius of $28\mathrm{\mu }\text{m}$, operating at $800\text{nm}$. At this wavelength, the differential AC Stark shift between the ground and excited states is negligible, facilitating the simultaneous use of trapping and Doppler cooling. The FORT was loaded by first forming a single-frequency nMOT operating at $\mathrm{\Delta }=-200\text{kHz}$ with a total beam intensity of $S=18$. Whilst the nMOT was running, the FORT beams were then applied for a total time of $35\text{ms}$ before the nMOT beams were switched off. The temperature of the atoms in the FORT was subsequently measured using time-of-flight expansion.

We included the effect of the FORT beams in our model by neglecting the small differential AC Stark shift of the cooling transition, and considering only the conservative optical dipole force experienced by atoms in the ground state. Thus, an extra acceleration is included in the Newtonian dynamics part of the model, given by:(12) $$\begin{array}{c}\hfill a_{\mathrm{DT}}(x,y,z)=\frac{1}{M}\mathrm{\nabla }U\left(x,y,z\right),\end{array}$$(12)

where(13) $$\begin{array}{cc}\hfill \mathbf{U}\left(x,y,z\right)& =U_0({\mathrm{e}}^{-2\left[{\left(x-x_0\right)}^2+{\left(z-z_0\right)}^2\right]/w^2}\hfill \\ \hfill & +{\mathrm{e}}^{-2\left[{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2\right]/w^2}).\hfill \end{array}$$(13) $U_0$ is the trap depth, $x_0$, $y_0$ and $z_0$ are linear offsets in the $\widehat{\mathit{x}}$, $\widehat{\mathit{y}}$ and $\widehat{\mathit{z}}$ directions, respectively and w is the $1/e^2$ radius of the FORT beams. Figure 8 shows the simulated effect of applying the crossed FORT beams to the nMOT as a function of time. The atoms clearly move into the high-intensity region where the FORT beams intersect. We also observe a small number of atoms leaking into each individual FORT beam which is in qualitative agreement with experimental observations. Ido et al. typically capture $\approx 20\%$ of the atoms from the nMOT into the FORT. However, the model predicts a value of approximately double this. We attribute this difference to the lack of collisional losses in the model.

To make a quantitative comparison with experiment, we simulate the temperature of atoms trapped in the FORT as a function of $U_0$. The experimental measurements shown in Figure 9 exhibit a linear dependence, with the temperature varying in the range $0.1-0.2U_0$. Also shown is the simulated temperature of the atoms captured in the FORT. The error bars on the simulated temperature result from the Maxwell–Boltzmann fit to the velocity distribution.

We clearly observe excellent agreement between theory and experiment, once again in the absence of any adjustable fitting parameters. Along with the images in Figure 8 these results show that it is the interplay between the optical dipole force and laser cooling that sets the temperature, rather than truncation of the velocity distribution or evaporative cooling.

Looking forward, these results illustrate that our model could be a useful tool for optimizing the loading parameters of FORTs and optical lattices, eliminating the time-consuming trial and error approach often used to explore the available parameter space. By adding in a differential AC Stark shift, the model could be easily extended to other trapping wavelengths.

7. Conclusions

In summary, we have constructed a semi-classical Monte-Carlo simulation in order to model the dynamics of a nMOT, in particular the ${}^{88}\mathrm{Sr}$$5{\mathrm{s}}^2{}^1{\mathrm{S}}_0\to 5\mathrm{s}5\mathrm{p}{}^3{\mathrm{P}}_1$ nMOT. We observed excellent quantitative agreement between theory and experiment without fitting parameters, replicating the spatial, thermal and temporal dynamics of the system. We have also shown that we can quantitatively produce accurate results which simulate the loading of a crossed FORT from an operational three-dimensional nMOT.

In future, we aim to implement atom–atom interactions into the model. It will, therefore, be possible to fully simulate loading into an optical trap in high-density regimes where atom–atom interactions become significant. This will further mitigate the trial-and-error approach of the best parameters for loading atoms into optical traps from nMOTs. We may also be able to simulate the dynamics of a system with large atom–atom interactions such as those displayed by Rydberg atoms, leading to a greater understanding of strongly interacting many-body systems.

Acknowledgements

We would like to thank I. Hughes and J. Keaveney for useful discussions. We would also like to thank T. Ido and H. Katori for supplying experimental data.

Additional information

Funding

This work was supported by EPSRC [grant number EP/J007021/], Royal Society [grant number RG140546] and European Union EU grants FP7-ICT-2013-612862-HAIRS, H2020-FETPROACT-2014-640378-RYSQ and H2020-MSCA-IF-2014-660028. The data presented in this paper are available for download at doi:10.15128/r24x51hj01g.

Notes

No potential conflict of interest was reported by the authors.

1 This is only true for the bosonic isotopes.

2 The simulation is written in Python and makes use of parallelization packages (psutil) in order to reduce the computation time. A typical running time to simulate a nMOT using a standard desktop computer with an Intel Core i5-4690 Processor with 5000 atoms up to a time of $15\text{ms}$ is $\approx$ 1 h.

3 The angular dependence of the spontaneous emission rate is neglected. Although the atoms scatter predominantly on a ${\mathit{\sigma }}^{-}$ transition (in regime II and II), which favours emission in the horizontal plane (by a factor of 2), the agreement we observe with the model suggests that this effect is not significant.