Optimally focused cold atom systems obtained using density-density correlations

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,⁷ 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.¹⁴ 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) = ∫dz ρ(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⁴ 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.

Monochromatic light of free-space wavelength λ and wavenumber k₀ = 2π/λ propagating through an object with complex relative permittivity ε(r) = ε/ε₀ and relative susceptibility χ(r) = ε(r) − 1, where ε is the permittivity, and ε₀ 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.

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¹⁵ 

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

$Ẽ(k2D,z)=∫d2r2Dexp(−ik2D·r2D)E(r)$⁠; and the transfer function for propagating a distance Δz in free space

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

$k2D2<k02$⁠,

$\tilde{\bf P}$

$P̃$ is oscillatory (propagating regime), and for

$k^2_{\rm 2D}>k_0^{2}$

$k2D2>k02$⁠, it is exponentially decaying (evanescent regime).

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.

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).¹⁶ 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₀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

$∂z2E′(r)$⁠. Like above, the spatial evolution of an initial E^′(r) can be partitioned into a spectral part

$\tilde{\bf P}^\prime ({\bf k}_{\rm 2D},\Delta z)$

$P̃′(k2D,Δz)$ and a coordinate 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³.

The paraxial equations allow us to introduce the depth of field

where

$k max 2$ is the largest 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₀].

We obtain the BL law by assuming that the system is thin along e_z, i.e., both δz ≪ d_dof, and

$P̃′(k2D,Δz)$ may be neglected. For purely absorbing materials where χ(r) ∝ iσ₀ρ(r), this gives the usual BL law

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

$∇⊥2$ term in Eq. (2) (valid when 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.²² 

In experiment, the BL law is generally applied by comparing the intensities I(r_2D) and I₀(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.

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;

$σ0=6π/k02$ is the resonant scattering cross-section; Γ is the atomic linewidth; and I_sat is the saturation intensity.¹⁴ The standard BL law is valid for dilute (⁠

$\rho \ll k_0^3$

$ρ≪k03$⁠, see Ref. 23), spatially thin systems (k₀δz ≪ 1), illuminated by low intensity (I₀ ≪ I_sat) probe beams.

The I₀ ≪ I_sat requirement can be lifted by introducing the intensity-corrected optical depth

which is related to the column density

of dilute (⁠

$ρ≪k03$⁠), spatially thin systems (k₀δz ≪ 1). Due to the limited dynamic range of the camera's pixels¹⁹ and the presence of background light, it is technically difficult to reliably detect uncorrected optical depths, larger than ≈4. Thus, we deliberately select I₀ > I_sat, saturating the transition with I₀ 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

$OD cor$ to behave asymmetrically when the focus is behind or in front of the cloud. In these images, there are no sharp features that identify the optimal focus at the micron level. Owing to the weak dependence of large-scale parameters such as peak-height or width on defocus, such precise focusing is not required in many experiments. As we see below, experiments that study correlations within such images are extremely sensitive to defocus and new methods are required. Our technique brings images such as these into focus, identifying an optimal focal plane at the ≈2 μm level.

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 ⁸⁷Rb Bose-Einstein condensates.²⁴ 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.²⁵ 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

$E0′(r2D,0−)=E0$⁠. Equation (8) shows that a thin absorber simply changes the amplitude of the field, leaving its phase untouched, and for simplicity, we assume this absorber has a Gaussian profile in the e_x-e_y plane with width w₀. Thus the electric field just following the absorber is changed by

$\delta E^\prime ({\bf r}_{\rm 2D},0^+) = -\delta E \exp [ -r^2/w_0^2]$

$δE′(r2D,0+)=−δEexp[−r2/w02]$⁠, with r² = x² + y². 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₀ in the e_x-e_y plane simply acquires an overall phase factor exp [−ik_2D · r₀]. 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₀)

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

$N$ with a random overall phase. Taking the magnitude squared gives the PSD

with

This quantity has zeros located at

$k zero [n]=2π(n+1/2)k0/d$ for integer 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

$k zero ′[m]=4πmk0/δz$ for non-zero integer 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

$k^\prime _{\rm zero}[1]=1.00\,\mu {\rm m}^{-1}$

$k zero ′[1]=1.00μm−1$⁠. Here again, we easily determine the optimal focus, 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.

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⁵ atom ⁸⁷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)

$f1̂=53.6 mm$⁠, collimated the light diffracted by the cloud and were separated by a distance

$D\break=\hat{f_1}+\hat{f_2}$

$D=f1̂+f2̂$ from a second pair of lenses with a

$\hat{f_2}=325\,{\rm mm}$

$f2̂=325 mm$ efl. The resulting 0.23 numerical aperture implies that a 10.6 μm diffraction-limited spot on our CCD sensor is larger than its 5.6 μm pixel size. The associated 1.7 μm spot-size on the cloud gives a d_dof = 18.6 μm depth of field in our imaging system.²⁰ 

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.²¹ 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₀)²/δz² + 1]^−1/4, a peaked function with the expected d^−1/2 behavior away from z₀. The fits give an optimal focus location of z₀[0] = 1836(2) μm using the zeroth order fringe or of z₀[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).

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.

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.