††thanks<EMAIL_ADDRESS>
# Inverse Transform Sampling for Efficient Doppler-Averaged Spectroscopy
Simulations

Andrew P. Rotunno  National Institute of Standards and Technology, Boulder, CO
80305, USA  Amy K. Robinson  Department of Electrical Engineering, University
of Colorado, Boulder, Colorado 80305, USA  Nikunjkumar Prajapati  Samuel
Berweger  Matthew T. Simons  Alexandra B. Artusio-Glimpse  Christopher L.
Holloway  National Institute of Standards and Technology, Boulder, CO 80305,
USA

(August 28, 2024)

###### Abstract

We present a thermal velocity sampling method for calculating Doppler-
broadened atomic spectra, which more efficiently reaches a smooth limit than
regular velocity weighted sampling. The method uses equal-population sampling
of the 1-D thermal distribution, sampling the ‘inverse transform’ of the
cumulative distribution function, and is broadly applicable to normal
distributions. We also discuss efficiencies from eliminating velocity classes
which don’t significantly contribute to observed atomic lines, and comment on
the application of this method in 2- and 3-dimensions.

##  I Introduction

Atomic vapors have recently grown in use as electromagnetic field sensors,
utilizing spectroscopic observation of atomic states for precise, transferable
measurement standards [1, 2, 3]. A major source of uncertainty in
electromagnetically-induced transparency (EIT) spectroscopy of warm Rydberg
atomic vapors comes from Doppler broadening [4], where the thermal
distribution of atomic velocities causes frequency shifts that result in
spectral broadening. Such effects are minimized in cold atomic systems [5], or
by using well-matched EIT schemes [6, 7]. However, Doppler broadening remains
a major effect in room-temperature and warmer vapor cells, broadening the D2
transition in ${}^{133}Cs$ from its decay linewidth of $2\pi\cdot 6$ MHz [8],
to a FWHM an order of magnitude larger at room temperature. To fit
experimental spectra, we must consider this Doppler broadening, and we aim to
increase efficiency when sampling thermal velocity distribution.

We use the illustrative case of two-photon Rydberg EIT to compute Doppler
broadening, where we calculate the EIT spectrum many times over a wide
sampling of velocities, and combine their spectra. The aim of this paper is
introduce a ‘natural’ sampling method for thermal velocity distributions, such
that fewer calculations are needed to realize convergence to the Doppler-
averaged curve, using a straightforward statistical transformation. For the
same number of divisions, we find a ‘population’ sampling method arrives at
smooth transmission curves much faster than ‘velocity’ sampling by being more
dense where the resonant EIT feature occurs, while also sampling into the
‘wings’ of the distribution. The general method of ‘inverse population
sampling’ has been used in various contexts. Here we apply it to thermal
sampling in atomic spectroscopy.

In Sec. II, we describe the relevant atom-photon theory, including the
Maxwell-Boltzmann distribution, Doppler shifts, and atomic EIT. In Sec. III,
we describe the proposed sampling method using the inverse error function. In
Sec. IV, we compare convergence rates between population and velocity sampling
methods. In the appendices, we address a few special cases and extensions of
the core paper. In Appx. A we examine using a velocity cut-off for efficient
line-shape calculations, in Appx. B, we discuss the application of this method
to 2D and 3D, and in Appx. C, we confirm the new method recovers additional
Doppler effects observed in Rydberg EIT. We review the optical master equation
model in Appx. D.

##  II Background Theory

This section touches on relevant background theory including Doppler shifts,
thermal distributions, and Doppler averaging transmission for Rydberg state
EIT, with more details on transmission given in Appx. D.

Figure 1:  Illustration of two-photon Doppler shifts for 133Cs Rydberg EIT.
Photon arrow lengths illustrate observed wavelength shifts.

###  II.1 Thermal Velocities

A population of atoms in thermal equilibrium (as in an alkali vapor cell),
have 3-D velocities v that are well-described by the Maxwell-Boltzmann
distribution $f_{MB}^{3D}(\textbf{v})$. This distribution gives the velocity
density distribution as a function of temperature $T$, particle mass $m$,
using the Boltzmann constant $k_{B}$ [9]:

$f_{MB}^{3D}(\textbf{v})dv^{3}=\left(\frac{m}{2\pi{}k_{B}T}\right)^{\frac{3}{2}}\exp\left(-\frac{m|v|^{2}}{2k_{B}T}\right)dv^{3}\,\,\,.$
(1)

While across three dimensions, atoms have most likely total speed
$\sqrt{\frac{2k_{B}T}{m}}$, we look at only one dimension, per the Doppler
shift’s directional selectivity. Using $v_{\sigma}=\sqrt{k_{B}T/m}$, we have
for any one dimension:

$f_{MB}(v)dv=\frac{1}{v_{\sigma}\sqrt{2\pi}}\exp\left(-\frac{v^{2}}{2v_{\sigma}^{2}}\right)dv\,\,\,.$
(2)

This distribution of velocities is shown in Fig. 2(a). Velocities in any one
dimension are normally distributed around 0 velocity, with standard deviation
$v_{\sigma}$.

###  II.2 Doppler Shifts

