1. Introduction

The physical picture of a gas of laser cooled and trapped atoms inside a magneto-optical trap (MOT) is very complex despite the fact that the MOT is commonplace in atomic physics experiments. Many models have been introduced to predict the behavior of the MOT, for examples see [1,2]. One important effect that has received attention is multiple photon scattering [1].

The basic model of multiple photon scattering [1], predicts that the MOT density becomes constant as a function of the number of atoms in the trap and uniform at large atom number, N_atoms ≫10⁶. A uniform density results in a flat-topped fluorescence intensity projection instead of a Gaussian one. The change of the density distribution as more atoms are added to the trap is a consequence of the net repulsive radiation pressure caused by multiple photon scattering.

Signatures of multiple scattering were observed in several studies, but at different atom number [2, 3, 4, 5, 6]. In some of these studies, the fluorescence intensity projection of the MOT was observed and it did not deviate from a Gaussian [2, 3]. Whether and how the density distribution of the MOT is affected by multiple scattering is a question that can give insight into light pressure forces.

In this paper, we measure the density distribution of atoms in a Cs MOT using computed tomography. We observe that the Cs MOT maintains a near Gaussian density distribution into the multiple scattering regime. The deviation from a pure Gaussian density distribution that is observed can be interpreted as a distortion of the linear restoring force that forms the MOT caused by the radiation trapping force [1]. Our experiment suggests that the differences between the observations of the MOT density distribution can be accounted for by the difficulty in interpreting 2-D projections of MOT fluorescence as a direct measure of the density distribution.

Measuring the density of laser atom traps is difficult but important to many experiments, particularly those involving collisions between cold atoms. Techniques for analyzing MOT density and volume rely on images taken of MOT fluorescence or absorption imaging. The images are 2-D projections of 3-D distributions, Fig. 1.

3-D distributions can be reconstructed from 2-D projections using computed tomography. The 3-D density distribution of the MOT can be recovered using an Abel transformation of CCD images of MOT fluorescence. The Abel transform is well suited to this task because the magnetic field used for the MOT is axially symmetric. The method is a non-destructive way to measure the density distribution of the trap. The technique can also be used in a similar fashion to investigate other axially symmetric traps and systems, such as optical dipole traps and ultracold plasmas. The Abel inversion has been used successfully to analyze ion images in physical chemistry and flames in combustion science [7, 8].

There are many different approaches to calculating the Abel transform. Each method varies in complexity, computation time, and sensitivity to noise. We chose to use two methods: the Fourier-Hankel method [8, 9], and the Gaussian basis-set expansion (BASEX) [10]. We find that both methods are useful and discuss the advantages and implementation of each one. If the CCD is calibrated, the Abel inversion can be done in real time, as fast as the acquisition rate of the CCD camera, to produce a density distribution.

 

Fig. 1. Mapping of a 3-D distribution onto a 2-D plane. Also depicted is the projection of a slice at constant z.

Download Full Size | PDF

2. Experimental setup and measurements

Our setup consists of a Cs vapor cell MOT. The chamber background pressure is 3 × 10^-10 Torr. The Cs vapor pressure is 3 - 5 × 10^-10 Torr when the trap functions. Shim coils are used to null the earth’s field. The MOT is formed between two plates, located 26.5 cm above a Z-stack multi-channel plate (MCP) detector. The trapping beams are provided by a single diode laser system [11] with a trapping beam radius of 1.2cm. Another diode laser system provides the repumping light.

During the experiments, the MOT fluorescence is imaged onto a CCD array. The CCD array contains 1004 × 1004 pixels each with dimensions of 7.4 μm × 7.4 μm. The number of atoms in the trap, N_atoms, is measured using a PMT.

In order to determine N_atoms, the scattering rate per atom is approximated as [2]

