Resonant absorption imaging is a common technique for detecting the two-dimensional column density of ultracold atom systems. In many cases, the system's thickness along the imaging direction greatly exceeds the imaging system's depth of field, making the identification of the optimally focused configuration difficult. Here we describe a systematic technique for bringing Bose-Einstein condensates (BEC) and other cold-atom systems into an optimal focus even when the ratio of the thickness to the depth of field is large: a factor of 8 in this demonstration with a BEC. This technique relies on defocus-induced artifacts in the Fourier-transformed density-density correlation function (the power spectral density, PSD). The spatial frequency at which these artifacts first appear in the PSD is maximized on focus; the focusing process therefore both identifies and maximizes the range of spatial frequencies over which the PSD is uncontaminated by finite-thickness effects.

## I. INTRODUCTION

Since the most important technique for obtaining properties of ultracold atoms is direct imaging, a well-designed, and well-aligned imaging system is crucial for obtaining high quality data which is valid at all length scales. While large scale properties such as the system's width or peak density can be obtained with little effort, significant care must be taken for experiments requiring very good spatial resolution,^{1,2} or those studying correlations.^{3,4} It is difficult to bring objects extended along the imaging axis, such as degenerate Fermi gases,^{5,6} 3D Mott insulators,^{7} and Bose-Einstein condensates (BECs),^{8,9} into focus particularly after time-of-flight (TOF) expansion because their spatial thickness often exceeds the imaging system's depth of field. Even for such objects, a high degree of accuracy in focusing is required to minimize imaging artifacts. Understanding and minimizing these artifacts is particularly important when studying density-density correlations, where the artifacts can be confused with the correlation signal under study.^{4,10–13} Here we describe a fairly generic technique for focusing on these extended objects which is far more precise than simply optimizing the “sharpness” of imaged atom clouds.

Absorption imaging is a ubiquitous approach for measuring the density distribution of ultracold atom systems.^{14} A probe beam illuminates the atomic system, and the resulting shadow is imaged onto a scientific camera, typically a charge-coupled device (CCD) or complementary metal-oxide-semiconductor (CMOS) detector. Ideally, the fraction of light absorbed would be directly related to the two-dimensional column density ρ_{2D}(*x*, *y*) = ∫d*z* ρ(*x*, *y*, *z*) of the atoms along the imaging direction **e**_{z}, where ρ(*x*, *y*, *z*) is the density of atoms. If the thickness δ*z* along **e**_{z} exceeds the imaging system's depth of field, then some of the atomic distribution must necessarily be out of focus, invalidating any simple relationship between absorption and column density. Given this, it is a challenge to obtain the optimal focal plane of the extended system that minimizes the artifacts resulting from this defocus, e.g., at the center of a distribution symmetric along **e**_{z}.

Typically a system is brought into a focus by minimizing the size or apparent diffraction effects from a compact object such as a trapped BEC; in many cases, no such compact reference at the desired image plane is available. In this paper, we present a technique for determining the optimal focus of absorption-imaged extended objects. Using this technique, we identify the focal plane within an accuracy of 2 μm for a δ*z* = 150 μm thick object. Specifically, given an object with density-density correlations^{4} with a spatial correlation length ℓ, we show that observations of correlations in the optical absorption as a function of camera position allow us to bring the object into focus to within a fraction of the depth of field associated with ℓ, even without knowing the details of the correlation function. This optimal focus is the camera position where the imaged auto-correlation function (ACF) most accurately reflects the atomic density-density correlations, minimizing both defocus-induced artifacts and the resolution limiting effect of the system's finite thickness.^{12,13}

In this paper, we review the basic theoretical formulation required to understand light propagating through an absorbing dielectric medium. We then consider several example images created by different idealized objects, in each case noting how to determine their optimal focus. Lastly, we experimentally apply this technique to images of BECs after TOF.

## II. THEORY

