Large molasses-like cooling forces for molecules using polychromatic optical fields: A theoretical description

Recent theoretical investigations have indicated that rapid optical cycling should be feasible in complex polyatomic molecules with diverse constituents, geometries and symmetries. However, as a composite molecular mass grows, so does the required number of photon scattering events necessary to decelerate and confine molecular beams using laser light. Utilizing coherent momentum exchange between light fields and molecules can suppress spontaneous emission and significantly reduce experimental complexity for slowing and trapping. Working with BaH as a test species, we have identified a robust, experimentally viable configuration to achieve large molasses-like cooling forces for molecules using polychromatic optical fields addressing both $X-A$ and $X-B$ electronic transitions, simultaneously. Using numerical solutions of the time-dependent density matrix as well as Monte Carlo simulations, we demonstrate that creation of Suppressed Emission Rate (SupER) molasses with large capture velocities ($\sim 40$ m/s) is generically feasible for polyatomic molecules of increasing complexity that have an optical cycling center. Proposed SupER molasses are anticipated to not only extend quantum control to novel molecular species with abundant vibrational decay channels, but also significantly increase trapped densities for previously laser-cooled diatomic and triatomic species.


A. Direct Molecular Laser Cooling
Optical control over atomic spatial degrees of freedom is one of the cornerstones of modern atomic physics [1,2] and quantum technologies [3,4]. In recent years, laser cooling and trapping methods have been successfully extended to a handful of molecular species [5,6].
Yet despite more than a decade of active research efforts, only three diatomic species (SrF [7], CaF [8,9] and YO [10]) have been trapped in three dimensions (3D) at microkelvin temperatures. While one-dimensional (1D) laser cooling of cryogenic molecular beams of diatomic BaH [11] and YbF [12], triatomic SrOH [13], YbOH [14] and CaOH [15], and even hexatomic CaOCH 3 [16] molecules has been achieved, the number of scattered photons demonstrated (∼ 100 − 1, 000) is still at least an order of magnitude below what is needed to achieve radiative slowing and 3D trapping. Therefore, the question of general prospects of utilizing laser slowing, 3D cooling and trapping for molecular species with new internal structures (e.g. BaH) or increased vibrational complexity (e.g. triatomics) remains largely unanswered.
Traditional Doppler slowing and cooling relies on a repeated process of directional photon absorption (resulting in k momentum transfer) followed by spontaneous emission to the initial set of states for the cycle to repeat [17]. While in many atoms, the use of specific angular momentum configurations for the ground and excited states together with the appropriate laser polarization can lead to an effective "two-level" system, the absence of strict vibrational selection rules for molecular electronic decays necessitates novel approaches to molecular laser cooling [5]. The probability of decay into a given vibrational level is described by the square of the overlap integral between the excited and ground vibrational wavefunctions 1 F v v , also known as a Frank-Condon Factor (FCF) for that transition [6].
For certain diatomic species with small off-diagonal FCFs (i.e. v = v ), one or two additional lasers can be used to repump molecules from excited vibrational levels v > 0 back to the ground vibrational state v = 0, enabling scattering of 10 4 photons [18][19][20] needed to slow molecular beams to below the capture velocity of a 3D molecular magnetooptical trap with v cap ≈ 5 − 10 m/s [21,22]. However, even for light triatomic species with 1 Following the established convention in the field of molecular spectroscopy, we mark quantum numbers with double (single) primes to refer to the electronic ground (excited) state.
relatively diagonal FCFs like CaOH (F 00 = 0.954), eight additional repumping lasers are needed to scatter ∼ 10 4 photons [15,23,24], thus, presenting a significant technical challenge for extending Doppler slowing and trapping methods to heavier (e.g. YbOH) or more complex (e.g. CaOCH 3 ) molecules. Towards this end, various alternative techniques have been developed for efficient momentum transfer from the laser light to atoms or molecules, while minimizing spontaneous emissions [25]. To date, the emphasis has been on developing novel experimental methods to achieve molecular slowing to v cap with a small number of spontaneously emitted photons [26][27][28], thus reducing the number of required repumping lasers. Here we present a novel cooling scheme that uses multifrequency light to rapidly dampen molecular motion in a wide range of velocity classes v eff v cap , while minimizing the number of spontaneous decay cycles. The proposed Suppressed Emission Rate (SupER) molasses could be either combined with coherent slowing techniques or used with previously magneto-optically trapped species to capture and cool molecules with v > v cap,MOT .

B. Coherent Optical Forces
Widely utilized optical slowing and cooling methods for atomic gases usually use a single optical frequency to address a specific "two-level" transition. Such radiative methods are characterized by the maximum force F rad = kΓ sp /2 affecting velocity classes within ∆v rad ≈ Γ sp /k, with the intrinsic spontaneous decay rate Γ sp limiting the maximum force as well as the capture range [17]. While a conservative dipole force arising from the gradient of the light shift can lead to strong confining forces, its utility in cooling atomic or molecular motion is severely limited since it averages out to zero over a spatial scale larger than light wavelength λ [25]. However, already more than thirty years ago it has been theorized that the dipole force can be "rectified" to maintain a constant sign over position scales much larger than λ by adding a second light field to spatially modulate the atomic energy levels and, therefore, the sign of the detuning for the initial dipole force laser field [29,30]. Shortly afterwards, Grimm and co-workers have conclusively demonstrated the effect of the rectified dipole force (RDF) on a sodium atomic beam achieving F RDF ≈ 4F rad [31], a factor of 2.5 lower than the initial prediction due to the presence of transverse atomic velocities larger than v ⊥ ∼ Γ sp /k [32]. In order to remedy the issue of a small velocity capture range of the RDF, other methods for generating large coherent optical forces have been proposed that realize coherent control of light -atom momentum exchange by tailoring the inversion of the atomic populations [33].
The use of counterpropagating amplitude-modulated light waves, leading to a stimulated light pressure on atoms, has been proposed as a viable method to achieve large force magnitudes F F rad over a wide velocity range v ⊥ Γ sp /k [34,35]. While the magnitude of the stimulated bichromatic force (BCF) can be explained using an intuitive resonant optical π-pulse interpretation [35,36], an understanding of the large velocity capture range requires a doubly-dressed atom picture [37]. In the simplest case, a two-level system interacts with collinearly superimposed bichromatic standing waves with equal intensities I BCF and symmetrically detuned by ±|δ Γ sp | from atomic resonance. By imposing a relative phase offset χ between the counterpropagating laser fields, the directionality of the force can be controlled by fixing the relative timing between the resulting beat pulse trains ( Fig.   1) [36]. Choosing I BCF and δ such that the Rabi frequency integrated over a single beat pulse area satisfies Ωt π ≈ π, efficient transfer of atomic population between the ground and excited states can be achieved at a rate of δ/π Γ sp [17]. Since each directional π pulse transfers k momentum to the atom, the order of magnitude for the bichromatic force 2 is F BCF ∼ kδ/π F rad . Even though a spontaneous emission rate can be significantly suppressed by reducing the excited state fraction with a properly designed pulse sequence [26,38], any spontaneously emitted photons will lead to quantum state decoherence and potential reversal of the momentum transfer direction. For an asymmetric χ phase choice required for a directional momentum transfer this leads to only order of unity reduction in the net magnitude of photon transfers averaged over a time greater than 1/Γ sp [36]. However, such processes can have important consequences for the limiting temperature of the ensemble when coherent optical forces are employed for cooling of lab-frame velocity. Careful understanding of limiting temperatures in stimulated transfer cooling methods has proven challenging [39,40], and here we develop a novel method of doing so using a continuous-time Markov Chain model detailed in Apps. A and B.
When considering BCF and other stimulated light forces like RDF, it is important to properly account for the atom's or molecule's finite velocity that could significantly affect the magnitude of the experimentally achievable decelerations (seen in the initial RDF ex- 2 As shown in App. A, the expression presented here is in fact exact for a two-level system: F BCF,2−level = kδ/π. periments, for example [32]). In order to create force profiles accurately depicting velocity dependence necessitates solving Liouville-von Neumann equations for density matrix evolution in the rotating wave, fixed-velocity approximations, followed by obtaining the force averaged over the ensemble F = −tr (ρ∇H), where ρ is the density matrix [25,37]. However, intuitively, the velocity capture range for the bichromatic force can be interpreted as arising from the relative dephasing between consecutive beat notes. Once the Doppler shift kv becomes comparable to the Rabi frequency Ω ∼ δ, the π pulses no longer lead to efficient population transfer, thus limiting the affected velocity range to v ∼ δ/k Γ sp /k [41] under conditions of large BCF detuning δ Γ sp .

