# Imaging moving atoms by holographically reconstructing the dragged slow
light

Yuzhuo Wang1<EMAIL_ADDRESS>Jian Zhao1 Current address: KLA-Tencor
Semiconductor Equipment Technology (Shanghai) Co., Ltd. Xing Huang1 Liyang
Qiu1 Lingjing Ji1 Yudi Ma1 Yizun He1 James P. Sobol2 Saijun Wu1
<EMAIL_ADDRESS>1Department of Physics, State Key Laboratory of Surface
Physics and Key Laboratory of Micro and Nano Photonic Structures (Ministry of
Education), Fudan University, Shanghai 200433, China.
2Keit Spectrometers, Harwell Campus, Oxfordshire UK.

###### Abstract

The propagation of light in moving media is dragged by atomic motion. The
light-drag effect can be dramatically enhanced by reducing the group velocity
with electro-magnetically induced transparency (EIT). We develop a systematic
procedure to accurately reconstruct the complex wavefront of the slow light
with single-shot measurements, enabling precise, photon shot-noise limited
spectroscopic measurements of atomic response across EIT even in presence of
generic atomic number fluctuations. Applying the technique to an expanding
cloud of cold atoms, we demonstrate simultaneous inference of the atomic
density distribution and the velocity field from the complex imaging data.
This inline imaging technique may assist a wide range of cold atom experiments
to access spectroscopic and phase space information with in situ and minimally
destructive measurements.

††preprint: APS/123-QED

##  I Introduction

Investigations of light-drag effects in moving media have historically
inspired people to better understand the nature of light, and for the
establishment of special relativity [1, 2]. More recently, with the advance of
quantum optical and cold atomic technology, there is resurgent research
interest in light-drag effects [3, 4, 5, 6, 7], particularly in coherently
prepared optical media featuring electro-magnetically induced transparency
(EIT) [8, 9], where the light-drag effect is dramatically enhanced by the
highly dispersive refractive index and enormously reduced group velocity [10].
The non-uniformly dragged dynamics of slow light in such medium have inspired
exciting concepts and developments on optical analogy of gravitational effects
[11, 12]. Practically, since the dragged optical phase senses the relative
motion between atoms and light, the effect is useful for inertial sensing and
to support quantum enhanced performance of light pulse atomic interferometers
[13, 14, 15, 5, 6].

Previous studies of light-drag effects in atomic medium usually resort to
optical interferometry for the phase readout [13, 3, 14, 5, 6]. The atomic
sample is typically placed within one arm of the interferometer. The
transmitted wavefront of the probe arm is compared with light in the reference
arm at the output beamsplitter to measure the dragged phase shift. Although a
single mode optical interferometer is ideal for sensing the average motion,
spatial dependent information is disregarded in such setups. Visualization of
the spatially dependent light-drag effect can form a powerful tool to retrieve
phase-space information of dilute atomic samples with generic fluctuation of
density and coupling strength to light. The phase space information may help,
for example, to suppress classical noises in atom interferometers [13, 16] or
to assist quantum feedback in degenerate gases [17, 18, 19]. However, such
visualization requires one to recover the complex wavefront of the slow light
within single shot camera exposure, at sensitivity close to the photon shot-
noise limit [20] so as to minimize the photon-recoil induced back-actions. The
required wavefront sensing capacity can in principle be achieved with
holographic imaging by recovering optical phase information from intensity
measurements [21, 22, 23, 24, 25, 26]. To this end, efforts have been made
over years for holographic imaging of cold atoms [27, 28, 20, 29, 30].
However, the achievable imaging accuracies in these work are limited by
various approximations, typically designed to resolve the “twin image problem”
associated with the phase ambiguity [31], such as by pre-assuming optical
properties of atoms [28, 30] or idealizing probe/reference wavefronts [27, 28,
20, 29, 30]. To suppress the associated errors is prerequisite to spatially
resolving the kHz-level 2-photon Doppler shift of a dilute gas. The resulting
atomic density-error immune spectroscopic imaging capacity has not been
demonstrated previously, to the best of our knowledge.

This work reports holographic reconstruction of slow light and spectroscopic
imaging of moving atoms under the EIT condition. For the purpose, we develop
an accurate procedure to faithfully recover optical phases from single-shot
intensity measurements. To resolve the phase ambiguity without typical
approximations [27, 28, 20, 29, 30], we apply a general iterative algorithm
[20] to reconstruct coherent forward scattering from atomic samples
illuminated by a precisely pre-characterized wavefront [32, 33, 34]. By
comparing the reconstructed coherent scattering with the unperturbed probe,
the complex optical response is obtained for the simultaneous inference of the
atomic density distribution and the velocity field. We also demonstrate
single-shot spectroscopic imaging across EIT resonance with photon-shot-noise
limited accuracy for the phase angle of the transmitted slow light, even in
presence of large shot-to-shot atom number fluctuations. Based on Gabor’s
holographic microscopy (GHM) [31] with intrinsic phase stability, the inline
method can be conveniently set up to assist a wide range of cold atom
experiments to access phase space information with in situ and minimally
destructive measurements. Our imaging technique also paves a practical pathway
toward precise spectroscopic imaging in presence of density fluctuations
generic to most ultra-cold atomic samples [35, 36, 37, 38].

Figure 1: Holographic sensing of slow light wavefront. (a) Experimental setup
and holographic reconstruction procedure: A sample of 87Rb atoms centered at
$z=z_{\rm A}$ is illuminated by a defocused probe light $E_{\rm p}$ collected
by an $NA=0.4$ imaging optics. By mounting the CCD camera on a linear stage,
the probe beam $E_{\rm p}$ can be pre-characterized by multi-plane intensity
measurements. The EIT condition (with the level diagram in (b)) is established
by first illuminating the atoms with the $E_{\rm c}$ control beam and then the
$E_{\rm p}$ probe beam. A controllable velocity distribution is achieved by
subjecting atoms to $E_{\rm SW}$ standing wave pulses along $y$. The camera
data $I_{1,2}$ with and without atoms are matched to the precisely pre-
characterized $E_{\rm p}=|E_{\rm p}|e^{i\phi_{\rm p}}$, including aberrations
and speckles, to retrieve the atomic coherent scattering $E_{\rm s}$. The slow
light wavefront and complex transmission $e^{i\varphi}=1+E_{\rm s}/E_{\rm p}$
can then be derived at $z_{\rm A}$ to retrieve the distribution of optical
depth ${\rm OD}=2{\rm Im}[\varphi]$ and phase shift $\phi={\rm Re}[\varphi]$,
and to infer the atomic column density $\rho_{\rm c}$ and velocity field $u$
simultaneously. See caption of Fig. 4 for the color coding of the $\rho_{\rm
c}-u$ plot. For clarity of display the $\tilde{\phi}_{\rm p}\equiv\phi_{\rm
p}-\phi_{\rm sp}$ plot gives the phase difference between $\phi_{\rm p}$ and
the optical phase of a best-matched spherical wave $\phi_{\rm sp}$. (c) The
timing sequence.

##  II Transmission of dragged slow light

We consider the experimental setup illustrated in Fig. 1a. A cold 87Rb atomic
sample centered at $z=z_{\rm A}$ plane is illuminated by a $\lambda_{p}=780$
nm probe beam $E_{\rm p}$ with frequency $\omega\approx\omega_{eg}$ resonant
to the $F=2-F^{\prime}=2$ hyperfine D2 transition. $E_{\rm p}$ also serves as
a reference beam to interfere with the coherently scattered light from the
atomic sample, $E_{\rm s}$, to be recorded on the $z=z_{\rm H}$ plane. The EIT
condition is established for $F=2$ ground state atoms by illuminating them
with a control beam $E_{\rm c}$ with frequency $\omega_{c}=\omega_{ea}$
resonant with the $F=1-F^{\prime}=2$ hyperfine transition.

The transmitted wavefront of the weak probe light through the atomic medium
can be expressed as $E_{\rm out}=E_{\rm p}+E_{\rm s}$. At the exit side of the
atomic sample, $E_{\rm out}(z_{\rm A})$ is related to $E_{\rm p}(z_{\rm A})$
through a simple relation, $E_{\rm out}(z_{\rm
A},\omega)=e^{i\varphi(\omega)}E_{\rm p}(z_{\rm A},\omega)$, for samples with
moderate thickness $L$ in the Raman-Nath regime. Here and in the following we
may omit the coordinate variables $x,y$ in the $E(z,\omega)$,
$\varphi(\omega)$ expressions for conciseness. For a dilute gas, the complex
phase $\varphi(\omega)$ is given by

$\varphi(x,y,\omega)=\int_{z_{\rm A}-L}^{z_{\rm A}}\frac{\chi({\bf
r},\omega)}{2}\frac{\omega}{c}\differential z,$ (1)

which is divided into $\varphi(x,y,\omega)=\phi+i{\rm OD}/2$ to represent the
phase shift $\phi(x,y)$ and optical depth ${\rm OD}(x,y)$. We consider the
atomic phase space distribution $f({\bf r},{\bf v})$ and the associated
spatial density $\rho({\bf r})=\int f({\bf r},{\bf v})\differential^{3}{\bf
v}$. Recall that EIT is associated with coherent spin-wave excitation with
wavevector ${\bf k}_{\rm eff}={\bf k}_{\rm p}-{\bf k}_{\rm c}$, with ${\bf
k}_{\rm p}$, ${\bf k}_{\rm c}$ the wavevectors of the probe and control fields
locally seen by the atoms respectively [39]. The Doppler-shift leads to an
inhomogeneously broadened polarizability, $\chi({\bf
r},\omega)=\int\differential^{3}{\bf v}f({\bf r},{\bf v})\alpha(\omega-{\bf
k}_{\rm eff}\cdot{\bf v})$, where $\alpha(\omega)$ is the atomic
polarizability with the imaginary part $\alpha_{\rm i}(\omega_{0})$ minimized
at the 2-photon resonance frequency $\omega_{0}=\omega_{c}-\omega_{ga}$. Here
$\omega_{ga}=2\pi\times 6.835$ GHz is the hyperfine splitting between the
$F=1,2$ ground states [40]. We consider $\tilde{\omega}_{0}=\omega_{0}+{\bf
k}_{\rm eff}\cdot{\bf v_{0}}$ to be a Doppler-shifted EIT resonant frequency
and Taylor-expand $\alpha(\omega-{\bf k}_{\rm eff}\cdot{\bf v})$ to have

