Spectroscopic localization of atomic sample plane for precise digital holography

In digital holography, the coherent scattered light fields can be reconstructed volumetrically. By refocusing the fields to the sample planes, absorption and phase-shift profiles of sparsely distributed samples can be simultaneously inferred in 3D. This holographic advantage is highly useful for spectroscopic imaging of cold atomic samples. However, unlike (e.g., biological samples or solid particles), the quasi-thermal atomic gases under laser-cooling are typically featureless without sharp boundaries, invalidating a class of standard numerical refocusing methods. Here, we extend the refocusing protocol based on the Gouy phase anomaly for small phase objects to free atomic samples. With a prior knowledge on a coherent spectral phase angle relation for cold atoms that is robust against probe condition variations, an ``out-of-phase'' response of the atomic sample can be reliably identified, which flips the sign during the numeric back-propagation across the sample plane to serve as the refocus criterion. Experimentally, we determine the sample plane of a laser-cooled $^{39}$K gas released from a microscopic dipole trap, with a $\delta z\approx 1~{\rm \mu m}$$\ll 2\lambda_p/{\rm NA}^2$ axial resolution, with a NA=0.3 holographic microscope at $\lambda_p=770~$nm probe wavelength.


Introduction
In a generic absorption imaging setup, the optical forward scattering from the sample under study is imaged together with the co-propagating probe light onto the imaging sensor arrays. The attenuation of the total intensity = | + | 2 records the in-phase component of relative to . Information on the out-of-phase component is lost. Similarly, in phase contrast imaging setups [1,2] where is phase-shifted by /2, the information loss occurs to the in-phase quadrature. Digital holography [3] recovers the full information by reconstructing the phase of relative to using holograms. For the case of inline holography [4,5], the holograms are simply out-of-focus interference fringes between and . With the full wavefront knowledge at hand, both the and fields can be volumetrically reconstructed around the sample planes via digital back-propagation. Furthermore, with sufficient knowledge of the samples, the reconstruction support self-consistent characterization of sparse samples for precise 3D microscopy [6][7][8]. During the process, to refocus each reconstructed sample image to its respective plane [9][10][11][12][13][14][15][16][17] is crucially important. For the purpose, various refocus schemes are developed based on priori knowledge of the samples and their interaction with light. Examples include the methods based on the edge sharpness and sparsity [9,12,14,15], the Gouy phase shift [13], by requiring imaging consistencies under multiple wavelength [10] and structured illumination [11], or even by deep learning of complex features [16,17].
The holographic advantages associated with 3D complex imaging of sparse samples can be highly useful for applications across fields [4,[18][19][20][21][22]. For atomic physics research, over the years efforts have been made for holographic imaging of cold atoms [23][24][25][26][27]. In a recent work, we show that an improved holographic technique with suppressed aberration and speckle noises supports simultaneous retrieval of atomic absorption and phase shift profiles with diffraction-limited spatial resolution and photon shot-noise limited sensitivity [28]. As illustrated in Fig. 1a here, the technique uses a precisely pre-characterized probe wavefront to recover the coherent atomic forward scattering with the hologram data ( Fig. 1a(ii)). Then, both and are numerically propagated from the camera plane = back to the sample plane = where the 2D optical depth OD( , ) = −2Re[log(1 + / )] and phase shift ( , ) = Im[log(1 + / )] (Fig. 1a(i)) are evaluated. Here, similar to the applications in other fields [9][10][11][12][13][14][15][16][17], to localize the atomic plane is crucially important for faithfully retrieving the generic atomic absorption and phase shift properties. However, unlike typical biological or solid samples with sharp boundaries, cold atoms in optical traps typically follow quasi-thermal distributions [29], without much distinct features as a priori criterion to perfect the sample-plane refocus. Nevertheless, in early efforts for holographic imaging of cold atoms [23][24][25], the sample-plane localization is still largely based on optimizing certain characteristic spatial features of atomic distribution, with moderate accuracies.
In fact, while cold atomic samples usually lack distinctive spatial structures, there are unique features constrained by atomic physics available for calibrating coherent imaging. For example, in aberration-free in-focus imaging, the power spectrum density of atomic density correlations is expected to be flat for non-correlated atoms [31]. The criterion is applied in ref. [27] to achieve precise refocus of aberrated phase-contrast atomic images. Other than exploiting spatial structures or correlations, in ref. [26] the authors suggest that for far off-resonant imaging (i.e., in Fig. 1b the probe detuning Δ is much larger than the atomic transition linewidth Γ), that atomic samples appear as phase objects with | ( , )| |OD( , )| across the sample becomes a refocus criterion to precisely locate the sample plane. This idea shares the same underlying physics with the refocusing method based on Gouy phase anomaly for 3D localization of small transparent particles [13]. As to be detailed shortly, relative to the uniform probe wavefront , the propagation of spatially confined picks up an extra Gouy phase G in its center in the far field. Therefore, the known relative phase relation between and for phase objects holds sensitively in the near field, naturally serving as a refocus criterion for locating the sample plane. Obviously, this diffraction phase criterion can be generalized for localizing atomic samples probed at arbitrary detuning Δ, if the phase angle 0 = arg( + OD/2) can be measured and compared with known values precisely.
In this work, we show that a recently demonstrated phase-angle spectroscopy [28] leads to a robust criterion for locating the sample plane in holographic microscopy of cold atoms, with an achievable axial resolution well below the diffraction limit. The prior spectral phase angle knowledge exploited in this method, the Eq. (11) relation to be discussed shortly, is easily understood in the linear optics regime where both OD and for a dilute gas can be evaluated analytically [24,32]. Practically, to achieve sufficient signals with short and nearly resonant exposures, saturation of atomic transitions can hardly be completely avoided. It is also known that the linear optical response of cold gases is prone to resonant dipole interactions [33][34][35]. Interestingly, we find that far beyond the linear response regime, the Eq. (11) relation can hold precisely for isolated atomic transitions probed with smooth and long pulses. Our method thus supports robust localization of atomic sample plane with flexible probe condition managements, for working with denser samples [34] and to achieve photon-shot-noise limited performances [25,28]. Experimentally (Fig. 1), we demonstrate the spectral refocus method by repeatedly probing an open hyperfine transition of 39 K D1 line with = 1 s pulses, to locate the atomic sample plane with sub-micron axial resolution.
Previously, the best refocus criterion for imaging cold atomic samples appears to be that based on the atomic shot-noise correlations [27]. In comparison, our method provides similar accuracy with a much stronger signal for rapid applications. In addition, instead of relying on regularizing contrast transfer functions [24,27,36], our holographic method directly supports a large depth of view, with diffraction-limited resolution [25], for future 3D spectroscopic imaging of sparse   The imaging optics is with magnification = 1 and a numerical aperture NA = 0.3. The (a,i) subplots are typical reconstructed optical depth OD and phase shift images at the = plane (experimental data, with peak OD≈ 0.2 and ≈ −0.15). The (a,ii) subplot gives the corresponding reduced hologram recorded by the camera at = (peak-to-peak ≈ 100 in terms of counts, with photon shot noise at the 20-level [28].). The scale bars are with 1 = 30 m and 2 = 1 mm respectively. (b): A 4-level diagram to represent the light-atom interaction at the 39 K D1 line ( = 770 nm) nearly resonant to the | (4S 1/2 , = 1) and | (4P 1/2 , = 2) transition. (c): Simulated OD and profiles associated with the D1 interaction under realistic experimental condition according to optical Bloch equations [30]. As detailed in Sec. 2.3, the spectroscopic phase angle 0 within the shadowed |Δ| ≤ Γ follows Eq. atomic samples.

Measurement principles
As schematically illustrated in the Fig. 1 setup, we consider holographic imaging of an atomic sample subjected to a spherical probe light illumination at wavelength . The atomic sample is centered at r = (0, 0, ) with spatial width { , } , so that the probe light propagates along through the sample with negligible wavefront curvature itself. We assume thin atomic samples. The length along the light propagation direction satisfies (1) when compared to the Rayleigh distance associated with the smallest spatial feature of interest of the sample characterized by an effective Gaussian width . The spherical in this work is derived from a defocused Gaussian beam. The spherical wave illumination [3] enhances the pixel-resolution and dynamic range of the camera sensors during holographic imaging [24,25]. Our method can be straightforwardly generalized to plane-wave illumination, as well as structured, complex illuminations.
The total field tot = + after interacting with the thin sample can be expressed as tot = . The complex phase shift [28] ( , , is related to the optical depth and phase shift as ( ) = + OD/2. Notice here and in the following we may omit the ( , ) variable in the 3D distribution functions if no ambiguity is induced. For a dilute sample of free atoms, the complex numbers are rotated from the real axis by the phase angle 0 = arg( ( )) according to the single-atom response. We re-write the complex phase during the backward propagation as Clearly, the out-of-phase component ⊥ vanishes at = , just like that the attenuation by thin phase objects is zeroed in-focus [13,26]. The effect is illustrated in Fig. 1d(i-ii) for the D1 line of 39 K in this work as to be detailed shortly, but is general for arbitrary atomic transitions as long as 0 is known for evaluating ⊥ with Eq. (5).
To simplify the following discussion on the propagation effect, we now model the atomic sample by the Gaussian profile with = . Propagating away from = , the defocused near the imaging center (with | | ≤ in Fig. 1c) picks up an additional phase, the Gouy phase G = −arctan − (Fig. 1d(iii)), relative to [13]. within a region of interest (ROI) defined by | | ≤ , we expect As to be detailed in the following, the Eqs. (3-6) relation can be exploited to locate the = plane with a given , data set, by minimizing the ⊥ components according to a prior 0 -knowledge.

Holographic , reconstruction
To experimentally evaluate ( , , ) according to Eq. (3), we need to volumetrically reconstruct , fields from experimental data first. The procedures to infer the probe wavefront from the pre-experimental characterizations, and from the single-shot atomic sample holograms, are detailed in our previous work [28]. Briefly, the probe wavefront at the camera sensor plane ( ) = √︁ ( ) ( ) is obtained first by a multi-plane Gerchberg-Saxton algorithm [37] with multiple { ( )} probe intensity measurements as inputs, using a −translating camera. The camera position is then fixed at = for experimentally recording the holograms with and without the cold atomic samples under study. With careful numerical adjustments and subtractions, the digital images are reduced to represent = | + | 2 and 0 = | | 2 respectively. An iterative twin-image removal algorithm [25] is then applied to retrieve from the reduced hologram = * + * + | | 2 ( Fig. 1a(ii)). With the full wavefront knowledge for and at hand, we numerically propagate both fields from the camera plane to locations around the sample plane, via the angular spectrum method [28,38], Hereˆ,ˆ− 1 represents the 2D Fourier transform and the inverse transform respectively: ( , , ) =ˆ( , , ) and ( , , ) =ˆ− 1 ( , , ). = 2 / is the wavenumber of the probe light.

A robust spectroscopic phase angle relation for cold atoms
Clearly, to construct a = refocus criterion with Eqs. (4-6), prior knowledge of the phase angle 0 is essential. With the thin sample condition by Eq. (1), the complex phase shift at = for a dilute sample is expected to follow the Beer-Lambert law as [28] Here is atomic density distribution and is the complex atomic polarizability. We therefore expect 0 = arg( ), which is precisely known in the linear optics regime [32]. However, as suggested in the Introduction, the validity of the linear analysis is practically prone to saturation and resonant dipole interaction effects.
Beyond the linear analysis, here we generally consider a sample of multi-level atoms interacting with a probe pulse (Fig. 1b). We assume the probe is strong enough so that inter-atomic resonant dipole interaction can be ignored [34]. In light of the fact that the holographic data recorded by the camera is given by = ∫ 0 ( ) , the effective polarizability for Eq. (8) is evaluated as Furthermore, when the probe frequency is nearly resonant to an isolated | − | transition, the induced complex dipole moment can be approximated as Here d is the dipole matrix element of the | − | transition. The coherence ( ) obeys the multi-level master equation for the atomic density matrix [39], as detailed in Appendix A.
With the key Eqs. (9)(10) assumptions, in Appendix A we show that for a smooth pulse with duration 1/Γ at detuning |Δ| ≤ Γ , 0 is approximated by with very good accuracy. For example, for the D1 line probed in this work as to be detailed in Sec. 3 (Fig. 1c), deviation from Eq. (11) within |Δ| ≤ Γ is less than 1% in terms of minimal ⊥ /| | by Eq. (5). Here = − Γ Γ+ and = − Γ+ parametrize the average Stark shift and optical pumping rate induced by all the | − | and | − | couplings that off-resonantly mix the atomic states during the time. The spectral isolation here requires the transition frequency to be far away from these ground-or excited-state sharing transitions, as well as all the other | − | transitions: The probe pulse can be fairly short and strong, as long as these other transitions are only excited perturbatively. The robustness of the phase-angle relation by Eq. (11) makes it particularly convenient to constrain the sample plane of dilute atomic gases in digital holography.

Spectroscopic sample plane localization
Experimentally, as to be detailed in Sec. 3, a set of holograms for a standard, free-space sample of dilute atoms is recorded in repeated preparation-measurement cycles at various detuning {Δ } with |Δ | ≤ Γ. The complex phase shifts { } are then evaluated according to Eq. (3). To locate the atomic sample plane, we numerically propagate , ( ) according to Eq. (7) to minimize the cost function which implicitly depends on the lineshape parameters , through Eq. (11), as well as the sample plane through Eq. (6). Each out-of-phase component ⊥ ( ) by Eq. (5) is averaged within the ROI defined by | | < (Fig. 1d). With the -minimization, we expect a linear relation between OD and within |Δ| ≤ Γ: To understand the accuracy of the −minimization by Eq. (12), we notice the ⊥ data entering the analysis is fundamentally limited by the uncertainty of the phase angle ( ) = arg( ( )) retrieved from the hologram data. In particular, with with , ∝ ROI | ( ) | 2 summing the elastically scattered photons from the ROI to the camera.
The photon shot noise affects the predictions to all the , , parameters. However, with Δ uniformly sampling |Δ| ≤ Γ so that , can typically be fixed quite precisely, we may simplify the analysis by ignoring the correlations so that = 0 suggests ⊥ + ⊥ = 0. Together with Eq. (6), we arrive at photon-shot-noise limited axial resolution Here = , is the number of all the elastically scattered photons entering the data analysis. The is a sample shape dependent factor, with = 1 for the Gaussian shaped sample.
We now discuss the choice of ROI for the coherent signal averaging in Eqs. (12) (13). As shown numerically in Fig. 1d(ii) (Appendix B), the rapid sign inversion for ⊥ along the light propagation direction is most pronounced near the center with | | < in the plot. With a large Δ = − to be comparable to (Eq. (2)), oscillatory ⊥ ( , ) is developed along due to the curvature mismatch between and . Therefore, for the -sized atomic sample, | | ≤ is a natural ROI choice to isolate the uniform ⊥ ( ) center for the ⊥ ( ) average. A too large ROI would reduce the -dependence of ⊥ ( ) comparing to Eqs. (5)(6). Conversely, a too small ROI would reduce the total number of elastically scattered photons entering the data analysis, compromising the photon shot noise limit (Eq. (15)). Practically, the atomic samples cannot always be approximated by Gaussian profiles. The size may not always be assumed as prior knowledge either. For the small samples to be discussed in this work (Fig. 1a), the ROI should in principle be refined toward | − 0 | < with a suitable central position 0 and width to balance the −sensitivity of the -function in Eq. (12) with the amount of ROI-photons . Practically, with the sample plane emphatically determined, we find that simply by thresholding the approximately refocused | | 2 intensity, e.g. ROI=1 for | | 2 > | | 2 max , ∼ 0.1, then nearly optimal sensitivity for the Eq. (12) minimization can be achieved.

Correcting high-order aberrations
From Eqs. (12)(15), the high quality atomic sample plane localization relies on high quality minimization of the −function toward the photon shot noise limit. The fit quality is affected by the imperfect optical system itself. For lensless holographic imaging [25], the optical transfer function according to free-space propagation is easily modeled. However, in most cold atom experiments, the imaging system usually requires an optical train to relay the coherent wavefronts from the samples in vacuum to the camera outside the vacuum [8,27,28], as schematically illustrated in Fig. 1a. Even for perfect optics, high-order aberration correction is required for volumetric imaging across a large depth of view.
Here, the high-order aberration correction can be achieved simultaneously with the = sample plane localization, through the minimization of −function by Eq. (6). For example, this can be achieved by setting up numerical Zernike plates [27] at = with coefficients entering the −function for the optimization. This is a step left for future work.

Methods
Our experiment demonstration is based on a 39 K holographic microscope on the D1 line. To facilitate the following discussions, in Fig. 1b we refer the 4 1/2 , = 1, 2 hyperfine ground states as | , | , and 4 1/2 , = 1, 2 excited states as | , | respectively [40]. As schematically illustrated in Fig. 1a, an NA = 0.3 optical train with magnification = 1 relays the probe wavefront and the forward scattering by the cold atomic sample to the digital camera. The probe wavelength = 770 nm is nearly resonant to the | − | hyperfine transition with a natural linewidth Γ = 2 × 5.96 MHz. The transition is spectrally isolated from | − | and | − | transitions by , = 2 × 55.5 MHz and , = 2 × 461.7 MHz respectively. Up to 10 3 atoms are laser-cooled to a temperature of tens of micro-kelvin [30] and loaded into a microscopic optical dipole trap (ODT) composed by a focused = 780 nm laser along . The approximately Gaussian-shaped atomic sample is with = 15 m along . The width , estimated to be less than 1 m is below the ( ) res = /NA = 2.6 diffraction limit [41] of the NA=0.3 holographic microscope [42]. In absence of imaging aberrations, we expect the diffraction-limited sample images to have an apparent Gaussian width ≈ ( ) res /2 √ 2 along   2)) is close to the diffraction-limited imaging axial resolution ( ) depth = 2 /NA 2 .
To spectroscopically locate the atomic sample plane, a sequence of two images ( ) 1,2 are recorded at each probe detuning Δ with and without the atomic sample in repeated measurements. The CCD camera with 1040 × 1392 pixels, each 6.45 × 6.45 m 2 in size, is effectively placed at = 10.4 mm to record holograms of the atomic sample at = 0.7 mm. The camera exposure is set as = 1 ms. As detailed in Appendix C, to avoid inhomogeneous light shifts that tend to invalidate Eq. (13), these "standard samples" are released from ODT before the holographic imaging. Notice the | − | transition is open: during the probe excitation, spontaneous | → | or | → | Raman scattering tends to quench the atomic population into the dark | . Therefore, instead of probing the atoms continuously, we interleave a train of = 1 s probe pulses with 4 s of trapping+cooling pulses composed of the ODT beam with a D1 molasses [30] blue detuned from | − | transition, which not only help to maintain the samples' shape, location, and temperature, but also depump the internal states back to | for the nearly resonant imaging. The probe pulse is set at a moderate ≈ 3 mW/cm 2 intensity in this work so that the near-resonant OD and are significant for the microscopic and dilute samples. Since the completion of this work, we have verified that our method works well at higher intensity, both numerically (Appendix A) and experimentally [28].
We follow the procedure outlined in Sec. 2.2 to reconstruct ( ) , at each probe detuning Δ , using the ( ) 1,2 hologram data taken from repeated measurements. Depending on the desired signal to noise ratio, ave = 2 ∼ 15 holograms obtained at a same measurement condition are averaged to ( ) 1,2 before proceeding the ( ) reconstructions. We then propagate both fields according to Eq. (7) to retrieve ( ) across the sample plane , for the -function minimization (Eq. (12)). The performance of the spectroscopic refocus method is evaluated by repeating the measurements, typically within an hour of measurement time, to check the consistency of the predicted = values. In addition, the results are compared with the more traditional method based on minimizing the apparent sample width, given by Here , 2 are the 2D average of , 2 in the − plane, weighted by the 2D complex phase magnitude | ( , , )| within a large enough ROI , during the numerical −scan of , . Experimentally we retrieve ( ) by fitting − averaged | ( , , )| with 1D Gaussian profiles. The factor of 4 is chosen so that = 2 √ 2 for the Gaussian beam model.