C. Cooling Properties of the Bichromatic Force
Since the bichromatic force does not vanish for atoms at rest and involves mostly coherent state transfers, it may seem surprising that rapid cooling of atomic beams has been achieved using BCF configurations 3 . However, the sharp edges of the BCF profiles can be used for compressing velocity distribution and achieving cooling of atomic motion. The frequencies of counter-propagating dual-frequency beams can be offset by opposite amounts, creating a situation where the atom or molecule undergoes efficient π-pulse transfers from ω 0 ± δ beat notes only at non-zero velocities (Fig. 1). The Doppler shift experienced by an atom or molecule is ±kv, so by letting ∆ = kv 0 we can center the force profile around a chosen nonzero velocity v 0 . By shifting the bichromatic force profile to be centered around a non-zero velocity, efficient longitudinal cooling of atomic beams has been demonstrated [36,41]. In at the profile's center. 3 In fact, transverse cooling of the metastable helium beam using the bichromatic force has been demonstrated in the regime of 1.5 emitted photons [42,43].

FIG. 1:
Interaction of a moving two-level system with a bichromatic light field. Pulse trains are created by a pair of counter-propagating two-color (±δ) beams offset from each other by a phase χ. To center the force profile around a non-zero velocity, pulses on one side have to be red-detuned, while pulses approaching from the other side have to be blue-detuned.
The use of large stimulated optical forces for achieving 3D cooling and kelvin-deep trapping of atoms and molecules was one of the primary motivations for the initial extensive development of such methods [29,34]. However, despite thirty years of research, the application of coherent stimulated forces to zero-velocity cooling (i.e. compression of velocity distribution towards zero lab-frame velocity) and confinement of atomic and molecular samples has been limited [27,44,45]. Partlow and co-workers have performed a landmark experiment on a helium beam to use spatially separated, shifted bichromatic force profiles with opposite phase χ and frequency shift ∆ [39]. Such a 1D collimation scheme required two sequential interaction regions acting on atoms with positive and negative initial velocities, respectively, and thus leading to experimental results emulating the effects of optical molasses. However, as pointed out by Partlow and co-workers [39], the underlying physical process was not resulting in a true damping force for velocities of interest and led to a different physical behavior than optical molasses cooling.
Here, we use the inherent multilevel structure of molecular radicals that limits cooling efficiency of traditional Doppler molasses to propose a novel laser cooling scheme with significantly higher velocity range and damping coefficient. By addressing two separate, yet radiatively coupled, two-level systems with polychromatic optical fields, we discover that it is possible to achieve a large velocity damping force with Suppressed Emission Rate (Su-pER). Furthermore, in the experimentally accessible regime, we demonstrate the feasibility of damping molecular motion to millikelvin temperatures on microsecond timescales. We show how SupER molasses force profiles can be created in a simple 4-level system, develop a new mathematical model for estimating the final temperature, and perform Monte Carlo simulations of the cooling dynamics to confirm the analytical estimates. Throughout the paper we use the barium monohydride (BaH) molecule as a test species for our time-dependent density matrix calculations, but also suggest a more general level scheme common to many diatomic and polyatomic radicals that could be utilized to create large molasses-like forces.
Therefore, we identify a way to use the internal complexity of molecular systems to enable their efficient quantum control for a wide variety of proposed applications [46]. Our work makes an important step towards experimental realization of kelvin-deep macroscopic (r λ) optical traps for molecules proposed more than thirty years ago [29,34].

II. SUPPRESSED EMISSION RATE MOLASSES
To realize cooling force profiles arising from a rapid coherent momentum exchange between light fields and molecules, we identify an appropriate multi-level quantum system that would allow us to combine two asymmetric (shifted) force-velocity profiles (e.g. shown in Fig. 2) without significant reduction in force vs velocity characteristics. As initially pointed out by Partlow and co-workers [39], it is impossible to simultaneously apply both force profiles on one single 2-level system without drastically perturbing individual force profiles. In order to circumvent this limitation, instead we consider two almost independent 2-level systems that are coupled only by spontaneous decay as shown in Fig. 3

FIG. 2:
A shifted bichromatic force profile for δ = 100 Γ, Ω = 3/2δ, χ = 45 • with opposing beams detuned by ∆ = −40 Γ. The profile is centered at v = ∆/k = −40 Γ/k and was obtained by numerically solving Liouville-von Neumann equations for density matrix evolution in the rotating wave, fixed-velocity approximations. The maximum force is equal to kδ/π ≈ 63 kΓ/2 and the width of the profile can be estimated to be equal to δ/2k = 50 Γ/k. Sharp vertical spikes are Doppleron resonances arising when integer numbers of red-detuned absorptions and blue-detuned emissions of photons (or vice versa) are resonant with the transition [47]. This will occur when (δ + kv) / (δ − kv) is rational. However, in previous measurements of the bichromatic force profiles [48,49], these narrow Doppleron resonances were not observed and are not expected to have any significant effect on real physical systems. component shown in Fig. 1. In order to obtain two asymmetric profiles those 2-level systems have to have opposite signs of the detunings and phase differences: ∆ 1 = −∆ 2 and χ 1 = −χ 2 .
While in the toy model δ has to only be much larger than Γ, in real systems we need it to be larger than naturally occurring energy splittings, such as hyperfine splitting, which quite often are on the order of couple Γ. Our specific choice, while arbitrary, should be applicable in many situations and be realizable using off-the-shelf acousto-optic modulators (AOMs). |E2 We found the most optimal profile for ∆ 1 = −15 Γ and depict it in Fig. 4. We should note that in such symmetrized system given no additional selection rules the light from both 2-level subsystems would couple to the other subsystem. In a real system, like that in BaH, light frequencies are vastly different and don't couple to both transitions simultaneously, yet still enable realization of SupER molasses as discussed below. As shown in Fig. 4, the resulting force profile shows remarkably strong forces and high capture velocities -the force peaks at around 35 kΓ/2 at velocity of about ±20 Γ/k, while at the same time the slope around zero velocity, representing velocity damping coefficient, is quite linear and steep. We can benchmark this force vs velocity curve against normal optical molasses realized with radiative forces. The comparison is shown in Fig. 5 after smoothing the BCF-induced force profile with a moving average filter in order to smooth out sharp spikes arising from multiphoton reasonances. Qualitatively, the obtained force profile perfectly resembles the radiative Doppler molasses, but on a much bigger scale: the slope near zero velocity is 16 times higher compared to radiative force in Fig. 5, peak forces are substantially bigger, and so are capture velocities. 4-level toy model and a regular optical molasses force profile (red graph and axes). BCF profile was obtained for ∆ ± = ∓15 Γ, δ = 100 Γ, Ω 0 = 3/2δ and χ ± = ±45 • and was smoothed using moving average filter, while the radiative force profile was drawn for I = I sat and detunings δ ± = ±1 Γ. Both scales differ by a factor of 20. While in the case of radiative forces, the slope is about 0.44 k 2 /2 for given parameters, in the BCF force profile it is close to 7 k 2 /2, a factor of 16 higher.
Rapid adjustment of the shape and magnitude of the cooling profile is easily achievable experimentally. Capture velocity and the damping force are modified when changing profileshifting detuning ∆ (see Fig. 1), which can be controlled by an AOM, by changing its RF drive frequency. The magnitude of the force depends on the choice of BCF detuning δ and the Rabi rate Ω for each two-level system shown in Fig. 3. In Fig. 6 we present comparison of smoothed force profiles for ∆ ranging from 10 Γ to 30 Γ. These resulting force-velocity profiles are effectively created by a sum of two opposing, shifted and re-scaled 2-level BCF profiles, as shown in Fig. 7. The scaling effect appears there, because in our 4-level model the atom or molecule spends on average less time in either of the 2-level subsystem interacting with their own respective bichromatic fields than in a simple situation of a BCF-driven 2level system. In fact, using the π-pulse approach (described in detail in Appendix B), we can predict that for our toy model such factor will be exactly equal to 4/7, and so the peak forces expected for any detuning δ are F BCF,mol = 4/7 kδ/π. However, this scaling factor has some limitations. In case of fields with more than 2-colors, the direct solution of 4-level system yields forces higher than one would obtain by re-scaling a 2-level system solution, thus reaffirming the necessity to perform a full quantum calculation with all levels and laser frequencies as done here rather than using scaled analytical results from isolated two-level systems.    The force profile can also be obtained for any 2n-color forces, though it is not immediately obvious what the benefits of adding additional frequencies are. If we assume that we always operate with a certain power per frequency component, in a simple 2-level system moving from 2-to 4-color laser fields increases the force and velocity range quite substantially while decreasing the time spent in the excited state [38]. Indeed, as shown in Fig. 8, the 4-color profile does not look that much different from the BCF profile. However, the peak force is much larger and remains such for higher velocities, which could result in higher capture velocity. In Fig. 9 we compare 4-color force profiles for various shifts ∆. We found the highest peak forces for |χ| = 25 • and Ω = 1.16 δ, where Ω is rate of every component.  The 4-color force profiles can be quite strongly adjusted by appropriately changing the phase χ, which experimentally is controlled by a relative length of the optical delay line between the counter-propagating multi-frequency laser beams of the same color (see Fig. 1), and Rabi rate Ω. In Fig. 10 we show a much wider profile with capture velocities as high as 60 Γ/k (which for molecule like BaH is equivalent to ≈ 72 m/s) and forces of the order of 50 kΓ/2 that was obtained for |χ| = 30 • and Ω = δ. Finally, we present comparison of 2-color molasses profile, narrow 4-color force profile and wide 4-color force profile in Fig. 11.
All these SupER molasses have high capture velocities, high peak forces and steep slopes ranging from 7 k 2 /2 for the BCF to 9.5 k 2 /2 for the narrow 4-color profile, enabling rapid damping of molecular motion.

