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

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\rm{s}^2~^1\rm{S}_0 \rightarrow 5\rm{s}5\rm{p}~^3\rm{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 (FORT), highlighting its use as a tool for optimising complex cold atom experiments.


Introduction
The advent of laser cooling and trapping [1][2][3][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][7][8][9]. While 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 threedimensional 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]. For more complex systems where optical pumping must be included, such as molecular MOTs, this model breaks CONTACT: Matthew P. A. Jones, Email: m.p.a.jones@durham.ac.uk down. A more accurate model is produced using the optical Bloch equations [16]. By performing an adiabatic elimination of the density matrix coherences, one is left with a series of rate equations. Atutov [17] has shown this model to be accurate at modelling a sodium MOT involving optical pumping whilst both Comparat [18] and Tarbutt [19] have utilised 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-µK temperatures [20]. 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 [21]. The narrow transition also facilitates the creation of high density atom clouds [22,23], since the radiation trapping that limits the density of conventional MOTs is suppressed. This has led to the a variety of studies in the high density regime, such as the study of quantum degenerate gases [22,[24][25][26][27], multiple scattering [28][29][30] and Rydberg blockade [31] to name but a few.
As the transition linewidth in an nMOT is so small, the transition is often power broadened, which prevents the conventional adiabatic elimination of the density matrix coherences used to a rate equation model [18,19]. 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 the concept of an 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.

Narrow-line MOTs (nMOTs)
The experimental configuration for an 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 λ) with the appropriate circular polarization. What makes nMOTs distinctive is the ratio η = Γ/ω R , where Γ is the natural linewidth of the cooling transition, and ω R = (4π 2 )/(2mλ 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 η 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 an nMOT η ≈ 1. This condition typically only occurs when cooling on narrow dipole-forbidden transitions. For example, the 88 Sr 5s 2 1 S 0 → 5s5p 3 P 1 transition that we consider in this paper has η = 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 [20].
Loftus et. al [23] showed that the behaviour of atoms in an nMOT is governed by the scaled detuning δ = |∆|/Γ (S), where Γ (S) = Γ √ 1 + S is the power-broadened linewidth, and ∆ = ω − ω 0 the laser detuning with ω 0 and ω the angular frequencies of the atomic transition and the cooling laser respectively. Here the parameter S = I/I Sat is the intensity of the cooling light I normalised by the saturation intensity I Sat . Three regimes can be identified. The regime that is closest to a conventional MOT occurs when δ ≈ 1, and S >> 1. This is illustrated in figures 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 ∆ is increased such that δ 1 but S > 1, then the trap no longer forms at the quadrupole centre, but is displaced to where the Zeeman shift ∆ω z = ∆. 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 δ 1, and S ≤ 1. As in the power-broadened regime, the MOT position is determined by ∆ 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 Sr is 460 nK [23].

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 → J = 1 transitions are completely closed 1 , 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 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

Mathematical formalism
We simulate the 88 Sr 5s 2 1 S 0 → 5s5p 3 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 r = (x, y, z) where the unit vector directionsx,ŷ andẑ 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 ω i and wave-vector k i and the helicity of each beam is illustrated in figure 2(c). The nMOT quadrupole field B is defined as B = γ (xx + yŷ − 2zẑ) where γ is the gradient of the magnetic field. The magnetic field splits the 3 P 1 state into three Zeeman levels j with splittings ∆ω z = gµ B |B| / where g = 3/2 is the g-factor and µ B is the Bohr magneton.
To reproduce the macroscopic dynamics of the nMOT, the simulation is performed for ∼ 5000 atoms. Initially, the atoms are uniformly, randomly placed into an userdefined ellipsoid in the laboratory frame with position r. The atoms are also assigned random velocity vectors 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 2 . Typically, an initial temperature of T = 1 µK is used.
The total simulation time t tot is broken into time-steps of length δt. At each time step, the probability that each atom scatters a photon from laser i via a transition j is given by where ρ ij ee is the steady state excited state population derived from standard two-level optical Bloch equations [16] and δt = 0.1/Γ e such that P ij 1. The coupling strength W j is dependent on the local magnetic field and the polarisation 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ẑ is directed along the local magnetic field vector. Firstly, the total electric field for each laser beam E is defined in the laboratory frame in the spherical basisˆ q [32,33] as where E = E 1 , E 0 , E −1 and where 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 where U is the transformation from the cartesian basis to the spherical basis and R (θ) is the rotation matrix which mapsk i onto B by an angle θ. This leads to a new polarisation vector Once P ij is known, random sampling from a 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 where g = g (0, 0, −1) is the acceleration due to gravity and the prime notation represents the final atom position or velocity after a time step δt. If a scattering event does occur, the atom evolves as where k i is the wavevector of the laser from which the atom initially absorbed a photon withk i its associated unit vector, andk s is a random unit vector representing the direction of spontaneous emission 3 . 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 thex −ẑ plane and calculating the column density alongŷ. 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 v. This allows us to obtain the spatial, thermal and temporal dynamics of the atom cloud.

Experimental Configuration
The experiment is described in detail elsewhere [34][35][36][37] and so it is briefly summarised here. Initially, atoms from a strontium oven were slowed using a Zeeman slower before a 'blue-MOT' was formed on the 5s 2 1 S 0 → 5s5p 1 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 nm which addresses the 5s 2 1 S 0 → 5s5p 3 P 1 transition was applied. This light was artificially broadened to match the Dopplerbroadened profile of the atoms in the blue-MOT. After sufficient cooling, the broadening of the 689 nm light was removed, leaving single-frequency light, and a cold, dense nMOT. The nMOT was imaged using absorption imaging on the 5s 2 1 S 0 → 5s5p 1 P 1 transition. Images are captured on a Pixelfly QE camera with a post-magnification pixel-size of 8 µ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 x and z directions in the conventional way.
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/e 2 radii of 30 to 300 µm, densities up to 1 × 10 12 cm −3 and  temperatures as low as 460 nK. The vertical magnetic field gradient is held at 8 Gcm −1 during the final stage of the nMOT.

Testing the Model
As discussed in section 2, the properties of an nMOT are significantly dependent on ∆ 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 ∆. The top row of figure 3 shows experimental absorption images of the nMOT at four different values of ∆. It is clear that the MOT 'sags' under gravity and forms at lower positions as ∆ 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 positionz 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 ∆ in figures 4 and 5 respectively. The normalised residuals R ν [38] are shown below each figure. The vertical FWHM saturates as a function of ∆ 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 ∆. 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.
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 normalised detuning δ varies, the nMOT crosses between the different regimes identified in section 2.
The dependence of the nMOT temperature on ∆ for two different values of S is shown in 6(a). Firstly we consider an nMOT operating close to the quantum regime with S = 1.9. The temperature is essentially independent of ∆, since as shown in figure  4 the position of the MOT 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. 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 |∆| > 2π × 140 kHz, but begins to deviate significantly from experiment close to resonance. Here, the powerbroadened linewidth begins to approach the Zeeman splitting in the excited state. Thus the MOT 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 . 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 [38] were χ 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 Γ (S) ≈ ∆, as expected.

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). Optimising 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.
To compare the model to experimental data, we simulated the experiment performed by Ido et al. [39]. Their FORT consisted of two crossed laser beams with 1/e 2 radius of 28 µm, operating at 800 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 ∆ = −200 kHz with a total beam intensity of S = 18. While the nMOT was running, the FORT beams were then applied for a total time of 35 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

U0/kB (µK)
Tv (µK) 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 [39]. where U 0 is the trap depth, x 0 , y 0 and z 0 are linear offsets in thex,ŷ andẑ 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 ≈ 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.2 U 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 optimising the loading parameters of FORTs and optical lattices, eliminating the timeconsuming 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.

Conclusions
In summary, we have constructed a semi-classical Monte Carlo simulation in order to model the dynamics of an nMOT, in particular the 88 Sr 5s 2 1 S 0 → 5s5p 3 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 trialand-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 atomatom interactions such as those displayed by Rydberg atoms, leading to a greater understanding of strongly interacting many-body systems.