Near source fluorescence spectroscopy for miniaturized thermal atomic beams

Miniature atomic beams can provide new functionalities for atom based sensing instruments such as atomic clocks and interferometers. We recently demonstrated a planar silicon device for generating well-collimated thermal atomic beams [Nat Commun 10, 1831 (2019)]. Here, we present a near-source fluorescence spectroscopy (NSFS) technique that can fully characterize such miniature beams even when measured only a few millimeters from the nozzle exit. We also present a recipe for predicting the fluorescence spectrum, and therefore, the source angular distribution, even under conditions of strong laser saturation of the probing transition. Monte Carlo simulations together with multi-level master equation calculations fully account for the influence of optical pumping and spatial extension of the Gaussian laser beam. A notable consequence of this work is the agreement between theory and experimental data that has allowed fine details of the angular distribution of the collimator to be resolved over 3 decades of dynamic range of atomic beam output flux.


I. Introduction
Microfabrication techniques have a remarkable opportunity to transform atomic sensors, normally laboratory scale devices, into portable instruments. Such instruments are urgently needed for precise navigation and timing, electromagnetic field sensing and gravimetry, all applications where atoms provide a basic reference standard * lichao@gatech.edu that is traceable to fundamental constants [2][3][4][5]. Prototypical examples of atomic platforms currently targeted for these applications are ultracold atoms on a chip and MEMS-based microfabricated alkali vapor cells [6,7].
In previous work, we demonstrated a continuous, miniature thermal atomic beam source that propagated along a silicon surface in microchannel arrays defined by photolithography [1]. Such well-collimated sources could find application in miniature atom interferometry [8] and enable precise studies of atom-surface interactions [9]. Moreover, unlike vapor cells, they do not suffer from Doppler broadening or wall collisions [10]. Thus, they combine the simplicity of thermal vapors with the spectroscopic purity of cold atom systems. Since these sources are necessarily miniature, however, one needs to characterize the atomic beam close to the source itself, where the atoms' spatial and velocity distributions are mixed together. This is in contrast to typical laboratory-scale atomic beam experiments operating in the far-field, where the two distributions can be separated [8,[11][12][13][14]. In this work, we use fluorescence spectroscopy to perform a four-dimensional tomography, reconstructing the spatial and speed distribution for the chip scale atomic beam. We term this near-source fluorescence spectroscopy (NSFS).