$\varphi(x,y,\tilde{\omega}_{0})\approx\varphi^{(0)}(x,y)+\delta\phi(x,y),$
(2)

with

$\varphi^{(0)}(x,y)=\int_{z_{\rm A}-L}^{z_{\rm A}}\rho({\bf
r})\frac{\alpha(\omega_{0})}{2}\frac{\omega_{0}}{c}\differential z,$ (3)

and

$\delta\phi(x,y)=\int_{z_{\rm A}-L}^{z_{\rm A}}\rho({\bf
r})\frac{\omega_{0}}{2}\frac{d\alpha}{d\omega}\frac{{\bf k}_{\rm
eff}\cdot\delta{\bf u}({\bf r})}{c}\differential z.$ (4)

From the complex transmission $e^{i\varphi}=E_{\rm out}/E_{\rm p}$ one shall
try to retrieve information on the velocity field $\delta{\bf u}({\bf r})={\bf
v}_{0}-{\bf u}({\bf r})$, with

${\bf u}({\bf r})=\frac{1}{\rho({\bf r})}\int\differential^{3}{\bf v}f({\bf
r},{\bf v}){\bf v},$ (5)

relative to a pre-estimated ${\bf v}_{0}$ frame. While $\varphi_{0}$ vanishes
at the 2-photon resonance for ideal EIT, practically $\varphi^{(0)}=i{\rm
OD}^{(0)}/2+\phi^{(0)}$ may arise with residual absorption and phase shifts
due to various mechanisms including non-ideal EIT preparation and ground state
relaxation [9]. By setting the probe frequency $\omega$ near
$\tilde{\omega}_{0}$, photon recoil heating during imaging can be managed to a
suitable level for minimally destructive measurement of $\varphi(\omega)$ at
within the EIT window.

We now remark on the relation between Eq. (4) and the light-drag effect. In
particular, the velocity-dependent phase shift $\delta\phi(x,y)$ is real, and
is a manifestation of the light-drag effect that modifies the light phase
velocity ${\bf v}_{\rm p}$, $\delta{\bf v}_{\rm p}\approx n_{\rm g}\delta{\bf
u}({\bf r})$, in the highly dispersive medium with group index $n_{\rm g}({\bf
r})\approx\rho({\bf
r})\frac{\omega}{2}\frac{\differential\alpha}{\differential\omega}$ [4, 5, 6,
7]. If the velocity field $\delta{\bf u}({\bf r})$ is $z$-independent, then
$\delta\phi(x,y)\approx{\bf k}_{\rm eff}\cdot\delta{\bf u}(x,y)\tau(x,y)$ with
$\tau(x,y)=\int_{z_{\rm A}-L}^{z_{\rm A}}\differential z/v_{\rm g}({\bf r})$
the optical delay introduced by the atomic medium. One may try to estimate
${\bf u}(x,y)$ along the ${\bf k}_{\rm eff}$ direction with the optical
$\delta\phi(x,y)$ readouts and additional $\tau(x,y)$ measurements. However,
in typical cold atom experiments the local atomic density $\rho({\bf r})$ is
not precisely known. On the other hand, GHM demonstrated in this work is able
to reconstruct $E_{\rm out}(z_{\rm A})$ and $E_{\rm p}(z_{\rm A})$ directly to
spectroscopically determine the 2-photon Doppler shift ${\bf k}_{\rm
eff}\cdot{\bf u}(x,y)$, so as to recover the velocity field ${\bf
u}(x,y)=-\delta{\bf u}(x,y)+{\bf v}_{0}$ along ${\bf k}_{\rm eff}$.

##  III Holographic reconstruction

Our holographic imaging method follows the philosophy of replacing standard
approximations [21, 22, 23, 24, 25, 26, 27, 28, 20, 29, 30] with measurements
whenever possible. The measurement procedure to be detailed in ref. [41] is
composed of four steps (Fig. 1a). First, before the holographic measurements,
the probe wavefront $E_{\rm p}$ is precisely pre-characterized with multiple
z-plane $I_{\rm p}(z)=|E_{\rm p}(z)|^{2}$ measurements. The 2D phase
information of $E_{\rm p}(z_{\rm H})=\sqrt{I_{\rm p}(z_{\rm H})}e^{i\phi_{\rm
p}(z_{\rm H})}$, $\phi_{\rm p}$, is recovered with an iterative Gerchberg-
Saxton algorithm [32, 33]. Here, for the nearly spherical $E_{\rm p}$, the
deviation ($\bar{\phi}_{\rm p}$ in Fig. 1) is mainly caused by aberration
through the imaging system and sources of speckles in the imaging path. Next,
$I_{1,2}$ are recorded with and without the atomic sample respectively.
Numerical translation of $E_{\rm p}$ is applied to minimize the difference
between $I_{2}$ and $|E_{\rm p}|^{2}$, thereby establishing the relation
$I_{1}=|E_{\rm p}+E_{\rm s}|^{2}$ and $I_{2}=|E_{\rm p}|^{2}$ even in presence
of slow alignment drifts. In the third step, $E_{\rm s}$ is retrieved from
$\Delta I=I_{1}-I_{2}$ with an iterative algorithm [20], applicable to general
probe beam wavefronts, to progressively removes both the twin and dc images
[42, 20, 43] based on the precise knowledge of $E_{\rm p}$ and known location
of the atomic sample. Also in this step, the actual imaging resolution ${\rm
NA}<{\rm NA}_{0}$ is flexibly chosen by digitally reducing the exposure area
of the camera for the reduced hologram $\Delta I$. In the final step, we
numerically propagate $E_{\rm p}$ and $E_{\rm s}$ to the $z_{\rm A}$ plane to
retrieve the complex phase $\varphi(x,y,\omega)=-i{\rm log}(\frac{E_{\rm
p}+E_{\rm s}}{E_{\rm p}})$. Regular ${\rm OD}(x,y)=2{\rm Im}[\varphi]$ and
$\phi(x,y)={\rm Re}[\varphi]$, as well as “dark ground” $|E_{\rm s}(x,y)|^{2}$
imaging data [44] are simultaneously retrieved. In addition, we obtain the
less often discussed phase angle $\beta(x,y)={\rm arg}[\varphi]$ and modulus
$|\varphi(x,y)|$ images.

In the GHM setup sketched in Fig. 1a, an ${\rm NA}_{0}$=0.4 imaging system
[45] relays both the probe light $E_{\rm p}$ and the scattered wavefront from
the atomic sample $E_{\rm s}$ to the $z=z_{\rm A}$ plane. The nearly spherical
$E_{\rm p}$ [31] leads to a reduction in the probe intensity at the hologram
plane $I_{\rm p}(z_{\rm H})=|E_{\rm p}(z_{\rm H})|^{2}$ relative to the
intensity $I_{\rm p}(z_{\rm A})=|E_{\rm p}(z_{\rm A})|^{2}$ effectively seen
by the atoms, a feature that has been utilized to enhance the detector dynamic
range limited atomic sensing capacity [20]. Here, to ensure plenty of camera
counts $I_{1,2}$ within a short probe interval $\tau_{\rm p}$, we choose
$z_{\rm A}=1.6$ mm and $z_{\rm H}=3.2$ mm with a moderate enhancement factor
of $(z_{H}/z_{A})^{2}=4$.

##  IV Complex atomic spectroscopy

We first demonstrate shot-noise limited spectroscopic imaging of stationary
atomic sample in presence of fluctuating atom number. Here the probe beam is
set with an intensity parameter $s\approx 2.2$ ($s=I/I_{\rm s}$ with
saturation intensity $I_{\rm s}=1.67~{}{\rm mW/cm}^{2}$) at the location of
the atoms. Holograms are recorded with an exposure time of $\tau_{\rm
p}=20~{}\mu$s. The atomic sample is composed of $N_{A}\approx 1\times 10^{4}$
atoms, nearly spherical with a diameter $L\approx 13~{}\mu$m and temperature
$T\approx 15~{}\mu$K. The control beam has an intensity of $\sim 7$ mW/cm2 to
maintain an EIT window width $\delta\omega_{\rm EIT}$ of 1 MHz. Holograms are
recorded in repeated experiments with the probe frequency $\omega$ scanning
across the 2-photon resonance $\omega_{0}$. Typical reduced holograms $\Delta
I=I_{1}-I_{2}$, with $\bar{I}_{1}\approx 300$ in terms of camera counts on an
$A_{\rm p}=6.45\times 6.45~{}\mu$m2 pixel area are given in Fig. 2a (top)
(Intensity of the probe for recording $I_{2}$ is increased by a factor of 25
to reduce the shot-noise contribution.). In addition to the expected rings of
interference signals, one clearly see the granular photon shot noise
background, which is at a 20-counts root-mean-square level. Following the
described procedure, we match $I_{2}$ to $|E_{\rm p}|^{2}$, and use $E_{\rm
p}$ to recover $E_{\rm s}$ from single holograms to retrieve the absorption
${\rm OD}(x,y)$ and phase shift $\phi(x,y)$. By averaging the simultaneously
reconstructed ${\rm OD}$ and $\phi$ images as in Fig. 2a over a region of
interest (ROI) that covers most of atomic signals, we obtain EIT curves for
average absorption $\overline{\rm OD}$ (Fig. 2b), phase shift $\bar{\phi}$
(Fig. 2e), phase angle $\bar{\beta}={\rm arg}(\bar{\varphi})$ (Fig. 2d) and
phasor diagram (Fig. 2c) of $\\{\bar{\phi},\overline{\rm OD}/2\\}$ plots [46].
From the phase shift data in Fig. 2(e), in particular, we can estimate an
average group index of $n_{g}\sim 3\times 10^{6}$ and a reduced group velocity
of light $v_{\rm g}$ at the 100 m/s level [10] across the sample within the
$\sim$1 MHZ EIT window.