where ω_L is the trapping laser frequency. The natural line width of the 6P _3/2 state is Γ = 5.22MHz [2] and the detuning is δ = -1.4Γ. The coefficients $C_1^2$ = $C_2^2$ = .73 are taken from [2]. Ω_tot is the Rabi frequency ${\Omega }_{\mathrm{tot}}=\Gamma \sqrt{\frac{I}{2I_{\mathrm{sat}}}}$ . I is the total intensity of the trapping beams and the saturation intensity is I_sat=1.1mW/cm². I_sat is chosen to be the σ ⁺ - σ ^- F=4→F′=5 saturation intensity because the C ₁ and C ₂ coefficients from [2] were used. The coefficients were obtained by fitting a two level atom model that used I_sat=1.1 mW/cm² to experimental data.

N_atoms is found by measuring the fluorescence power P emitted from the MOT, N_atoms = P/R. The fluorescence from the MOT is collected onto the PMT. Measurements of N_atoms using a PMT and Abel inverted CCD images allow the MOT density to be determined.

To study the onset of multiple scattering, we used a combination of laser intensities and magnetic field gradients (10–15 G/cm) to investigate the density distribution of the trap as a function of N_atoms. N_atoms was limited by our setup, primarily by the trap loading rate, trapping laser intensity and magnetic field gradient. We varied N_atoms to investigate a region that spanned the temperature limited and multiple scattering regimes, Fig. 2 (a).

One signature of multiple scattering is a linear scaling of volume with N_atoms. In the multiple scattering regime, the density as a function of N_atoms becomes constant. Our measurements indicate that multiple scattering occurs at a density of n_ms ~ 2.2 × 10¹⁰ cm^-3 . As can be seen in Fig. 2 (a), the density increases linearly with N_atoms until n_ms is achieved. The abrupt change from the temperature limited to the multiple scattering regime is clearly visible in Fig. 2 (a).

 

Fig. 2. (a.) MOT density versus N_atoms trapped. (b.) MOT volume versus N_atoms trapped. (c.) ${\mathrm{N}}_{\text{atoms}}^{1/3}$ versus MOT radius. (d.) Temperature as a function of detuning as measured from time-of-flight velocity distributions for I = 30mW/cm² and magnetic field gradient dB/dz= 15G/cm. Open triangles (△) represent data analyzed using the BASEX transform, closed circles (●) were analyzed using the Fourier-Hankel transform.

Download Full Size | PDF

A plot of volume versus N_atoms is shown in Fig. 2 (b). The volume is constant until N_atoms =N_ms ~ 2.0 × 10⁶. At N_ms, the volume starts to grow linearly. The linear growth of the volume as a function of N_atoms is a signature of multiple scattering. The density and N_atoms can be combined to verify the minimum value of the multiple scattering radius r_ms = .2 mm [2]. This measurement agrees with prior work [1].

The MOT radius scales as ${\mathrm{N}}_{\text{atoms}}^{1/3}$ in the multiple scattering regime [12]. The ${\mathrm{N}}_{\text{atoms}}^{1/3}$ scaling was observed experimentally for Rb [6]. Our data is consistent with a ${\mathrm{N}}_{\text{atoms}}^{1/3}$ dependence on MOT radius as shown in Fig. 2 (c). Multiple scattering occurs for ${\mathrm{N}}_{\text{atoms}}^{1/3}$ > 125 in our Cs MOT. This value agrees with other experiments with Cs [4].

The temperature of the MOT was measured by observing the expansion of the 89D Rydberg state. High n Rydberg states (n ~ 100) are excited with a Coherent 699–21 ring laser at ~ 508 nm. The ring laser is focused through the MOT with a beam waist of 110 μm and gated by an acousto-optic modulator. Upon excitation, the Rydberg atoms expand at their thermal velocity. After a delay, the Rydberg atoms are pulse-field ionized and projected onto the MCP detector. The ion signal is fed into a multi-channel analyzer and a time-of-flight velocity distribution is accumulated. The experiment is repeated at a rate of 1kHz. The ion count rate is ~ 10 - 100ions/s. This rate assures that there is no significant broadening of the velocity distribution due to Coulomb repulsion.