II. Experimentation
The fabrication procedure for the chip scale atomic collimator as well as the fluorescence measurement protocols have been described in detail in the previous work [1]. We briefly review these here, showing the major components of the NSFS measurement in Fig. 1a-1c. Thermal atomic beams generated by 20 silicon microcapillaries attached to an effusive oven were probed a few millimeters A microfabricated silicon device with 20 collimating channels sits inside a slit at the center of the copper head. The red arrow along +ŷ represents the laser beam and the gray arrow along +ẑ represents the atomic beams.
The shaded region in blue shows the field of view of our fluorescence collecting system. (b) shows the cartoon top-view of such a collimator (termed an ordinary collimator, in contrast to the cascaded collimator of (c)) before bonding the sealing wafer on top. The chip is 3 mm × 5 mm in size with 0.5 mm thickness. For an individual channel, l/d = 30 (l = 3 mm and d = 100 µm). (c) shows the top-view of a 3-stage cascaded collimator (∼ 660 µm/stage and ∼ 500 µm/gap) with the same overall length and diameter for an individual channel before wafer bonding.
(d) Estimated spatial distribution of the normalized resonant photon scattering rate, s/(1 + s) in terms of probe saturation parameter s = I( r)/I sat , on planes defined by y = 0 (side view) and x = 0 (top view), for three different laser intensities corresponding to s max (or s 0 ) = 1, 10, 100. Each panel is 4 × 12 mm. The laser beam has Gaussian radius (for 1/e 2 intensity) w x = 0.5 mm and w z = 1.4 mm. Bars at the top represent the array of collimator outputs. Two lines forming an angle β as a guide depicts atoms emitted within the collimator full-width at half-maximum (FWHM) angle β = 2θ 1/2 = 3.6 • [15].
after the nozzle exit. 20 atomic beams spaced by 150 µm center-to-center distance, co-propagate along the +ẑ direction on the y − O − z plane as defined in Fig. 1. An external cavity diode laser at 780 nm used as a fluorescence probe is scanned over a 1 GHz range across the 87 Rb D 2 F = 2 to F ′ = 3 transition at a rate of 5 Hz. The laser beam is propagating perpendicular to the traveling direction of the atom beam along the +ŷ direction centered at z 0 ≈ 6 mm and it is linearly polarized along the z−axis to maximize the collected fluorescence. Fluorescence emitted from the volume shown in Fig. 1a at the intersection of the laser and atomic beams is collected through two 2-inch lenses (not shown) located ≃3 inches above this volume. Light is collected onto a photodiode and the photocurrent is amplified by a low noise current amplifier. Because of the transverse Doppler effect, the fluorescence collected at different laser detunings is sensitive to the transverse speed distribution of the atomic beam, which contains information about the collimating performance of our microfabricated silicon microcapillary array-better collimation means narrower transverse speed distribution. Part of the laser output is injected into a rubidium vapor cell at room temperature for saturated absorption spectroscopy, calibrating its operating frequency and assisting the scan control.
We performed two sets of experiments for the two different types of collimators shown in Fig. 1b and 1c. The cascaded collimator (1c) generates 40 times purer atomic beams compared to the ordinary type single-stage collimator (1b). This is because the two gaps efficiently release atoms whose trajectory deviates from the central axis, thus, behaving as a transverse velocity filter, as discussed in Ref. [1]. For each collimator, we recorded fluorescence spectra over 10 different laser intensities adjusted by varying the probe laser power. All other experimental parameters such as laser beam alignment, beam width and propagation direction were kept identical from one collimator to another, so that the spectra could be directly compared with one another.