Figure 2:  Spectroscopic imaging across the EIT resonance. (a) From top to
bottom: Reduced holograms $\Delta I$ and the derived ${\rm OD}$ and $\phi$
images with probe detunings $\omega-\omega_{0}$ from $-6$ MHz to $6$ MHz with
$\sim$2 MHz interval. An elliptical ROI with $\mathcal{A}_{\rm ROI}\approx$
130 $\mu$m2 is marked with white dotted lines. The images are from a sub-set
of single-shot data presented in Figs. (b-e). Figs. (b)(d)(e) are ROI-averaged
spectra for optical depth ($\overline{\mathrm{OD}}$), phase angle
($\overline{\beta}$) and phase shift ($\overline{\phi}$) respectively. Single-
shot data are plotted with semi-transparent symbols to illustrate the
statistical noise level, while the mean value for $N_{m}=9$ are plotted with
dark dots. The data is further plotted as a phasor diagram in (c), with 2D
errorbars giving the standard deviations along the $|\bar{\varphi}|$ and
$\beta$ directions. The solid curves are fit to a phenomenological EIT model
[41].

Unless working with confined single atoms [37], a standard challenge in cold
atom laser spectroscopy [35, 36, 38] is associated with fluctuating atom
number and coupling strengths during the measurements. Here, the spectroscopy
is recorded by repeatedly preparing the cold atomic samples. The shot-to-shot
atom number fluctuation translates into a spreading of data points in both
$\overline{\rm OD}$ (Fig. 2b) and $\overline{\phi}$(Fig. 2e) spectroscopy.
Looking into individual $\overline{\rm OD}$ or $\overline{\phi}$ data, it is
impossible to isolate such fluctuation from frequency-dependent atomic
response of interest. To efficiently suppress the noise, one may have to
resort to normalizing the atom number with double measurements as carefully as
possible [38]. Here, a unique strength of the complex spectroscopy is to
unveil the correlation of such fluctuation, as demonstrated in the phasor
diagram Fig. 2(c) where the EIT resonance leads to a circle of complex atomic
response. Unlike the photon shot noise, the spreading of the spectroscopic
data in the Fig. 2(c) phasor diagram is highly annisotropic since the shot-to-
shot atom number fluctuation merely shifts the data points along the
$|\overline{\varphi}|$ direction.

Figure 3:  Statistical analysis of the image noise. (a) Illustrated with a
$|E_{\rm s}|$ image, the noise level $n_{\overline{\rm OD}}(\mathcal{A})$ and
$n_{\overline{\phi}}(\mathcal{A})$ are evaluated from 100 complex-$\varphi$
images as $n_{\overline{\varphi}}(\mathcal{A})=\sqrt{n_{\overline{\rm
OD}}^{2}/4+n_{\bar{\phi}}^{2}}$, within area $\mathcal{A}$ excluding the atom
sample. (b): Dependence of $n_{\bar{\varphi}}(\mathcal{A})$ on radius $r_{M}$:
Blue square ([A]) and green square ([B]) symbols are based on experimental
$\varphi$-data measurements in presence and absence of the atom sample
respectively. Black disks ([C]) are based on simulated shot-noise limited
data. The black line gives $n=1/\sqrt{N_{\mathrm{p}}(\mathcal{A})}$, with
$N_{\mathrm{p}}(\mathcal{A})$ to be the total number of photon counts through
area $\mathcal{A}$ detected by the camera. Notice here for $r_{M}$ smaller
than half the NA=0.13-limited $\delta x=6~{}\mu$m, the imaging noise are
numerically smoothed down, leading to deviation of the data points from the
power law, or even below the shot noise limit. (c): Noise of the phase angle
$\overline{\beta}$ for atomic polarizability $\alpha$, estimated by
$\varphi$-data within $\rm ROI$. The red diamond symbol ([D]) gives root-mean-
square value of the shot-to-shot $\overline{\beta}$ fluctuation estimated from
groups of 36 measurements. Blue squares ([E]) give
$n_{\bar{\beta}}(\mathcal{A}_{\rm
ROI})=n_{\overline{\varphi}}/\sqrt{2}|\bar{\varphi}|$ estimated according to
the [A] method in Fig. (b). The black line ([F]) gives
$n_{\bar{\beta}}(\mathcal{A}_{\rm ROI})=\xi/\sqrt{N_{\rm s}(\mathcal{A}_{\rm
ROI})}$, with $\xi=|(e^{-i\overline{\varphi}}-1)/\overline{\varphi}|$ to be an
optical depth dependent correction factor [41].

We now focus on the phase angle $\overline{\beta}={\rm
arg}(\overline{\varphi})$, which can be resolved by the complex spectroscopy
from single shot measurements. For the dilute and optically thin samples,
$\bar{\beta}$ is directly decided by the phase angle of the atomic
polarizability $\alpha(\omega)$ (Eq. (1)) and is insensitive to the atom
number fluctuation. Indeed, the shot-to-shot noise level for
$\overline{\beta}$ in Fig. 2(d) within the EIT window is close to being
photon-shot-noise limited. The shot-noise-limited performance is demonstrated
in Fig. 3. Here, with $|\varphi|\ll 1$ and therefore $\varphi\approx-iE_{\rm
s}/E_{\rm p}$ at the EIT resonance, the phase angle measurement is equivalent
to interferometrically estimating the optical phase of the forward scattering
$E_{\rm s}$ relative to the reference field $E_{\rm p}$ [47] through an
aperture defined by the ROI. It is known that the holographic microscopy is a
“spatial heterodying” measurement [27]. Therefore, the phase angle measurement
here has a resolution bounded by $n_{\bar{\beta}}>1/\sqrt{N_{\rm
s}(\mathcal{A}_{\rm ROI})}$, with $N_{\rm s}(\mathcal{A}_{\rm ROI})$ to be the
amount of elastically scattered photons by the atoms through ROI to be
detected by the camera in term of interference with $E_{\rm p}$. More
generally, for optically thick samples we have
$n_{\bar{\beta}}>\xi/\sqrt{N_{\rm s}(\mathcal{A}_{\rm ROI})}$, with
$\xi=|(e^{-i\overline{\varphi}}-1)/\overline{\varphi}|$ to be an optical depth
dependent correction factor [41]. We verify the shot-noise-limited phase angle
resolution by directly evaluating the total number of scattered photons from
the reconstructed $E_{\rm s}$ field at $z=z_{\rm A}$ as $N_{\rm
s}(\mathcal{A}_{\rm ROI})=\sum_{\rm ROI}|E_{\rm s}|^{2}$ after proper
normalization (Fig. 3c, the line plot) [41]. For $\omega\approx\omega_{0}$
with optimal EIT transparency we typically find $N_{\rm s}(\mathcal{A}_{\rm
ROI})\approx 2000$, predicting shot-noise-limited $n_{\overline{\beta}}\approx
0.02$. The shot-to-shot $\bar{\beta}$ fluctuation as displayed by the diamond
symbols in Fig. 3(c) is twice as large at the EIT resonance. The lightly
noisier $\bar{\beta}$ is likely caused by power fluctuation of EIT control
field ${\bf E}_{c}$ in repeated measurements which affects $\alpha(\omega)$
itself [41]. With increased light scattering when $\omega-\omega_{0}$ is off
the EIT resonance, the $\bar{\beta}$ fluctuation deviates further from the
$\xi/\sqrt{N_{\rm s}(\mathcal{A}_{\rm ROI})}$ prediction, due to the technical
noise.

The precision of the phase-angle measurement is further characterized by
estimating $n_{\bar{\beta}}=n_{\overline{\varphi}}/\sqrt{2}|\bar{\varphi}|$
(Fig. 3c the square symbols) in an atomic sample independent manner using
individually reconstructed complex data of $\varphi=\phi+i{\rm OD}/2$. The
method assumes the noise level in each $\rm OD$, $\phi$ images is
approximately constant at the vicinity of the sample location ROI. As
illustrated in Fig. 3(a), with the probe photon shot noise dominating the
holography noise channels and with $|z_{H}-z_{A}|\gg\mathcal{A}/\lambda$ for
the shot noise to digitally propagate [20], the standard deviation
$\\{n_{\bar{\phi}}(\mathcal{A}),n_{\overline{\rm OD}}(\mathcal{A})/2\\}$
evaluated as root-mean-square values from repeated measurements is found to be
quite isotropic and can thus be simultaneously characterized by the noise
level to the complex phase
$n_{\overline{\varphi}}\equiv\sqrt{n^{2}_{\bar{\phi}}+n^{2}_{\overline{\rm
OD}}/4}\approx\sqrt{2}n_{\overline{\phi}}$. From Figs. 3(a)(b) we see
$n_{\overline{\varphi}}$ itself is close to the probe photon shot limit
$1/\sqrt{N_{\rm p}}$ [20], decided by the amount of probe photons $N_{\rm
p}(\mathcal{A})$ through the imaging area $\mathcal{A}$ detected by the
camera. Here, with $n_{\overline{\varphi}}(\mathcal{A})\approx
1.3/\sqrt{N_{\rm p}(\mathcal{A})}$, the slightly larger
$n_{\overline{\varphi}}(\mathcal{A})$ is associated with additional technical
noises, and may also be related to imperfect digital holographic processing
such as those during $I_{1,2}$ subtraction [41]. With the shot-noise-limited
$n_{\bar{\varphi}}$, the phase angle noise level in Fig. 3c (the square plot)
evaluated as $n_{\bar{\beta}}(\mathcal{A}_{\rm
ROI})=n_{\overline{\varphi}}(\mathcal{A}_{\rm
ROI})/\sqrt{2}|\bar{\varphi}(\mathcal{A}_{\rm ROI})|$ (with modulus
$|\bar{\varphi}(\mathcal{A}_{\rm
ROI})|=\sqrt{\overline{\phi}^{2}+\overline{\rm OD}^{2}/4}$ evaluated within
ROI) is close to the $1/\sqrt{N_{\rm s}(\mathcal{A}_{\rm ROI})}$ shot noise
limit, as expected.