The Doppler effect is an apparent shift in frequency of a wave when an
observer or source are moving relative to one another [10]. The observed
frequency shift $\Delta{}f\equiv{}f^{\prime}-f_{0}$ from a stationary source
with frequency $f_{0}$ is
$\Delta{}f=-\frac{v}{c}f_{0}=-\frac{v}{\lambda_{0}}$, where $v$ is the speed
of the observer in the wave propagation direction, $c$ the wave velocity
(here, the speed of light), and wavelength $\lambda_{0}=c/f_{0}$. We also
write this Doppler shift as $\Delta\omega=-\textbf{v}\cdot\textbf{k}$ using
angular frequency $\omega=2\pi f$, the photon’s wavevector
$|\textbf{k}|=2\pi/\lambda_{0}$, and the dot product to determine alignment
with an atom’s 3-D velocity v. The factor $v/c$ for typical atoms at room
temperature is of order $10^{-6}$, but optical transitions of order $10^{14}$
Hz leave shifts at the $10^{8}$ Hz, nearly of order GHz, far wider than laser
or decay linewidths.

###  II.3 Transmission Spectrum

In single-photon absorption, this broadening is well-characterized by a Voigt
profile, the combination of Lorentzian (atomic) and Gaussian (thermal)
linewidths. Some limits such as the ‘weak probe’ approximation allow closed-
form expressions for multi-photon EIT [11, 12], so we discuss here the general
case with strong probe and coupling lasers, when ‘brute force’ parameter scan
calculations are required. We examine a particular absorption Doppler feature
of the weak probe case in Appx. C. In Appx. D, we detail of the calculation of
two-photon spectroscopy in Rydberg EIT [1, 13]. Here, we discuss the parts
relevant to the thermal focus of this paper.

We consider a three-level, two-photon system, with the ‘probe’ laser driving
the ground state to an intermediate state at 852
nm$=\lambda_{p}=2\pi/\textbf{k}_{p}$ (detuned by $\Delta_{p}$), and the
‘coupling’ laser which connects the intermediate state with a Rydberg state at
509 nm$=\lambda_{c}=2\pi/\textbf{k}_{c}$ (detuned from the atomic resonance by
$\Delta_{c}$). Given an atomic velocity v, the observed detuning
$\Delta^{\prime}$ of the probe and coupling are shifted from the lab-frame
values $\Delta$:
$\Delta_{p}^{\prime}=\Delta_{p}-\textbf{v}\cdot{}\textbf{k}_{p}$ and
$\Delta_{c}^{\prime}=\Delta_{c}-\textbf{v}\cdot{}\textbf{k}_{c}$. When
$\textbf{k}_{p}$ and $\textbf{k}_{c}$ counter-propagate, we can write the
observed detunings in terms of one-dimensional velocity:
$\Delta_{p}^{\prime}=\Delta_{p}-v\frac{2\pi}{\lambda_{p}}$ and
$\Delta_{c}^{\prime}=\Delta_{c}+v\frac{2\pi}{\lambda_{c}}$. Extension to
higher dimensions is discussed in Appx. B. Holding $\Delta_{p}=0$, any $v\neq
0$ will require a shift in experimental detuning to
$\Delta_{c}=v(\textbf{k}_{p}-\textbf{k}_{c})$ to bring the EIT back to
resonance. Note the vast frequency/wavelength difference in this two-photon
scheme is the source of apparent detuning when atoms have $v\neq 0$. In
effect, by scanning $\Delta_{c}$ and holding $\Delta_{p}=0$, we observe many
velocity classes by the detuning of their EIT peak when brought into
resonance. These constituent velocity-class transmission curves are given
later in Sec. IV. Given the large mis-match between $\lambda_{p}$ and
$\lambda_{c}$, only a narrow range of velocities contribute to the Doppler-
broadened EIT lineshape.

As laid out in Appx. D, the atomic coherence $\rho_{12}(v,\Delta_{c})$
determines probe transmission spectrum across $\Delta_{c}$ for a particular
$v$. Calculation of steady-state $\rho_{12}$ is the most time-intensive step,
as it must be computed multiple times over a range of $v$ in order to combine
them. We integrate $\rho_{12}$ over velocities, weighted by the Maxwell-
Boltzmann distribution, yielding the Doppler-broadened spectrum:

$\rho_{12,D}(\Delta_{c})=\int_{-\infty}^{\infty}\rho_{12}\left(v,\Delta_{c}\right)f_{MB}(v)\leavevmode\nobreak\
dv\,\,\,.$ (3)

We seek to minimize the number of calculations required to determine this
Doppler-averaged transmission curve, which is is the primary topic of this
paper.

##  III Velocity Sampling

The most straightforward approach to evaluate the velocity-space integral of
Eq. 3 is by constructing a discrete equally-spaced ‘scan’ of velocities over a
few standard deviations, to cover nearly all of the atomic ensemble’s
exhibited velocities. We refer to this as the “velocity sampling” method, and
we plot a representative sampling in Fig. 2(a). Each bin of velocity is
weighted by the Maxwell-Boltzmann distribution to determine what fraction of
the population is between $v$ and $v+dv$. This renders Eq. 3 discretized for
computation as:

