Multiple scattering is studied in a Cs magneto-optical trap (MOT). We use two Abel inversion algorithms to recover density distributions of the MOT from fluorescence images. Deviations of the density distribution from a Gaussian are attributed to multiple scattering.
©2005 Optical Society of America
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 .
The basic model of multiple photon scattering , predicts that the MOT density becomes constant as a function of the number of atoms in the trap and uniform at large atom number, Natoms ≫106. 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 . 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) . 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.
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  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, Natoms, is measured using a PMT.
In order to determine Natoms, the scattering rate per atom is approximated as 
where ωL is the trapping laser frequency. The natural line width of the 6P 3/2 state is Γ = 5.22MHz  and the detuning is δ = -1.4Γ. The coefficients = = .73 are taken from . Ωtot is the Rabi frequency . I is the total intensity of the trapping beams and the saturation intensity is Isat=1.1mW/cm2. Isat is chosen to be the σ + - σ - F=4→F′=5 saturation intensity because the C 1 and C 2 coefficients from  were used. The coefficients were obtained by fitting a two level atom model that used Isat=1.1 mW/cm2 to experimental data.
Natoms is found by measuring the fluorescence power P emitted from the MOT, Natoms = P/R. The fluorescence from the MOT is collected onto the PMT. Measurements of Natoms 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 Natoms. Natoms was limited by our setup, primarily by the trap loading rate, trapping laser intensity and magnetic field gradient. We varied Natoms 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 Natoms. In the multiple scattering regime, the density as a function of Natoms becomes constant. Our measurements indicate that multiple scattering occurs at a density of nms ~ 2.2 × 1010 cm-3 . As can be seen in Fig. 2 (a), the density increases linearly with Natoms until nms is achieved. The abrupt change from the temperature limited to the multiple scattering regime is clearly visible in Fig. 2 (a).
A plot of volume versus Natoms is shown in Fig. 2 (b). The volume is constant until Natoms =Nms ~ 2.0 × 106. At Nms, the volume starts to grow linearly. The linear growth of the volume as a function of Natoms is a signature of multiple scattering. The density and Natoms can be combined to verify the minimum value of the multiple scattering radius rms = .2 mm . This measurement agrees with prior work .
The MOT radius scales as in the multiple scattering regime . The scaling was observed experimentally for Rb . Our data is consistent with a dependence on MOT radius as shown in Fig. 2 (c). Multiple scattering occurs for > 125 in our Cs MOT. This value agrees with other experiments with Cs .
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, kB 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
We measured the temperature as a function of detuning for a total trap laser intensity of I= 30 mW/cm2. 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 . The data in  was measured using absorption time-of-flight spectroscopy.
To summarize, our measurements of Natoms, 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 . 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 ,
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
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 . 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 fk (r, z), where the Ck are expansion coefficients, the projection is
The Gk (x,z) are the Abel transform of the basis functions fk (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 Ck 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 Ck . Finally, the expansion coefficients are used to reconstruct the density distribution.
If we let i = 1…Nx and j = 1…Nz where Nx×Nz 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 (xi ,zj ) are the pixel coordinates in the projection space . The function h(x - xi ,z - zj ) 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 - xi ,z - zj ) = δ(x - xi )δ(z - zj ), where (xi ,zj ) is the pixel coordinate on the CCD image.
For axially symmetric distributions, fk (r,z) is separable, fk (r,z) = ρm (r)ξn (z). The separability of fk (r,z) makes it convenient to define
The definition of Z nj 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 . 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 . 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 . The images and density distributions that we measured were all predominantly Gaussian. Since Natoms < 107, 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 further confirms that the MOT is in this regime for the experiment [2, 21].
The Butterworth filter is a standard filter to apply for digital image processing . The Butterworth filter is
where k is the spatial frequency in units of inverse length. k 0 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 0 determines how many spatial frequencies will be used to reconstruct the distribution. For the Fourier-Hankel and BASEX method, we set k 0 = 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 R2 was <1 part in 103 (0.04%) while for Fourier-Hankel the error was <1 part in 104 (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 . 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 . This is one of the reasons that the BASEX method can be superior for data with high resolution and large dynamic range.
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 2 radius to calculate the trap volume.
It has been shown that, when large numbers of atoms, Natoms > 106, 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 , causing the density to become constant as greater numbers of atoms are trapped. A second consequence is that, for Natoms ≫ 106, 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 Natoms < 106.
For Natoms < 106, 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 Natoms where we expect to observe perturbations of a Gaussian density distribution. As Natoms 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 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 R2 values of a Gaussian fit to the images and density distributions are plotted against Natoms 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 Natoms >2×106=Nms. Note that the scatter in the BASEX data at low Natoms 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 Natoms 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 Natoms 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 Natoms 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 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.
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  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).
References and links
1 . D. W. Sesko , T. G. Walker , and C. E. Wieman , “ Behavior of neutral atoms in a spontaneous force trap ,” J. Opt. Soc. Am. B 8 , 946 ( 1991 ). [CrossRef]
2 . C. G. Townsend , N. H. Edwards , C. J. Cooper , K. P. Zetie , C. J. Foot , A. M. Steane , P. Szriftgiser , H. Perrin , and J. Dalibard , “ Phase-space density in the magneto-optical trap ,” Phys. Rev. A 52 , 1423 ( 1995 ). [CrossRef] [PubMed]
3 . M. Drewsen , Laurent P. , A. Nadir , G. Santarelli , A. Clairon , Y. Castin , D. Grison , and C. Salomon , “ Investigation of sub-Doppler cooling effects in a cesium magneto-optical trap ,” Appl. Phys. B 59 , 283 ( 1994 ). [CrossRef]
4 . S. Grego , M. Colla , A. Fioretti , J. H. Muller , P. Verkerk , and E. Arimondo , “ A cesium magneto-optical trap for cold collision studies ,” Opt. Commun. 132 , 519 ( 1996 ). [CrossRef]
5 . C. Gabbanini , A. Evangelista , S. Gozzini , A. Lucchesini , A. Fioretti , J. H. Muller , M. Colla , and E. Arimondo , “ Scaling laws in magneto-optical traps ,” Europhys. Lett. 37 , 251 ( 1997 ). [CrossRef]
6 . A. Vorozcovs , M. Weel , S. Beattie , S. Cauchi , and A. Kumarakrishnan , “ Measurements of temperature scaling laws in an optically dense magneto-optical trap ,” J. Opt. Soc. Am. B 22 , 943 ( 2005 ). [CrossRef]
9 . L. M. Smith , D. R. Keefer , and S. I. Sudharsanan , “ Abel Inversion Using Transform Techniques ,” J. Quant. Spect. Radiat. Transfer 39 , 367 ( 1988 ). [CrossRef]
10 . V. Dribinski , A. Osssadtchi , V. A. Mandelshtam , and H. Reisler , “ Reconstruction of Abel-transformed images: The Gaussian basis-set expansion Abel transform method ,” Rev. Sci. Instrum. 73 , 2634 ( 2002 ). [CrossRef]
11 . K. R. Overstreet , J. Franklin , and J. P. Shaffer , “ Zeeman effect spectroscopically locked Cs diode laser system for atomic physics ,” Rev. Sci. Instrum. 75 , 4749 ( 2004 ). [CrossRef]
13 . T.M. Brzozowski , M. Maczynska , M. Zawada , J. Zachorowski , and W. Gawalik , “ Time-of-flight measurement of the temperature of cold atoms for short trap-probe beam distances ,” J. Opt. B: Quant. Semiclass. Opt. , 4 , 62 ( 2002 ). [CrossRef]
14 . C.D. Wallace , T.P. Dinneen , K.Y.N. Tan , A. Kumarakrishnan , P.L. Gould , and J. Javanainen , “ Measurements of temperature and spring constant in a magneto-optical trap ,” J. Opt. Soc. Am. B 11 , 703 ( 1994 ). [CrossRef]
15 . J.P. Shaffer and Ph.D. Thesis , University of Rochester ( 1999 ).
16 . D. Hoffmann , P. Feng , and T. Walker , “ Measurements of Rb trap-loss collision spectra ,” J. Opt. Soc. Am. B 11 , 712 ( 1994 ). [CrossRef]
17 . R. N. Bracewell , The Fourier Transform and Its Applications ( McGraw-Hill , 2000 ).
18 . S. Kawata and Y. Ichioka , “ Iterative image restoration for linearly degraded images. II. Reblurring procedure ,” J. Opt. Soc. Am. 70 , 768 ( 1980 ). [CrossRef]
19 . C. R. Gebhardt , T. P. Rakitzis , P. C. Samartzis , V. Ladopoulos , and T. N. Kitsopoulos , “ Slice imaging: A new approach to ion imaging and velocity mapping ,” Rev. Sci. Instrum. 72 , 2001 ( 2001 ). [CrossRef]
20 . G. H. Golub , P. C. Hansen , and D. P. O’Leary , “ Tikhonov Regularization and Total Least Squares ,” SIAM J. Matrix Anal. Appl. 21 , 185 ( 1999 ). [CrossRef]
22 . K. Iizuka , Engineering Optics ( Springer-Verlag , 1987 ).