##  V Velocity field sensing

We create atomic samples with spatially dependent velocity field ${\bf u}({\bf
r})$ by subjecting the dipole trapped atoms with standing wave diffraction,
followed by a $\tau_{\rm tof}$ expansion along $y$ in a quasi-2D trap with
$x-z$ confinements (Fig. 1). The standing wave is formed along $y$ by
$\lambda_{\rm a}=795~{}$nm counter-propagating D1 light (Fig. 1a) blue-detuned
from the $F=1-F^{\prime}=2$ hyperfine resonance by 1 GHz. The pulse is
powerful ($I_{\rm a}=10{{}^{3}}$W/cm2) and short ($\tau_{\rm a}=50~{}$ns) for
precise modeling in this investigation. Microscopically, a periodic repulsive
light shift potential is formed by the optical standing wave for the
bidirectional acceleration of atoms along $y$. Macroscopically, after the
optical impact the phase space distribution of atomic sample along $y$ is
characterized by $f(y,v_{y})=\rho(y)g(v_{y})$ with a velocity distribution
$g(v)$ determined by the sinusoidal standing wave shift. The control field
${\bf k}_{\rm c}$ is aligned in the $y-z$ plane with a small angle
$\iota=11.1^{\circ}$ relative to ${\bf e}_{y}$ [48]. The nearly free expansion
leads to $f(y,v_{y},t)=f(y-v_{y}t,v_{y})$, resulting in a velocity field ${\bf
u}(y,t)$ (Eq. (5)) approximately along $y$.

We probe the expanding atomic sample after $\tau_{\rm tof}=20~{}\mu$s for a
$\tau_{\rm p}=10~{}\mu$s exposure time. The exposure is short enough that most
atoms hardly move a distance beyond the numerically-reduced NA$\approx
0.13$-limited imaging resolution $\delta x\approx 6~{}\mu$m [41]. The probe
field $E_{\rm p}$ transmitted through the expanding sample is holographically
recorded by the camera. We retrieve the phase shift holographically as the
real part of the complex phase $\varphi$ (Fig. 4(a,iii)), together with the
absorption (Fig. 4(a,ii)). For a certain probe frequency $\omega$, a
particular velocity class of atoms satisfying $\omega-\omega_{0}={\bf k}_{\rm
eff}\cdot{\bf v}_{0}$ falls into the EIT window with enhanced transmission.
Around the EIT window, the phase shift $\phi$ always increases linearly with
$y$ for all the probe frequency $\omega$, suggesting the same $u(y)-y$
dependence as expected from the point-source expansion. With a reduced peak
atomic density $\rho$ for the expanding cloud at a $5\times 10^{11}/{\rm
cm}^{3}$ level, the group velocity of the slow light is increased to $v_{\rm
g}\approx 1000$ m/s, and it takes $\tau=L/v_{\rm g}\approx 20$ ns for the
probe pulse to trespass the atomic sample at the EIT center. The optical phase
of the probe wavefront locally follows the atomic motion, leading to
$\delta\phi(x,y)=k_{u}u(x,y)\tau$ phase shift in the transmitted wavefront for
${\bf v}_{0}=0$ in particular. The ${\bf u}$-component of ${\bf k}_{\rm eff}$
assuming ${\bf u}=u{\bf e}_{y}$ is $k_{u}=k_{\rm c}{\rm cos}(\iota)+\epsilon
k_{\rm p}$, with $\epsilon={\bf k}_{\rm p}\cdot{\bf e}_{y}/k_{\rm p}$
numerically evaluated from the $E_{\rm p}$ wavefront at the atomic location.
For $|u|<0.3$ m/s within the linear dispersion regime, the maximum shift
$\phi_{\rm M}$ is at a 50 mrad level, suggesting a dragged transverse
wavefront displacement by $(\delta y)_{\rm M}=\phi_{\rm M}/2\pi\lambda_{\rm
p}\approx 5$ nm during the slow light travel. However, without a prior
knowledge of the atomic density distribution, precise characterization of the
atomic motion cannot be realized with the $\phi(x,y)$ image alone.

Here, with the imaging data expanded from real to complex numbers by
holography (See the $|\varphi|-\beta$ color domain plot in Fig. 4(b,iii-iv).),
we have the opportunity to estimate atomic density $\rho$ and the velocity
component $u$ simultaneously from single shot measurements. In particular, for
the thin-lens regime of light scattering by the dilute sample, the atomic
column density $\rho_{\rm c}(x,y)$ and velocity field $u(y)$ can be inferred
from the $\varphi(x,y)$ data by minimizing
$\mathcal{L}=\sum_{x,y}|2\varphi(x,y)-k_{\rm p}\rho_{\rm c}\alpha(\omega-
k_{u}u)|^{2}$ [41]. The atom polarizability $\alpha(\omega)$ is obtained
beforehand by fitting a simple EIT model to the spectroscopic measurements as
those in Fig. 2. Typical inferred $\rho_{\rm c}-u$ fields are given in Figs.
4(v,vi) probed at $\omega=\omega_{0}$ for the stationary and expanding atomic
clouds respectively. For comparison, we plot velocity field $u=u_{y}(y)$ in
Fig. 4(vii) estimated directly as $u_{y}=y/(\tau_{\rm tof}+1/2\tau_{\rm p})$
following the point source expansion [41]. Notice a curved $\rho_{\rm c}$
distribution suggests that the pulsed standing wave light shift induces
certain “slosh” motion of atoms along $x$ during the expansion in the dipole
trap. We expect similar motion along $z$ which is close to ${\bf k}_{\rm p}$
to be EIT-sensitive. The small discrepancies between the Fig. 4(b,vi-vii)
velocity field estimations could be due to the additional $z$-motion sensed by
the EIT spectroscopy not to be recorded by the displacement along $y$ in the
same images.

Figure 4: Holographic reconstruction of spatial-dependent light drag effect
for velocity field sensing. (a) From top to bottom: Reduced holograms $\Delta
I$ and the derived ${\rm OD}$ and $\phi$ images for the standing-wave-
accelerated atomic samples. The holograms are taken near the EIT condition
with the two-photon detuning $\delta$ scanned across the EIT resonance. Each
shot-noise-limited hologram is averaged 80 times to enhance the signal to
noise for the display. (b) Reduced holograms for atomic sample at rest (i) and
expanding along $y$ (ii), recorded with zero 2-photon detuning. The
reconstructed complex-$\varphi$ data are displayed with color-domain images in
(iii)(iv). From the $\varphi$-data, smooth $\rho_{c}-u$ distributions are
inferred in (v,vi). Direct $\rho_{c}-u$ retrieval assuming $u_{y}(y)$
following the one-source expansion is given in Fig. (vii) [41]. Here the
velocity is coded in color and atomic column density is coded in brightness.

##  VI Summary and outlook

In this work, we have demonstrated a technique to precisely reconstruct the
full wavefront of slow light transmitted through cold atomic samples with
single-shot intensity measurements. We have shown that as in other coherent
imaging schemes [49, 20], the quality of the reconstructed spectroscopic
imaging data is close to that limited by the shot noise of the elastically
scattered photons from atoms, and can thus be minimally destructive for
velocity field sensing. The key advantage of the holographic imaging method is
the ability to retrieve absorption and phase shift simultaneously. By
expanding the imaging data from real to complex numbers, the mixed information
on the velocity field and density distributions can be disentangled and
separately inferred.

Comparing with previous works [21, 22, 23, 24, 25, 20, 26], our method is a
simplest application of Gabor’s original holography [31], but with various
pre-assumptions of unknowns removed for better accuracy. Future development of
this holographic spectroscopic imaging approach, particularly toward imaging
through optics with larger numerical aperture [50, 51], may offer a platform
that combines high spatial resolution [50, 51, 52, 53] with interferometric
precision for spectroscopic imaging [37] and velocity sensing [54, 13, 55, 56]
with ultracold atoms.

## Funding Information

National Key Research Program of China (2017YFA0304204, 2016YFA0302000);
National Science Foundation of China (11574053). Shanghai Scientific Research
Program (15ZR1403200).

## Appendix

###  VI.1 Holographic retrieval of probe transmission and atomic response

A key technique in this work is to reconstruct the transmitted wavefront
$E_{\rm out}=E_{\rm p}+E_{\rm s}$ to be compared with $E_{\rm p}$ so as to
simultaneously obtain the atomic absorption ${\rm OD}(x,y)$ and phase shift
$\phi(x,y)$ from single-shot holographic measurements. The holographic
technique is composed of four steps as following.

First, before the atomic imaging measurements, we pre-characterize the
wavefront of the probe field, $E_{\rm p}$, through multi-plane intensity
measurements followed by a multi-plane Gerchberg-Saxton (G-S) algorithm to
reconstruct the phase information [32, 33].