The thermal distribution for the atoms in the overlap region of the MOT and laser focus is

σ_i represents the standard deviation of the focal spot size in the respective dimension, t is the delay time, k_B is Boltzmann’s constant and T is the temperature. Since the time-of-flight distribution is measured in one dimension, z, Eq. (2) has to be integrated over all velocities and the x and y directions. After integration, the spatial distribution is

Eq. (3) is a Gaussian with a standard deviation ${\sigma }_z\left(t\right)=\sqrt{\frac{k_{B}T}{m}{t}^2+{\sigma }_z^2}$ [13]. The spot size of 110μm was measured with a CCD and verified by fitting the data to σ_z (t).

We measured the temperature as a function of detuning for a total trap laser intensity of I= 30 mW/cm². The results are shown in Fig. 2 (d). The measurement shows the same temperature dependence as a function of trap laser detuning that was observed in [1]. The data in [1] was measured using absorption time-of-flight spectroscopy.

To summarize, our measurements of N_atoms, volume and temperature are found to be in good agreement with prior work [1, 2, 3, 4, 5, 6, 14, 15, 16]. The signature of multiple scattering is clearly observed in our data. The measurements give us confidence that our experimental setup produces results consistent with other experiments.

3. Theory of the Abel inversion algorithms

The density distribution of the MOT can be recovered from CCD images of trap fluorescence. A CCD image is the 2-D projection of a 3-D density distribution, Fig. 3. In cases with axial symmetry, the 3-D distribution can be reduced to 2-D by the replacement $r=\sqrt{x^2+y^2}$. The 3-D density distribution is represented by a function I(r,z). The 2-D projection of I(r,z) is

where x is the distance from the symmetry axis in projection space. For the Abel transform, each value of z in the projection is treated independently. One line of pixels for constant z is referred to as a “slice” of the image. One slice of the distribution and its projection are illustrated in Fig. 1. The 3-D distribution I(r, z) can be recovered from the 2-D mapping P(x,z) by using the inverse Abel transform [17],

To implement the inverse Abel transform to recover I(r, z), a numerical algorithm is required. The primary difficulties in evaluating Eq. (5) are the derivative and the singularity. The derivative dP(x,z)/dx is problematic when the image has sharp features or contains excessive noise. Any practical application of the Abel inversion requires smoothing of the data because of the derivative in Eq. (5). Smoothing is most easily accomplished by filtering. In our work we applied a low pass Butterworth filter.

The Abel inversion requires axial symmetry, but it can be a challenge to properly identify the symmetry axis. The symmetry axis of the projected distribution can be difficult to identify numerically due to noise and systematic errors in the experiment. For a MOT, beam misalignment, beam shape and beam intensity imbalances predominantly contribute to asymmetry. Great care was taken to ensure a symmetric MOT.

Routines for symmetrizing experimental data for inversion have been developed [7, 8, 9, 18]. The method that we used relies on exploiting symmetry. Axial symmetry implies that the Fourier components of the 2-D projection should be real and even. Removing the imaginary components in reciprocal space leaves a symmetric projection.

The Fourier-Hankel approach [8, 19] starts by performing a fast-Fourier transform (FFT) on each slice. In Fourier-space, the CCD image is filtered and recentered. Symmetrization and smoothing are followed by a discrete, inverse Hankel-transform.

We can recover the distribution I(r) of a slice from its projection P(x) by first considering the Fourier transform of P(x),

Transforming to polar coordinates and substituting the zero-order Bessel function allows Eq. (6) to be written

Eq. (7) is the zero-order Hankel transform of I(r). Because the zero-order Hankel transform is identical to its inverse [17], we can recover I(r) using

The procedure for numerically implementing Eq. (8) begins with integrating over the spatial coordinate by performing a FFT. The real part of the transformed projection is then symmetrized and smoothed. Finally, the inverse Hankel transform is numerically applied to recover the distribution I(r). The procedure is repeated for each slice.