Monochromatic light of free-space wavelength λ and wavenumber *k*_{0} = 2π/λ propagating through an object with complex relative permittivity ε(**r**) = ε/ε_{0} and relative susceptibility χ(**r**) = ε(**r**) − 1, where ε is the permittivity, and ε_{0} is the electric constant, is described by the vectorial wave equation for the electric field **E(r)**:

In a medium where ε(**r**) is slowly varying, the right-hand side (rhs) of Eq. (1) can be neglected, reducing Eq. (1) to separate scalar wave equations

for each vector component of **E(r)**, e.g., we might have **E(r)** = *E***(r)****e**_{x} for light linearly polarized along **e**_{x}.

### A. Wavefield propagation

Here, we cast the above scalar wave equation into the form

suitable for light predominantly traveling along **e**_{z}. For a known field configuration at *E*(**r**) (such as the probe laser before it interacts with the atoms), Eq. (2) has the formal solution

describing the field propagated a distance Δ*z* along **e**_{z} [plus minus sign].

Wave propagation in free space [i.e., χ(**r**) = 0 in Eq. (2)] is solved exactly in the angular spectrum representation^{15}

for a forward going wave, with the 2D position **r**_{2D} = (*x*, *y*) and wavevector **k**_{2D} = (*k*_{x}, *k*_{y}); the Fourier-transformed wavefield

*z*in free space

The transfer function behaves differently in two regions of spatial frequencies: for

Meanwhile, considering only χ(**r**) [neglecting the first term in the rhs of Eq. (2)], the absorption and refraction of light traveling a distance Δ*z* is described by

Unlike the usual Beer-Lambert (BL, discussed in Sec. II B), this expression alone does not reflect a good approximation to beam propagation for systems of any significant thickness.

### B. Beer-Lambert law and the paraxial approximation

To better understand the independent influence of the beam's propagation and its interaction with matter, we apply the paraxial approximation to Eq. (2), allowing us to draw an analogue between the paraxial wave equation and the Schrödinger equation, which can be solved numerically using a split-step Fourier method (SSFM).^{16} To understand the difference between Eq. (5) and the usual BL law, we again turn to Eq. (2), now assuming that the electric field can be written as *E*(**r**) = exp (*ik*_{0}*z*)*E*^{′}(**r**), where *E*^{′}(**r**) is a slowly varying envelope along **e**_{z}. Inserting this form into Eq. (2) gives the paraxial wave equation

where the assumed weak *z* dependence of *E*^{′}(**r**) allowed us to drop

*E*

^{′}(

**r**) can be partitioned into a spectral part

**Q**

^{′}(Δ

*z*), with

For the paraxial approximation to be valid, the condition |χ(**r**)| ≪ 1 must also hold: otherwise the **Q**^{′}(Δ*z*) evolution would lead *E*^{′}(**r**) to depend strongly on *z*.

We numerically evolve the paraxial wave equation [Eq. (6)] along **e**_{z} using a split-step Fourier method (SSFM),^{17,18} where the operators in the rhs of Eq. (6) are split into two: one operator represents wave propagation in a uniform medium using Eq. (7) and the other operator takes into account the effect of refractive index variation using Eq. (8). In the SSFM, we alternately apply the two evolution operators with steps of size Δ*z*. For each step, the complex amplitude *E*′(**r**) is propagated first by **P**′(Δ*z*/2), then by **Q**′(Δ*z*), and then again by **P**′(Δ*z*/2). The resulting symmetrized split evolution

has its first correction at order Δ*z*^{3}.

The paraxial equations allow us to introduce the *depth of field*

where

*k*

_{2D}of interest and

*l*

_{min}= 2π/

*k*

_{max}is the corresponding minimum length scale [these might be specified by: the maximum significant wavevector in χ(

**k**

_{2D},

*z*); the resolution of the physical imaging system; or at most by

*k*

_{0}].

We obtain the BL law by assuming that the system is thin along **e**_{z}, i.e., both δ*z* ≪ *d*_{dof}, and

**r**) ∝

*i*σ

_{0}ρ(

**r**), this gives the usual BL law