Next, for the holographic imaging, defocused images $I_{1}=|E_{\rm p}+E_{\rm
s}|^{2}$ and $I_{2}=|E_{\rm p}|^{2}$ are recorded at a “holography” plane
$z_{\rm H}$, with and without the atomic sample respectively, as outlined in
the main text (Fig. 1). Taking into account ambient noise that affect both
$I_{1}$ and $I_{2}$, in this step a numerical procedure is applied to
optimally match $I_{2}$ to the pre-characterized $E_{\rm p}$ and to derive the
reduced hologram $\Delta I=I_{1}-I_{2}$.

Thirdly, $E_{\rm s}$ is numerically reconstructed from the reduced hologram
$\Delta I=I_{1}-I_{2}$ with an iterative algorithm to remove the twin and dc
noises [43].

Finally, $E_{\rm s}$ and $E_{\rm p}$ are numerically propagated to atomic
plane $z=z_{\rm A}$ to evaluate the complex phase $\varphi(x,y)=\phi+i{\rm
OD}/2=-i{\rm log}(\frac{E_{\rm p}(z_{\rm A})+E_{\rm s}(z_{\rm A})}{E_{\rm
p}(z_{\rm A})})$.

As in the main text, we may omit the $x,y$ (and $k_{x},k_{y}$) coordinate
variables in the notions of wavefronts and intensities if ambiguity would not
be introduced.

To propagate a wavefront $E(x,y,z)$ from $z=z_{\rm A}$ to $z=z_{\rm B}$, we
use the angular spectrum method with the following notations:

$\begin{array}[]{l}E(z_{\rm B})=\hat{U}(z_{\rm B}-z_{\rm A})E(z_{\rm A}),{\rm
with}\\\
\hat{U}(L)={\hat{F}}^{-1}e^{i\sqrt{k^{2}-k_{x}^{2}-k_{y}^{2}}L}\hat{F}\end{array}$
(6)

Here $k=2\pi/\lambda$ is the wavenumber of the light field, and $\hat{F}$
represent the 2D Fourier transform [57], $i.e.$,
$\hat{F}E(x,y,z)=E(k_{x},k_{y},z)$ and
$\hat{F}^{-1}E(k_{x},k_{y},z)=E(x,y,z)$.

In the following we detail the four-step holographic reconstruction procedure.

####  VI.1.1 Pre-characterization of $E_{\rm p}$

Figure 5: Schematic setup for the measurement (a) and numerical reconstruction
(b) of the probe light wavefront $E_{\rm p}$. The intensity profile $I_{\rm
p}(z_{j})=|E_{\rm p}(z_{j})|^{2}$ is recorded at multiple $z-$planes with a
digital camera mounted on a translation stage. The flow diagram in (b)
describes a single 3-plane G-S iteration for the complex $E_{\rm p}$
reconstruction. The full reconstruction is composed of multiple 3-plane G-S
iterations with $\\{I_{\rm p,j}\\}$ data randomly chosen from the pre-recorded
data-set.

As in Fig. 5(a), we use a CCD camera mounted on a stepper-motor-driven linear
stage to sample the probe field intensity distribution $I_{\rm
p}(x,y,z)==|E_{\rm p}(x,y,z)|^{2}$ in 3D, by translating the camera along $z$
and record a set of 2D intensity profiles $I_{\rm p,j}=I_{\rm p}(z_{j})$. Here
the $x-y$ plane is defined by the camera sensor chip. The direction normal to
the plane is thus referred to as the $z-$direction, which is adjusted to be
parallel to the direction of translation to within $\sim$10 mrad angular
precision. At each $z_{j}$-plane, the intensity $I_{j}$ is adjusted with an
acoustic-optical modulator (AOM) to ensure sufficient counts on the camera
without saturation. The distance $z_{j}$ is digitally recorded, which is
subjected to numerical adjustments as detailed in the following. The imaging
data is referred to as a $\\{I_{\rm p}\\}$-set.

The 3D intensity sampling is subjected to systematic noises such as multiple
reflections involving the camera sensor chip itself. To recover $E_{\rm p}$
from the imperfect $\\{I_{\rm p}\\}$ measurement, it is important to sample
many $z_{j}-$ planes to suppress the position-dependent noise so as to
determine $E_{\rm p}$ through the average intensity field. To achieve the
purpose without demanding computer memory, we apply a multi-plane G-S
algorithm that samples the $\\{I_{\rm p}\\}$-set in small groups, detailed as
following. It should be noted that although we take a strongly focused probe
beam (NA=0.15 after all the imaging optics) as the example in this work, the
method is applicable to reconstructing general form of probe beam wavefronts.

  1. (1)

Split pixels: each image of the camera-recorded $I_{\rm p}$-set, as
represented by a $M\times N$ 2D matrix, is split by a factor of $n_{\rm
s}\times n_{\rm s}$ to enhance the spatial resolution $\Delta x=\Delta
x^{(0)}/n_{\rm s}$. Here $\Delta x^{(0)}$ is the physical camera sensor
resolution. For the CCD camera in this work, with $M=1240$, $N=1392$, and
$\Delta x^{(0)}=6.45~{}\mu$m, we find $n_{\rm s}=4$ is large enough to support
the complex $E_{\rm p}$ reconstruction at $NA=0.15$. We choose $n_{\rm s}=6$
with large enough margin in this work.

  2. (2)

Set a “bank” plane: We choose a “bank” plane $z_{\rm b}$ and the associated
probe beam intensity profile $I_{\rm b}$ from the $\\{I_{\rm p}\\}$-set. The
phase of the wavefront $E_{\rm b}$ is initialized as a spherical wave phase
$\phi_{\rm p}^{\rm sw}=k\sqrt{x^{2}+y^{2}+z_{\rm b}^{2}}$. To do this, we have
assumed $z=0$ for the $E_{\rm p}$ focal plane, which is estimated by
additional camera images near $z=0$. An initial guess of $E_{\rm p}$ as
$E_{\rm p}^{GS}=\sqrt{I_{\rm b}}e^{i\phi_{\rm p}^{\rm sw}}$ starts the G-S
iteration in the next step. With $z_{\rm b}$ chosen as the “bank” plane,
$I_{\rm b}$ is subsequently removed from the $\\{I_{\rm p}\\}$ data-set.

  3. (3)

Three images from the $\\{I_{\rm p}\\}$-set are randomly chosen with intensity
distribution $I_{\rm p,i},I_{\rm p,j},I_{\rm p,k}$ recorded at
$z_{i},z_{j},z_{k}$. With the angular spectral method by Eq. (6), $E_{\rm
p}^{\rm GS}$ is first propagated to $z_{i}$ as $E_{\rm p}^{\rm
GS}\leftarrow\hat{U}(z_{i}-z_{\rm b})E_{\rm p}^{\rm GS}$.

  4. (4)

The following G-S iteration among the $i,j,k$ planes are applied for $n_{\rm
A}$ times as following (Fig. 5b): Cyclically propagate $E_{\rm p}^{\rm GS}$
from $z_{i}$ to $z_{j}$, $z_{k}$ and back to $z_{i}$. In each step and before
the next numerical propagation, the amplitude of $E_{\rm p}^{\rm GS}$ is
replaced with the camera-recorded $\sqrt{I_{\rm p}}$ while keeping the phase
information. We usually set $n_{\rm A}=30$.

  5. (5)

With the $z_{i,j,k}$ values initially estimated from the linear stage
readouts, the $z_{i,j,k}$ optimization procedures are sparsely inserted into
the G-S iteration in step 4. For example, when $E_{\rm p}^{\rm GS}$ is
propagated from $z_{i}$ to $z_{j}$, the value $z_{j}$ is adjusted so as to
optimize the inner product between the $|E_{\rm p}^{\rm GS}(z_{j})|^{2}$ and
$I_{\rm p,j}$ images. The corrections to $z_{i,j,k}$ are updated in the
following iterations.

  6. (6)

Add to the “bank” wavefront: Propagate $E_{\rm p}^{\rm GS}$ back to the “bank”
plane $z_{\rm b}$. Add $E_{\rm p}^{\rm GS}$ to the “bank” waveform $E_{\rm
p}^{\rm Bank}$ with an exponential weight: $E_{\rm p}^{\rm
Bank}\leftarrow\frac{1}{1+\eta}(E_{\rm p}^{\rm Bank}+\eta E_{\rm p}^{\rm
GS})$. Here $\eta$ is a weight factor that determine the convergence speed and
stability. After the bank-addition procedure, the G-S wavefront for the
iteration is updated with the best estimated value at the “bank” location:
$E_{\rm p}^{\rm GS}\leftarrow E_{\rm p}^{\rm Bank}$.

  7. (7)

Repeat the (3)-(6) loop $n_{\rm B}$ times until the reconstruction similarity
$s=\langle\Delta I_{\rm t}|\delta|E_{\rm p}^{\rm GS}(z_{\rm t})|^{2}\rangle$,
an “inner product” between the reconstructed $|E_{\rm p}^{\rm GS}(z_{\rm
t})|^{2}$ and $I_{\rm t}$ profiles converge. Here $I_{\rm t}$ the intensity
distribution at a target plane far away from the planes for the $\\{I_{\rm
p}\\}$ set. The reconstructed reference wavefront is finally obtained at
$z=z_{\rm b}$ as $E_{\rm p}=E_{\rm p}^{\rm GS}$, and is now free to propagate
to any other $z$-plane.

Here, with $\Delta I_{\rm A}=I_{\rm A}-\bar{I}_{\rm A}$ and $\Delta I_{\rm
B}=I_{\rm B}-\bar{I}_{\rm B}$ the similarity is defined as

$\langle\Delta I_{\rm A}|\Delta I_{\rm B}\rangle=\frac{\overline{\Delta I_{\rm
A}\Delta I_{\rm B}}}{\sqrt{\overline{\Delta I_{\rm
A}^{2}}}\sqrt{\overline{\Delta I_{\rm B}^{2}}}}.$ (7)