The BASEX method is an alternative to the Fourier-Hankel method [10]. In cases where the projection has high resolution or a large dynamic range, the BASEX method can be more efficient. The essential feature of the BASEX method is an expansion of the image in a set of functions that have analytic Abel inversions.

The BASEX method can be understood if we recognize that the Abel transform is linear. If we expand the distribution I(r, z) in terms of N basis functions f_k (r, z), where the C_k are expansion coefficients, the projection is

The G_k (x,z) are the Abel transform of the basis functions f_k (r, z). Choosing basis functions that have analytic Abel inversions avoids the problems with taking derivatives of noisy data. The basis can also be chosen to obtain rapid convergence of the expansion.

Calculation of the expansion coefficients C_k is most easily accomplished using linear algebra. The method begins by filtering and symmetrizing the image in frequency space. The CCD image is then returned to projection space and expanded in a basis of Abel transformed functions one slice at a time to determine the C_k . Finally, the expansion coefficients are used to reconstruct the density distribution.

If we let i = 1…N_x and j = 1…N_z where N_x×N_z is the number of pixels in the projection, we can cast the projection as a matrix multiplication

or more concisely P = CG. The Abel transform of the basis functions is

where (x_i ,z_j ) are the pixel coordinates in the projection space [10]. The function h(x - x_i ,z - z_j ) is determined by experiment. It is an instrumental weight that can be set equal to a delta function if no instrumental weighting is needed. We used h(x - x_i ,z - z_j ) = δ(x - x_i )δ(z - z_j ), where (x_i ,z_j ) is the pixel coordinate on the CCD image.

For axially symmetric distributions, f_k (r,z) is separable, f_k (r,z) = ρ_m (r)ξ_n (z). The separability of f_k (r,z) makes it convenient to define

and

The definition of Z n_j allows us to determine X _mi for a fixed z. This is analogous to treating a slice in the Fourier-Hankel method. The 2-D projection is

or in matrix form P = X ^T CZ. Solving for C gives C = APB, where A = (XX ^T + λ_x I)^-1 X and B = Z ^T(ZZ ^T + λ_z I)^-1. λ_x and λ_z are the regularization parameters for X and Z respectively. Regularization is necessary if the basis functions ρ_m (r) and ξ_n (z) under or over-determine the matrices X and Z [20]. Regularization keeps A and B from becoming singular in a numerical calculation.

The basis function ρ_m (r) needs to have an analytic Abel inversion. For this reason, we chose

where σ is equal to the pixel size [10]. Similarly, ξ_n (z) = ρ_n (z) simplifies the problem. For the BASEX method, the resolution of the CCD influences the size of the basis. Basis set computation can take considerable time. Once computed, though, the basis set can be reused to analyze images taken with the same CCD imaging system. We chose the number of basis functions to be half the number of pixels in each dimension of the image so that our matrices were well conditioned. The images were filtered and symmetrized in Fourier-space as in the Fourier-Hankel method.

4. Comparison of methods

We did not observe irregularities in the images or inversions that indicated preferential absorption within the MOT or that the MOT was optically thick enough to affect the inversions. It would be expected that absorption of photons from the far side of the MOT by atoms on the near side would be observed as the flat-topped distribution suggested in [1]. The images and density distributions that we measured were all predominantly Gaussian. Since N_atoms < 10⁷, the optical thickness will not be large enough for a photon to scatter more than twice before escaping the trap cloud. The fact that the radius of the trap scales like ${\mathrm{N}}_{\text{atoms}}^{1/3}$ further confirms that the MOT is in this regime for the experiment [2, 21].

 

Fig. 3. Comparison of raw image to the Fourier-Hankel and BASEX transform methods. Additional structure is visible in both transformed images. The raw image was taken for N_atoms = 2.0×10⁶, a laser intensity of 30mW/cm², and a magnetic field gradient dB/dz = 15G/cm.

Download Full Size | PDF

The Butterworth filter is a standard filter to apply for digital image processing [22]. The Butterworth filter is