Results
We first present typical 2D ( ) profiles in Fig. 2 around the sample location = . The probe detuning is chosen as Δ = Γ in this example with substantial optical peak depth OD ≈ 0.2 and peak ≈ −0.15 respectively (see Fig. 1a(i)). To improve the display, the holographic data is averaged over ave = 14 images (see Fig. 1a(ii)). We rotate the complex with the known 0 (Δ) according to Eq. (4), which, according to the -minimization to be described shortly, is adjusted to be arctan(2 + ) with = −0.84, = 0.08 (fit from experiment data, see Table. 1 "average" column). In Fig. 2a we see the ⊥ almost vanishes at = while some weak fringes are still seen, due to the uncompensated high-order aberrations (Sec. 2.5). On the other hand, substantial ⊥ ( ) are developed at Δ = ±20 m. In Fig. 2d the ⊥ ( = 0, , ) similar to Fig. 1d(ii) is given, where we see the experimental data matched very well with the theoretical expectation. The complex phase shift ( ) is given in Fig. 2c with the color-domain plots. At the precisely refocused = , the ( ) becomes "monomorphous" with a uniform = arg( ( )) distribution [24], as expected.
Next, the highly -sensitive Eqs. (12)(13) criterion is illustrated in Fig. 3 with the ROI-averaged ( ) − Δ curves. Here we still have ave = 14. The 2D ( ) distribution around = as those in Fig. 2 are reproduced in the inset plots of Fig. 3(a,b,d,e). We evaluate ( ) as described in Sec. 2.4, within the ROI that are marked with dashed circles in the insets. The Fig. 3(a-c) data are according to the experimental geometry, but for linear response of ideal 2-level atoms ( = −1, = 0). The solid curves in Fig. 3(d-f) are instead numerically generated by adjusting the saturation of the 2-level model (Appendix B) together with = −0.84, = 0.08 parameters to fit the experimental data. In Fig. 3(c,f) we see the linearity of the phase angle is strongly impacted by the deviation from the sample plane by a distance as small as Δ ≈ ±10 m to be comparable to 2 , in agreement with Eq. (6).
Having introduced the general spectroscopic features of the reconstructed ( ) near = , we now present details of the -plane localization by minimizing the − function in Eq. (12). Specifically, for a set of 25 { , Δ } data with Δ scanning from -6 MHz to 6 MHz by 0.5 MHz steps, we globally minimize to obtain the opt , opt , parameters. We then plot the normalized cost function,  13)).
in Fig. 4a vs . To check the consistency of the localization, the same procedure is repeated with five data sets, with holograms in each set averaged by a moderate ave = 2 ∼ 3. The min ( ) curves are compared with the apparent width ( ) according to Eq. (16) evaluated with the same data set. The detailed numbers are given in Table 1. For a clear comparison, in both Fig. 4a and Table 1 the distance is evaluated relative to the average-from the five min −based predictions. Experimentally, from Fig. 4a we see min ( ) on the displayed vertical scale almost vanishes. Indeed, min ( ) ≈ 0.2% (Table 1) is less than a tenth of min ( ± Δ ) by a slight defocusing distance Δ = 5 m. For comparison, ( ) hardly changes by 30% by the same defocus. As by Line 2 of Table 1, the = localization in repeated measurements show remarkable consistency with an estimated standard deviation of = 0.3 m. For comparison, the sample width method in Fig. 4a (orange lines) and Line (5-6) of Table 1 show a substantially larger = 2 m. The much better performance by the spectroscopic method is associated with the aforementioned strong −dependence of the cost function min ( ), due to the Gouy phase anomaly (Eq. (6)), which makes the spectroscopic method substantially more resilient to imaging noises.
Ideally, according to Eq. (14), we expect min ( ) = 1/ in the photon-shot-noise limit. for the Fig. 4a and Table 1 results, the total number of elastically scattered photons in each data set of holograms at the 25 Δ detuning is estimated to be ∼ 10 6 , taking into account the ∼ 20% quantum efficiency of our camera (Pico Pixelfly). However, as in Fig. 2a our imaging system is not ideal so that the observed min ( ) ≈ 2 × 10 −3 1/ is instead limited by high-order aberrations. The observed min ( ) is also substantially larger than those caused by deviation from the Eq. (11) relation, due to transient and multi-level effects (Appendix A) which is numerically estimated at a 10 −5 level [43]. In other words, there is an intrinsic uncertainty to the phase angle 0 ≈ ( min ( )) 1/2 ≈ 0.05, due to the imaging system smearing itself. With min ( + Δ ) = min ( ) + Δ 2 /2 in quadratic form, we attribute the observed min ( ) as unreliable modeling that limits the absolute axial resolution to = √︃ min ( )/ ≈ 1.0 m. Finally, it is useful to remark that although the data in Table 1 suggests any optical drifts between the atomic sample and camera is small in repeated measurements, practically any drifts during the Δ -scan measurements effectively increase the sizes of the "average samples" entering the data analysis. In that case, our min method should still operates to find the average planes.