$\rho_{12,D}(\Delta_{c})=\frac{\Delta
v_{i}}{v_{\sigma}\sqrt{2\pi}}\sum_{v_{i}=-3v_{\sigma}}^{3v_{\sigma}}\rho_{12}\left(v_{i},\Delta_{c}\right)\exp\left(-\frac{v_{i}^{2}}{2v_{\sigma}^{2}}\right)$
(4)

where common factors are brought outside the sum. In Eq. 4, the sample $v_{i}$
spans three standard deviations, with equal spacing $\Delta{}v_{i}$. For a
symmetric summation over $N_{s}$-$\sigma$ (where $N_{s}=3$ in the above
example) using $N_{v}$ points, we have velocity resolution
$\Delta{}v_{i}=2N_{s}v_{\sigma}/N_{v}$ across the entire range. An even
velocity sampling means that in order to sample densely near $v_{i}=0$, one
must also densely sample the distribution’s ‘wings,’ where the contribution to
a resonant interaction is minimal, both by being off-resonance and diminished
by the Maxwell-Boltzmann density.

###  III.1 Local Refinement

To significant effect in computational efficiency, we can break the summation
of Eq. 4 into parts: one with ‘fine’ velocity spacing at the region of
interest, the rest remains ‘coarse.’ This fine segment can be centered around
arbitrary velocities ($v\neq 0$), in the event that detuning in either field
requires an off-resonant calculation. Imposing bounds at $v_{\downarrow}$ and
$v_{\uparrow}$, and using ‘fine’ and ‘coarse’ velocity spacing
$\Delta{}v_{f}\ll\Delta{}v_{c}$, we can split Eq. 4 into partial sums with
different sample spacing and bounds:

$\begin{split}\rho_{12,D}(\Delta_{c})=&\\\ \leavevmode\nobreak\ \frac{\Delta
v_{c}}{v_{\sigma}\sqrt{2\pi}}&\sum_{v_{i}=v_{\uparrow}}^{3v_{\sigma}}\rho_{12}\left(v_{i},\Delta_{c}\right)e^{\frac{-v_{i}^{2}}{2v_{\sigma}^{2}}}\\\
+\frac{\Delta
v_{f}}{v_{\sigma}\sqrt{2\pi}}&\sum_{v_{i}=v_{\downarrow}}^{v_{\uparrow}}\rho_{12}\left(v_{i},\Delta_{c}\right)e^{\frac{-v_{i}^{2}}{2v_{\sigma}^{2}}}\\\
+\frac{\Delta
v_{c}}{v_{\sigma}\sqrt{2\pi}}&\sum_{v_{i}=-3v_{\sigma}}^{v_{\downarrow}}\rho_{12}\left(v_{i},\Delta_{c}\right)e^{\frac{-v_{i}^{2}}{2v_{\sigma}^{2}}}\end{split}\,\,\,\,\,.$
(5)

In practice, the bounds $v_{\uparrow}$ and $v_{\downarrow}$ are brought in
until their cutoff is observed in the ‘wings’ of the calculated EIT peak,
typically a fraction of $1\sigma$, and the coarse velocity spacing accounts
for the rest of the overall absorption. We must re-weight these for
combination by $\Delta{}v_{f}$ and $\Delta{}v_{c}$, as velocity bin size
scales the Maxwell-Boltzmann density. This splitting into coarse and fine
velocity ranges allows one to focus computational resources on near-resonant
velocity classes, which heavily depend on small changes in $v$, while
economizing on far off-resonant velocity classes for which the absorption
profile is nearly flat.

\begin{overpic}[width=433.62pt]{fig2abc.pdf} \put(14.0,80.0){(a)}
\put(60.0,80.0){(b)} \put(13.0,37.0){(c)} \end{overpic} Figure 2: Contrasting
sampling methods, using $N_{v}=51$ for 133Cs at 300 K, where
$v_{\sigma}\approx 137$ m/s. (a) An equal-velocity sampling uses Gaussian
statistics (3-$\sigma$, Eq. 2), varying population in each sample. (b) An
equal-population partition uses the inverse transform (Eq. 7), varying
velocity density. (c) Plotting both on the same probability scale, showing the
dynamic velocity spacing of the population sampling method.

###  III.2 Population Sampling Method

We are motivated by this _ad hoc_ sampling in the region of interest to
construct a ‘natural’ velocity sampling that over-samples near resonance, and
under-samples in the wings, for computational efficiency in the resonant case.
To this end, we sample velocities that are equally spaced in the Gaussian-
distributed population or percentile space, such that each step represents an
equal population of atoms, $1/N_{v}$ for $N_{v}$ velocity bins. In Fig.
2(b,c), we show this equal-population method’s dynamic velocity sampling,
using $N_{v}=51$ for illustration. This method of percentile sampling is
broadly applicable, often given the name “inverse transform sampling” [14], or
probit analysis [15] when applied to population statistics.