where k is the spatial frequency in units of inverse length. k ₀ sets the cutoff frequency and n is the order of the filter. One advantage of the Butterworth filter is that it maintains it’s shape for higher orders. The order sets the sharpness of the cutoff.

The Butterworth filter was chosen because it has a linear response over the pass band, but does not have a sharp cutoff compared to other standard filters. Despite the rather broad transition band, the stop band of a Butterworth filter has large, essentially infinite, attenuation. The phase response is more linear than a Chebyshev or elliptic filter and as a result does a reasonable job of minimizing distortion. The application of the low pass Butterworth filter is similar to Gaussian smoothing in the image domain.

We took great care not to over filter the data in our experiments. Because the MOT is ~ 1mm in size, k ₀ determines how many spatial frequencies will be used to reconstruct the distribution. For the Fourier-Hankel and BASEX method, we set k ₀ = 7mm^-1 with order n = 15. The cutoff for the spatial frequency is less than the Nyquist frequency to avoid aliasing but large enough to accurately reproduce the features of the image. Consistency of the results of the two methods was checked for each CCD image.

The differences between the inversion techniques can be seen in Fig. 3. We estimated the methodological error for the Fourier-Hankel and BASEX algorithms by running a computer generated Gaussian through each routine. The BASEX error in R² was <1 part in 10³ (0.04%) while for Fourier-Hankel the error was <1 part in 10⁴ (0.004%). Fig. 3 shows a typical data set viewed perpendicular to the symmetry axis. Both inversion methods show the existence of structure that is masked in the 2-D projection. The density distribution is more clearly defined for the Abel inversion. The Fourier-Hankel transform is smoother because the filter cutoff is low enough to prevent oversampling. Oversampling can cause an amplification of noise because when the difference between neighboring points of the image decreases towards the noise level the inversion becomes inaccurate [7]. The biggest problem with noise when using the Fourier-Hankel method is oversampling. The BASEX method shows more structure because the basis functions are narrow enough to reproduce sharp features [10]. This is one of the reasons that the BASEX method can be superior for data with high resolution and large dynamic range.

 

Fig. 4. Correlation of the (a.) raw images, (b.) Fourier-Hankel transform (●) and BASEX (⊕) to a Gaussian fit. The inset in (b.) shows an adjacent averaging of the R² data.

Download Full Size | PDF

There is a slight difference in the volume calculated using the two different inversion methods. It can be seen in Fig. 2 (b)–(c) that the Fourier-Hankel method returns density distributions that indicate larger volumes. In the BASEX method the basis functions are able to produce more structure that can reproduce abrupt changes in the distribution. We also found it much more difficult to implement the centering algorithm for the BASEX method. Misidentifying the center by as much as a single pixel caused errors in the BASEX inversion around the symmetry axis, r ~ 0. The centering errors also contributed to the differences in volume. Both methods avoid having to make a choice between the 1/e and the 1/e ² radius to calculate the trap volume.

5. Discussion

It has been shown that, when large numbers of atoms, N_atoms > 10⁶, are confined in a MOT, multiple photon scattering plays an important role in the trap [1,12]. Multiple photon scattering leads, under most conditions, to a repulsive force that perturbs the MOT. The multiple scattering force is a sum of the restoring force due to attenuation of the trap lasers by the atoms and the repulsive radiation trapping force.

One signature of multiple photon scattering is that the size of the MOT will scale as ${\mathrm{N}}_{\text{atoms}}^{1/3}$, causing the density to become constant as greater numbers of atoms are trapped. A second consequence is that, for N_atoms ≫ 10⁶, the density of the trap becomes uniform. The resulting CCD image of MOT fluorescence becomes flat-topped consistent with the Abel inversion of a uniform density distribution.

The uniform density effect is a result of the balance between the MOT trapping force and the multiple scattering force. The high atom number regime is to be contrasted with N_atoms < 10⁶.

For N_atoms < 10⁶, the multiple scattering force is weak compared to the trapping force. The trapping force is effectively linear. The MOT behaves like a Boltzmann gas in a harmonic potential and the density distribution is Gaussian.