describing the attenuation of the free space optical intensity *I*(**r**) = *c*ε_{0}|*E*(**r**)|^{2}/2 by absorbers of density ρ(**r**) and scattering cross-section σ_{0}. This BL result can also be obtained without the paraxial approximation by first neglecting the

*k*

_{max}δ

*z*≪ 1: a more strict requirement than in the paraxial approximation where we had δ

*z*≪

*d*

_{dof}) and again assuming |χ(

**r**)| ≪ 1, a small relative susceptibility.

^{22}

In experiment, the BL law is generally applied by comparing the intensities *I*(**r**_{2D}) and *I*_{0}(**r**_{2D}) measured with and without atoms present, respectively. This relates the optical depth

to the 2D column density. In cold atom experiments, this column density is the primary observable in experiment.

### C. Absorption imaging

Here we consider systems of ultracold atoms illuminated by laser light on a cycling transition, where the atom-light interaction is described by a complex relative susceptibility

ρ(**r**) is the atomic density; δ is the laser's detuning from atomic resonance;

*I*

_{sat}is the saturation intensity.

^{14}The standard BL law is valid for dilute (

*k*

_{0}δ

*z*≪ 1), illuminated by low intensity (

*I*

_{0}≪

*I*

_{sat}) probe beams.

The *I*_{0} ≪ *I*_{sat} requirement can be lifted by introducing the intensity-corrected optical depth

which is related to the column density

of dilute (

*k*

_{0}δ

*z*≪ 1). Due to the limited dynamic range of the camera's pixels

^{19}and the presence of background light, it is technically difficult to reliably detect uncorrected optical depths, larger than ≈4. Thus, we deliberately select

*I*

_{0}>

*I*

_{sat}, saturating the transition with

*I*

_{0}such that OD

_{cor}< 3.

In addition, the spatial thickness of many cold atom systems exceed the depth of field leaving parts of its distribution along imaging direction inevitably out of focus, thereby invalidating Eq. (12). Even for dilute clouds (after sufficient TOF), images taken an equal distance above and below the focal plane can differ. This lack of symmetry makes a straightforward determination of the optimal focus difficult (lensing effects from even slightly off-resonance imaging beams and aberrations in the imaging system can complicate the situation further.)

To illustrate this difficulty, we consider images of BECs with the focal plane displaced a distance *d* = −54 μm, 0 μm, and 54 μm from the BECs' center (see Fig. 1). Because the BEC is thick compared to the depth of field, Eq. (12) does not hold; in addition lensing effects cause the cloud's peak

### D. Modeling

To obtain a basic understanding of our approach, we first consider the defocused image of a 1 μm thick absorbing medium, inhomogenous **e**_{x}-**e**_{y} plane, bounded above and below by vacuum, with, χ(**r**) = *ig*(*x*, *y*) for *z* ∈ (−0.5 μm, 0.5 μm), where *g*(*x*, *y*) ⩾ 0 is a Poisson distributed random variable. Like atoms illuminated on resonance, this medium has a purely imaginary susceptibility. The illuminating light is modeled by a plane wave with wavelength λ = 780 nm suitable for imaging our ^{87}Rb Bose-Einstein condensates.^{24} While this object has no visible structure, by virtue of its spectrally flat density-density correlation function, it can be brought into focus.

The imaged intensity pattern *I*(*x*, *y*) from this 1 μm layer appears random at various distances from focus, but its correlations become oscillatory. To reveal this information, we turn to its spatial power spectral density: the magnitude squared of *I*(*x*, *y*)'s Fourier transform.^{25} The power spectral density (PSD) is circularly symmetric in the spatial frequency **k**_{2D} = (*k*_{x}, *k*_{y}) plane. Fig. 2(a) shows the PSD in this *k* = |**k**_{2D}| “radial” direction as a function of distance from focus *d*. This PSD has a fringe pattern; the wavevector of the first minimum exceeds the maximum imaged wavevector only near the image's focus at *d* = 0 μm.