A. Experimental spectra
In general, laser spectroscopy is most precise at low illumination (low saturation parameter s = I/I sat ), in order for the measurement to be minimally perturbative of the system under study. Here I sat is the saturation intensity for the particular optical transition being probed. However, fluorescence detection benefits from a higher probe intensity in order to overcome background noise caused by environmental light and detector/electronic noise sources. Therefore, in practice one needs to work at an appreciable value of s. In our case we would like to resolve the fluorescence spectrum at large detunings to probe the atomic emission at large angles to the collimator axis. Therefore, we have developed a method to extract this velocity distribution even at finite values of s where saturation cannot be ignored. Fig. 2a shows our main experimental result. We show experimentally measured spectra over a 200 MHz range of detunings that agree very precisely with a full threedimensional numerical calculation according to Eq. (9). By contrast, the simple one-dimensional theory used in our earlier work [1] shows clear deviations at the level of 30% that are particularly pronounced for the ordinary collimator of Fig. 1b. While the best agreement occurs, as expected, for low saturation parameters, we have also systematically investigated in Fig. 2b and 2c the influence of laser intensity on the fluorescence spectra up to saturation parameter s 0 = I 0 /I sat = 250. Here I 0 is the laser peak intensity and I sat = 3.05 mW/cm 2 is for linearly polarized light coupling the 87 Rb D 2 F = 2 to F ′ = 3 transition. Good agreement is found throughout the range of parameters explored.
In Fig. 2a we have plotted data on the blue side of the atomic resonance, as it avoids contamination from the F = 2 to F ′ = 2, 1 transitions occurring at -267 and -424 MHz, respectively, as well as contributions from the 85 Rb isotope's hyperfine transitions [16]. Therefore, we can regard our atoms as a pure two-level system for the F = 2 to F ′ = 3 transition, with the spectral wings truly representative of the atoms' transverse Doppler velocity distribution. The transition natural linewidth is Γ = 2π × 6.1 MHz, but other mechanisms contribute to broadening the spectra. These include power broadening, transit-time broadening, Doppler broadening, Zeeman broadening, collision broadening, as well as the finite laser linewidth [17]. The power broadening dominates over these other mechanisms for the parameters of our experiment, and scales as Γ √ 1 + s with respect to the saturation parameter. For example, at 100 C, the transit time for atoms going though the laser beam is estimated to be 2w z /v ≈ 9 µs contributing a transit-time broadening about 0.1 MHz. The nozzle heater coil wrapped onto the copper head (shown in Fig. 1a) has a DC current running through it that generates a magnetic field with a strength less than 1 Gauss, corresponding to less than ∼ 0.23 MHz line broadening [16]. No significant Zeeman broadening is also verified by momentarily turning off the nozzle heater and monitoring if the spectra get narrower. Since the mean free path is larger than the probe beam diameter, the collision broadening can be neglected because of the low atom number density of 10 8 atoms/cm 3 6 mm away from the nozzle exit at T = 100 C. The laser linewidth is < 1 MHz. Hereafter, we discuss mainly power broadening and the Doppler effect, since all other mechanisms together contribute a total broadening at the ∼ 1 MHz level that can be neglected.
It is interesting to compare the different theoretical curves shown in Fig. 2a. Our 1D theory assumes that the atom number density only varies with respect to the transverse coordinate y, and not the transverse coordinate x. It is a useful approach if the laser interrogation occurs far from the source, i.e. z 0 ≫ 2w, where w is the Gaussian beam waist of the probe. In this limit, one can neglect off-axis atomic trajectories along the xdirections, as well as laser intensity variations. In reality, however, our probe is very close to the source, as seen in Fig. 1a and 1d, and one must account for the full three-dimensional nature of both the laser beam profile as well as atomic trajectories through the beam. For our data, we see that the 1D theory accurately reproduces the spectral full-width at half-maximum (FWHM) but not the wings of the data. Atomic trajectories at large angles to the central axis contribute most to these wings.
To gain a further qualitative understanding, we can divide the population into two components: atoms traveling inside θ 1/2 (beam component) and outside θ 1/2 (diffuse component). The half width of the transverse speed distribution for the beam component can be reasonably estimated by 2vθ 1/2 , while for the diffuse component it is much larger since atoms emitted into the diffuse component originated from atom-tube reflections within the source. We can think of the laser beam as a bundle of light rays traveling along y with different x and z coordinates. All rays that lie in the x = 0 plane intersect both the beam component and the diffuse component. However, rays in the z = z 0 plane with |x| > z 0 θ 1/2 will interact only with the diffuse component. Thus, the latter will measure a broader spectra. We quantify this result further in Fig. 4c. The final spectrum will be a sum over rays with different coordinates, weighted by the contributions of different atomic densities, trajectories and laser intensities. For our simulations, we assumed that the atomic speed distribution emitted into any direction obeys the Maxwell-Boltzmann distribution (in traditional atomic physics texts this is known as setting the deformation function to 1 [18,19]).
Armed with this theoretical approach, we apply it to a variety of different saturation parameters used in the  experiment. In Fig. 2b and 2c, we show the peak height and half width at half maximum (HWHM) vs. saturation parameter s 0 of values swept from 1 to 250 for two sets of measured fluorescence spectra for the ordinary collimator and the cascaded collimator, respectively. Our theory can successfully reproduce the measured peak fluorescence intensity as well as the HWHM of the measured spectra. This agreement was achieved without any adjustable parameters, assuming the known saturation parameter of I sat = 3.05 mW/cm 2 for linearly polarized light coupling the 87 Rb D 2 F = 2 to F ′ = 3 transition. What appears to be peculiar about this result is that there is no apparent saturation effect at all. The measured curves in Fig. 2b continue to increase with laser intensity, although the rate of growth at high values of s is slower for the cascaded collimator than for the ordinary collimator. This apparently counterintuitive behavior is a consequence of the laser interrogation region being close to the source, a condition not satisfied in most experiments, where a clear saturation can be observed [11]. As the laser intensity increases, the scattering rate in the center of the beam saturates, while in the wings, it can continue to grow. The ordinary colli-mator contains a diffuse component that encounters the laser beam wings. By contrast, the cascaded collimator does not contain this diffuse component and therefore shows a less pronounced rate of growth. This effect is well-captured by our three-dimensional theory, although the measured spectral HWHM was slightly lower than predicted, as will be discussed later.