III. TEMPERATURE IN SUPER MOLASSES
Having developed a mathematical approach to the π-pulse model of the polychromatic force dynamics using continuous-time Markov chains (CTMC) (App. A), we are now able to estimate the final temperature of atoms or molecules under the influence of SupER molasses.
The diffusion coefficient, which determines the limiting temperature, should not only include terms related to spontaneous emission from the excited states, but also the term specific to this model, which is related to the uncertainty in "position" in the excitation-stimulated emission cycle. This term is effectively the variance in momentum transfer from the optical field to the atom or molecule. In our model, the system quickly relaxes to a steady state that's achieved at long times t 1 Γ 12 +Γ 21 and the variance term can be written as (App. B): where σ 2 (ε) is dependent on fraction of time ε spent in the excited state per cycle (it is not the ensemble average excited state population ρ ee appearing in the density matrix and it is discussed in Appendix A as well as in Ref. [50]). Apart from the contribution from PCF, total variance of the momentum transfer will have a contribution arising from the spontaneous decay as well: with ρ ee being the time-averaged excited state population (in case of the 4-level system we have shown, ρ ee = ρ E1 + ρ E2 ). We should also note that the variance shown here is variance of momentum transfer in one dimension and is the reason behind the 1/3 factor in the second term on the right hand side in the equation above.
Obtained result shows a quadratic dependence on kδ which is consistent with results experimentally shown in Ref. [39] and first estimated in Ref. [51] for dipole forces. In case of the bichromatic force used in the helium collimation experiment, authors obtained diffusion , while the value we obtain from our model is ≈ 0.6( kδ) 2 /Γ.
From the calculated variance we can obtain the diffusion coefficient D: which is then related to the limiting temperature T L . The PCF molasses force around v = 0 is linear and can simply be written as F = −βv, where −β is the slope. Using the definition of the diffusion coefficient and by assuming mass of the atom or molecule is M , at equilibrium one can write: By associating the limiting temperature T L with kinetic energy, i.e. M v 2 = k B T L and using the definition of Doppler temperature T D = Γ/2k B , we get: In the equation above the parameter β is measured in natural units for this problem - is on the order of 10, in Eq.(2) the spontaneous emission term is negligible given detuning δ typical for this problem. Slope β that appears in this equation can be estimated (following [39]) to be β ∼ δ/4πΓ k 2 /2 for presented profile, which leads to T L ∼ 8 σ 2 (ε) δ/π T D . For example, for BCF molasses presented here σ 2 (ε BCF ) ≈ 20 and δ = 100 Γ, so T L ≈ 5.09 × 10 3 T D , which for BaH is about 137 mK. For a more realistic and asymmetric system σ 2 (ε BCF ) ≈ 21 and so, given the slope of similar magnitude, the final temperature would be almost exactly the same. The linear dependence of T L on the detuning shows that it might be beneficial to keep δ as small as possible, while still keeping δ Γ condition fulfilled.
While the estimated limiting temperature appears high compared to the Doppler limit, it can actually be made lower. The temperature at equilibrium in such system is not only determined by the slope of the profile, but also by how strongly the polychromatic forces act on the atom or molecule around zero velocity. In the derivation of the formula for variance of momentum transfer (App. B) we assumed that in every state considered the force acting on an atom or molecule is 2 kδ/π. However, around zero velocity that does not need to be the case. Assuming that the force around v = 0 is equal to F 0 , we can simply substitute kδ/π ∼ F 0 . If we removed all the constants appearing naturally in both β and F 0 , we'd simply obtain: clearly showing that the final temperature mainly depends on the ratio F 2 0 /β. This shows us that effectively, depending on how the basic 2-level system profiles ( Fig. 7) are aligned to create a full molasses-like forces, we obtain different values of β and F 0 . For small shifts ∆ both force and slope are high as shown in Fig. 7(a). If we make our detuning too high, like in Fig. 7(c), the force around zero velocity becomes small, but because the profiles are far from the center, the slope of the effective profile is very small as well.
In between there should exist an optimal configuration, where slopes of single PCF profiles are each other's continuation. There, the slope should remain high, while the force F 0 should be relatively small and Fig. 7(b) depicts such situation. In this configuration the smallest temperature ought to be achieved. Given that width of BCF profiles is approximately δ/2k, we should expect that the most optimal cooling forces will appear around ∆ ± = ∓δ/4k. To confirm these estimates, We analyze the effective forces, damping coefficients and limiting temperature T L using Monte Carlo simulations of an exact realization of the 4-state model described in App. B.