We begin to calculate the velocity of each partition of population using the
error function $\textrm{erf}(v)$ [16]. A Gaussian’s cumulative distribution
function is
$\frac{1}{2}\left(1+\textrm{erf}\left(\frac{v-\mu}{\sqrt{2}\sigma}\right)\right)$,
for mean $\mu$ and standard deviation $\sigma$. The error function spans
$v=(-\infty,\infty)$, and yields a value which integrates signed population
away from the mean, from $0$ to $v$:

$\textrm{erf}(v)\equiv\frac{2}{\sqrt{\pi}}\int_{0}^{v}e^{-t^{2}}dt$ (6)

To perform the inverse transform sampling, we use the inverse of the error
function $\textrm{erf}^{-1}(\eta)$ or $\textrm{inverf}(\eta)$, where $\eta$
spans the signed population space $(-1,1)$ with $N_{v}$ points. The inverse
error function transforms a linearly-spaced (or randomized) population into an
inhomogeneously-spaced distribution of velocity values. The function
$\textrm{erf}^{-1}(\eta)$ does not have a closed form expression, and is
implemented numerically [17].

In pseudo-code, the velocity distribution vector $v_{vec}$ of a Gaussian-
distributed, normalized population is:

$v_{vec}=v_{\sigma}\sqrt{2}\cdot{}\textrm{erf}^{-1}(\eta_{vec})$ (7)

where $\eta_{vec}$ is the discretized equal-spaced population sampling vector
from $-1$ to $1$ with $N_{v}$ elements. We plot $v_{vec}$ against $\eta_{vec}$
in Fig. 2(b). The endpoints $-1,1$ yield infinities, so the full vector
$\eta_{vec}$ is generated with $N_{v}+2$ points, then we use elements $[2]$ to
$[end-1]$, for starting index 1. If $N_{v}$ is odd, the center point 0 is
included. In the limit of high $N_{v}\gtrsim 100$, the gap $dv$ near 0
approaches $v_{\sigma}\cdot\sqrt{2/\pi}/N_{v}$. Since $\sqrt{2/\pi}\approx
0.8$, this sampling is approximately as dense as equal-spaced velocities
spanning only a symmetric $0.4\sigma$, while actually sampling asymptotically
into the distribution wings. While this inverse-transform sampling method
works well for resonant EIT, it is less efficient when significant detuning is
employed, whereas the ‘coarse/fine’ method works well for arbitrary $v$.

The calculation includes a new normalization factor, dividing by $N_{v}$,
which we can do with two different heuristic interpretations: First, we can
consider that each velocity class is only a portion ($1/N_{v}$) of the total
atom number. In this case, we sum over the velocities to get total
susceptibility from all of them. In a second interpretation, each $\rho_{12}$
represents the transmission spectrum if the entire atomic population had that
velocity class. In this case, we average over the velocities, to see how the
total population transmits light. Note in either case, the sum of $\rho_{12}$
or susceptibility over velocities is divided by $N_{v}$, and each are equally
weighed.

Having combined susceptibilities, we can compute total transmission using
Beer’s law, Eq. 8 in Appx. D. We note that adding susceptibilities (inside the
exponential) and taking one exponential is nominally fewer computations than
taking an exponential of each susceptibility and multiplying those
transmission values, as the two are mathematically equivalent. Additionally,
this method skips the Maxwell-Boltzmann weighting multiplication, which also
grows with $\mathcal{O}(N_{v})$, although these operations are insignificant
with respect to finding the steady-state atomic coherence matrix $\rho$, which
remains the largest computational load.

##  IV Results

Figure 3: Comparison of simulation results using population (left) vs.
velocity (right) sampling methods, using parameters listed in Sec. IV. Top
row: Transmission spectra (offset) over detuning for many constituent velocity
samples. Middle row: Transmission for a two-dimensional parameter scan over
detuning and velocity. Bottom row: Doppler-averaged total transmittance
spectra.  Figure 4:  Comparison of calculated transmission curves over
coupling detuning for the population sampling method (red) and the velocity
sampling method (blue), as $N_{v}$ is varied (plot titles).

We illustrate the new “population” sampling method with the “velocity”
sampling method in Fig. 3. We compare computed EIT transmission spectra of
Cesium’s 55$S_{1/2}$ state, holding the same number of sampling points
($N_{v}=81,N_{\Delta}=201$) and all other parameters held constant
($\lambda_{p}\approx 852.3$ nm, $\lambda_{c}\approx 509.4$ nm,
$\Omega_{p}/2\pi\approx 18$ MHz, $\Omega_{c}/2\pi\approx 3.4$ MHz,
$\Delta_{p}=\Delta_{c}=0$, temperature of 293 K, cell length 75 mm), using Eq.
7 to generate $v_{vec}$ for population sampling on the left, and velocity
sampling over 3$\sigma$ on the right side. This comparison shows that
population sampling converges to a smooth transmission curve with the same
parameters, whereas velocity sampling exhibits sampling artifacts.

The first row in Fig. 3 plots individual velocity sample transmission curves
(i.e., the EIT signal), where we have offset the transmission of each curve to
illustrate how they combine into the total transmission curve. The middle row
of Fig. 3 represents the same information in a surface plot over detuning
horizontally, and velocity sample point vertically, using the color to
represent transmission value. The bottom plot of Fig. 3 shows the resultant
curve when combining velocities for the two cases. The population sampling
method provides denser sampling in the resonant region for equal $N_{v}$,
while ‘bumpy’ artifacts remain in velocity sampling case, due to insufficient
sampling density.