B. Optical pumping
In addition to the saturation behavior explored earlier, we also deduced and quantified the influence of high probe laser intensity optical pumping on the hyperfine transitions. In Fig. 2d, we show broader spectra inclusive of all relevant hyperfine transitions as we scanned the laser frequency. The laser scan rate was approximately 12 GHz/second, and hence the laser frequency stayed within the natural linewidth Γ = 6.1 MHz of any resonance for approximately 0.5 ms, longer than both the 9 µs transit time as well as the optical pumping time. Optical pumping to the dark F = 1 ground state can occur, whence less fluorescence is collected on the F = 2 → F ′ = 1, 2 tran-sitions. Once atoms are driven into the excited states labeled by F ′ = 1, 2, they can spontaneously decay back to ground states label by F = 1, 2 since they are allowed by the selection rule (indicted by the dashed lines). Atoms in F = 2 state can by excited by the light field into the excited state again, but those pumped into F = 1 are 6.8 GHz detuned away, hence they will not fluoresce anymore. The optical pumping rate to the dark state F = 1 increases with laser intensity and thus relatively less fluorescence can be collected in comparison to the cycling F = 2 to F ′ = 3 transition. This effect is more pronounced at high saturation parameter, as shown in the spectra plotted in Fig. 2e. This effect was also confirmed by multi-level master equation simulations that computed the excited state populations (and therefore the emission line strengths) for each hyperfine transition. We discuss these simulations in detail in section IV C.

C. Long term collimator output
While the collimating device is relatively simple, for practical applications such as clocks and other precision instruments, it needs to be run for years without servicing. Traditional beam clocks have had to deal with the issue of nozzle clogging [19].  limited lifetime tests under continuous operation for a period of 6 months at various fixed temperatures. Our tests could both i) probe the overall flux emitted from the source using a photodetector to capture the fluorescence as well as ii) measure the output of every single individual microchannel using a CCD camera focused directly on the output of the chip. Such channel-by-channel testing is not easily accomplished using 3-dimensional capillary array sources [19,20], and our results are therefore a novel test of atomic beam behavior. For our tests, the Rb oven body was kept at a constant temperature, while the chip temperature was between 10 and 30 degrees higher. For the total flux measurement, we probed the atomic beam 60 mm downstream from the chip after the beam had passed through a stainless steel plate with a 9 mm hole that was 36 mm away from the chip. This aperture prevented any Rb vapor from accumulating in the probe region and contaminating the measurement. As before, the laser probe was perpendicular to the atomic beam with a similar fluorescence collection, detection and amplification setup. The optical collection efficiency was 1.15% and the power of the probe laser was 247 µw, with a saturation parameter I/I sat = 3.81. The spectrum of the fluorescence was deconvolved to obtain the transverse velocity distribution, which was used to calculate an average scattering rate R sc . We used R sc to calibrate the total throughput of the collimator [13]. The result is shown in Fig. 3. We set the oven temperature to 100 degrees at the beginning, 125 degrees at 1148 hours, 150 degrees at 1319 hours and 125 degrees again at 2387 hours. At the time of writing, this oven had run for over 4200 hours and none of the microchannels showed any sign of clogging. Pictures of the atomic beams near the nozzle were taken several times during the test, and showed that all 29 channels remain unclogged. It can be noticed in Fig. 3 that the flux is decaying slowly, which might result from the migrating of Rb inside the oven from hotter spots to colder spots. Coating of Rb on the vacuum windows over time and the drift of laser alignment may also contribute to this decaying effect. Since the flux is around 11 times higher at 150 degrees and 3 times higher at 125 degrees, our result implies a continuous operation time at 100 degrees of over 19000 hours without failure. This test proved that our microfabricated atomic beams are reliable and robust at different temperatures and can have a very long lifespan.