IV. MONTE CARLO SIMULATIONS
To prove that the estimated temperature follows the model we have created we present results of a Monte Carlo simulation of an ideal 4-level toy model for a 2-color light field.
In it we assumed that the system consists of four states (described in App. B). While in the π-pulse model we assumed that the force in every state is 2 kδ/π, here we assume it is velocity-dependent and follows a typical BCF force profile, re-scaled in such a way that at maximum it is equal to the mentioned 2 kδ/π value. Using such values results in effective force profiles seen in figures provided before. It is important to note that the temperature model and obtained formulas presented in App. A and B work in the regime where the interaction can actually be described in the framework of a π-pulse model. This occurs for specific parameters (such as χ = 45 • for 2-color fields) and for velocities, for which the force is at its maximum. The model might not hold at the edges of force-velocity profiles (|v − v 0 | > δ/2k), which is the regime we expect to be in in the cooling process (v ∼ 0). Therefore, the proportionality constant in Eq. 3 could be difficult to predict, given that σ 2 (ε) stems from the near-perfect π-pulse behavior. Additionally, the slope β that appears in the formula might not be the effective slope that can be read directly from provided effective force-velocity profiles.
Occurrence of transitions between different states in the simulations was assumed to follow Poissonian statistic with rates given in the appendix A. We have also included recoils from spontaneous emission events even though they play a very limited role. The force profiles used were re-scaled versions of those seen in Fig. 7(a) and Fig. 7(b), which should show a typical (the former) and a perfect (the latter) configuration for the molasses. In the simulation these profiles were separate and, like in mentioned figures, already smoothed with a moving average filter. We have also assumed that the molecule undergoing the cooling process is BaH with Γ = 2π × 1.15 MHz, M = 139 AMU and k = 2π/1060.7867 nm. Experimentally, to obtain these profiles for our test species the lasers would have to be detuned from the resonance by δ ≈ 115 MHz, and, assuming beams with uniform power density and diameter of 5 mm, have a total power of ∼ 4.4 W.
The effective damping coefficients (slope) that we expect for ∆ = ±15 Γ are β ≈ 4.9 k 2 /2 and β ≈ 3.1 k 2 /2 for ∆ = ±20 Γ, while the force around v = 0 is expected to be F 0 ≈ 103.4 kΓ/2 and F 0 ≈ 42.4 kΓ/2 respectively. If the developed model holds in both described regimes, we expect that the limiting temperature will be (in the natural units): which should lead to T L ≈ 147 mK for ∆ = ±15 Γ and T L ≈ 39 mK for ∆ = ±20 Γ. Eq. 4 can also be written using the notation used in the appendix as: where F IV is already the effective and re-scaled force around v = 0 that is obtained in a 4-state model (App. B) and that can be read directly from the presented force profiles, σ 2 IV (ε) = σ 2 (ε)/2 and is equal to 10.06 for BCF in a symmetric system, and µ IV (ε) is the re-scaling factor equal to 4/7 for BCF in such system.
In the simulation we started with molecules distributed with σ v = 5 m/s , which for BaH corresponds to temperature of approximately 0.4 K. The simulation was run in steps of 5 ns for a total time of 200 µs, which was more than enough to reach the final limiting temperature for both considered SupER molasses configurations. Figure 12 shows molecular temperature at different times averaged over multiple simulations for ∆ = ±15 Γ molasses (top frame) and for ∆ = ±20 Γ molasses (bottom frame). The limiting temperatures obtained were 125.84 mK and 37.91 mK respectively, relatively close to the CTMC model estimates. The model was also proven to be correct by the steady state populations which quickly relaxed to predicted levels of η C 1 = η C 2 = 11/28 and η W 1 = η W 2 = 3/28 (App. B). Assuming that the temperature follows an exponential decay curve, i.e.
To investigate the capture velocity we have also performed simulations for ∆ = ±20 Γ molasses with molecules starting with σ v = 30 m/s corresponding to T 0 ≈ 14.7 K for BaH.
Simulation was performed to reach final time of 200 µs, which was enough to show the approximate capture velocity. Figure 13 shows comparison between velocity distributions at t = 0 and at times of t = 10, 20, 40, 60 µs. We can see that the molecules with |v| 40 Γ/k (∼ 47 m/s for BaH) accumulate around v = 0 showing that we can consider this to be the effective capture velocity of the SupER molasses for our parameters. In general, we could expect the capture velocity to be equal to the typical width of the force vs velocity profile of ∼ δ/2k. Note that the simulated cooling time of SupER molasses for BaH agrees with the characteristic BCF timescale of τ BCF = π/ (4ω r ) ≈ 100 µs [25] where ω r ≡ k 2 / (2M ) ≈ 2π × 1.3 × 10 3 s −1 is the recoil frequency for BaH. The value of τ BCF represents the timescale over which a molecules is accelerated across the full velocity range ∼ δ/k and gives an approximate upper bound on the cooling time in the molasses configuration.

V. APPLICATION TO REAL MOLECULES
With the theoretical results established for an ideal 4-level system (the minimum number of states required for the proposed scheme), we relax the simplifying assumptions in order to apply our cooling method to real molecular systems. We develop a general scheme to obtain rotational closure (while satisfying the requirements dictated by the 4 level toy model presented above), at the same time accounting for spin-rotation, fine and nuclear hyperfine structure for molecules with an unpaired electron spin. As we show below, additional internal substructure present in molecular radicals enables generic experimental realization of the SupER molasses cooling scheme.
14: Energy correlation diagram for alkaline earth metal (M) atoms bonded to ionic ligands (L). As can be seen from the diagram, multiple closely-lying electronic states generically arise for ML molecules with diverse constituents and symmetries. Using the lowest two non-degenerate excited electronic states (e.g. 2 Π and 2 Σ) for molecules of C ∞v or 2 A and 2 A for C S molecules will allow application of the molasses methods described here. While we indicated selective ligands in the diagram, other possible constituents include M = Ca, Sr, Ba and L = F, OH, CCH, OCH 3 , CH 3 , SH. For a comprehensive list of ML monovalent derivativies of the alkaline earth metals that have an electronic structure of the type presented here refer to Refs. [52,53].