In Fig. 6 we plot the evolution of the similarity $s$ during a typical $E_{\rm
p}$ reconstruction, together with the reconstructed wavefront phase $\phi_{\rm
p}$ at typical steps of the G-S iteration. The G-S process typically take
$n_{\rm b}=1000$ steps toward full convergence. The final intensity similarity
$s\approx 0.99$ suggests a similar level of fidelity for the reconstructed
wavefront $E_{\rm p}=\sqrt{I_{\rm p}}e^{i\phi_{\rm p}}$. The slow iteration is
a result of low-efficiency correction to the large scale wavefront
aberrations. In ref. [58] an optimization procedure that substantially speeds
up the process, as well as a full characterization of the wavefront
reconstruction fidelity will be discussed.

Figure 6: Evolution of the reconstruction fidelity $f$. Typical phase
difference between iterated wavefront and the spherical wave $\delta\phi_{\rm
p}=\phi_{\rm p}-\phi_{\rm p}^{\rm sw}$ (Fig. a-e) are also presented. The
final wavefront

####  VI.1.2 Hologram recording, reduction and $E_{\rm p}$ matching

In this step, holograms of atomic sample are recorded as $I_{1}$ and $I_{2}$,
with and without the atomic sample respectively as outlined in the main text
(Fig. 1). To suppress the impact of vibration noises during the $I_{1,2}$
recording, we minimize the inter-frame time to be merely 10 ms just to allow
the atomic sample to escape the probe area. We then optimally subtract away
the ambient background recorded by additional background measurements,
following a procedure detailed in ref. [20]. The procedure allows us to obtain
the reduced hologram $\Delta I=|E_{\rm p}+E_{\rm s}|^{2}-|E_{\rm p}|^{2}$ for
the $E_{\rm s}$ retrieval in the next step. To account for drifts of optics
alignments, We perform a 3D numerical shift of $E_{\rm p}$ to optimally match
the pre-characterized $|E_{\rm p}(z_{\rm H})|^{2}$ with $I_{2}$.

####  VI.1.3 $E_{\rm s}$ retrieval

We now divide the reduced hologram $\Delta I$ with the optimally-matched
$E_{\rm p}^{*}$ to obtain our initial guess of $E_{\rm s}$,

$E_{\rm diff}=\frac{\Delta I}{{E_{\rm p}}^{*}}=E_{\rm s}+\frac{{E_{\rm
s}}^{*}E_{\rm p}}{{E_{\rm p}}^{*}}+\frac{|E_{\rm s}|^{2}}{{E_{\rm p}}^{*}}$
(8)

Figure 7: Schematic illustration of the $E_{\rm s}$, twin image, DC term. Here
$E_{\rm s}$ is focused at atomic plane $z=z_{\rm A}$. DC and twin images are
focused approximately at $z=0$ and $z=-z_{\rm A}$ respectively. To remove the
twin and dc noise, a $\hat{P}$ operation can be introduced by picking up the
$E_{\rm s}$ signal with concrete support.

As illustrated in Fig. 7, $E_{\rm diff}$ contains information of $E_{\rm s}$
that properly focuses at the atomic sample plane $z=z_{\rm A}$. In addition,
the twin-image term $E_{\rm s}^{*}E_{\rm p}/E_{\rm p}^{*}$ focuses
approximately at $z=-z_{\rm A}$, while the dc noise $|E_{\rm s}|^{2}/E_{\rm
p}^{*}$ focuses at $z=0$. Due to the inline geometry, these additional noise
terms overlap with $E_{\rm s}$ and there is no simple method to isolate
$E_{\rm s}$ without a prior knowledge. Here we employ minimal amount of a
prior knowledge on the known location of the atomic sample, and use an
iterative algorithm [20] to remove the twin and dc noises and to faithfully
retrieve $E_{\rm s}$ from $E_{\rm diff}$, as following:

First, we define an aperture operator $\hat{P}$, which sets $\hat{P}E(z_{\rm
A})=E(z_{\rm A})$ inside the aperture and $\hat{P}E(z_{\rm A})=0$ outside the
aperture (Fig. 7). With the concrete support of the $E_{\rm s}$ signal by the
a prior location information, $\hat{P}E_{\rm s}(z_{\rm A})=E_{\rm s}(z_{\rm
A})$, $\hat{P}E_{\rm diff}(z_{\rm A})$ picks up all the $E_{\rm s}$ but
removes part of the twin and dc noises. To iteratively improve the isolation
of the $E_{\rm s}$, we take advantage of the symmetry between $E_{\rm s}$ and
the twin-image terms to introduce a complex conjugation operator $\hat{C}$ at
$z=z_{\rm H}$ plane as [20],

$\hat{C}E(z_{\rm H})=\frac{E^{*}(z_{\rm H})E_{\rm p}(z_{\rm H})}{E_{\rm
p}^{*}(z_{\rm H})}.$ (9)

A iteration procedure then follows

$E_{\rm s}^{(n+1)}=\hat{U}\hat{C}(E_{\rm diff}-\hat{U}^{-1}\hat{P}E_{\rm
s}^{(n)})$ (10)

to remove the twin and dc noise.

The initial value of of the $E_{\rm s}^{(n)}$ is set as $E_{\rm
s}^{(n=0)}=\hat{U}E_{\rm diff}$. As discussed in ref. [20, 43], $E_{\rm
s}^{(n)}$ converges to $E_{\rm s}$ with a speed determined by the fraction of
light field energy removed by $\hat{P}$ at each step. Unlike the work in ref
[20] where a complex $E_{\rm p}$ composed of two point sources are introduced,
here we simply use a nearly spherical $E_{\rm p}$. Fast convergence of the
iteration is guaranteed by a small atomic sample with size
$\sigma^{2}\ll\lambda z_{\rm A}$ to support a relatively small $\hat{P}$
$E_{\rm s}$ aperture. On the other hand, to image large atomic sample with
$\sigma^{2}\sim\lambda z_{\rm A}$ we should increase the $E_{\rm p}$
complexity accordingly as those in ref. [20].

Finally, after the convergence of iteration by Eq. (10), the removal of the dc
noise is achieved by updating the iteration relation as [20]

$E_{\rm s}^{(n+1)}=\hat{U}\hat{C}(E_{\rm diff}-\hat{U}^{-1}\hat{P}E_{\rm
s}^{(n)}-\frac{|\hat{U}^{-1}\hat{P}E_{\rm s}^{(n)}|^{2}}{{E_{\rm p}}^{*}}).$
(11)

Figure 8: Convergence of the residual energy $\varepsilon^{(n)}$ during the
twin-DC removal iteration. Here $\varepsilon^{(n)}$ is summed over the
displayed imaging area on the top, with the $\hat{P}$ area excluded.

We characterize the convergence of iteration with a residual energy,
$\varepsilon^{(n)}=\sum|(1-\hat{P})E_{\rm s}^{(n)}|^{2}$ summing over an area
substantially larger than the sample size at $z_{\rm A}$. Typical convergence
of the residual is shown in Fig. 8. Notably, the phase-recovery algorithm here
does not assume any special form of reference wavefront $E_{\rm p}$, nor para-
axial approximations that limit the spatial frequency of the wavefront. The
method is therefore naturally compatible with the full wavefront
reconstruction process based on precise $E_{\rm p}$ characterization.

####  VI.1.4 ${\rm OD}$ and $\phi$ retrieval

In this last step, we propagate the reconstructed $E_{\rm p}$ and $E_{\rm s}$
to the plane of atomic sample, $z=z_{\rm A}$, to obtain $\rm OD$ and $\phi$
data as the imaginary and real part of complex phase $\varphi(x,y)=-i{\rm
log}(\frac{E_{\rm p}+E_{\rm s}}{E_{\rm p}})$. Here, with full knowledge of
$E_{\rm p}$ and $E_{\rm s}$, various coherent imaging-processing steps can be
taken to suppress imaging aberration and speckle noises. For example, to
reduce the impact of speckle noises in $E_{\rm p}$ which are introduced by the
imaging system itself (not seen by the atoms), we first propagate $E_{\rm p}$
to $z=0$ to limit its value within the spread of $E_{\rm s}$, before
propagating the filtered $E_{\rm p}$ back to the $z=z_{\rm A}$ plane for the
atomic response retrieval. It is worth noting these speckle noises should
equally affect signals in regular absorption images where a similar correction
is difficult to make.

###  VI.2 Details on the experimental setup and data analysis

####  VI.2.1 Fitting the complex EIT spectrum

Figure 9: Schematic of a 3-level system system to obtain the phenomenological
EIT response.

The complex spectrum obtained from the holographic measurements can be fit to
a phenomenological EIT model to assist the data analysis. The fit process is
detailed in this section. Here, instead of a full analysis of EIT for the
multi-level atom in presence of residual magnetic field, we use a simple
3-level model (Fig. 9) and perform a linear analysis of the atomic response to
$E_{\rm p}$, leading to the standard expression of the atomic polarizability,

$\alpha(\omega_{\mathrm{p}},\omega_{\mathrm{c}})=\alpha_{0}+\frac{C}{\frac{|\Omega_{\mathrm{c}}|^{2}}{4(\Delta_{\mathrm{p}}-\Delta_{\mathrm{c}}-i\kappa/2)}-(\Delta_{\rm
p}+i\Gamma/2)}.$ (12)