IV. Theory
Our theoretical calculations combined i) Monte-Carlo simulations of the atomic flux in the molecular flow regime based upon the actual experimental geometry of the collimator imported from a CAD model and ii) an atom-by-atom computation of the fluorescence spectrum based upon each atom's interaction with the laser beam along its particular trajectory using the geometry defined   in Fig. 1a. For the latter, we used a master equation simulation to compute the population of each excited hyperfine level to deduce the fluorescence rate. Details of these calculations are provided in the subsequent sections.

A. Monte Carlo simulations
We used Molflow+, a test-particle Monte Carlo simulator dealing with molecular flow, to simulate atomic trajectories that pass through the laser beam [22]. We assume that the collimators operate largely in the transparent regime [15], where the mean free path λ for atomic collisions is much larger than the length of an individual collimating channel l c . For our parameters, l c = 3mm, while λ can be estimated using the following equation where P is the pressure and d the atomic diameter. This condition begins to break down at temperatures of 120 C (150 C) for the single-stage (cascaded) collimators, as shown in Table I. However, uncertainties in the exact temperature and atomic cross-section make this an approximate assessment. In practice, we have not seen any discernable deviation from the transparent regime for oven temperatures up to 150 C for both collimators. The output beam can be fully characterized by an an-gular distribution function κf (θ): where we use the notation of [18]. For a single orifice, f (θ) = cos θ and therefore κ = 1. More generally for a tube, κ ∼ 1/W is the "peaking factor" of the collimator, depending inversely on its transmission probability W = (4d/3l c )(1 + 4d/3l c ) −1 , where l c , d are the tube length and diameter, respectively [18]. One can import CAD geometries for their collimators into Molflow+ and specify the inlet surface as an cosine emitter with sticking factor=1 mimicking the physics that atoms can return to the source. Sticking factors for all internal surfaces of the collimator are set to be zero since no Rb chemical absorption/physical condensation is assumed for crystalline silicon maintained at a temperature that is higher than the rubidium vapor temperature. At the channel exit, another facet with sticking factor=1 captures all atoms emitted and records their angles while hitting the surface. The atomic gas is considered to be isothermal neglecting the thermalization over the whole simulation when atoms enter the microchannels kept at higher temperature than the oven. After releasing enough test particles since the statistical error varies with 1/ √ N [25], the angular distribution function can be recovered using the following relation without distinguishing their speed distribution where with C(θ n ) represents the total number of particle falling into the n th bin spanned by the differential solid angle integrated over the specified range and ∆θ represents the θ angle sampling step. We did not sample φ angle assuming the angular distribution function has a built-in rotational symmetry with respect to its center line (z−axis), given that the cross-sectional shape of collimating tubes has negligible influences to its angular distribution function based on Ref. [21]. κf (θ) can then get appropriately normalized through Using Molflow+'s advanced facet parameters, we have finer samplings imposed for small angles θ to achieve better precision depicting the peaking behavior of the angular distribution function: ∆θ = 3.14E-4 for θ ∈ [0, 3.14E − 2] for the single-stage collimator while ∆θ = 1.57E-4 for θ ∈ [0, 1.57E − 2] for the cascaded collimator.
Another 500 values are equally sampled for θ angle outside these ranges up to π/2, and only 1 value for φ angle from 0 to 2π. Results for the single-stage collimator and the cascaded collimator are plotted in Fig. 4a, where the angular distribution functions have been multiplied by its corresponding transmission probability for comparison with the same peak value 1 noticing the fact that W · κf (θ = 0) = 1. For simple geometries as the ordinary type of collimators, an analytical expression for the angular distribution function κf (θ) exists, hence we plot the results given by Molflow+ together with the analytical expression. The good agreement between the two confirms the validity of Molflow+, which can be applied to any nontrivial structures, e.g. the cascaded collimators.
At the center axis of the Gaussian laser beam as shown in Fig. 1a, we then set up two square targets facing towards the collimator output. The two targets are sitting 6 mm away from the nozzle exit and detecting the number of hits per unit area per unit time. The impingement rate distributions for two different types of collimators are shown in Fig. 4b. More than 99% of atoms are actually emitted into the halos for the traditional type single-stage collimator. This can smear out the fine spectra features, contaminate nano structures or micro cavities [12,26,27], and reduce transparency of optical accesses. In contrast, the cascaded collimator with the identical physical size can produce purer atomic beams with 40 times fewer atoms emitted into large angles while maintaining the same axial beam intensity of nAv/(4π) [19].