Discussions
The last twenty years witness rapid developments of quantitative imaging techniques in digital holographic microscopy, with applications across fields [4,[18][19][20][21][22]. In comparison, holographic imaging for atomic physics research has been underdeveloped. A list of unique technical challenges needs to be addressed [28,41], before the holographic method can be applied with sufficient accuracy for imaging the highly fragile ultra-cold atomic samples. This work aims to resolve a particular challenge: the precise localization of the sample plane for retrieving the generic optical response of the atoms. The difficulty arises from the fact that typical atomic samples are spatially featureless. Previously, the only effort to address the problem appears to be exploiting atomic shot-noise correlations in phase-contrast imaging [27]. In this work, instead of relying on spatial information to form refocus criterion, we propose to utilize characteristic spectroscopic features of atomic transitions for precise refocus in holographic microscopy. The underlying principle is to exploit the additional diffraction phase in the forward direction picked up by small objects, known as Gouy phase anomaly [13], that leads to deviation of apparent spectroscopic responses from those predicted by theory. The idea has already been demonstrated for localizing transparent objects [13,26]. We combine the diffraction phase idea with the unique ability of holographic microscopy for resolving the complex phase shift [28], and propose a spectroscopic criterion to robustly localize the atomic sample plane. The proposal not only utilizes the fact that for dilute, thin samples the spectral phase angle is insensitive to atomic density fluctuation [28], but also exploit an interesting phase-angle relation (Eq. (11)) which holds precisely for multi-level atom driven by strong optical pulses (Appendix A). Experimentally, this work demonstrates super-resolved sample plane localization during digital holography of a diffraction-limited, laser-cooled 39 K sample with sub-micron repeatability. This axial resolution is improved from the traditional method based on fitting the sample widths (Fig. 4a) by nearly an order of magnitude, in presence of imaging noises in our system. The absolute axial resolution of ≈ 1 m is yet limited by high-order aberrations of the imaging system itself (Sec. 3), that can be minimized in future work (Sec. 2.5). With the improvements, we expect the absolute axial resolution to reach the sub-micron level too, to be even smaller than the sample size itself. Our method can be applied to larger samples, where the atomic density fluctuations [27,31] lead to the required diffraction phase shifts. By properly choosing a set of ROIs (Eq. (12)) (Sec. 2.4), spectroscopic signatures of density-fluctuating features at various length scales of interest can be exploited to efficiently locate the central planes of the samples with the holographic microscope. Finally, it is important to note that in our experiment, the peak atomic density of about 10 13 cm −3 1/ 3 is quite dilute, while the peak optical depth OD max < 0.5 is still small. To exploit our method for localizing samples with higher OD and density, stronger pulses should help to suppress contribution of resonant dipole interactions that would otherwise modify the line shape [34,35] to compromise the Eq. (11) criterion.
With the precise knowledge of the sample plane location, the complex spectroscopy method in this work can be uniquely powerful for resolving the phase angle information of ultra-cold samples next [28]. In particular, to infer nontrivial, correlated optical responses of high OD, high density gases [34,35], the atomic sample plane can be spectroscopically located by probing a strong, isolated atomic transition with strong enough pulses first, as in this work. The phase-angle spectroscopy of the cooperative responses of the denser samples can then be reliably retrieved, in presence of density fluctuations generic to cold atomic samples which typically prevent regular imaging methods from obtaining accurate spectroscopic information in single-shots [44,45].
Our spectroscopic method can also be extended to locate multiple samples in digital holography, where the highly precise sample plane localization forms an excellent starting point for complex spectroscopic imaging [28,43] of sparsely distributed cold atomic samples in 3D [46,47].

Funding information
The authors acknowledge support from National Key Research Program of China (2022YFA1404200; 2017YFA0304204); NSFC under Grant No. 12074083; and the Original Research Initiative at Fudan University.

Disclosures
The authors declare no conflicts of interest.

Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Acknowledgments
We thank Dr. Liyang Qiu for helping on numerical calculations in this project.

A. Derivation of Eq. (11)
In the following we show that for spectrally isolated | − | transition as those in Fig. 1a and  [30,48]. The results [43] for = 1 s pulses suggest the Eq. (11) relation supports min < 10 −4 (to compare with min ( ) in Table (1)) for the open D1 transition of 39 K in this work, improves further for resolved D2 cooling transition lines such as those for 87 Rb [28], and becomes nearly perfect for 2-level atoms with a similar linewidth Γ. These are for weak pulses, as well as for intense pulses that strongly saturate the | − | transition.
As illustrated in Fig. 5, we generally consider the interaction between a probe pulse and a multi-level atom. The probe frequency is nearly resonant to the | − | transition with detuning Δ = − . In addition to the Ω coupling, off-resonant Ω , Ω couplings, as well as separate Ω couplings are all induced. All the couplings contribute to the induced dipole moment d and therefore the atomic polarizability ( ) to attenuate the intensity and shift the phase of . Here, with the small detuning Δ of order Γ and further with the aforementioned spectral isolation condition, the atomic response to is dominantly due to the 2-level coherence. Writing the density matrix element in the frame co-rotating with the probe field, we have the single-frequency dipole moment given by Eq. (10) in the main text. The equation of motion for is according to the master equation formalism [39] = Here the Rabi frequency is defined as Ω = ·d ℏ , with other transitions follow similarly. Importantly, the Ω and Ω couplings in the 2nd line of Eq. (18) are off-resonant. For "not too strong" pulses so that the level splittings still dominate the time scales, the dynamics of these coherences largely follow and can thus be adiabatically eliminated: We therefore rewrite Eq. (18) as Equation (20) is coupled to ( ), ( ) dynamics that we have not written out explicitly. However, for understanding 0 = arg( ), it suffices to integrate Eq. (20) to have HereΔ = (Δ + ) + (Γ + )/2. We then evaluate the nominator of Eq. (9) in the main text, Here Λ( ) = |Ω ( )|( ( ) − ( )). To arrive at Eqs. (22)(23) with smooth ( ), we have ignored the time-dependence of ( ), ( ) so that Δ = 0. Accordingly, we regard , as average values for Eq. (11). Furthermore, assuming smooth pulse with Λ ( ) (0) = 0, Eq. (23) can be evaluated with integration by parts, For the smooth pulse with 1/Γ, then the leading term dominates in the integration. Assuming > , has a phase angle 0 = −arg(1/Δ) that obeys Eq. (11). Here, it is interesting to note that for symmetric |Ω ( )| ≈ |Ω ( − )| pulse profiles, the second term Ω ( ) Λ ( ) Δ 2 tends to average itself out with the integration. The average effect works for both weak and strong, off and resonant pulses, as long as ( ) − ( ) approximately follows the same symmetry at the scale. Clearly, this average effect works best for closed transitions, and is compromised for open transitions as in this work, due to the population quenching (atoms leaving {| , | } to other states at the 1/Γ time scale).
We finally remark that the Eq. (10) approximation requires the | − | transition is strong, with large enough dipole element d comparing to nearby off-resonant transitions. For atomic states with Zeeman degeneracy, the multi-level interaction with can be decomposed into degenerate 2-level dynamics [49] to "monomorphously" contribute to the optical response. The Eq. (11) relation is therefore still valid. Finally, the Eq. (11) relation is robust against inhomogeneous broadening effects, such as those due to Zeeman shifts, Doppler shifts, or a finite laser linewidth (as in this work [50]), as long as these broadenings are small comparing to the transition linewidth Γ. These conclusions are supported by numerical simulations [30,48], with details to be presented in a separate paper [43].

B. Simulation of complex phase shift by 2-level samples
In this section we provide details on the simulation of propagation through a weakly-driven 2-level atomic sample. In all the simulations leading to Fig. 1c, Fig. 3(a-c), the atomic sample assumes a Gaussian distribution (r) with a spatial width matching the experimental situation for the 39 Here an −parameter is introduced to effectively account for the saturation effect in the simulation [39], with = 0.5 for the Fig. 3(d-f) simulation. Taking ( ) from experimental pre-characterization (Sec. 2.2), numerical propagation of ( ) across the sample follows a split-operator method [51], by interleaving the freespace propagation (Eq. (7)) tot ( + ) = ( ) tot ( ) with spatial-dependent phase shift tot ( + d ) = tot ( ) ( ( )−1) d in small steps d . Starting from tot = ( − ), we obtain tot = ( + ) across the full sample length in the model. The full field is then decomposed into tot = + to find the "true" complex phase shift ( ) at the atomic central plane = according to Eq. (3), as well as to numerically generate 1,2 at the camera plane = for verifying our holographic reconstruction and sample plane localization algorithms. In this section we provide additional experimental details for holographic imaging of 39 K atoms. As mentioned in Sec. 3, the { ( ) 1,2 } data at detunings {Δ } are obtained with repeated cycles of atomic sample preparation and measurements. A single cycle takes about 2 seconds time. The timing sequence for a single cycle is shown in Fig. 6. Firstly, a standard magneto-optical trap (MOT) [29] with D2 line cooling captures ∼ 10 6 atoms with a temperature ≈ 0.5 mK. We then apply D1 molasses cooling [30], with ∼ 100 mW of 770 nm light delivered to the atoms with the same MOT beams. During this cooling stage, the D1 molasses contain two sidebands to address both hyperfine ground states 4 1/2 , = 1, 2 at Raman 2-photon resonance [30]. The single-photon detuning to = 2 is set as Δ = 2 × 16 MHz. To create the microscopic sample for the imaging study, an optical dipole trap (ODT) at wavelength 780 nm is switched on simultaneously. The ODT is fairly strongly focused by another NA=0.3 lens array to reach a Gaussian waist of ≈1.5 m at the MOT center, with a trap depth ≈ 15 MHz. This combined D1 in-trap cooling of about 3 ms is able to prepare the = 10 3 atomic sample at a temperature of ≈ 10 K. Next, in the holographic imaging step, we perform "switching detection" by interleaving = 1 s D1 probe pulses with = 4 s D1 single-sideband cooling + 780 trapping pulse, for a total exposure time of 1 ms. The level diagrams for the light-atom interaction during probing and cooling are summarized in Fig.7. In particular, during the D1 in-trap cooling that only addresses the = 2 hyperfine level, atoms are not only cooled and trapped, but also "depumped" back to = 1 for the next probe.