The primary objective of this effort is to reduce the number times we have to
calculate the steady-state solution for $\rho_{12}$ before reaching a ‘smooth’
simulated spectroscopy absorption profile. We show a qualitative comparison of
transmission curves for a range of $N_{v}$ in Fig. 4. The low density velocity
sampling is apparent in both cases for low $N_{v}$, but the population method
is seen to converge to a smooth curve for moderate values of $N_{v}$.

In Fig. 5, we quantitatively compare results from both methods while varying
$N_{v}$ to see when each method converges to a smooth transmission curve. To
measure the computational convergence rates as $N_{v}$ is increased, we use
two metrics: the extrema transmission values of the spectrum in Fig. 5(a), and
the ‘error’ quantified as RMS of residuals from the population method at
$N_{v}=101$ in Fig. 5(b). The actual values of $N_{v}$ required for
convergence will depend heavily on other parameters in any particular
application. Plotted in Fig. 5(c), we show computation time required between
the two methods, which is nearly equal, dominated by the steady-state $\rho$
calculation, which occurs $N_{v}\times N_{\Delta}$ times in either method. The
highest $N_{v}$ were calculated first, and ran slightly faster than the linear
trend of $\mathcal{O}(N_{v})$. The advantage in computation speedup is
illustrated in Fig. 5(d), where we plot the RMS residual ‘error’ against
computation time. The population sampling method reaches 1% RMS residuals in
less than half the computation time as the velocity sampling method. By over-
sampling velocities near zero, far fewer samples are needed to obtain accurate
spectra.

Figure 5: Calculation time using odd $N_{v}$ for the population method in red
red$\times$’s and the velocity method in blue blue$\circ$’s. (a) Total
transmission maxima and minima. (b) RMS of differences between each spectrum
and the population curve at $N_{v}=101$. (c) Time per spectrum calculation.
(d) Residual RMS against calculation time.

##  V Conclusion

We have presented a method for calculating Doppler-broadened transmission
spectra for Rydberg EIT measurements with more efficient sampling over
population than velocity. Rather than using a Maxwell-Boltzmann weighting over
an equal-spaced velocity distribution, we use an equal-population partition to
‘naturally’ sample densely near zero velocity, and sparsely sample into the
thermal distribution’s wings. This significantly lowers the required number of
velocity-sampling curves to be calculated before ‘smooth’ convergence of the
transmission curves is obtained. In this work, we reach 1% error in less than
half the time of the typical method for the chosen parameters. Conversion of
Doppler-sampling scripts to the new method is as simple as adjusting the
velocity vector that is sampled over, and changing the Maxwell-Boltzmann
weighting operation of Eq. 3 into an average over velocities.

##  Appendix A Cheap (Un-physical) Calculations

Figure 6: Transmission spectra from partial velocity populations $0<\xi<1$
around $v=0$.

There is one further simplification that can be implemented for rough line-
shape calculations for the resonant case (including for our typical
application of Autler-townes splitting [1]), but does not give physical total
transmission values. Looking at the velocity dependence in the middle row of
Fig. 3, we notice that only the central portion near $v=0$ experiences
resonant EIT, and the rest is absorbed independent of $\Delta_{c}$. Instead of
calculating the entire velocity distribution, we can restrict the bounds of
$v_{vec}$ to just the lowest velocities, which contribute most to resonant
spectral features. This is synonymous with taking the above ‘fine’ partition
alone, and tossing the ‘coarse’ parts. Implementing this is as easy as using a
tighter bound, $0<\xi<1$ to generate $\eta_{vec}$ from $-\xi$ to $\xi$.
Maintaining sampling density now requires smaller $N_{v}$ by a factor of
$\xi$, reducing computation time. This sub-sampling reveals the approximate
line shape of EIT features, without giving an appropriate total transmission,
lacking the remaining probe absorption.

We illustrate the effect of varying $\xi$ in Fig. 6, where we adjust $N_{v}$
to maintain sampling density, such that the computation time lowers by a
nominal factor of $\xi$ while giving approximate line shapes. As a result of
this ‘velocity band’ narrowing, we use a scaling factor of
$N_{0}\rightarrow\xi{}N_{0}/N_{v}$ in Eq. 9, to account for the sub-sampling
of the population space. One can lower $\xi$ until it ‘cuts into the wings’ of
the EIT feature, which is more evident for smaller probe powers than the one
in our simulations, which broadens the line to order $\Omega_{p}\approx$18
MHz. Lowering $\xi$ crucially allows us to reduce $N_{v}$ at the same time to
achieve the same sampling density. This method is most useful for resonant,
poorly Doppler-matched beams, where only the near-zero velocities exhibit EIT.

##  Appendix B Higher Dimension Extension

This sampling method can be easily extended to apply to EIT configurations
with three or more lasers which are non counter-propagating, as in planar (2D)
three-photon configuration [18, 19, 11], or in general 3D. These schemes use
angle-tuning to better remove residual Doppler mis-match, depending on the
choice of wavelengths used [20].