The physical origin of this structure can be understood by turning to the paraxial wave equations [Eqs. (7) and (8)], and by first studying a single absorber at **r** = 0 illuminated by a plane wave

**e**

_{x}-

**e**

_{y}plane with width

*w*

_{0}. Thus the electric field just following the absorber is changed by

*r*

^{2}=

*x*

^{2}+

*y*

^{2}. The propagation of such a gaussian mode by a distance

*d*along

**e**

_{z}can be solved exactly in the paraxial approximation, and in the spectral basis this is

The total field from an absorber located at a different location **r**_{0} in the **e**_{x}-**e**_{y} plane simply acquires an overall phase factor exp [−*i***k**_{2D} · **r**_{0}]. We now compute the experimentally relevant optical depth by taking the reverse Fourier transform of the full electric field, computing the intensity, then the optical depth, and taking the Fourier transform to obtain (retaining terms of order δ*E*/*E*_{0})

with the same overall phase factor depending on the initial position. Averaging over *N* randomly placed absorbers therefore gives an overall signal scaling as

with

This quantity has zeros located at

*n*. In our numerical simulation, the minima follow the functional form

*k*

_{zero}[

*n*] =

*A*[

*n*]|

*d*|

^{−1/2}as shown by the dotted lines in Fig. 2, with

*A*[0] ≈ 5.08 and

*A*[1] ≈ 8.76 for the first and second zeros: the expected values for

*A*[

*n*]. Thus, for this thin apparently structureless system, fringes in the PSD allow us to identify the focal plane.

To demonstrate the technique of finding optimal focus of an extended object, we now consider a second disordered scattering potential with a columnar structure, now 100 μm thick, i.e., χ(**r**) = *ig*(*x*, *y*) for *z* ∈ (−50 μm, 50 μm), where again *g*(*x*, *y*) ⩾ 0 is a Poisson distributed random variable. This object's PSD is plotted as a function of distance *d* from its center in Fig. 2(b); in addition to the same fringe pattern as for the 1 μm thick case, the PSD now vanishes at specific spatial frequencies independent of *d*. To model this, we note that the absorbers can now be located at a distance *z* from the symmetry plane, so in Eq. (13), we replace *d* → *d* − *z* and integrate *z* from −δ*z*/2 to δ*z*/2, which ultimately gives the PSD

This predicts the appearance of additional zeros located at

*m*(this is an artifact of the box-like density distribution of atoms, and would be greatly softened in real systems where the density drops smoothly to zero). In our example, the lowest order horizontal fringes is located at

*d*= 0 μm, from the diverging curved-fringes.

Next, we consider a scattering potential fully disordered in 3D, again with a 100 μm thickness, i.e., χ(**r**) = *ig*(*x*, *y*, *z*) for *z* ∈ (−50 μm, 50 μm), where *g*(*x*, *y*, *z*) ⩾ 0 is a Poisson distributed random variable. In this case, the independent random scatterers along imaging direction causes the PSD to rapidly loose structure with increasing *k*_{2D} (see Fig. 2(c)). Here too, our random scatter model can be applied, giving

This reduces to our earlier result when δ*z* → 0 for a thin system and shows that, while the same fringes exist, they are rapidly attenuated for larger spatial frequencies, where the signal approaches a constant background value. However, in principle the curved-fringes still allow the optimal focus to be identified.

## III. OPTIMAL FOCUSING OF ELONGATED BOSE-EINSTEIN CONDENSATES

Using on our model, we now consider absorption imaged BECs and implement the technique presented in Sec. II to find the optimal focus.

We prepared *N* = 7 × 10^{5} atom ^{87}Rb Bose-Einstein condensates in the |5S_{1/2}, *F* = 1, *m*_{F} = 0⟩ electronic ground state in a crossed-dipole trap with frequencies ω_{x, y, z} = 2π × (3.1, 135, 135) Hz. *In situ*, the BECs were javelin shaped owing to the extremely anisotropic confining potential. After a 17 ms to 21 ms TOF, we repumped into the *f* = 2 manifold, and resonantly imaged on the |5S_{1/2}, *f* = 2, *m*_{F} = 2⟩ to |5P_{3/2}, *f* = 3, *m*_{F} = 3⟩ cycling transition with a λ ≈ 780.2 nm probe laser.

The imaging system consisted of a CCD camera and two pairs of lenses functioning as a compound microscope, magnifying the intensity pattern at the object by a factor of ≈6 at the image plane. The first pair of objective lenses, with effective focal length (efl)

*d*

_{dof}= 18.6 μm depth of field in our imaging system.

^{20}

Instead of varying the distance from focus by physically moving imaging lenses or the CCD, we changed the time during which the BEC fell along **e**_{z} and obtained absorption images with TOF times *t*_{TOF} from 17.0 ms to 21.0 ms. At these TOFs, the condensates' radii were *R*_{y, z} ≈ 75(5) μm and *R*_{x} ≈ 210(10) μm. Initially, the cloud was elongated in the harmonic trap with aspect ratio 43 : 1. The initial 43 : 1 aspect ratio was reduced to 2.65 : 1 after TOF, and the transverse size of the cloud exceeded the imaging depth of field by a factor of 8.

Figure 3(b) shows the 1D PSD of the atoms' corrected optical depths along **e**_{z}, which is directly related to the absorption intensity through Eq. (11). The fluctuations in the BEC's density distribution behave like the randomly modulated χ(**r**) in our example systems, creating a recurring fringe pattern in the PSD spectrum as obtained in Fig. 3(c). The fringes are quite pronounced for quasi one-dimensional BECs, where initial phase fluctuations map into pancake-shaped density fluctuations arrayed along the initially long axis after TOF.^{21} Despite the decreased contrast at high spatial frequencies due to the BEC extent along **e**_{z}, we clearly observe fringes curving as a function of *t*_{TOF} in Fig. 3(c). This allows us to determine the optimal focus of the system.

From the above experimental data, we fit the two lowest order fringes to *k*_{m}[(*d* − *z*_{0})^{2}/δ*z*^{2} + 1]^{−1/4}, a peaked function with the expected *d*^{−1/2} behavior away from *z*_{0}. The fits give an optimal focus location of *z*_{0}[0] = 1836(2) μm using the zeroth order fringe or of *z*_{0}[1] = 1837(2) μm using the first order fringe. These values correspond to a TOF of 19.36(1) ms. We are thus able to determine the optimal focus within ≈2 μm or equivalently ≈10 μs in TOF. Comparing the experimental data to the theoretical forms, we notice that the fringes are slightly asymmetrical with their locations slightly below theoretical ones for larger TOF. Based on our simulations, this likely results from the *z* dependent magnification of our imaging system, which changes by about 10% as the atoms fall from 1420 μm to 2150 μm (17 ms to 21 ms TOF).

## IV. SUMMARY

We presented a systematic method to bring clouds of ultracold atoms, particularly initially elongated BECs, into an optimal focus. The density fluctuations in the BECs after TOF acted like random scatterers, creating diffraction pattern which changed predictably as a function of distance from the optimal focus. Using TOF absorption imaging, we demonstrated this method, pinpointing the optimal focus of the BEC to within 2 μm for a 150 μm thick BEC. This robust technique is easily implemented, requires no hardware changes, and uses a minimum of computation.

## ACKNOWLEDGMENTS

We thank F. E. Becerra, A. Hu, and W. D. Phillips for a careful reading of the article. We acknowledge the financial support from the NSF through the Physics Frontier Center at JQI, and the ARO with funds from both the Atomtronics MURI and DARPA's OLE Program.

## REFERENCES

**r**) in Eq. (1) cannot be safely neglected for systems where |χ(

**r**)| is large or sufficiently rapidly varying, although this would generally imply a breakdown of the paraxial approximation as well.

*z*= 1 μm step size.