B. Fluorescence spectral computation
Under the transparent regime, the on-axis atom number density at a certain location z with respect to the nozzle exit is given by where n 0 is the atom number density inside the oven, and A is the cross-section area of a single tube. For N t =20 tubes with diameter d, the optical depth at T =100 • C can be estimated using where σ = 3λ 2 /2π is the photon scattering cross section [16]. OD is about 6% at maximum for the 87 Rb D 2 F = 2 to F ′ = 3 transition. Neglecting the absorption, we can assume the laser beam intensity only varies with respect to x and z while it does not change much along its propagation axis y. The Gaussian laser beam intensity distribution can be written as For multiple atomic beams, the measured output voltage is where δ is the laser detuning, ω 0 gives the energy per photon, R resp is the responsivity of the photodiode, G is the current amplifier gain, η is the overall photon collecting efficiency, A b is the natural abundance of 87 Rb, f p is the population fraction for a certain hyperfine level, N t is the total number of microcapillaries, f c is the correction factor [1] for the theoretical total throughputṄ and B is determined by the fluorescence collecting volume. Numerical values of these scaling factors are given in Table II. While seemingly complicated, Eq. (9) can be formulated into a simple linear algebra problem after discretizing the integral (see the appendix) where A(s 0 ) is a matrix mapping the angular distribution function f into the predicted fluorescence spectra V f for a certain saturation parameter s 0 . Theoretical results for Fig. 2(a)-(c) are readily obtained once A(s 0 ) is computed for the ten different experimental saturation parameters. One can either predict the angular distribution function using Monte Carlo simulations and then predict the fluorescence spectra treating Eq. (10) as a simple forward problem, or measure the fluorescence spectrum and inversely find the angular distribution function [28,29].
In the low saturation parameter regime, a broader spectrum can be measured for the single-stage collimator while probing the fluorescence spectrum off the y − O − z plane. We quantitatively verify this by computing the fluorescence spectra using Eq. (9) along the dashed lines and the center lines, and normalize their peak value to be 1 for convenience (see Fig. 4c). That's why the 1-D approximation predicts the narrower spectrum compared to the measured data as reported in our previous work [1]. For the cascaded collimator, spectra measured at different x−locations are similar, since the number of atoms emitted from tube walls has already been suppressed by 40 times. Atoms passing through the region near the dotted lines for the cascaded collimator mainly come from the beam component rather than the diffuse component.
In the high saturation parameter regime, a narrow spectrum can be measured in the Gaussian wings of the laser beam compared to the center when probing the fluorescence spectrum on the y − O − z plane. Following the same procedure above, the spectra are calculated along lines at different z−locations (see Fig. 4d). The spectra computed at z = z 0 is broader than the spectra computed at z = z 0 − w z where s = s 0 /e 2 solely because the power broadening is larger. These spectra share the same non-trivial transverse speed distribution, and thus we can compare the power broadening directly. At the Gaussian wings of the laser beam, the spectral half width converges towards their intrinsic Doppler broadening as shown in Fig. 4c. Analytical estimates are not available for the half width of the measured fluorescence spectra when the Doppler and power broadening are comparable to each other, considering the atomic/laser beam inhomogeneous effects if probing close to the source.