The Gaussian velocity distribution applies independently over each direction.
Therefore in 2-D, we can loop over both velocity directions (i.e. for each
$v_{x}$ calculate all $v_{y}$ combinations), sampling the same $v_{vec}$
values, independently. For each combination of $v_{x}$ and $v_{y}$, we then
calculate the Doppler-adjusted detuning
$\Delta^{\prime}=\Delta-\textbf{v}\cdot{}\textbf{k}$ for each laser, where for
convenience, at least one laser can be along $\hat{x}$, whose effective
detuning is only changed in the outer loop. An off-axis beam’s Doppler
adjustment $\textbf{v}\cdot{}\textbf{k}$ can be inferred from the atom’s
trajectory $\theta=\tan^{-1}\left(\frac{v_{y}}{v_{x}}\right)$ and total speed
$|v|=\sqrt{v_{x}^{2}+v_{y}^{2}}$. An extra normalization factor of
$N_{0}/N_{v}^{2}$ must be included before the sum over both dimensions, as
each sample represents a further division of the population space from the 1-D
case. This approach is easily extended from two to three dimensions, nominally
increasing calculation time from $\mathcal{O}(N_{v}^{2})$ to
$\mathcal{O}(N_{v}^{3})$ from 2-D to 3-D respectively, making reductions in
$N_{v}$ all the more significant.

##  Appendix C Low-Probe Case

In addition to Doppler broadening of lines, velocity-dependent absorption
curves give other Doppler effects, such as enhanced _absorption_ on either
side of an EIT line [11]. To ensure our method captures these features as
well, we compare the population and velocity sampling methods in Fig. 7. This
calculation uses $\Omega_{p}/2\pi\approx 0.40$ MHz, and
$\Omega_{c}/2\pi\approx 3.1$ MHz. Since this feature is near the resonance, it
requires significantly higher sampling density than cases with significant
probe Rabi broadening, as in the rest of this paper. Again, the population
sampling method appears to converge closer to a smooth curve for the same
$N_{v}$ compared to the velocity sampling method.

Figure 7: Doppler-based absorption enhancement for a small probe power,
comparing results from either sampling method.

##  Appendix D Master-equation model

We find the time-independent solutions of a master-equation model for
calculating EIT signals. The model presented here is for 133Cs atoms, although
a similar model for 87Rb is presented in Ref. [1]. The power of the probe beam
measured on the detector (the EIT signal, i.e., the probe transmission through
the vapor cell) is given by [21]

$P_{out}=P_{in}\exp\left(-\frac{2\pi L\,\,{\rm
Im}\left[\chi\right]}{\lambda_{p}}\right)=P_{in}\exp\left(-\alpha
L\right)\,\,\,,$ (8)

where $P_{in}$ is the power of the probe beam at the input of the cell, $L$ is
the length of the cell, $\lambda_{p}$ is the wavelength of the probe laser,
$\chi$ is the susceptibility of the medium seen by the probe laser, and
$\alpha=2\pi{\rm Im}\left[\chi\right]/\lambda_{p}$ is Beer’s absorption
coefficient for the probe laser. The susceptibility for the probe laser is
related to the density matrix component ($\rho_{21}$) by the following [22]

$\chi=\frac{2\,{\cal{N}}_{0}\wp_{12}}{E_{p}\epsilon_{0}}\rho_{21_{D}}=\frac{2\,{\cal{N}}_{0}}{\epsilon_{0}\hbar}\frac{(d\,e\,a_{0})^{2}}{\Omega_{p}}\rho_{21_{D}}\,\,\,,$
(9)

where $d=2.02$ is the normalized transition-dipole moment [8] for the probe
laser and $\Omega_{p}$ is the Rabi frequency for the probe laser in units of
rad/s. The subscript $D$ on $\rho_{21}$ presents a Doppler averaged value, per
this paper’s theme. ${\cal{N}}_{0}$ is the total density of atoms in the cell
and is given by

${\cal{N}}_{0}=\frac{p}{k_{B}T}\,\,,$ (10)

where $k_{B}$ is the Boltzmann constant, $T$ is temperature in Kelvin, and the
pressure $p$ (in units of Pa) is given by [8]

$p=10^{-217.3571+\frac{1088.676}{T}-0.08336185\cdot{}T+94.88752\log_{10}(T)}$
(11)

In Eq. (9), $\wp_{12}$ is the transition-dipole moment for the
$\ket{1}$-$\ket{2}$ transition, $\epsilon_{0}$ is the vacuum permittivity, and
$E_{p}$ is the amplitude of the probe laser E-field.

The density matrix component ($\rho_{21}$) is obtained from the master
equation [22]

$\dot{\boldsymbol{\rho}}=\frac{\partial\boldsymbol{\rho}}{\partial
t}=-\frac{i}{\hbar}\left[\mathbf{H},\boldsymbol{\rho}\right]+\boldsymbol{\cal{L}}\,\,\,,$
(12)