Here the parameter $C$ determines the $|g\rangle\leftrightarrow|e\rangle$
effective transition strength. The
$\Delta_{\mathrm{p}}=\omega_{\mathrm{p}}-\omega_{\mathrm{eg}}$ and
$\Delta_{\mathrm{c}}=\omega_{\mathrm{c}}-\omega_{\mathrm{ea}}$ are the probe
and control laser detunings respectively. $\kappa$ is a phenomenological
damping constant to account for relaxation between the ground state
$|g\rangle$ and $|a\rangle$. A non-zero $\alpha_{0}$ is expected from real
atomic response due to non-ideal optical pumping, additional magnetic
coupling, as well as off-resonant coupling of $E_{\mathrm{p}}$ to additional
hyperfine transitions.

For the dilute and optically thin atomic sample in this work, the complex
phase shift $\varphi(\omega_{\mathrm{p}},\omega_{\mathrm{c}})$ obtained from
the holography is expected to be proportional to
$\alpha(\omega_{\mathrm{p}},\omega_{\mathrm{c}})$. We thus apply Eq. (12) to
fit the experimental complex phase from the $\mathrm{OD}$ and $\phi$ data in
Figs. 2, 3, 4 of the main text for the case of stationary atoms. The free
parameters
$\alpha_{0},C,|\Omega_{\mathrm{c}}|^{2},\kappa,\Gamma,\Delta_{\mathrm{c}}$ for
the Fig. 2 fit are given by: $\alpha_{0}=-0.026+0.38i$, $C=8.7$,
$\Delta_{\mathrm{c}}/(2\pi)=-1.1$MHz, $\Gamma=13.7/(2\pi)$ MHz,
$|\Omega_{\mathrm{c}}|^{2}/(2\pi)^{2}=11.7$ MHz2, $\kappa/(2\pi)=-0.79$ MHz.

####  VI.2.2 Photon shot noise limit to the $\overline{\beta}$ measurement

In the main text we considered optically thin samples which induce complex
phase shifts $\varphi\ll 1$, so that a measurement of phase angle $\beta={\rm
arg}(\varphi)$ is equivalent to measuring the relative phase between $E_{\rm
s}$ and $E_{\rm p}$, leading to $n_{\bar{\beta}}>1/\sqrt{N_{\rm s}}$ shot
noise limit through an imaging aperture $\mathcal{A}$ in the atomic sample
plane.

Here we derive the more general bound $n_{\bar{\beta}}>\xi/\sqrt{N_{\rm s}}$
as in the caption of Fig. 3 in the main text. In particular, from the
definition of the complex phase $\varphi=-i{\rm log}(1+E_{\rm s}/E_{\rm p})$,
we have $\delta\overline{\varphi}=\delta\overline{E_{\rm s}}/(\overline{E_{\rm
s}}+\overline{E_{\rm p}})$, with $\delta\overline{\varphi}$,
$\delta\overline{E_{\rm s}}$ evaluated as the rms value of the complex number
within an imaging area as those in Fig. 3a of the main text. When the imaging
noise is dominated by the probe photon shot noise, we have
$n_{\overline{\beta}}=|\frac{\delta\overline{\varphi}}{\sqrt{2}\overline{\varphi}}|$,
and further $|\frac{\delta\overline{E_{\rm s}}}{\sqrt{2}\overline{E_{\rm
s}}}|>1/\sqrt{N_{\rm s}}$. Following
$\frac{\delta\varphi}{\varphi}=\frac{\delta\overline{E_{\rm
s}}}{\overline{E_{\rm
s}}}\frac{1-e^{-i\overline{\varphi}}}{\overline{\varphi}}$ we arrive at
$n_{\overline{\beta}}>\xi/\sqrt{N_{\rm s}}$ with a optical-depth dependent
factor $\xi=|\frac{1-e^{-i\overline{\varphi}}}{\overline{\varphi}}|<1$.

####  VI.2.3 Inference of $\rho_{c}-u$ distribution

This section describes the procedure to infer atomic density $\rho_{c}$ and
velocity field $u$ from the holographically reconstructed complex $\varphi$
data such as those for Fig. 4 in the main text.

For the model based $\rho_{\rm c}-u$ inference leading to Fig. 4(vi), the
numerical method parametrizes $\rho_{c}$ and $u$ fields to efficiently
minimize the cost function $\mathcal{L}=\sum_{x,y}|2\varphi(x,y)-k_{\rm
p}\rho_{c}(x,y)\alpha(\omega-k_{u}u(x,y))|^{2}$ where $\alpha(\omega)$ is
obtained by fitting the EIT-data as described in Sec. VI.2.1. Generally,
assumptions are tailored into the parameters for specific type of sample
characterizations. Here, taking advantage of a prior knowledge that the atomic
density distribution $\rho_{\rm c}$ is localized and smooth, we decompose
$\rho_{\rm c}$ into $N_{\rm G}$ Gaussian packets as $\rho_{\rm
c}=\sum_{j=1}^{N_{\rm
G}}C_{j}e^{-(x-X_{j})^{2}/A_{j}^{2}-(y-Y_{j})^{2}/B_{j}^{2}}$ parametrized by
parameters $\\{C_{j},X_{j},Y_{j},A_{j},B_{j}\\}$. On the other hand, the
velocity field $u(x,y)$ follows a polynomial expansion with
$u(x,y)=\sum_{m,n=0}^{M,N}c_{m,n}x^{m}y^{n}$ parametrized by $\\{c_{m,n}\\}$.
With in mind the limited signal/noise in single $\varphi(x,y)$ data, we
restrict the parameters to $m=0$ and $n=0,1,2$ during the numerical
optimization as following.

The optimization procedure starts with $N_{\rm G}=1$ with quasi-Newton methods
provided by Matlab (fminunc). After the $\mathcal{L}$ minimization converge,
$N_{\rm G}$ is increased by one, with the new Gaussian packet inheriting the
optimal parameters of the last Gaussian, except for a randomly shifted
$X_{j},Y_{j}$ combined with halved $A_{j},B_{j}$ for a next round of nonlinear
optimization. The process of “splitting the last Guassian” is continued up to
$N_{\rm G}=7$ to obtain the model fit in Fig. 4(vi) in the main text, beyond
which we find over-fitting of artifacts from apparent noises occurs.

For the time-of-flight method leading to Fig. 4(vi), we assume that the atomic
position recorded during the $\tau_{\rm tof}<t<\tau_{\rm tof}+t_{\rm p}$
camera exposure time are perfectly correlated to their velocity through the
point-source expansion. By setting $y=0$ for the unperturbed atomic position
(Fig. 4(iii)), the velocity field is estimated as $u(y)=y/\tau_{\rm mean}$,
with $\tau_{\rm mean}=\tau_{\rm tof}+1/2\tau_{\rm p}$ to be the average free-
flight time. We then derive $\rho_{\rm c}(x,y)=2|\varphi(x,y)|/|k_{\rm
p}\alpha(\omega-k_{u}u(x,y))|$ from the $\varphi(x,y)$ data.