It is expected that the density distribution of the MOT makes a continuous transformation from Gaussian to flat-topped as the multiple scattering force grows larger. Our data spans a region of N_atoms where we expect to observe perturbations of a Gaussian density distribution. As N_atoms increases and the multiple scattering force competes with the MOT restoring force, we should observe less correlation with a Gaussian. This regime is interesting because many experiments have observed the ${\mathrm{N}}_{\text{atoms}}^{1/3}$ scaling of trap size but not detected any change of the density distribution [2, 4, 14, 15].

We tested the hypothesis that the density distribution is perturbed by multiple photon scattering by fitting a Gaussian to CCD images and densities obtained with the Abel inversion algorithms. To fit the data, we used a Levenberg-Marquardt non-linear least squares method. The R² values of a Gaussian fit to the images and density distributions are plotted against N_atoms in Fig. 4. Although other parameterizations of the density distribution have been proposed [4, 16], we did not use them. We explored deviations of the density distribution from a Gaussian that are expected for a Boltzmann gas in the harmonic potential of the trap.

The two Abel inversion methods show a departure from a Gaussian density distribution, as expected in the multiple scattering regime, of ~0.5%. The decrease in correlation shown in Fig. 4 is evidence for a redistribution of density. The departure is slight and is not evident in the CCD images. It shows up as a reduced correlation to a Gaussian at larger atom number, Fig. 4 inset, where N_atoms >2×10⁶=N_ms. Note that the scatter in the BASEX data at low N_atoms is a result of noise in the CCD images. In general, the centering difficulties with the BASEX method caused the correlations obtained using that method to show more scatter than the Fourier-Hankel data. The spread in the data at large N_atoms may partially be caused by the sensitivity of the rescattering force to laser detuning and intensity (through the Doppler effect and optical pumping) as well as the magnetic field gradient. The changes in these parameters that were necessary to vary N_atoms also affect the rescattering force. We checked for direct correlations with each of the variables but did not observe any correlation with a single parameter. Although the deviations are small in Fig. 4, it is anticipated because N_atoms is at the low end of the multiple scattering regime. The MOT restoring force still dominates the density distribution.

Our results clarify earlier experiments that observed a trap size that scaled as ${\mathrm{N}}_{\text{atoms}}^{1/3}$ but did not detect a perturbation of the MOT density distribution [2, 4, 14, 15]. These results can be explained by the size of the perturbations at the onset of the multiple scattering regime. The deviations from a Gaussian distribution are small and not easily measured, particularly when the density distribution is not determined but the 2-D projection of fluorescence is analyzed instead.

6. Conclusions

We have shown that the density distribution of a MOT can be obtained reliably through Abel inversion of trap fluorescence images. The method is non-destructive so it can be used to monitor changes in trap density as an experiment is performed. For atom traps, the Fourier-Hankel method is probably better than the BASEX method because it is less sensitive to the centering algorithm, the MOT image does not have a large dynamic range and the MOT has a smooth distribution. Comparing both methods is a useful measure of accuracy.

We used the Abel inversion method to show that the density distribution of a Cs MOT does become perturbed at the onset of the multiple scattering regime. The differences in previous observations have most likely resulted from the fact that there is a continuous transition from the pure harmonic MOT, where the density distribution is Gaussian, to the extreme case of multiple scattering, where the density is uniform.

Further theoretical work on the multiple scattering regime is complicated. A modification of the simple model in [1] that accounts for optical pumping, the Doppler effect, the magnetic field and polarization may be able to describe the multiple scattering force. A new theory could yield insight into atom field interactions and light pressure forces.

Finally, Abel transform techniques are not limited to studies of the MOT. They can be applied to other types of axially symmetric traps, like optical dipole traps, and images of cold plasmas. The programs used for this work are available from the authors.

We would like to acknowledge funding from the Research Corporation, the OSRHE and the Air Force Office of Scientific Research (FA9550-05-0328).