where $\mathbf{H}$ is the Hamiltonian of the atomic system under consideration
and ${\boldsymbol{\cal{L}}}$ is the Lindblad operator that accounts for the
decay processes in the atom. The $\mathbf{H}$ and $\boldsymbol{\cal{L}}$
matrices for the three different tuning schemes are given below.

We numerically solve these equations to find the steady-state solution for
$\rho_{21}$ for various values of Rabi frequency ($\Omega_{i}$) and detunings
($\Delta_{i}$). This is done by forming a matrix with the system of equations
for $\dot{\rho}_{ij}=0$. The null-space of the resulting system matrix is the
steady-state solution. The steady-state solution for $\rho_{21}$ is then
Doppler averaged [22]

$\rho_{21_{D}}=\frac{1}{\sqrt{\pi}\,\,v_{\sigma}}\int_{-3v_{\sigma}}^{3v_{\sigma}}\rho_{21}\left(\Delta^{\prime}_{p},\Delta^{\prime}_{c}\right)\,\,e^{\frac{-v^{2}}{v_{\sigma}^{2}}}\,\,dv\,\,\,,$
(13)

where $v_{\sigma}=\sqrt{k_{B}T/m}$ and $m$ is the mass of the atom. We use the
case where the probe and coupling laser are counter-propagating. Thus, the
frequency seen by the atom moving toward the probe beam is upshifted by $2\pi
v/\lambda_{p}$ (where $v$ is the velocity of the atoms), while the frequency
of the coupling beam seen by the same atom is downshifted by $2\pi
v/\lambda_{c}$. The probe and coupling beam detuning is modified by the
following

$\Delta^{\prime}_{p}=\Delta_{p}-\frac{2\pi}{\lambda_{p}}v\,\,\,{\rm
and}\,\,\,\Delta^{\prime}_{c}=\Delta_{c}+\frac{2\pi}{\lambda_{c}}v\,\,\,.$
(14)

For the three level system, the Hamiltonian can be expressed as:

$\footnotesize H=\frac{\hbar}{2}\left[\begin{array}[]{cccc}0&\Omega_{p}&0\\\
\Omega_{p}&-2\Delta_{p}^{\prime}&\Omega_{c}\\\
0&\Omega_{c}&-2(\Delta_{p}^{\prime}+\Delta_{c}^{\prime})\\\
\end{array}\right]\,\,,$ (15)

where $\Delta_{p}^{\prime}$ and $\Delta_{c}^{\prime}$ are the Doppler-shifted
detunings, and $\Omega_{p}$ and$\Omega_{c}$ are the Rabi frequencies
associated with the probe laser and coupling laser respectively. The detuning
for each field is defined as

$\Delta_{p,c}=\omega_{p,c}-\omega_{o_{p,c}}\,\,\,,$ (16)

where $\omega_{p,c}$ are the on-resonance angular frequencies of transitions
$\ket{1}$-$\ket{2}$, $\ket{2}$-$\ket{3}$, and $\ket{3}$-$\ket{4}$,
respectively; and $\omega_{p,c}$ are the angular frequencies of the probe
laser and coupling laser, respectively. The Rabi frequencies are defined as
$\Omega_{p,c}=|E_{p,c}|\frac{\wp_{p,c}}{\hbar}$, where $|E_{p,c}|$ are the
magnitudes of the E-field of the probe laser, the coupling laser, and the RF
source, respectively. Finally, $\wp_{p}$ and $\wp_{c}$ are the atomic dipole
moments corresponding to the probe, coupling, and RF transitions.

For the three-level system, the ${\cal{L}}$ matrix is given by

$\footnotesize{\cal{L}}=\left[\begin{array}[]{cccc}\Gamma_{2}\rho_{22}&-\gamma_{12}\rho_{12}&-\gamma_{13}\rho_{13}\\\
-\gamma_{21}\rho_{21}&\Gamma_{3}\rho_{33}-\Gamma_{2}\rho_{22}&-\gamma_{23}\rho_{23}\\\
-\gamma_{31}\rho_{31}&-\gamma_{32}\rho_{32}&-\Gamma_{3}\rho_{33}\end{array}\right]\,\,\,,$
(17)

where $\Gamma_{i}$ are the state-to-state spontaneous decay rates, and
$\gamma_{ij}=(\Gamma_{i}+\Gamma_{j})/2$ use the total decoherence rates of
either state. We have ground state $\Gamma_{1}\leavevmode\nobreak\
=\leavevmode\nobreak\ 0$, D2 line [8]
$\Gamma_{2}=2\pi\cdot\leavevmode\nobreak\ 5.222$ MHz, and state-to-state decay
from $\ket{55S_{1/2}}\rightarrow\ket{6P_{3/2}}$ of
$\Gamma_{3}=2\pi\cdot\leavevmode\nobreak\ 308$ Hz. The decoherence terms
$\gamma_{3,i}$ use total Rydberg decay rate at $T=300$ K of
$2\pi\cdot\leavevmode\nobreak\ 83.867$ kHz [23].