A. Rotational Level Schemes
As detailed above, in order to realize large cooling forces primarily due to coherent momentum exchange between multi-frequency laser beams and molecules we must work with two separate excited states that each decay into both ground states. One way this can be accomplished for molecular radicals is when the |G1 and |G2 states in Fig. 3 are the rotational ground (N = 0) and second (N = 2) excited levels in the ground vibrational manifold, while the excited states consist of the ground vibrational, first rotational (N = 1) level of an excited state manifold with two different sufficiently separated sub-manifolds (e.g. spin-rotation components). Angular momentum selection rules dictate that these excited states will decay to both ground states we have selected. This scheme is shown, using as general of notation as possible, in Fig. 15. This scheme provides for the general requirements outlined by the toy 4-level model, but there are additional nuances introduced by potential hyperfine and spin-rotation splittings in various states. Figure 14 demonstrates energy correlation diagrams for molecular radicals of different structural symmetries. In the presented scheme, the B 2 Σ + electronic state (second excited electronic state arising from the mixing between the (n − 1) dσ and npσ orbitals), which exists in all diatomic and polyatomic molecules that have been considered for optical cycling applications thus far [14,[54][55][56][57], is chosen as the excited state. Angular momentum selection rules dictate that only the included hyperfine and rotational states take part in the cycle, although there are some basic criteria for this scheme to work. Firstly, because we would like to use different J manifolds as |E1 and |E2 levels in Fig. 3 of our 2-level subsystems, the spin-rotation splitting in the excited electronic state ∆ BΣ J has to be much larger than the either of the PCF detunings δ BΣ 1/2 or δ BΣ 3/2 . Otherwise, the N = 0 state might couple to J = 3/2 or N = 2 to J = 1/2, reducing the velocity damping coefficients. Secondly, for the same reason we have to make sure that both transitions are separated in frequency by more than the PCF detunings, which means that rotational splitting ∆ N and spin-rotation splitting ∆ BΣ J have to be quite different.
While the branching ratios in this scheme are not balanced, polychromatic forces acting on both subsystems can be balanced by an appropriate choice of PCF detunings and Rabi rates (App. B). However, the FCF for the X 2 Σ + (ν = 0) ↔ B 2 Σ + (ν = 0) transition have to be high enough to allow multiple scattering events to occur. Albeit, they don't have to be perfect -one of the reasons polychromatic fields are extremely promising in the context of molecules is that they generate high forces while suppressing spontaneous emission [28].
FCFs of F BΣ 00 0.95 should be sufficiently high to allow the forces to create observable effects. Additionally, a repump laser might be added and population recycled, as was shown in multiple diatomic and polyatomic systems [5].
Finally, hyperfine splittings and potential creation of dark states in the N = 2 manifold have to be discussed. Ideally, we would like the hyperfine splittings in a sub-level to be smaller than the PCF detunings, e.g. ∆ XΣ F 1/2 ≈ ∆ BΣ coupling. In the case of N = 2 ↔ J = 3/2 transitions we also have the ground state spin-rotation splitting ∆ XΣ J to take into account. If it is smaller than the PCF detunings, then all hyperfine states will be coupled resulting in creation of 12 dark states. If the splitting is larger than the detunings, we can couple only the J = 3/2 manifold in the ground state without creating any dark states, but with population accumulating in the J = 5/2 manifold.
In both cases we could use an auxiliary transition, driven by Ω aux , to a different excited state, like A 2 Π 1/2 , which would bring the population back into the cycle. Alternatively, dark state remixing method could be used such as adding a magnetic field, polarization switching or microwave-induced AC Stark shift. If the spin-rotation splitting in the B 2 Σ + state is not sufficiently large, a different scheme can be utilized, shown in Fig. 16. Here, excited states |E1 and |E2 are two A 2 Π Ω states with a different total angular momentum projection of Ω = 1/2 and Ω = 3/2. In fact, in some molecules (like BaH) such scheme might work better, due to more beneficial FCF's compared to the B 2 Σ + state [58]. In this configuration, we exploit the fact that the excited states are parity doublets and given positive parity of the ground states, we can choose excited states of negative parity to create our light-coupled 2-level systems. In all other aspects this scheme is analogical to the one described before, and so the same constraints and criteria apply, with the spin-rotation splitting in the excited state being replaced by Ω-splitting ∆ Ω . 16: Alternative realization of polychromatic molasses-like forces in a real molecular level structure using two parity-doublet A 2 Π electronic state with different total angular momentum projections Ω.

B. Analysis of Bichromatic Forces in Barium Monohydride
Motivated by the prospects of ultracold hydrogen production via molecular laser cooling followed by photo-dissociation [59], recently there has been an increased experimental [11,60] and theoretical [58,61,62] interest in direct laser cooling and trapping of alkaline-earthmetal monohydrides. High mass imbalance, low Doppler cooling limit, and small photon recoil velocity make BaH an extremely attractive candidate for producing ultracold atomic hydrogen via zero-energy photo-fragmentation [59]. Unfortunately, the very same inherent molecular characteristics make laser cooing of BaH experimentally challenging [11]. However, fine and hyperfine structure of BaH in the rotational states involved in the optical cycling process [63] together with technically accessible transition wavelengths make it an ideal candidate for stimulated slowing and cooling using polychromatic optical forces.
In the electronic ground state X 2 Σ + , the ν = 0 N = 1 rovibrational state has an exceptionally large spin-orbit splitting of 8.6 GHz, while the hyperfine splitting is unresolvable in the J = 1/2 state [63]. This allows the transition from X 2 Σ + J = 1/2 to A 2 Π 1/2 J = 1/2, where the hyperfine splitting is similarly unresolvable, to be addressed simultaneously on all transitions with equal detuning between transition and carrier frequency, while leaving the J = 3/2 states unperturbed even when the bichromatic detuning is significantly larger than the decay lifetime (∆ J /Γ AΠ = 1180, leaving plenty of range over which Γ AΠ δ ∆ J holds), and thus enabling realization of the stimulated force significantly larger than any possible radiative force.
When this transition is driven with π-polarized light, a set of four separate, radiatively cross-coupled two-level transitions are all driven at equal strength and equal transition frequency. As a result, applying a set of BCF optical fields to this transition results in a force which (when neglecting any off-transition decays from the A state) is nearly identical to that achieved in simple two-level BCF. This includes that there is no formation of dark states, avoiding the requirement of dark-state destabilization as has been needed in prior application of BCF to molecules [64]. This was verified by numerical simulation of the BCF force on the four ground state, four excited state system. Our numerical simulations solved the Liouville-von Neumann equations for density matrix evolution in the rotating wave, fixed-velocity approximations, as in previous simulations of the BCF on molecules [65].
An obvious weakness of this model is that, with the J = 3/2 states unaddressed, the off-diagonal decays from the A state are far from negligible. Even considering only decays to X 2 Σ + v = 0 N = 1, one third of spontaneous decays of the excited states should end in J = 3/2. If this is allowed to continue undisturbed, the system will quickly go dark and stop feeling force after all population is pumped out of J = 1/2 into J = 3/2 and a small fraction into other rovibrational states.
In the context of creating a sustained force, this can be remedied by using additional optical fields to drive transitions from J = 3/2 in such a way that population eventually returns to the BCF-driven transition. In a simple BCF scheme, this can be done by addition of a CW optical field which drives J = 3/2 to B 2 Σ + N = 0 J = 1/2. This state also decays primarily to X 2 Σ + v = 0 N = 1, with one third of these decays ending in J = 1/2, back in the BCF cycle. This comprises an indirect repumping scheme for BCF as discussed in [65] and previously implemented for SrOH BCF deflection [66]. This scheme is illustrated in Fig. 18 but with weak repumping beams to address the J = 3/2 to B 2 Σ + N = 0 J = 1/2 transition. Notably there are four sets of states, ground and excited states in both the BCF and repump transitions as depicted in Fig. 3. These sets of states will be referred to as |G1 , |E1 , |G2 , |E2 and the labels g, e, g R and e R will refer to the time-averaged ensemble population in each set, correspondingly.
With sufficiently strong optical fields, so that Ω Γ in each case, population returning to the BCF or repump cycle from the opposite cycle can be taken as a small perturbation to the population dynamics, and overall population can be estimated accurately by assuming the relative populations in each cycle will be identical to that which they would have absent the other cycle and that relative population between cycles is determined by equilibrium rate equations. In particular, eΓ 12 = e R Γ 21 , where Γ 12 is the decay rate from |E1 to |G2 and Γ 21 is the decay rate from |E2 to |G1 as shown in Fig. 3.
Optical dark states will exist in the J = 3/2 to B 2 Σ + transition, for any fixed choice of repump polarization. Assuming a remixing magnetic field and a saturated CW repump, the populations in the repump transition will equilibrate to have a proportionality equal to that of the number of states: e R /g R = N |E2 /N |G2 . In this case, with four excited states and eight ground states, 1/3 of the population in the repump cycle will be in the excited B 2 Σ + state at any given time.
On the BCF transition manifold |G1 → |E1 , the proportion of population between ground and excited states depends on the BCF (or polychromatic force) driving. This will produce a characteristic time-average excited state population P e = e/(e + g). The J = 1/2 to A 2 Π 1/2 transition in BaH behaves nearly identically to a two-level system in response to BCF, as discussed above. The optimal BCF or 4-color PCF optical fields in a two-level system are known from previous works, along with the time-average excited state populations that results [38]. After the decay of any transient behavior, the equilibrium populations between the two cycles will be reached when cross-decay occurs at equal rates.
Taken together, these considerations are sufficient to determine the participating fraction, i.e. the fraction of molecular population which is in the BCF cycle at any given time: The effective time-averaged force at a given velocity can then be taken to be equal to the time-averaged force that would be achieved if the BCF cycle were closed, multiplied by this participating fraction which is a function of the time-averaged excited state fraction in the closed cycle at that velocity. This can be compared to the maximum radiative force that would result were both the X 2 Σ + J = 1/2 ↔ A 2 Π 1/2 and X 2 Σ + J = 3/2 ↔ B 2 Σ + transition manifolds driven with resonant, saturated CW optical fields. In that case, P e would equal one half, and the expected radiative force due to both transitions would be, where k and Γ are the wavenumber of and decay rate along the transition in question, This force can significantly exceeded the optimal radiative force with experimentally achievable BCF irradiances. Figure 17 shows simulated effective force profiles for a bichromatic detuning of 189 MHz, which would require a per-beam irradiance of 22.5 W/cm 2 , or in other words a total summed irradiance of 90 W/cm 2 across all four BCF beams, to have the optimal BCF Rabi frequency at this detuning. Given the experimentally realized laser power at 1060 nm [67], which is the wavelength for the X − A optical cycling transition in BaH, our calculations indicate strong feasibility of using the described BCF-driving scheme to achieve rapid slowing of cryogenic BaH molecular beam in a short (∼ few cm) distance.
Using realistic experimental parameters (5 W and 1 mm radius beam), we anticipate that achieving a total summed irradiance of 500 W/cm 2 is feasible, leading to potential for even larger force enhancements.

C. SupER Molasses in BaH
Given a high potential for an effective realization of BCF in BaH on the X 2 Σ + ↔ A 2 Π 1/2 transition as described in Sec.V B, a different optical scheme can be used to realize the SupER molasses cooling configuration. As presented in Fig. 18, instead of using two states with different rotational quantum numbers N , we choose to use two J states in N = 1 rotational state, which are separated by about 8.6 GHz [63]. As the excited states we use two different electronic states: A 2 Π 1/2 with positive parity and B 2 Σ + in its ground rotational level N = 0.
In this system, many of the criteria listed before are fulfilled -both FCFs are greater than 0.95 [58] and both transitions are far apart in the frequency space (λ AΠ ≈ 1060.7868 nm and λ BΣ ≈ 905.3197 nm). As was detailed in Sec. V B, polychromatic forces can be created using the excitation to the A 2 Π 1/2 state. The hyperfine splitting ∆ AΠ F in the A 2 Π 1/2 electronic state and ∆ XΣ F 1/2 in J = 1/2 manifold of the X 2 Σ + state are both very small (less than 4 Γ AΠ ∼ 2π × 4 MHz [63]), and so by using a π-polarized light fields with δ AΠ detuning having any reasonable value much larger than Γ AΠ these forces will be created.
The X 2 Σ + ↔ B 2 Σ + transition considered will also create dark states in the J = 3/2 manifold. Fortunately, the g-factors are large enough [63] to provide efficient remixing with the help of a magnetic field of modest strength. Given that we already need high detunings δ BΣ , the Zeeman splittings should not cause any additional problems. Finally, while the branching ratios in this case are symmetric, the decay rates are different for both electronic states (Γ AΠ ≈ 2π × 1.15 MHz and Γ BΣ ≈ 2π × 1.21 MHz); although, as was mentioned before, this asymmetry can be easily adjusted for by an appropriate choice of detunings and Rabi rates. To obtain force profile in a realistic BaH quantum level system, we have again solved the Liouville-von Neumann equations for density matrix evolution for a system that included Zeeman sublevels in states depicted in Fig. 18. For simplicity, we have solved it by assuming that F AΠ 00 = F BΣ 00 = 1. We have chosen to center the X 2 Σ + ↔ A 2 Π 1/2 transition on the F = 1 ↔ F = 0 frequency, while also assuming that ∆ AΠ F = 2π × 2 MHz and ∆ XΣ F 1/2 = −2π × 2 MHz. The X 2 Σ + ↔ B 2 Σ + transition was centered at frequency placed symmetrically between F = 2 ↔ F = 1 and F = 1 ↔ F = 1 transition frequencies.
We have assumed a presence of 12 G ambient magnetic field that defined the quantization axis. Zeeman splitting was obtained using experimentally obtained effective linear g-factors [63]. The X 2 Σ + ↔ A 2 Π 1/2 transition light field was polarized along the quantization axis, while the other light field was perpendicular to it. Finally, we have Doppler-shifted the frequencies by appropriate amounts, that is ∆ 1 = −δ AΠ /4 ≈ −2π × 67.5 MHz and ∆ 2 = δ BΣ /4 ≈ 2π × 63.1 MHz, and chose χ 1 = −χ 2 = 45 • . The SupER molasses force profile obtained is shown in Fig. 19.
The obtained force profile is quite symmetric, despite the fact that both bichromatic fields act on quite different level structures. It is also linear around zero velocity. We can estimate the slope to be β ≈ 3.24 F rad /(m/s) ≈ 3.47 k 2 AΠ /2, and while we do not know the exact value of F 0 , we can place an upper bound equal to the maximum force seen in the profile, i.e. F 0 45 F rad = 39.45 k AΠ Γ AΠ /2. Using Eq.(5) with σ 2 (ε) ≈ 21 (asymmetric 4-level system) we can place an upper bound on the temperature to be T L 193.8 mK.
We also see that the capture velocity gets as high as 120 m/s, which is consistent with the δ/2k estimate provided earlier. The temperature equivalent to capture velocity in these molasses is T cap = M δ 2 /4k B k 2 , and for v cap = 120 m/s in BaH is equal to T cap ≈ 60 K, providing further evidence that the SupER molasses method does indeed realize a way to cool and confine molecules in kelvin-deep optical potentials as thought after for more than thirty years following initial speculations by Kazantzev [29] and Voitsekhovich [34]. The unusual structure at small negative velocities is due to summation of a close-toregular bichromatic force profile obtained via the X 2 Σ + ↔ A 2 Π 1/2 transition, and the offcenter resonant (Doppleron) peak always appearing at k(v − v 0 ) = ±δ/3 in all bichromatic force profiles, and here created via the X 2 Σ + ↔ B 2 Σ + transition. The latter transition's profile is centered around v 0 = δ BΣ /4k BΣ ≈ 57.1 m/s, so we would expect its off-center peak to appear at v = v 0 − δ BΣ /3k BΣ ≈ −12.5 m/s, which is where we observe the dip.

VI. CONCLUSIONS AND FUTURE PROSPECTS
We have presented a novel experimentally viable method for achieving large optical molasses-like cooling forces for molecules using polychromatic optical fields driving coherent dynamics in a four-level system. Using direct numerical solutions of the time-dependent density matrix as well as Monte Carlo simulations of the cooling dynamics, we provide evidence that achieving rapid damping of a wide velocity capture range towards zero velocity should be possible for diatomic and polyatomic molecules with various constituents and geometries. Proposed Suppressed Emission Rate (SupER) molasses method relies on spontaneous emission coupling between two coherently-driven two-level systems and should be realizable with many complex nonlinear molecules for which scattering ∼ 100 − 1, 000 photons has been previously proposed [14,56] or already experimentally demonstrated [16,24].
We anticipate that large velocity damping coefficients together with a broad velocity capture range will enable extension of laser-based cooling and coherent quantum control to novel The average fraction of time an atom 4 spends in the excited state on the correct cycle can be associated with an average fraction of particles in the wrong cycle. In the π-pulse model, this fraction can be obtained from the optical pulse shape and interval between consecutive beatnote pulses, and it is also what determines the average excited state population in an ensemble, as well as the photon scattering rate. Assuming that the fraction of time the atom spends in the excited state in the correct cycle is ε and the natural decay rate of the excited state is Γ, the scattering rate for these atoms is simply εΓ, which is to say that on average ε fraction of them undergoes a decay. On the opposite cycle, the atom spends 1 − ε fraction of time in the excited state, so then the scattering rate is (1 − ε)Γ.
The photon scattering process is a random process described by a Poisson distribution.
Therefore, the process of changing cycles is a random process as well, and the waiting time between events can be modelled as an exponential distribution with rate λ 1 = εΓ on the correct cycle, and λ 2 = (1 − ε)Γ on the opposite one. Hence, we consider a CTMC with two different transition rates -the correct-cycle state C transitions to the wrong-cycle state W with rate λ 1 , while the wrong-cycle state W transitions to state C with rate λ 2 . Graph shown in Fig. 20 depicts this configuration: FIG. 20: Simple schematic of BCF π-pulse model. The atom moves from state C to W with rate λ 1 and from W to C with rate λ 2 .
In this picture, to obtain both mean force and variance of momentum transfer we need to know probability P(t) of occupying one of the two states as a function of time. We can find it by using the Kolmogorov forward equation P (t) = P(t)Q, where Q is the generator matrix. In such CTMC the generator matrix is simply [70]: Solution to this equation is: By finding eigenvalues and eigenvectors of the generator matrix, one can calculate the ex-ponent of the Qt matrix. In the end, by using P(0) = I we obtain: Having a certain initial state α = [α 1 , α 2 ] with α 1 + α 2 = 1, the probability of occupying state C or W is simply αP(t). However, the steady-state probabilities at t 1 λ 1 +λ 2 are independent of the initial state: We can now evaluate the expected momentum transfer in time T . Before a spontaneous emission occurs, momentum of 2 k is exchanged between the light field and the atom with approximate constant frequency of δ/π, where δ is the detuning of one of the laser components of the 2-color force. Therefore, if we define random variables: where 1 C/W is an indicator variable for atom's state 5 X being C/W at time t , and random variable Θ(t) is state's so-called occupancy time, we can find the expected value of momentum exchange in a given state given the initial condition α: The last part of that formula is the probability of occupying state X at time t given initial state α. For example, for α = [1, 0] and X = C: (αP) C (t) = P 11 (t) = λ 2 λ 1 + λ 2 + λ 1 λ 1 + λ 2 e −(λ 1 +λ 2 )t . 5 The generic label X to indicate the cycle in which atom resides during the dynamics under the influence of the coherent stimulated forces should not be confused with the ground electronic state for molecules, which is also customarily denoted as X.
In our model, the atom moves in one direction in state C and in the opposite direction in state W . Therefore, the expected total momentum transfer at time T is: The above equation can also be thought of as part of time-averaged expected value of the force: Integration leads to: which for T (λ 1 + λ 2 ) −1 simplifies to: Average force is then: and it also simplifies to: where the last term can be thought of as proportion of time the atom or molecules spends in state C minus proportion of time it spends in state W .
Notice that the result above is independent of the initial condition we have chosen, validating our initial assumption. That independence, of course, is related to the independence of the stationary state of CTMC on the initial conditions. The result for large times T can be also obtained simply from P(t → ∞) -one can show that a time-averaged function of the states (here, state occupancy) is simply the expected value of the function with respect to the stationary distribution. Other way of deriving the average force in this model, would be to calculate average reward (force) per cycle.
To estimate variance of the momentum transfer we first note a few general things about the occupancy time random variables. First, it should be obvious that Θ C (T ) + Θ W (T ) = T , and thus: We also see that: = Var Θ C (T ) + Var(Θ W (T )) + 2Cov (Θ C (T ), Θ W (T )) = 2Var Θ C (T ) + 2Cov (Θ C (T ), Θ W (T )), which shows us that Cov (Θ C (T ), Θ W (T )) = −Var Θ C (T ). Because the momentum transfer depends on the difference between occupancy times, we need to find: We therefore need to only find the variance in occupancy time of one of the CTMC states.
To calculate variance in momentum transfer for a state up to time T we start with writing the definition of variance: where the last term was part of the previous calculation. We concentrate on the first term: The expectation value of indicator random variables has to be calculated with care. We can write it in the following way assuming s < t: Given previously chosen initial conditions and remembering that Markovian process is memoryless: After multiplication of above terms and re-defining s ≡ min (s, t) and t ≡ max (s, t), we obtain (for α = [1, 0]): Integrating Eq.(A6) and using previously found expectation value of momentum transfer, one finds variance of the occupancy time: The above result has a constant term, term linear in time, exponentially decaying terms and a mixed term. In general, it can be found that variance of a reward in CTMC can only have specific terms 6 [71], and terms in Eq.(A7) fall into that category. For large times T → ∞ only one term survives, and so the overall variance in momentum transfer can be found to be: This also allows us to find the diffusion coefficient D: Like before, we can define a new variable, σ 2 II (ε) = 16 ε(1 − ε) and use the average force we found previously in Eq.(A5) to obtain: There are several interesting aspects we can notice about the expression for calculated variance in momentum transfer presented in Eq.(A8). Firstly, it grows linearly with time, which is the same as for the radiative force. Secondly, it follows ∝ δ 2 proportionality, which for polychromatic forces having detunings of the order of 100 Γ means that the variance will be quite substantial. However, if we consider ratio of the standard deviation to the mean: we realize that the distribution becomes narrower the longer the process. In Fig. 21 we show histograms of simulated polychromatic forces obtained from our model. The distribution of the force follows the distribution of the occupancy times. Occupancy time in single state in our model is a sum of independent exponential random variables and is therefore distributed with an Erlang distribution (special case of Gamma distribution), with mean and variance determined by the CTMC. Force, being proportional to difference in occupancy times, is distributed as difference of two Gamma distributions and closed form of its momentgenerating function can be found [72]. Fortunately, for longer interaction times, due to the Central Limit Theorem, the distributions approach a Gaussian distribution, which we included in our figure. Its mean is given by Eq.(A4), while its variance is related to Eq.(A8) -if σ 2 p is variance in the momentum transfer distribution, σ 2 f = σ 2 p /T 2 will be the variance in the force distribution.
The value of the aforementioned ratio grows for ε approaching 1/2, but it, as well as the variance itself, can be brought arbitrarily close to zero by making ε small. For instance, adding additional colors reduces value of this parameter. At optimum force, ε ≈ χ/π, where χ is phase difference between counter-propagating beams in the polychromatic force (refer to Fig. 1). For example, while for a 2-color force χ = π/4, for a 4-color force χ ≈ π/6, and so ε ≈ 0.167. Effects of decreasing ε can be seen in Fig. 21.  It is also worth noting that the velocity diffusion (known as beam "pluming") that could be observed when performing slowing or deflection using bichromatic forces should be pretty small comparing to the overall effect observed due to what was just mentioned.
Finally, for all polychromatic forces, the Rabi rate Ω is of the order of the detuning δ, so the variance Var p ∝ Ω 2 . Already Cohen-Tannoudji divided optical forces into two categories depending on their origin: dissipative and reactive [73]. Polychromatic forces are reactive according to that definition, just like dipole forces. He showed that for such forces the momentum dissipation tensor D, which is associated with variance in momentum transfer and was shown in Eq.(A10), should scale as square of the Rabi rate, which is consistent with our result. It is also consistent with value obtained in Refs. [39,51].
Additionally, we can connect the average excited state population in an ensemble ρ ee with the average time ε a single atom spends in the excited state, which simultaneously is the fraction of atoms currently in the wrong cycle. These can be tied together in a very simple way: ε fraction of atoms spends 1 − ε time on average in the excited state, while 1 − ε of them spends ε fraction of time in the excited state. We can then write: Because 0 ≤ ε ≤ 0.5 we obtain: Finally, we can re-write the expected value of force and time-averaged variance of momentum transfer in terms of ρ ee :

Appendix B: Variance Estimation in PCF Molasses
To estimate variance in momentum transfer and, from there, the limiting temperature in SupER molasses, we will use the CTMC model introduced in App. A. In all generality, the diagram of our system is depicted in Fig. 22. We consider 4 states: C 1 , W 1 , C 2 and W 2 .
The first two correspond to states an atom can be in when it is feeling PCF acting on one of two 2-level systems (Fig. 3). Accordingly, states C 2 and W 2 correspond to the other 2-level system. States marked with letter C are states, where, like in the simple PCF model in the previous section, the atom spends most of its time ("correct" cycle). These states in PCF molasses will create force in opposing directions. Similarly, cycles marked with letter W are the ones, where an atom spends less time ("wrong" cycles). In those cycles the momentum is transferred in direction opposite to the direction in their respective C states.
We assign the average time spent by atoms in states C in the excited states |E1 and |E2 (Fig. 3) with respect to the total time spent in one 2-level system as ε 1 , ε 2 and 1 − ε 1 , 1 − ε 2 for average proportion of time spent in states |G1 and |G2 respectively. Having defined these variables we can find rates for all of our states: These rates are set up in a similar fashion as in the simple PCF model -we assume that when an atom is in one of the 2-level systems the situation is just like the model analyzed in App. A. Then, the rate at which the atom leaves the state due to spontaneous emission is just the natural decay rate Γ times the proportion of time it spends in that excited state considering only the 2-level system in which the cycle of excitation and stimulated emission occurs for this atom or molecule, which is simply ε or 1 − ε.
In our 4-state model the decay can always go to all three other states. A spontaneous decay that doesn't cause an atom or molecule to change a 2-level system (so when it switches between states C 1 and W 1 or C 2 and W 2 ) is assumed to have a branching ratio r 1 and r 2 for the first and second two level system respectively. Switching 2-level system occurs with branching ratios of 1 − r 1 and 1 − r 2 . Comparing that to Fig. 3 gives: When such a switch of a 2-level system happens, the decay can lead to either cycle C or W , but the probabilities are not equal. In Fig. 22 probability of entering state C i is labeled q i . However, in our situation they are simply the proportions of time an atom spends in respective states within each 2-level system, so q i = 1 − ε i 7 . We now can write the generator matrix for this CTMC: Such system is, however, difficult to solve analytically. Even finding the expected occupancy time or the stationary state η of this CTMC, which corresponds to a left eigenvector associated with the zero eigenvalue, is very challenging. Fortunately, the symmetrized version of the system simplifies the situation. Before we move forward, we note a few characteristics of the expected value of the force in this more general system that emulates the population dynamics during the SupER molasses cooling process.
If both 2-level systems are kept at optimum, for example with optimal parameters of Ω 1,2 = 3/2 δ 1,2 and |χ| = 45 • for BCF, in the stationary state the proportion of time spent in state C 1 or W 1 with respect to total time spent in the first 2-level system (i.e. in either of those states) is the same as corresponding proportions in the second 2-level system: where we assume that stationary state ratios are not necessarily the average times spent in energy levels of the system (ε i and 1 − ε i ). In the end, the expected time-averaged value of the force for both 2-level systems should be: In the formulas above we already assumed that both 2-level systems would generate opposing forces (due to χ 1 = −χ 2 ). In molasses we'd like F 1 = F 2 , so that there's no net force at zero velocity. Taking that condition and by multiplying both sides by 1 we obtain: we assume no phase coherence and, therefore, non-zero probability of ending in either of the 4 states.
Because the ratios are the same for both 2-level systems, these terms will drop out: We're now left with total proportions of time spent in the first and second 2-level systems.
Atom decays from the first 2-level system to the second with a total rate (1 − r 1 )Γ 1 and from the second back to the first with a rate (1 − r 2 )Γ 2 . This creates its own two-state CTMC, which we have already solved, so: Using the above, we finally arrive at a criterion for detunings δ i that has to be met to properly balance power in an asymmetric system (like in the BaH molecule considered in Sec. V B and V C): Criterion shown in Eq.(B1) can be understood intuitively: detuning, which determines rate of the cycle of spontaneous and stimulated emission, has to be higher, if energy of the scattered photon is smaller or if spontaneous emission causing switching of the 2-level system occurs more often, which is determined be either the decay rate or the branching ratio. Now, we can move forward with simplification of the model to obtain analytical estimates.
As in the main section, we symmetrize the CTMC and assume that all relevant parameters are identical in both 2-level systems: The simplified state graph for the CTMC is shown in Fig. 23. We also used values for branching ratios in BaH: r 1 = r 2 = 2/3. This model was used in Monte Carlo simulations, results for which are shown in Sec. IV.
FIG. 23: Schematic for a symmetrized and simplified π-pulse model PCF molasses.
This system has a generator matrix: For such simplified system we can identify the stationary distribution: As expected, in symmetrized system proportion of time spent in both 2-level systems is the same, that is η C 1 + η W 1 = η C 2 + η W 2 = 1/2. And so the time-averaged expected value of the force created by both 2-level systems is: Just like in the case of a 2-level system (see Eq. A5), we can write the force exerted through one of the subsystems as: with: Calculating variance of the occupancy time, and therefore the momentum transfer, is more challenging. We first note that by using similar tricks as in the previous section, we can show that variance of the momentum transfer, which is proportional to = 4 2 k 2 δ 2 π 2 × 4Var (Θ C 1 (T ) + Θ W 2 (T )) = 16 2 k 2 δ 2 π 2 Var (Θ C 1 (T ) + Θ W 2 (T )) where Var Θ(T ) is variance of any of the occupancy times. To calculate it, we first take a step back and look at general solutions of the Kolmogorov forward equation. Because the generator matrix of CTMC is negative semi-definite, its eigenvalues are non-positive.
The zero eigenvalue is related to the stationary distribution, while others add exponentially decaying parts to the matrix P. In all generality, we can write: where P ij (t) describes probability of being in state j at time t given the system in state i at time 0. Therefore, η j , the component of the stationary distribution for state j, is the same for any initial state i. Here, ν k is the k-th eigenvalue (ν 0 = 0) and u k ij is a function of eigenvectors multiplying the exponential part. For an initial state α and an n-state system we have: where at the end we used the fact that i α i = 1. From the above and the fact for all k > 0 eigenvalues ν k < 0, we easily see that at T → ∞ only the η j term survives. In variance calculations we need the square of the expectation value: [E α Θ j (T )] 2 = η 2 j T 2 − 2η j To obtain variance, we also require the expectation value of the square of the occupancy time. We first find that (for t > s): Plugging in appropriate values for the probabilities and, as before, by re-defining s ≡ min (s, t) and t ≡ max (s, t), we get: When calculating the variance for large T , both the quadratic term as well as the term that depends on the initial conditions will cancel out, leading to: which is a result that is, as expected, linear in time and independent of the initial conditions.
Analogically, one can show that for large T : In our simplified system we can find the eigenvalues: First, we should note that ν 3 eigenvalue becomes 0 at ε = 0. In fact, at ε = 0 the CTMC stops being recurrent and so there's no well-defined stationary distribution. Physically, an atom will be trapped in one of the C states, thus moving continuously in one direction. In such situation we simply obtain a deterministic continuous momentum transfer with zero variance.
In case of non-zero ε we should expect that variance will not behave as for BCF in 2-level system. For small ε atom or molecule will spend a lot of time in one 2-level system, before jumping to the other one, so the variance will be high. Indeed, using Eq.(B4) and Eq.(B5) for our generator matrix Q, we obtain the variance, which can be written as: where the numerator is non-zero at ε = 0. Therefore, as expected, the variance diverges, when ε becomes small, and is smallest at ε = 0.5 reaching value of: For experimentally achievable ε (that is ε 0.1), variance stays very close to the given limiting value (changes by at most a factor of 2). At ε = 0.5, we wouldn't generate any force at either of 2-level systems. For BCF ε = 0.25 and the variance is approximately: 1 T Var p BCF ≈ 20.11 2 k 2 δ 2 π 2 Γ .
Using the found eigenvalues, we can re-write the function µ IV found in Eq.(B2) as: and similarly, like in Eq.(A10), we can write the diffusion coefficient as: with: We can also notice a relationship between forces F II in the 2-level system and F IV in this model. Namely: This shows that we should expect forces acting on the two-level subsystem discussed here to be just re-scaled versions of the normal 2-level polychormatic forces with re-scaling factor of −1/2ν 1 , which is equal to exactly 4/7 for the bichromatic fields.
In a more general case we can simply use Eq.(B4) and Eq.(B5) directly on whatever combination of occupancy times is appropriate in the system. In general we can write: where sum is over all the states in the model. For example, in a more realistic BCF molasses model in BaH described at the beginning of this section with Γ 1 = Γ 2 , k 1 = k 2 and δ 1 = δ 2 , but with ε ≡ ε 1 = ε 2 , we would obtain: with γ ≡ k 2 δ 2 /k 1 δ 1 . Omitting the term preceding the variance of occupancy times, we have a 1 = 1, a 2 = −1, a 3 = −γ and a 4 = γ. Because the forces in asymmetric systems have to be balanced according to Eq.(B1) to create molasses centered at zero velocity, we know that γ = Γ 2 (1 − r 2 )/Γ 1 (1 − r 1 ).
Evaluating Eq.(B4) and Eq.(B5) algorithmically can be done with ease as long as we are able to find eigenvectors and diagonalize generator matrix Q. In general, assuming the eigenvectors of Q are columns in a matrix V and eigenvalues are diagonal elements of D, we have Q = V DV −1 and so: where exponential of eigenvalue matrix simply has exponents of eigenvalues on its diagonal.
After matrix multiplication one obtains values in cells of P as given in Eq.(B3). To easily get values for η j and u k ij we can instead create a matrix U k = V D k V −1 , where D k is defined as: Then, we simply obtain u k ij that we need in Eq.(B4) and Eq.(B5) as the ij-th cell of matrix U k , i.e. U k ij = u k ij . The same matrix gives us η j = U 0 ij , where the equality holds for any index 0 ≤ i ≤ n − 1 in a system with n states. In summary, evaluating variance and covariance of occupancy times, and therefore variance of momentum transfer in π-pulse models for PCF, boils down to finding eigenvalues and eigenvectors of the generator matrix of the appropriate CTMC. Using this method we can numerically find that in a more realistic, asymmetric and Finally, we should note that applicability of this model and all the formulas to the actual polychromatic forces is limited to situations when the interaction can be actually approxi-