C. Master equation simulations
Fig. 2e reveals that the red side of the fluorescence spectra manifests different line strength for the F = 2 to F ′ = 2, 1 transitions under different laser intensities due to optical pumping effects. A closed-form analytical expression for R sc , as required in Eq. (9), isn't available for these two transitions in contrast to the F = 2 to F ′ = 3 cycling transition. Therefore, to quantify these optical pumping effects, we need to run master equation simulations to track the dynamical evolution of the excited state populations contributing to the fluorescence.
The master equation reads as the following [1,30] where the first part is for the atom-light interaction and the second part is for the spontaneous decay. At low light level (peak intensity I 0 =3.8 mW/cm 2 ) and oven temperature T=100 C, an ensemble of atoms are assumed to be traveling with the most probable velocity v pB = 327 m/s for an atomic beam [31] going through the light field corresponding to Eq. (8), from (z 0 -1.5w z ) to (z 0 +1.5w z ) where the light intensity has dropped to a significantly low level . We initialize the density matrix according to the Boltzmann distribution. Fig. 5 shows the population dynamics for both F ′ = 3 and F ′ = 2 excitations. Our principal observation is that the F ′ = 3 excitation can be effectively treated as twolevel. Fig. 5a shows that the excited state population . The laser frequency is on-resonance with the 87 Rb D 2 F = 2 to F ′ = 3 cycling transition for (a), and is on-resonance with the F = 2 to F ′ = 2 non-cycling transition for (b), where we have a separate axis on the right for the population in excited states.
tracks its steady state value P e (F ′ = 3, r(t)) = 1 2 I( r(t))/I sat 1 + I( r(t))/I sat where I( r(t)) is the intensity at the instantaneous location r(t) of the atoms. The difference between Eq. (12) and the master equation simulation was negligible. From Eq. (12) we can get the on-resonance photon scattering rate R sc ( r(t)) = ΓP e (F ′ = 3, r(t)) for atoms at any specific locations. Fig. 5b shows the populations in |F = 2 , which can contribute to the collected fluorescence through temporarily occupying |F ′ = 2 and are eventually transferred into the dark state |F = 1 that has a 6.8 GHz detuning from the excitation field. An additional dark state is |F = 2, m F = 0 for which there is a zero coupling matrix element to |F ′ = 2, m F = 0 [17].
These states become substantially populated within ∼ 5 µs. The presence of a kink feature in the solid red curve at the beginning stage of the simulation is because the Rabi flopping takes a time about 1/Γ to reach steady state due to spontaneous emission, and during that time the system evolves coherently. The integrated number of photons scattered is computed from the excited state probabilities and yields good quantitative agreement with the measured spectra in Fig. 2e. These procedures were previously documented (see Ref. [1]) but not the master equation simulation data of Figure 6.

V. Conclusion
In summary, we have developed NSFS for scenarios where laser intensities, atoms' spatial distributions and speed distributions are all wound together. The cloggingfree continuous operation and the individual-channel addressability identify these on-chip elements as promising candidates, generating well-collimated atomic beams for clocks, interferometery or the atomic-optical device hybridization. [32][33][34][35][36]. Looking forward, combining Monte Carlo and master equation simulations allows the tracking of atom trajectories and their internal states simultaneously, which is indispensable for customizing the laser manipulation of atoms close to or within a chip-scale source itself [37]. The simple recipe phrased in linear algebra for reconstructing the speed/angular distribution of the atomic beams is important for precision spectroscopy [38] or atom scattering experiments [39], and useful for guiding collimator and oven designs [40][41][42]. where we've dropped all constants as coefficients and A(v, r, θ, φ, δ) represents A(v, r, θ, φ, δ, s 0 ) = sin θ · R sc (r, θ, φ, s 0 , δ, v) · v 2 e −(v/α) 2 .
The discrete version of Equation 18 will be Now, the theoretical prediction for the fluorescence spectrum has been completely formulated into a problem of the following type where f represents the angular distribution function f (θ), V f represents the fluorescence spectrum V f (δ) and matrix A(s 0 ) represents the mapping between the two under different saturation parameters. Interestingly, matrix A(s 0 ) does not depend on what type of collimator/angular distribution function we have at all. It intrinsically characterizes the mapping once the laser probe geometry and power are well defined. For absolute fluorescence signal calibration, we absorb previously dropped scaling coefficients back to Eq. (20) by the end of the computation. Scaling factors defined in Eq. (9) are summarized in Table II.