## References

  * Holloway  _et al._ [2017] C. L. Holloway, M. T. Simons, J. A. Gordon, A. Dienstfrey, D. A. Anderson, and G. Raithel, Electric field metrology for si traceability: Systematic measurement uncertainties in electromagnetically induced transparency in atomic vapor, Journal of Applied Physics 121, 233106 (2017).
  * Artusio-Glimpse  _et al._ [2022] A. Artusio-Glimpse, M. T. Simons, N. Prajapati, and C. L. Holloway, Modern rf measurements with hot atoms: A technology review of rydberg atom-based radio frequency field sensors, IEEE Microwave Magazine 23, 44 (2022).
  * Sedlacek  _et al._ [2012] J. A. Sedlacek, A. Schwettmann, H. Kübler, R. Löw, T. Pfau, and J. P. Shaffer, Microwave electrometry with rydberg atoms in a vapour cell using bright atomic resonances, Nature Physics 8, 819 (2012).
  * Leahy  _et al._ [1997] C. Leahy, J. T. Hastings, and P. Wilt, Temperature dependence of doppler-broadening in rubidium: An undergraduate experiment, American Journal of Physics 65, 367 (1997).
  * Peters  _et al._ [2012] T. Peters, B. Wittrock, F. Blatt, T. Halfmann, and L. P. Yatsenko, Thermometry of ultracold atoms by electromagnetically induced transparency, Physical Review A 85, 063416 (2012).
  * Javan  _et al._ [2002] A. Javan, O. Kocharovskaya, H. Lee, and M. O. Scully, Narrowing of electromagnetically induced transparency resonance in a doppler-broadened medium, Physical Review A 66, 013805 (2002).
  * Finkelstein  _et al._ [2019] R. Finkelstein, O. Lahad, O. Michel, O. Davidson, E. Poem, and O. Firstenberg, Power narrowing: counteracting doppler broadening in two-color transitions, New journal of physics 21, 103024 (2019).
  * Steck [1998] D. A. Steck, Cesium D line data (1998), revision 2.2.1, 21 November 2019.
  * Maxwell [1860] J. C. Maxwell, V. illustrations of the dynamical theory of gases.—part i. on the motions and collisions of perfectly elastic spheres, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 19, 19 (1860).
  * Doppler [1903] C. Doppler,  _Ueber das farbige Licht der Doppelsterne und einiger anderer Gestirne des Himmels: Versuch einer das Bradley’sche Aberrations-Theorem als integrirenden Theil in sich schliessenden allgemeineren Theorie_ (K. Böhm Gesellschaft der Wissenschaften, 1903).
  * Carr  _et al._ [2012] C. Carr, M. Tanasittikosol, A. Sargsyan, D. Sarkisyan, C. S. Adams, and K. J. Weatherill, Three-photon electromagnetically induced transparency using rydberg states, Optics letters 37, 3858 (2012).
  * Finkelstein  _et al._ [2023] R. Finkelstein, S. Bali, O. Firstenberg, and I. Novikova, A practical guide to electromagnetically induced transparency in atomic vapor, New Journal of Physics 25, 035001 (2023).
  * Tanasittikosol  _et al._ [2011] M. Tanasittikosol, J. Pritchard, D. Maxwell, A. Gauguet, K. Weatherill, R. Potvliege, and C. Adams, Microwave dressing of rydberg dark states, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 184020 (2011).
  * Olver and Townsend [2013] S. Olver and A. Townsend, Fast inverse transform sampling in one and two dimensions, arXiv preprint arXiv:1307.1223 10.48550/arXiv.1307.1223 (2013).
  * Finney [1952] D. J. Finney,  _Probit analysis: a statistical treatment of the sigmoid response curve_ (Cambridge university press, Cambridge, 1952).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1964).
  * Wichura [1988] M. J. Wichura, Algorithm as 241: The percentage points of the normal distribution, Journal of the Royal Statistical Society. Series C (Applied Statistics) 37, 477 (1988).
  * Grynberg and Cagnac [1977] G. Grynberg and B. Cagnac, Doppler-free multiphotonic spectroscopy, Reports on Progress in Physics 40, 791 (1977).
  * Thaicharoen  _et al._ [2019] N. Thaicharoen, K. Moore, D. Anderson, R. Powel, E. Peterson, and G. Raithel, Electromagnetically induced transparency, absorption, and microwave-field sensing in a rb vapor cell with a three-color all-infrared laser system, Physical Review A 100, 063427 (2019).
  * Shaffer and Kübler [2018] J. Shaffer and H. Kübler, A read-out enhancement for microwave electric field sensing with rydberg atoms, in  _Quantum Technologies 2018_, Vol. 10674 (SPIE, 2018) pp. 39–49.
  * Yarin [1995] A. Yarin,  _Optical Electronics, 4th ed._ (CBS College Publishing, New York, 1995).
  * Berman and Malinovsky [2011] P. R. Berman and V. S. Malinovsky,  _Principles of Laser Spectroscopy and Quantum Optics_ (Princeton University Press, Princeton, 2011).
  * Šibalić  _et al._ [2017] N. Šibalić, J. D. Pritchard, C. S. Adams, and K. J. Weatherill, Arc: An open-source library for calculating properties of alkali rydberg atoms, Comput. Phys. Commun. 220, 319 (2017).