## References

  * Einstein [1905] A. Einstein, On the Electrodynamics of Moving Bodies, Ann. Phys. 17, 891 (1905).
  * Gordon [1923] B. W. Gordon, On the propagation of light in the theory of relativity, Ann. Phys. 72, 421 (1923).
  * Strekalov  _et al._ [2004] D. Strekalov, A. B. Matsko, N. Yu, and L. Maleki, Observation of Light Dragging in a Rubidium Vapor Cell, Phys. Rev. Lett. 93, 023601 (2004).
  * Safari  _et al._ [2016] A. Safari, I. D. Leon, M. Mirhosseini, O. S. Magaña-loaiza, and R. W. Boyd, Light-Drag Enhancement by a Highly Dispersive Rubidium Vapor, Phys. Rev. Lett. 116, 013601 (2016).
  * Kuan  _et al._ [2016] P. C. Kuan, C. Huang, W. S. Chan, S. Kosen, and S. Y. Lan, Large Fizeau’s light-dragging effect in a moving electromagnetically induced transparent medium, Nat. Commun. 7, 13030 (2016).
  * Chen  _et al._ [2020] Z. Chen, H. M. Lim, C. Huang, R. Dumke, and S. Y. Lan, Quantum-Enhanced Velocimetry with Doppler-Broadened Atomic Vapor, Phys. Rev. Lett. 124, 93202 (2020), arXiv:2003.03491 .
  * Solomons  _et al._ [2020] Y. Solomons, C. Banerjee, S. Smartsev, J. Friedman, D. Eger, O. Firstenberg, and N. Davidson, Transverse drag of slow light in moving atomic vapor, Opt. Lett. 45, 3431 (2020), arXiv:2004.04414 .
  * Fleischhauer  _et al._ [2005] M. Fleischhauer, A. Imamoglu, and P. J. Marangos, Electromagnetically induced transparency, Rev. Mod. Phys. 77, 633 (2005).
  * Novikova  _et al._ [2012] I. Novikova, R. L. Walsworth, and Y. Xiao, Electromagnetically induced transparency-based slow and stored light in warm atoms, Laser Photonics Rev. 6, 333 (2012).
  * Vestergaard Hau  _et al._ [1999] L. Vestergaard Hau, S. E. Harris, Z. Dutton, and C. H. Behroozi, Light speed reduction to 17 metres per second in an ultracold atomic gas, Nature 397, 594 (1999).
  * Leonhardt [2002] U. Leonhardt, A laboratory analogue of the event horizon using slow light in an atomic medium, Nature 415, 406 (2002).
  * Rosenberg [2020] Y. Rosenberg, Optical analogues of black-hole horizons, Phil. Trans. R. Soc. A 378, 20190232 (2020), arXiv:2002.04216 .
  * Matsko  _et al._ [2003] A. B. Matsko, N. Yu, and L. Maleki, Gravity field measurements using cold atoms with direct optical readout, Phys. Rev. A 67, 043819 (2003).
  * Zimmer and Fleischhauer [2004] F. Zimmer and M. Fleischhauer, Sagnac interferometry based on ultraslow polaritons in cold atomic vapors, Phys. Rev. Lett. 92, 253201 (2004).
  * Haine  _et al._ [2015] S. A. Haine, S. S. Szigeti, M. D. Lang, and C. M. Caves, Heisenberg-limited metrology with information recycling, Phys. Rev. A 91, 041802 (2015), arXiv:1411.5111 .
  * Dickerson  _et al._ [2013] S. M. Dickerson, J. M. Hogan, A. Sugarbaker, D. M. Johnson, and M. A. Kasevich, Multiaxis inertial sensing with long-time point source atom interferometry, Phys. Rev. Lett. 111, 083001 (2013), arXiv:1305.1700 .
  * Engels  _et al._ [2003] P. Engels, I. Coddington, P. C. Haljan, V. Schweikhard, and E. A. Cornell, Observation of Long-Lived Vortex Aggregates in Rapidly Rotating Bose-Einstein Condensates, Phys. Rev. Lett. 90, 170405 (2003), arXiv:0301532 [cond-mat] .
  * Wright  _et al._ [2013] K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, Driving phase slips in a superfluid atom circuit with a rotating weak link, Phys. Rev. Lett. 110, 025302 (2013).
  * Galitski  _et al._ [2019] V. Galitski, G. Juzeliunas, and I. B. Spielman, Artificial gauge fields with ultracold atoms, Phys. Today 72, 38 (2019).
  * Sobol and Wu [2014] J. P. Sobol and S. Wu, Imaging cold atoms with shot-noise and diffraction limited holography, New J. Phys. 16, 093064 (2014).
  * Cuche  _et al._ [1999] E. Cuche, P. Marquet, and C. Depeursinge, Simultaneous amplitude-contrast and quantitative phase-contrast microscopy by numerical reconstruction of Fresnel off-axis holograms, Appl. Opt. 38, 6994 (1999).
  * Coppola  _et al._ [2004] G. Coppola, P. Ferraro, M. Iodice, S. D. Nicola, A. Finizio, and S. Grilli, A digital holographic microscope for complete characterization of microelectromechanical systems, Measurement Science and Technology 15, 529 (2004).
  * Kemper  _et al._ [2005] B. Kemper, D. Carl, J. Schnekenburger, I. Bredebusch, M. Schäfer, W. Domschke, and G. von Bally, Investigation of living pancreas tumor cells by digital holographic microscopy, Journal of Biomedical Optics 11, 34005 (2005).
  * Khmaladze  _et al._ [2008] A. Khmaladze, M. Kim, and C. M. Lo, Phase imaging of cells by simultaneous dual-wavelength reflection digital holography, Optics Express 16, 10900 (2008).
  * Greenbaum  _et al._ [2012] A. Greenbaum, W. Luo, T. W. Su, Z. Göröcs, L. Xue, S. O. Isikman, A. F. Coskun, O. Mudanyali, and A. Ozcan, Imaging without lenses: Achievements and remaining challenges of wide-field on-chip microscopy, Nat. Methods 9, 889 (2012).
  * De Haan  _et al._ [2020] K. De Haan, Y. Rivenson, Y. Wu, and A. Ozcan, Deep-Learning-Based Image Reconstruction and Enhancement in Optical Microscopy, Proc. IEEE 108, 2949575 (2020).
  * Kadlecek  _et al._ [2001] S. Kadlecek, J. Sebby, R. Newell, and T. G. Walker, Nondestructive spatial heterodyne imaging of cold atoms, Opt. Lett. 26, 137 (2001).
  * Turner  _et al._ [2005] L. D. Turner, K. F. E. M. Domen, and R. E. Scholten, Diffraction-contrast imaging of cold atoms, Phys. Rev. A 72, 031403 (2005).
  * Smits  _et al._ [2020] J. Smits, A. P. Mosk, and P. van der Straten, Imaging trapped quantum gases by off-axis holography, Opt. Lett. 45, 981 (2020), arXiv:1911.09590 .
  * Altuntas and Spielman [2021] E. Altuntas and I. B. Spielman, Self-Bayesian Aberration Removal via Constraints for Ultracold Atom Microscopy, Phys. Rev. Res 3, 043087 (2021), arXiv:2108.07106 .
  * Gabor [1972] D. Gabor, Holography, 1948-1971, Science 177, 299 (1972).
  * Gerchberg and Saxton [1972] R. W. Gerchberg and W. O. Saxton, A Practical Algorithm for the Determination of Phase from Image and Diffraction Plane Pictures, Optik 35, 237 (1972).
  * Ivanov  _et al._ [1992] V. Y. Ivanov, V. P. Sivokon, and M. A. Vorontsov, Phase retrieval from a set of intensity measurements : theory and experiment, J. Opt. Soc. Am. A 9, 1515 (1992).
  * Latychevskaia [2019] T. Latychevskaia, Phase retrieval for digital holography, J. Opt. Soc. Am. A 36, 31 (2019).
  * Gomez  _et al._ [2006] E. Gomez, L. A. Orozco, and G. D. Sprouse, Spectroscopy with trapped francium: Advances and perspectives for weak interaction studies, Reports Prog. Phys. 69, 79 (2006).
  * Lu  _et al._ [2013] Z. Lu, P. Mueller, G. W. F. Drake, and S. C. Pieper, Colloquium : Laser probing of neutron-rich nuclei in light atoms, Rev. Mod. Phys. 85, 1383 (2013).
  * Marti  _et al._ [2018] G. E. Marti, R. B. Hutson, A. Goban, S. L. Campbell, N. Poli, and J. Ye, Imaging Optical Frequencies with 100 muHz Precision and 1 . 1 mu m Resolution, Phys. Rev. Lett. 120, 103201 (2018).
  * Li  _et al._ [2020] R. Li, Y. Wu, Y. Rui, B. Li, Y. Jiang, L. Ma, and H. Wu, Absolute Frequency Measurement of Li 6 D Lines with khz-Level Uncertainty, Physical Review Letters 124, 063002 (2020).
  * Fleischhauer and Lukin [2002] M. Fleischhauer and M. D. Lukin, Quantum memory for photons: Dark-state polaritons, Phys. Rev. A 65, 022314 (2002).
  * [40] D. A. Steck, Rubidium 87 D Line Data, http://steck.us/alkalidata (revision 2.2.1, 21 November 2019). .
  * [41] See supplementary material on details of experimental setup, holographic imaging process, and data analysis .
  * Kreis and Juptner [1997] T. M. Kreis and W. P. O. Juptner, Suppression of the dc term in digital holography, Opt. Eng. 36, 2357 (1997).
  * Sobol [2014] J. P. Sobol,  _Imaging cold atoms with shot-noise and diffraction limited holography_ , Ph.D. thesis, Swansea University (2014).
  * Pappa  _et al._ [2011] M. Pappa, P. C. Condylis, G. O. Konstantinidis, V. Bolpasi, A. Lazoudis, O. Morizot, D. Sahagun, M. Baker, and W. von Klitzing, Ultra-sensitive atom imaging for matter-wave optics, New J. Phys. 13, 115012 (2011).
  * Li  _et al._ [2018] X. Li, F. Zhou, M. Ke, P. Xu, X.-D. He, J. Wang, and M.-S. Zhan, High-resolution ex vacuo objective for cold atom experiments, Appl. Opt. 57, 7584 (2018).
  * Cronin and Perreault [2004] A. D. Cronin and J. D. Perreault, Phasor analysis of atom diffraction from a rotated material grating, Phys. Rev. A 70, 043607 (2004).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * He  _et al._ [2020] Y. He, L. Ji, Y. Wang, L. Qiu, J. Zhao, Y. Ma, X. Huang, D. E. Chang, and S. Wu, Geometric control of collective spontaneous emission, Phys. Rev. Lett. 125, 213602 (2020), arXiv:1910.02289 .
  * Lye  _et al._ [2003] J. E. Lye, J. J. Hope, and J. D. Close, Nondestructive dynamic detectors for Bose-Einstein condensates, Phys. Rev. A 67, 043609 (2003).
  * Bakr  _et al._ [2009] W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice, Nature 462, 74 (2009), arXiv:0908.0174 .
  * Sherson  _et al._ [2010] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Single-atom-resolved fluorescence imaging of an atomic Mott insulator, Nature 467, 68 (2010), arXiv:1006.3799 .
  * Boll  _et al._ [2016] M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Spin- and density-resolved microscopy of antiferromagnetic correlations in Fermi-Hubbard chains, Science 353, 1257 (2016).
  * Yang  _et al._ [2020] B. Yang, H. Sun, R. Ott, H. Y. Wang, T. V. Zache, J. C. Halimeh, Z. S. Yuan, P. Hauke, and J. W. Pan, Observation of gauge invariance in a 71-site Bose–Hubbard quantum simulator, Nature 587, 392 (2020), arXiv:2003.08945 .
  * Mølmer  _et al._ [1992] K. Mølmer, D. S. Weiss, M. Kasevich, and S. Chu, Theoretical analysis of velocity-selective Raman transitions, Phys. Rev. A 45, 342 (1992).
  * Su  _et al._ [2010] E. Su, S. Wu, and M. Prentiss, Atom interferometry using wave packets with constant spatial displacements, Phys. Rev. A 81, 043631 (2010).
  * Carey  _et al._ [2019] M. Carey, J. Saywell, D. Elcock, M. Belal, and T. Freegarde, Velocimetry of cold atoms by matter-wave interferometry, Phys. Rev. A 99, 023631 (2019), arXiv:1802.02190 .
  * Kim [2010] M. K. Kim, Principles and techniques of digital holographic microscopy, Spie Reviews 1, 018005 (2010).
  * [58] Since the preparation of this work, we have made additional developments to improve the speed of the probe wavefront pre-characterization, on phase-recovery for imaging large atomic sample with a non-iterative method, and on holography-based aberration corrections. The results will be given in a future publication. .