Dynamics of a buffer-gas-loaded, deep optical trap for molecules

We describe an approach to optically trapping small, chemically stable molecules at cryogenic temperatures by buffer-gas loading a deep optical dipole trap. The ~10 K trap depth will be produced by a tightly-focused, 1064-nm cavity capable of reaching intensities of hundreds of GW/cm$^2$. Molecules will be directly buffer-gas loaded into the trap using a helium buffer gas at 1.5 K. The very far-off-resonant, quasielectrostatic trapping mechanism is insensitive to a molecule's internal state, energy level structure, and its electric and magnetic dipole moment. Here, we theoretically investigate the trapping and loading dynamics, as well as the heating and loss rates, and conclude that $10^4$-$10^6$ molecules are likely to be trapped. Our trap would open new possibilities in molecular spectroscopy, studies of cold chemical reactions, and precision measurement, amongst other fields of physics.

molecules. Buffer-gas cooling with 4 He in the trap volume is therefore sufficient to load molecules into the trap. After equilibration, the buffer gas is pumped out of the chamber, leaving a trapped sample in the laser beam. These methods rely only on a molecule's DC polarizability, which is nonzero for any species, and do not require a particular energy level structure or magnetic or polar molecules. This not only allows a large number of species that cannot be trapped by existing methods to be loaded into our dipole trap, but also for two or more different species to be trapped at the same time.
Although our methods should be applicable to molecules of any symmetry, we will limit ourselves to discussing linear molecules for simplicity. Table I summarizes our results. Molecules with typical mean DC polarizabilities (α s 2Å 3 , averaged over all orientations), high ionization energies (I 0 ∼ 12 eV), few atoms (≤ 3), and low DC polarizability anisotropies (∆α s 1Å 3 ) are good candidates for our trap, but other molecules could be trapped with small modifications to the trap design. With so many molecules to consider, for clarity we will frequently refer to a hypothetical molecule Q as a representative of molecules we wish to trap. Q is a small, chemically stable (SCS) molecule with characteristics similar to many of the molecules in Table I. We choose the polarizability and ionization energy of Q to be α s = 2Å 3 and I 0 = 12 eV, respectively, based on the species listed in Table I.
Furthermore, we assume Q's polarizability to be isotropic (i.e., zero ∆α s ), and take Q to have a boiling point of T B = 200 K and a molecular mass of m Q = 30 u. These are typical values for SCS molecules, and we note that the trapping is not sensitive to the exact values.

C. Outline
The experimental design, including the cavity and the buffer-gas cells, is discussed in Section II. Section III treats the very far-off-resonant, quasielectrostatic dipole trapping of molecules. Section IV considers the effects of the high-intensity light on the molecules in the trap. The dynamics of buffer-gas cooling and loading into the trap are the subject of Section V. Finally, the ionization-based detection of molecules is reviewed in Section VI, and a summary and outlook are given in Section VII. formed by 1064-nm light (red beam) inside a 100-mm-long, near-concentric build-up cavity, whose mirrors are held at room temperature. Helium buffer gas (He; shown in green) at 1.5 K and molecules (black and purple) at a higher temperature flow through fill lines into two opposing cryogenic cells, held at 1.5 K. The molecules thermalize with the buffer gas and flow, along with the buffer gas, through conical apertures into the loading region. Here, the molecules are loaded through buffer-gas collisions into, and subsequently trapped in, the distinct lattice sites formed by the antinodes of the standing wave of the cavity (inset on top right). Shutters in front of the cells can quickly interrupt the flow out of the cells to create an isolated sample of trapped molecules. Resonance-enhanced multiphoton ionization (REMPI) with a second laser beam (blue beam) ionizes the molecules, which are then detected with a microchannel plate (MCP). Also shown are radiation shields at 50 K and 4 K (with attached charcoal sorbs used for cryopumping of He), and non-reflective (NR) material to absorb scattered light. The drawing is to scale.

II. Experimental design
A. Near-concentric build-up cavity for 1064-nm light The high-intensity, 1064-nm trap light is produced by a near-concentric build-up cavity, modeled on a demonstrated cavity designed by our group [50]. The demonstrated cavity has a finesse of F ≈ 37000 and a power enhancement factor P circ /P in ≈ 9000 (accounting for coupling inefficiencies and technical losses). By coupling input powers P in ≈ 15 W into the fundamental TEM 00 mode of the cavity, circulating powers of P circ ≥ 125 kW have been achieved. With the 20-mm-long, symmetric cavity operating 10 µm from concentricity, the All molecules listed are linear symmetric tops. Entries are calculated assuming a laser intensity of I = 300 GW/cm 2 at a wavelength of λ = 1064 nm. The mean value α s and the anisotropy ∆α s of the DC polarizability determine the trap depth T trap , with ∆T trap being the contribution to T trap from molecular alignment by the intense trap light (see Appendix A). The Rayleigh and rotational Raman scattering rates, W R and W RR , are proportional to the square of α s and ∆α s , respectively (see Section IV A and Section IV B). The ionization energy I 0 is used to estimate the ionization rate W i (see Section IV C). Additionally, the mass m and boiling point T B are given. Resonance-enhanced multiphoton ionization (REMPI) schemes for detection of each molecule are shown (see Section VI), where n + l refers to a scheme absorbing a total of n + l photons at the given wavelength(s). Entries and relevant spectra are taken from a variety of public databases and individual publications [53][54][55][56][57][58].  30 200 - † Atomic Xe is included as its properties make it a good test species for experimental designs. ‡ Included for reference.
a Also: 2+1, 202 nm [71]. b Photodissociation of N 2 O into N 2 and O, followed by 2+1 REMPI of N 2 , using the same laser beam at 204 nm. mode has been focused to a waist of w 0 ≈ 8.7 µm (1/e 2 intensity radius), resulting in an intensity ≥ 400 GW/cm 2 at the antinodes of the cavity's standing intensity wave [50].
Compared to [50], we require a five-fold increase in the cavity length to 100 mm, which creates space for the cryogenic system surrounding the cavity focus. To achieve the same mode waist, the longer cavity needs to be aligned closer to concentricity, which makes it more sensitive to misalignments. In particular, thermal mirror deformation from laser-induced, local heating has a tendency to increase the mode waist and hence decrease the intensity in the cavity. While a detailed analysis is beyond the scope of this work, using mirrors made from ultra-low expansion glass (ULE), coated with an ultra-low-absorption reflective coating (with an absorption as low as 0.4 ppm at 1064 nm (unpublished measurement by our group)), will allow us to achieve an intensity of I = 300 GW/cm 2 with a mode waist of w 0 = 8 µm, based on extrapolating the data of [50,72].
Our cavity mirrors will be kept outside the cryogenic system to prevent buildup of frozen molecules on their surfaces. The high intensity and high finesse of the cavity limits the use of optical elements like windows in the cavity's optical path. Therefore, as shown in Fig. 1, our cryogenic system will have apertures that let in the cavity light. To prevent diffraction losses, the diameter of these apertures is here chosen to be four times the 1/e 2 intensity beam diameter.
B. Buffer-gas cells A 1.5-K buffer gas is used to both cool and load the molecules into the trap. Many aspects of the design of buffer-gas cells benefit from the extensive work done in other experiments [26,27,[73][74][75]. However, as opposed to creating a buffer-gas beam [74,75] or loading our trap within a buffer-gas cell [26,27], we here opt to populate a millimeter-scale loading region between two cells, centered on the dipole trap, with buffer gas and cold molecules.
This minimizes the amount of gas pumpout required and allows for optical access. Furthermore, since our experiment is critically dependent on highly reflective cavity mirrors whose sensitivity to ablation byproducts is unknown, we will not source molecules in the cells using laser ablation, and instead will flow molecules into the buffer-gas cells through heated fill lines. The resulting dual-buffer-gas-cell geometry is depicted in Fig. 1, and uses two closely spaced, cylindrically symmetric, opposing 1.5-K cells of ∼30-mm dimensions and conical apertures with a diameter of 5 mm. This design was validated to achieve the required densities and sufficient thermalization using the direct simulation Monte Carlo (DSMC) method [76]. Vacuum is maintained outside the loading region by differential pumping with activated charcoal sorbs on the 4-K shield at 6 L/s/cm 2 [77]. Pumping out the loading region is done by shutting off flow from the buffer-gas cells using rapidly actuating, cryogenic shutters.
Based on the demonstrated shutter of [25], we assume a shutter can actuate in 1 ms, which sets the pumpout timescale from the loading region in our experiment.

C. Heat load on cryogenic system
The cells will be maintained at 1.5 K by thermal contact to a 4 He-filled 1-K pot, which is pre-cooled and radiation-shielded by 4-K and 50-K stages of a pulse-tube cryocooler.
Commercially available cryocoolers are able to provide 100 mW of cooling power at 1.5 K, which sets the heat-load budget for the experiment.
The conductive heat load from the heated molecule fill lines has been shown to be manageable by previous experiments, including one using a fill line for water held at 280 K connected to a similar buffer-gas cell [78]. Convective heat loads will be managed by differential pumping with activated charcoal. Radiative heat loads on the buffer-gas cells through the apertures in the radiation shields are estimated to be less than 10 mW.
Another major heat load is scattering of the high-intensity cavity light off the mirrors through the apertures onto the buffer-gas cells. Polishing the cells to enhance their reflectivity, and carefully designing the radiation shields, including adding non-reflective (NR) material to the 50-K shields (see Fig. 1), manages this heat load (see Appendix B).

A. Trap depth
Most small, chemically stable (SCS) molecules have limited activity in the optical and near-infrared (NIR) spectrum. Therefore, a dipole trap using light at a wavelength of λ = 1064 nm is red-detuned by several harmonics from the first electronic transition of a typical SCS molecule. In this regime, the commonly used rotating wave approximation does not apply, and the light produces a quasielectrostatic trap by creating a dipole potential [79][80][81][82] where α s is the mean DC polarizability [83] and · denotes a time average over the optical period 2π/ω = λ/c (ω: optical angular frequency, c: speed of light). The oscillating electric field is given by E = E 0 cos(ωt)Ê, withÊ the normalized polarization vector and E 0 the electric field amplitude. The intensity is I = 0 c |E| 2 = 0 cE 2 0 /2 ( 0 : permittivity of free space). From Eq. (1) the intensity required for a 10 K trap depth depends simply on α s as The trap depth of several K is comparable to the rotational energy level spacing of many SCS molecules. The trap may therefore significantly hybridize the rotational levels and align the molecules' most polarizable axis with the optical polarization, leading to an increase in the trap depth that is not accounted for in Eq. (1) [49,84,85]. The trap depth Boltzmann's constant) as listed in Table I therefore includes a correction to Eq. (1) of ∆T trap , which depends on the polarizability anisotropy ∆α s and is discussed in Appendix A. This correction is usually only a few percent of T trap , but for highly anisotropic molecules like CS 2 it can be substantial.
As seen in Table I, most molecules have T trap 10 K at our intensity I = 300 GW/cm 2 , more than six times the buffer gas temperature. Species with particularly low polarizabilities, like H and H 2 , may still be trapped deeply if much higher intensities can be generated.
Helium buffer gas atoms are also attracted to the trap center, but with a smaller trap depth of T trap = 0.95 K due to their low polarizability. We ignore this small effect in the rest of this publication.
Molecules typically have a rich spectrum of pure rovibrational transitions. Although the trap light is blue-detuned from all pure rovibrational transitions, these transitions do not contribute significantly to the dynamic polarizability at 1064 nm [86][87][88], and hence do not lead to antitrapping. Typical SCS molecules are in their electronic ground state and predominantly in their rovibrational ground state at our chosen buffer gas temperature of T = 1.5 K.

B. Optical potential
The trap itself is characterised firstly by the trap depth T trap = max |U |/k B . Since our trap light is the Gaussian TEM 00 mode of an optical cavity, the spatial dependence of the trapping potential, in cylindrical coordinates measured relative to the focus of the cavity (placed at the origin), is where the approximation is valid near the focus of the trap. Here, k = 2π/λ, w(z) = Appreciable molecule loading will only occur in the region where the trap potential is large compared to the buffer gas temperature. The trap volume is therefore characterized approximately by the volume within which U (r, z) ≥ 3k B T , which is sensitive to the dimensionless trap depth parameter η 0 = T trap /T : We note that V thus scales as w 4 0 . This volume is split non-uniformly into 4 η 0 3 − 1z R /λ distinct lattice sites, spaced axially by λ/2.
We propose Q be loaded from a buffer gas at 1.5 K into our 1064-nm dipole trap with a peak intensity of 300 GW/cm 2 and an 8 µm waist. The resulting trap depth is 9.4 K, and η 0 = 6.3. The corresponding approximate trap volume is V = 1.8w 2 0 z R = 2.2 × 10 −8 cm −3 , and there are 743 distinct lattice sites with a site depth greater than 3k B T .

IV. Effects of high-intensity light on trapped molecules
In this section, we review Rayleigh and Raman scattering, photoionization, and photodissociation of molecules, and determine the class of molecules which is unaffected by these processes.

A. Rayleigh scattering
In a quasielectrostatic trap, photon scattering is dominated by Rayleigh scattering [79], with a cross section of σ R = 8π 3 α 2 s /(3 2 0 λ 4 ) [89]. Even at the high intensity of I = 300 GW/cm 2 , the Rayleigh scattering rate W R = Iσ R / ω is only ∼ 10 3 s −1 for the molecules in Table I. The maximum amount by which the Rayleigh scattering can be enhanced by the cavity is given by the Purcell factor F P = 6λ 2 F/(π 3 w 2 0 ) [90], which is 119 for our cavity. Taking this into account, the heating rate from photon scattering is <0.1 K/s for all molecules in Table I [82], which is negligible in our experiment. Similarly, the low scattering rate means the maximum density limit from reabsorption of scattered light is > 10 35 cm −3 , which can be ignored [91]. In reality, F P = 119 overestimates the cavity enhancement of the Rayleigh scattering, as effects such as the recoil shift and collisional broadening will shift the scattered photons off resonance with the cavity. The molecule absorbs a cavity photon and emits a photon of a different energy, the energy difference accounting for a change in rotational state of the molecule. In linear molecules with zero electronic angular momentum about the symmetry axis, the Raman selection rule is ∆J = 0, ±2, where J is the rotational angular momentum of the molecule [92]. RRS from the J = 0 ground state of molecules occurs at a rate [57,93] where λ ≈ λ is the wavelength of the scattered photon, ω = 2πc/λ is the incident angular frequency, and ρ = 3/4 is the depolarization ratio of the Raman transition [94] ( = h/2π: reduced Planck constant). In rare cases, e.g., for NO, the molecular ground state has nonzero electronic angular momentum about the symmetry axis, in which case Eq. (4) is modified by different selection rules and Placzek-Teller coefficients [95].
Molecules with W RR 70 s −1 , such as N 2 , CO, O 2 and HCl, are ideal first candidates for our proposed experiment, as they can be isolated from the buffer gas, and potentially evaporatively cooled, faster than they are rotationally heated (see Section V B 5). The RRS cross section scales with ∆α 2 s and λ −4 , so molecules with higher polarizability anisotropies can still be trapped and isolated as described in this work using a trap with λ > 1064 nm (see Table I).

C. Ionization
The molecules suitable for trapping in Table I  meaning non-resonance-enhanced ionization occurs via multiphoton ionization (MPI) [99].
The ionization rate W i scales as W i = σ n I n , where n = 11 is the number of photons needed to reach the continuum and σ n is the MPI cross section.
It is difficult to estimate σ n ab initio. As a starting point, we use Popruzhenko's formula [100] for hydrogen-like atoms, which is a Perelomov-Popov-Terent'ev (PPT) ionization formula [101]. The formula matches (within two orders of magnitude) experimental results for ionization rates of noble gases [102,103] and air [104]. PPT formulas also match ab initio estimates of ionization rates for polyatomic molecules [105]. Popruzhenko's formula predicts W i = 2.8 × 10 −9 s −1 for Q. W i scales down superexponentially with increasing I 0 , meaning molecules with slightly higher ionization energies have substantially smaller ionization rates, as seen in Table I. The ionization rate can be minimized by either using a molecule with a high ionization energy, or a molecule with a high polarizability which requires less laser intensity for trapping. To compare these approaches, we compute the ionization rate at an intensity for which the trap depth is 10 K, W i,10 K = W i (I| Ttrap=10 K ). This is plotted for different values of the mean DC polarizability α s and ionization energy I 0 in Fig. 2, ignoring the correction to the  Table I (except CS 2 , which is outside the shown region) are marked with red crosses on the plot. Evidently, W i,10 K is much more sensitive to the ionization energy than the polarizability. trap depth ∆T trap discussed in Appendix A. The near-vertical contours roughly represent lines where the ionization energy corresponds to an integer multiple of the photon energy, so the absorption of an extra photon becomes necessary to ionize the molecule, thereby greatly decreasing the ionization rate. By comparison, the higher intensity necessary for lower polarizabilities has a much smaller influence on the ionization rate. Thus, ionization energy is a far more important quantity than polarizability in the selection of suitable molecules, granted that a sufficiently high laser intensity can be achieved.
An alternate way to reduce ionization rates is to trap molecules at a longer wavelength.
The ionization rate depends primarily on the number of photons needed to reach the ionization continuum, so molecules with low ionization energies may be suitable for trapping with light at λ > 1064 nm.
Popruzhenko's formula provides a first estimate of the ionization rate, but can be inaccurate in general. There are many examples where multiphoton ionization occurs through an intermediate multiphoton resonance, which can enhance the ionization rate by several orders of magnitude when compared to Popruzhenko's formula [106][107][108]. Resonance-enhanced multiphoton ionization (REMPI) is now a routine experimental technique for ionization of different molecules, as will be discussed in Section VI. It is challenging to estimate REMPI cross sections a priori, particularly due to the lack of detailed spectra of general molecules that are free of spectral broadening. For this reason, we seek molecules with an ab initio ionization rate estimate of 10 −9 s −1 (I 0 12 eV) in Popruzhenko's formula as first candidates for our trap, so that even if we ignore the contribution of electronic resonances to the ionization rate, we are still unlikely to see ionization in the trap. We note that limits on the ionization rate can be placed with a room-temperature experiment using our trap light, which is especially helpful for determining the suitability of molecules with higher ab initio ionization rates such as NO or CS 2 (see Table I).
In a high-density gas, laser-induced breakdown becomes the dominant mechanism for ionization in the beam. We have studied breakdown in high-intensity, CW laser beams in detail, and the results are shown in Appendix G. We have determined that the buffer gas is too low-density for breakdown to be a concern during trap loading.

D. Dissociation
Photodissociation of molecules can occur by two distinctively different mechanisms depending on the character of the light that drives it. The first mechanism is direct electronic excitation to a dissociative state. In principle, any molecular state located above any bond's dissociation threshold is dissociative. In practice, though, excited states above the dissociation continuum can be long-lived, as the probability of direct tunneling out of a quasibound excited state into the dissociation continuum is extremely small due to the Franck-Condon principle [92]. Only certain electronic excitations lead to dissociation. In SCS molecules, the dissociative states are either highly excited quasibound states, which have been probed through absorption of single vacuum-ultraviolet photons [109,110], or excited states above the ionization continuum, as has been seen in multiphoton dissociation experiments using NIR and visible pulsed lasers [108,[111][112][113]. In particular, [108] finds that dissociation using NIR lasers scales highly nonlinearly with the peak intensity, as expected for multiphoton absorption. Thus, following the arguments on ionization in Section IV C showing that simultaneous absorption of enough photons to reach these highly excited states is unlikely, and provided Q's electronic transitions are far from resonant with the low harmonics of the 1064-nm light, we expect not to be limited by this form of dissociation.
The second mechanism, infrared multiphoton dissociation (IRMPD), is usually observed in polyatomic molecules or ions using mid-to long-wavelength infrared light, such as from a CO 2 laser at 10 µm [114], resonant with pure rovibrational transitions. IRMPD is a heating mechanism, so the dissociation probability depends on the total energy absorbed rather than peak intensity, provided the light is near a resonance and above a relatively low threshold intensity [115]. As the latter is easily exceeded with our trap light, we must ensure that the IRMPD will not occur at 1064 nm for our molecules of choice.
The three-step mechanism of IRMPD is explained in detail in [116][117][118]. Firstly, light tuned to a pure rovibrational transition of a molecule repeatedly excites a vibrational mode, until the anharmonicity of the internuclear potential shifts the transition frequency off resonance with the light. This step is impossible in homonuclear diatomic molecules, which do not admit dipole-allowed, pure rovibrational transitions. Secondly, once excited to a high vibrational energy, in a polyatomic molecule, the different molecular vibrational modes are split into a number of closely spaced energy levels, and a "quasicontinuum" of energy levels forms. In the quasicontinuum, heating and ergodic mixing of vibrational quanta occur. The quasicontinuum cannot form in any diatomic molecules, and is unlikely to form in triatomic molecules, because there are not enough different vibrational modes to form a closely spaced structure [117,119]. Finally, once the molecule's total vibrational energy is higher than the dissociation energy, the molecule soon dissociates through accumulation of excitations in a particular dissociative bond. The molecules in Table I all have ≤ 3 atoms, and therefore are not susceptible to IRMPD.
For molecules with many atoms, IRMPD may still be suppressed because 1064-nm light is several times more energetic than a single vibrational quantum in any molecule. In the dipole and harmonic approximations, rovibrational selection rules only allow molecules to absorb one vibrational quantum of energy at a time [120], so absorption in high harmonics of rovibrational transitions is unlikely, and the IRMPD mechanism cannot begin. In polyatomic molecules that do not contain H atoms, absorption at optical frequencies appears to be much smaller than Rayleigh scattering [88], so it is likely that the trap light is far enough off-resonant from the fundamental vibrational modes to suppress IRMPD. For molecules containing H atoms, NIR frequencies correspond roughly to the third harmonic of a fundamental vibrational mode, and absorption has been seen in room-temperature experiments [121,122]. However, we are not aware of any corresponding data at cold temperatures.
Thus, for molecules with more than three atoms the overall likelihood of IRMPD in our trap is unclear, and we therefore do not include them as candidate molecules in Table I.

V. Buffer gas dynamics
A. Buffer-gas cooling The first step toward loading the optical trap is cooling the molecules to cryogenic temperatures through buffer-gas cooling [23]. Molecules of all masses and initial temperatures appear to be amenable to buffer-gas cooling [23,123]. As seen in Table I, we focus on molecules with a boiling point T B around 300 K or lower, so they can be loaded into the cold buffer gas in gas phase through a heated fill line without significant thermal load on the cryogenic system. In the design shown in Fig. 1, buffer-gas cooling occurs in two opposing 1.5-K cells to thermalize the molecules to T = 1.5 K before buffer-gas loading the trap.
This section focuses only on the cooling in the cells, and the loading dynamics are left to Section V B.
A cold gas of helium at temperature T = 1.5 K is pumped into cells, also at T , at a moderate to low density which we set to be n He = 10 15 cm −3 . Q is pumped into the cells through heated fill lines to make up a small fraction of the number density, which we set to be 1/100. Through collisions with He atoms, Q comes from a warm temperature T i to T , however it does not solidify as long as it does not collide with the cell walls. To produce a cold gas of Q, we therefore need the mean free path l of Q in the buffer gas to be short compared to the cell dimensions. l is given by [124] where ν is the Q-He collision frequency, n He is the buffer gas number density, σ the Q-He collision cross section, and v Q and v rel the Q and mean Q-He relative velocities, respectively.
The last equality, labeled with "eq", holds in thermal equilibrium.
Based on [125][126][127], we estimate that SCS molecules have collision cross sections with He of order 10 −14 cm 2 at T = 1.5 K, and we take the Q-He collision cross section to therefore be 10 −14 cm 2 . With n He = 10 15 cm −3 , l is 0.34 mm at 1.5 K.
We simulate the thermalization dynamics, using the methods described in the Appendix of [76]. Q molecules are drawn by rejection sampling from an effusive initial speed distribution at T i , directed out of a heated fill line at the origin in theẑ direction. Each molecule of the molecule at that point. Trajectories continue until the effective temperature reaches the 1.5-K buffer gas temperature. The green ellipse has semi-minor and -major axes equal to the x and z standard deviations respectively, showing the 1-σ spread of the positions of molecules when thermalization occurs. (b) Mean z-distance traveled, z, by 1000 simulated Q molecules into the cold buffer gas as a function of their initial temperature T i and mass m Q . Not shown are the x and z standard deviations, which are both less than 10 mm at all T i , m Q points. The coordinates corresponding to all molecules of Table I (except H 2 , which is outside the shown region) are marked with a red x.
undergoes billiard-ball collisions with He atoms in the buffer gas until it cools to the buffer gas temperature T .
A 2-D projection of 50 molecular trajectories from a fully 3-D simulation is shown in Fig. 3 (a). Each point marks a collision, and the color axis represents the effective temper- The green ellipse shows the 1-σ spread of positions where the molecules come to the buffer gas temperature of 1.5 K. The z coordinates of the centers of the green ellipses z for different values of m Q and T i are shown in Fig. 3 (b). Based on this figure, many SCS molecules will thermalize in a 30-mm-tall cell.

B. Buffer-gas loading
With the molecules thermalized to 1.5 K, we next discuss loading them into the optical trap. Since the dipole force is conservative, some supplementary dissipation in the trap volume is required for trapping [1]. Since here Q-He collisions are the only universal source of dissipation, buffer-gas cooling must occur inside the trap volume to achieve loading.

Loss-free, equilibrium trapped molecule number
The energy-level splittings associated with the trap frequencies are three orders of magnitude smaller than the buffer gas temperature, so the loading dynamics are semiclassical.
Molecules arrive at the loading region after passing through one of the buffer-gas cells, so are already thermalized to 1.5 K. In our closely-spaced, opposing cell geometry (Fig. 1), the molecule and He number densities in the loading region are approximately the same as in the cells during loading, with n Q,LR = 10 13 cm −3 , and n He = 10 15 cm −3 , respectively. If trap losses are negligible, the number density of trapped molecules n(x) during trap loading must eventually follow a Boltzmann distribution [128]: where for the truncation of the distribution at the trap depth.
We now integrate Eq. (6) over space to determine the trapped molecule number. We will index distinct lattice sites with j, refer to their positions as z(j) = jλ/2, and describe the trap depth of each lattice site by the site depth parameter η(j) = η 0 w 2 0 /w 2 (z(j)) (η(0) = η 0 = T trap /T ). By integrating Eq. (6) one site at a time, we can determine the total number of trapped molecules N Q (j) in each lattice site, as well as the total number of trapped molecules N Q = j N Q (j). The sum to j = ±∞ diverges logarithmically, so we truncate the sum beyond sites where η(j) = η 0 w 2 0 /w 2 (z(j)) ≤ 3. The results are shown in Fig. 4, with N (j) plotted on the y axis and the sum N shown in the legend. The Boltzmann distribution predicts that several millions of molecules will be trapped for each species shown, at peak trapped densities of order 10 15 -10 16 cm −3 .
Eq. (6) does not, however, indicate the timescale of equilibration. This makes it hard to determine the effect of losses on the number of trapped molecules per site, and it fails to estimate the number of trapped molecules remaining in the trap after the buffer gas and untrapped molecules are pumped out (n Q,LR → 0). We therefore consider a microscopic approach to trap loading.  6), assuming no losses. Solid curves indicate molecules trapped in sites with site depth parameter η(j) = η 0 w 2 0 /w 2 (z(j)) > 3, dashed curves indicate molecules bound in the trap in weak lattice sites with η(j) ≤ 3. The legend gives the total number of trapped molecules N , given by N = j N (j), with the sum over all j such that η(j) > 3.

Microscopic, ergodic loading model
In our trap, the trap dimensions (set by w 0 = 8 µm and λ = 1064 nm) are much smaller than the mean free path l = 340 µm, so our loading model will differ from models used in buffer-gas-loaded magnetic trapping experiments, where l is small compared to the trap dimensions [23].
Starting with an empty trap, any free molecule passing through the trap potential can only be loaded if a collision with a He atom reduces its energy to below the local trap depth. Of the free molecules that collide with a He atom in the trap, a substantial fraction f 0 will not lose enough energy in the collision to become bound. Since l w 0 , no other collision will occur in the trap volume for these molecules, so they are lost. The remaining molecules do lose enough energy to initially be bound in the trap, however not all of these molecules will go on to thermalize at the bottom of the trap potential. A fraction f 1 of molecules will be initially bound, but never fall to a total energy consistent with a truncated Boltzmann distribution that would indicate thermalization. Meanwhile, a fraction f 2 (with f 0 + f 1 + f 2 = 1) will be initially bound and also go on to thermalize. Molecules which thermalize will live in the trap until they eventually are ejected after a mean number of collisions k Q-He with the buffer gas.  3. The model will compute N Q (j) without explicitly determining the local density distribution in the trap, and will therefore not completely account for out-of-equilibrium loss rates.
With these approximations in mind, we can write down an equation for dN Q (j)/dt.
Molecules will be loaded into an empty trap at a rate n Q,LR (t)V e (j)ν(n He (t))f 2 (j), dependent on the number of free molecules in the trap volume n Q,LR (t)V e (j), the Q-He collision frequency ν(n He (t)), and the fraction of collisions in V e (j) that lead to loading f 2 (j). Molecules will be collisionally ejected from the trap at a rate [ν(n He (t))/k Q-He (j)]N Q (j), and two-body collisions between trapped molecules lead to an additional loss term [β Q-Q /V e (j)]N 2 Q (j), where β Q-Q is the Q-Q two-body loss coefficient. Rotational Raman scattering (RRS) leads to an additional loss term Γ RR dependent on the RRS rate W RR . We neglect a non-collisional one-body loss term −ΛN Q (j), because ionization, dissociation, and recoil heating are small (see Section IV). Other sources of trap heating are also small, as discussed in Appendix D. Overall, as opposed to a commonly used V eff = π 2 w 2 0 λ/η 3/2 [129]. This effective volume V eff is only appropriate when the density distribution in the trap is close to a Boltzmann distribution.
Our approximate loading model, on the other hand, is only strictly valid when losses are high (β Q-Q 10 −11 cm 3 /s), in which case the trap density distribution is pinned near a constant value n Q (j) ≈ ν(n He (t))f 2 (j)n Q,LR (t)/β Q-Q . The molecules are therefore evenly distributed over the volume V e (j) in the high-loss limit. In the low-loss limit, on the other hand, the choice of the effective two-body loss volume is not critical.
Eq. (7) can be integrated once β Q-Q and Γ RR are specified. In magnetic traps, two-body loss is usually caused by spin-changing inelastic collisions between trapped particles, and typical values of β Q-Q range from 4 × 10 −11 cm 3 /s [28] to 9 × 10 −13 cm 3 /s and lower [130].
Our trap, however, is insensitive to the molecules' internal state, so we are likely not limited by this kind of loss. In Appendix E 2, we show that collision-induced absorption, seen in O 2 gases, also does not cause significant two-body loss when trapping O 2 at the densities expected in our experiment. In optical trapping experiments on molecules, particularly bialkali molecules, a high two-body loss rate is observed. In this so-called "universal loss", a close to unity fraction of collisions between molecules in the trap are lossy. The mechanism is discussed in Appendix E 1, but in short, we believe universal loss is unlikely in our experiment due to the weak interactions between, and the high excitation energies of, SCS molecules, and consider universal loss only as a possible "worst case scenario". It is more likely that β is small enough to be ignored, and the loss is dominated by the other loss terms in Eq. (7).
The loss due to RRS, Γ RR , has a complicated form, which is detailed in Appendix F.
Nevertheless, the effect of Γ RR is simple. During buffer-gas loading, rotational cooling is efficient [123], and RRS does not lead to substantial loss. Once the buffer gas is removed, rotational heating causes exponential decay of the ground state trapped population at a rate W RR .

Loading simulation results
Integrating Eq. (7) to find approximate values for N Q (t) is done in four stages, and the results are shown in Fig. 6. We start in Fig. 6 (a) with the simplest case of negligible two-body loss, and show the evolution of trapped Q as a function of time for various optical intensities. The pumpout begins at t = 0 ms, when the cryogenic shutter actuates in 1 ms to stop flow into the loading region. We conservatively assume molecules will cryopump to the surface of the shutter, rather than bounce off it, to underestimate the trap loading during the pumpout, and conservatively assume the He buffer gas will take longer to evacuate the loading region to overestimate the trap losses during pumpout. To this end, we set the He and molecule pumpout timescales from the loading region to be 2 ms and 0.5 ms, respectively.
After t = 0 ms, n He and n Q,LR decay exponentially with these pumpout timescales. To model the effect of He film formation on the outside of the buffer-gas cells, we let n He saturate at With increasing intensity, the number of molecules trapped during the loading phase increases, and also the fraction of these molecules retained in the trap after pumpout increases.
Our proposed experiment, with I = 300 GW/cm 2 , appears to trap about 10 5 Q molecules in our approximate model. The inset in Fig. 6 (a) shows the evolution of the number of molecules over time in each site for I = 300 GW/cm 2 . Just before pumpout, at t = 0 ms, the trapped molecule numbers predicted by Eq. (7) (solid curve) resemble a Boltzmann distribution (dotted curve), except that the approximations involved in Eq. (7) lead to an underestimation of the trapped molecule number by a factor of 6. After 50 ms, when the trap is isolated from background-gas collisions, we see that molecules near the focus of the trap are retained more than molecules far from the focus, since the trap depth η(j) is higher.
In Fig. 6 (b), we consider the effect of two-body losses. The details of the pumpout are the same as in Fig. 6 (a), but the integration time is increased so that the long-time behavior is visible. The blue curve, with β Q-Q = 10 −10 cm 3 /s, represents universal loss, while the red curve, with β Q-Q = 0 cm 3 /s, represents no two-body loss. We see that two-body loss leads to a reduction in the trapped molecule number during loading, the fraction of molecules that survive the pumpout, and the number of molecules in the long-time limit. Roughly speaking, in the high-β Q-Q limit, each order-of-magnitude reduction in β Q-Q leads to an order-of-magnitude increase in the trapped molecules 200 ms after pumpout.
The inset of Fig. 6 (b) shows how molecules are distributed among lattice sites in the universal loss limit (β Q-Q = 10 −10 cm 3 /s). Pinning of the trapped molecule number near the trap focus is clearly visible, which retroactively justifies the use of V e as the effective two-body loss volume instead of V eff . With universal loss, only a few molecules per lattice site survive after the buffer gas is pumped out, but because the molecules are spread over many lattice sites, the total number of trapped molecules remains large enough for sensitive detection schemes to detect them (see Section VI).
Third, we consider some real molecules (N 2 , CO, O 2 , HCl) in Fig. 6 (c). We ignore RRS in this figure for clarity, and treat it separately. For each molecule, we recompute the parameters f 2 (η(j)) and k Q-He (η(j)), which change with the molecular mass and η 0 , and then separately integrate the case of no two-body loss (solid curves) and universal loss (dashed curves) to demonstrate a best and worst case scenario for the number of each molecule we can expect to trap. Even our worst case estimates based on universal loss suggest thousands of molecules will be trapped at peak densities of 10 11 cm −3 (at 100 ms, averaged over the center lattice site volume V e (0)), which will be enough to demonstrate trapping. In the case of no universal loss, on the other hand, about 10 4 O 2 and N 2 molecules, 10 5 CO molecules, and Finally, we exemplify the effect of RRS in Fig. 6 (d) by studying its impact on N 2 (see Appendix F). During buffer-gas loading, the curves for the case without RRS (blue) and with RRS (orange) are indistinguishable, owing to efficient rotational cooling by the buffer gas.
After the buffer-gas pumpout, exponential decay of the rotational-ground-state, trapped population is observed at the rate W RR .

Evaporative cooling and other considerations
In our models, for the case when β Q-Q and W RR are small, we have so far ignored the evaporative effect of Q-Q elastic collisions after the buffer-gas pumpout. They will lead to additional Q loss, but also decrease the sample temperature and thus increase η 0 through evaporative cooling, making them fundamentally different in nature to the losses included  (7) for different scenarios. Loading occurs at negative t, with He density n He = 10 15 cm −3 and Q density n Q,LR = 10 13 cm −3 . Pumpout to produce an isolated sample begins at t = 0, after which n He and n Q,LR exponentially decay with a timescale of 2 and 0.5 ms, respectively. n He saturates at 10 11 cm −3 to model He film desorption dynamics. In (a), Q is treated assuming no two-body loss for various optical intensities I. Inset shows the distribution of molecules across lattice sites at t = 0 ms (solid curve) and t = 50 ms (dashed curve). The t = 0 ms curve can be compared to the analytic result from Eq. (6) (dotted curve), which shows that the approximate model Eq. (7) underestimates N Q (j) during loading in the low-loss limit by a factor of 6. In (b), Q is treated with fixed intensity I = 300 GW/cm 2 , but varying two-body loss coefficients β Q-Q . The highest value of β Q-Q = 10 10 cm 3 /s corresponds to the "worst case" of universal loss, which we believe to be unlikely in our trap (Appendix E 1). Inset shows the distribution of molecules across lattice sites at different times for the case of universal loss. In (c), we integrate Eq. (7) for N 2 , CO, O 2 , and HCl with I = 300 GW/cm 2 , ignoring rotational Raman scattering (RRS). For each molecule, we consider the case of no twobody loss (solid curves), as well as the unlikely "worst case" of universal loss (dashed curves). β Q-Q coefficients for universal loss are computed as in Appendix E 1. For the case of no two-body loss, the peak trapped densities at t = 100 ms (averaged over the central lattice site volume V e (0)) are 1.7 × 10 12 cm −3 (O 2 ), 3.1 × 10 12 cm −3 (N 2 ), 1.4 × 10 13 cm −3 (CO), and 2.3 × 10 14 cm −3 (HCl). In (d), we integrate Eq. (7) for N 2 , comparing the cases of including and excluding RRS to exemplify its effect. The effect of RRS is modeled in Appendix F.
in Eq. (7) [133]. We compute the initial evaporation timescale at constant η 0 [134,135] to be N Q (0)/Ṅ Q (0) = τ ev = 15 ms for the central lattice site of our trap, 1.5 times less than [134]. Therefore, for species with W RR < 70 s −1 , at least some amount of evaporation can be achieved initially, allowing for a reduction of the trap depth which in turn proportionally reduces the RRS rates.
For molecules with low losses, we expect similar evaporative cooling dynamics to [134], where 800 optically trapped Rb atoms are evaporated to 40 atoms with η ∼ 5, resulting in a 1000-fold reduction in temperature and increase in phase-space density. At the now reduced optical intensity required to maintain trapping of tens of molecules per lattice site at ∼1.5 mK (phase-space density ∼ 10 −4 ), the sites could be combined using a bichromatic light field [136]. The resulting sample of thousands of molecules at mK temperatures could be evaporated further towards the ultracold regime.
We have also thus far ignored the effect of trapped He atoms. Although a small amount of sympathetic evaporative cooling can be expected from the rapid evaporation of trapped He [137], the number of He atoms that will survive the buffer-gas pumpout is negligible due to their low trap depth of η 0 = 0.6.
One final consideration during He pumpout is the buffer-gas "wind" dragging molecules out of the trap, as observed in buffer-gas-loaded magnetic traps [132]. In our trap, the Q-He collision frequency is slow compared to the trap frequencies, so the Q position and velocity in the trap is randomized between collisions. The wind therefore does not provide a unidirectional drag force on trapped molecules, so need not be considered in our trap.

VI. Detection
To detect the molecules in the trap both during and after loading, a sensitive and background-free scheme is required. Absorption or fluorescence detection techniques, commonly used for cold and ultracold atoms, are challenging to implement for most SCS molecules due to the lack of optical cycling transitions.
Resonance-enhanced multiphoton ionization (REMPI) of molecules with an intense UV laser pulse, in combination with the detection of the charged products, is both sensitive and background-free [138]. The UV laser pulses can be derived from a frequency-converted, tunable dye laser, and a microchannel plate (MCP) can serve as a detector. Due to its resonant nature, REMPI only ionizes a given species, but not any other species in the background gas.
It can also resolve internal states, allowing a measurement of the rovibrational temperature of the molecules. Furthermore, differential AC Stark shifts from the trap light will likely lead to the resonant frequencies for trapped and untrapped molecules to be different on the order of the trap depth of ∼200 GHz (∼10 K). Thus, given a narrow enough resonance, REMPI can distinguish between trapped and untrapped molecules. For example, [139] has used REMPI to resolve different rotational levels in N 2 at 15 K, corresponding to a resolution of better than 120 GHz. This also opens up the intriguing possibility to selectively remove molecules with a certain kinetic energy from the trap, which could be of use in forced evaporative cooling schemes, or to forcibly remove rotationally excited molecules from the trap.
Although REMPI is only applicable to molecules which have selection-rule-allowed multiphoton transitions [140], the technique is widely applicable. In addition to all of the molecules listed in Table I, REMPI has been demonstrated on other symmetric tops like NH 3 [141] and benzene [142], asymmetric tops like SO 2 [143], H 2 O [144], methanol and ethanol [145], radicals like NH [140], SF 2 [146] and OH [147], aromatics and organic compounds [148], and a wide variety of other atomic and molecular species [149].
An alternative scheme is nonresonant ionization of molecules in a tightly focused, ultrashort laser pulse, and subsequent characterization of the products using time-of-flight mass spectrometry (TOFMS) [102,107,108]. This allows for the simultaneous detection of arbitrary molecular species in the trap and thus the monitoring of chemical populations as a function of time during cold chemical reactions. However, trapped and untrapped molecules cannot be easily distinguished with this technique, leading to a large background signal especially during the loading phase. Moreover, the long TOFMS path and ion optics required to distinguish different mass products will make the technique challenging in our experiment. We note that, in principle, an electron beam could be used for non-resonant ionization, but this will lead to a larger background signal compared to a laser beam as gas outside the focal volume will also be ionized.
Finally, the trap light scattered off the molecules could be used for detection. However, even for the high intensities assumed here, the Rayleigh scattering rate is only ∼ 10 3 s −1 (see Section IV A), which will be difficult to distinguish from other sources of scattered light in a realistic apparatus. This could in principle be overcome by using a second resonant, but otherwise empty, cavity to enhance the Rayleigh scattering rate [90].

VII. Summary and outlook
In this work, we have studied trapping small, chemically stable ( heavy molecules and radicals that are of interest in atmospheric and interstellar chemistry [18], and in tests of the Standard Model [15,150].
Moreover, the literature on high-intensity, continuous-wave laser-molecule interactions is sparse, so although we have determined that molecules with ≤ 3 atoms and I 0 ≥ 12 eV are unlikely to be destroyed by the trap light, this does not necessarily imply that other molecules will be destroyed. We believe that many molecules outside these constraints will still be suitable for trapping. In the case of ionization, because of the approximations made in Popruzhenko's formula used here to estimate the ionization rate W i , we here have placed a rather conservative constraint of W i < 10 −9 s −1 . By measuring the ionization rates of molecules in our high-intensity cavity, without even the need for buffer-gas cooling, the robustness of molecules against ionization can be directly determined. As for dissociation, the arguments made above do not obviously rule out molecules like the planar BF 3 (α s = 2.4Å 3 , I 0 = 15.7 eV [53]) or the spherical top SF 6 (α s = 4.5Å as potential trap candidates. Although they have more than three atoms, initial excitation of vibrational modes may be suppressed because the 1064-nm trap light is so far from the fundamental vibrational modes in these molecules. Similarly, although room-temperature data for hydrogen-containing molecules with more than three atoms, like CH 4 (α s = 2.5Å 3 , I 0 = 12.6 eV [53]) and CH 3 OH (α s = 3.2Å 3 , I 0 = 10.8 eV [53]), indicate some absorption at near-infrared frequencies [121,122], it is unclear that this will lead to infrared multiphoton dissociation (IRMPD) at cold temperatures. Deuteration and halogenation of these molecules may also reduce the risk of IRMPD by lowering their fundamental vibrational frequencies. Measurements of molecular ionization and dissociation rates in our cavity provide valuable information about the kinds of molecules we can trap in our proposed experiment, but also fill a gap in the literature surrounding the interaction of high-intensity, continuouswave lasers with molecules.
Nevertheless, even the cold chemical reactions of very simple molecules such as those in Table I are difficult to simulate, and therefore interesting to study [17,18]. At cold and ultracold temperatures, accurate descriptions of chemical reactions require a fully quantum mechanical treatment [17,18,151,152]. Despite this, state-of-the-art calculations, including that of the predissociative state lifetime of cold collisional complexes of diatomic bialkali molecules [153], still rely on semiclassical approximations, which are not valid for SCS molecules (see [152] and Appendix E 1). Our proposed trap is insensitive to most molecular properties and could be loaded with multiple different molecular species simultaneously, allowing for the study of a diverse range of cold chemical reactions. Our trap will therefore provide valuable information about the transition between classical and quantum mechanical descriptions of chemical reactions, and help benchmark new theoretical and numerical techniques to compute the dynamics of cold chemical reactions.
Spectroscopy on cold and ultracold, trapped molecules is another promising application of our trap. For example, radio searches for interstellar organic molecules are partly lim-ited by a lack of available experimental data to compare to observed spectra [19,154]. In many cases, computational models are being used to augment experimental observations of molecular spectra, leading to their more frequent use in molecule searches [155]. Laboratory measurements of cold molecular spectra in our trap will therefore not only directly assist interstellar molecule searches, but also provide valuable data to calibrate computational models and prove their general accuracy. Likewise, astrophysical studies of the variation of the proton-to-electron mass ratio µ from the observation of molecular spectra, such as of CH 3 OH, would benefit from improved laboratory measurements of the relevant transitions [14]. Atomic-or molecular-beam-based precision measurements [15,156,157] could be im-

A. Rotational state hybridization and extra trap depth
Eq. (1) is an approximate expression for the trap depth. Here, two corrections are discussed: firstly, the mean dynamic polarizability α s (λ) at λ = 1064 nm is usually slightly larger than the mean DC polarizability α s used in Eq. (1). A correction of order (λ 1 /λ) 2 ∼ 5%, where λ 1 is the wavelength of first electronic excitation of the given molecule, may be warranted for many of the molecules in Table I [ 79]. The correction is difficult to estimate accurately for most molecules, but is in any case of little consequence to the experiment.
The second correction to Eq. (1) is a result of the hybridization of the rotational states of the molecule, and can have large consequences for the trapping of molecules with large polarizability anisotropies ∆α s . In general, symmetric top molecules have a different polarizability along their symmetry axis (α ) and perpendicular to this axis (α ⊥ : ∆α s = α − α ⊥ , α s = (α +2α ⊥ )/3). In this analysis, these polarizabilities are assumed to be constant within any given rotational band. In a field-free setting, the rotational eigenstates of a molecule will be thermally populated such that there is no molecular alignment. However, optical fields of sufficient intensity will dress these rotational states and align molecules so that their maximally polarizable axis aligns with the optical polarization [49,162]. Linear molecules are described by the Hamiltonian [49,84] Here, θ is the angle between the molecule's symmetry axis and the electric field with amplitude E 0 in a molecule-fixed frame, J is the rotational angular momentum of the molecule, and B is the rotational constant. The time-independent Schrödinger equation takes the form of an oblate spheroidal wave equation [49,84] where z = cos θ, m is the projection of angular momentum onto the optical polarization axis,J is the quantum number which adiabatically turns into J in the limit of zero electric field, and ψJ ,m e imφ and uJ ,m are the eigenfunction and energy of the state with quantum numbersJ, m, respectively, where φ is the azimuthal angle about the optical polarization axis. The solutions for ψJ ,m in Eq. (A2) are the angular oblate spheroidal functions S |m|J as defined in section 21.6.4 of [163]. The eigenvalues of Eq. (A2), and hence the energies uJ ,m , are readily computed using standard libraries (e.g., scipy.special.obl cv of [164]). For ∆αE 2 0 /4B 1, we can write a power series expansion for uJ ,m (see Section 21.8.1 in [163]), from which the trap potential UJ ,m is given by For the special case ofJ = 0 and m = 0, we have where I = 0 cE 2 0 /2 is the optical intensity. In the main text, we use U 0,0 for the trap potential U , and define the trap depth as T trap = max|U 0,0 |/k B , and the extra trap depth resulting from molecular alignment as ∆T trap = T trap − α s I/2 0 ck B .
For some of the molecules in Table I, including those treated in the detail in the main text, ∆T trap is much smaller than 1 K, and the degree of alignment (computed using the Hellmann-Feynman theorem as described in [49]) is small. These molecules are trapped in barely hybridized rotational ground states. However, some of the molecules on the list, such as CO 2 , N 2 O, Cl 2 , and CS 2 , see a substantial increase in the trap depth due to rotational alignment, which scales nonlinearly with the intensity around 300 GW/cm 2 . For these molecules, the character of the rotational ground state is highly aligned with the optical polarization (note that this effect is here ignored in the calculation of rotational Raman scattering rates).
The buffer gas collision frequencies (∼ 10 5 s −1 ) and trap frequencies (<100 MHz) are small compared to the rotational constants of the molecules in Table I (>3 GHz). Thus, trapped molecules adiabatically follow the dressed rotational ground state as they traverse the trap. We note in passing that this is unlike a previously proposed experiment to trap polar molecules in microwave fields [34]. In particular, the polarizability at 1064 nm is independent of the rotational state, so avoided crossings do not open up avenues for rotational state changes during trap traversal in our experiment, as opposed to [34].
In precision spectroscopy of atoms, magic wavelength schemes [158] are used to cancel differential light shifts between two given states by making the dynamic polarizability α s (λ), and thus UJ ,m , of the states equal by choice of a suitable ("magic") trap light wavelength.
To first order in intensity, a magic wavelength can still be found for molecules, although the wavelength must be tuned to not simply cancel the difference in each state's α s (λ) (except in the special case ofJ = 0 [159,160]), but rather, to cancel the first-order term in Eq. (A3) [161]. In general, more magic parameters are required to cancel higher-order terms in Eq. (A3). Rotational state hybridization can represent a large additional contribution ∆T trap to T trap at I = 300 GW/cm 2 (up to 37 % of T trap in the case of CS 2 ) when compared to other corrections [162], such as: nonlinear corrections due to electronic and vibrational contributions to molecular hyperpolarizabilities, which typically contribute a few mK in additional trap depth at I = 300 GW/cm 2 [165][166][167][168][169], and; linear corrections due to higher order multipoles in the multipole expansion [160,161], which typically are orders of magnitude smaller than the leading linear term [170]. Small differences in the polarizability components α s (λ) and ∆α s (λ) within each rotational band [161] additionally modify the analysis starting with Eq. (A1).

B. Heat load from scattered trap light
A major heat load on the cryogenic system in our experiment is scattering of the highintensity cavity light. Our mirrors will be superpolished to ∼1Å RMS (root-mean-square) surface roughness, and the resultant scattered light will be a few hundred mW. From white light interferometry measurements of our previous mirrors' surface profiles, we have determined the angular distribution of scattered light [171][172][173]. About half of the scattered light is diffusely scattered, and is managed by placing the mirrors more than 10 mm from the apertures in the 50-K shields. The other half is scattered into a small cone around the cavity mode. From ray-tracing this scattered light through the near-concentric cavity, we know this light will be incident on either a buffer-gas cell, or the inside of one of the cryocooler's radiation shields, rather than the outside of the radiation shields. The buffergas cells will be polished and highly reflective, so light incident on them will be reflected towards the radiation shields (mostly the 4-K shield). We aim to absorb the scattered light that misses the buffer-gas cells on the 50-K shields rather than the 4-K shields, due to their larger cooling capacity. This is achieved by designing the 50-K shields with a smaller angular size than the 4-K shield apertures when viewed from the opposite mirror, so the scattered light will be absorbed on the non-reflective (NR) material on the 50-K shields shown in An alternative approach to handle the scattered light is to conically indent the radiation shields around the cavity mode, allowing the apertures to be closer to the cavity focus, and therefore smaller, so that almost no scattered light enters the cryogenic system in the first place. This approach has small consequences for the pumpout timescale, and has not been investigated thoroughly.
All other heat loads from scattering of the high-intensity light, including Rayleigh scattering from the buffer gas and trapped molecules (nW) and scattering of transmitted input light not matched into the cavity mode (<1 mW), are negligible for our experiment.

C. Details of ergodic loading simulations
At their core, the ergodic loading simulations work by storing the energy E of Q molecules, and assigning initial conditions for Q-He collisions stochastically based on this energy within the ergodic approximation, before enacting a billiard-ball elastic collision and tracking the change to the energy after the collision, which is again stored, and so on. This is numerically efficient, since no trajectories need to be integrated to simulate the dynamics. The simulations rely on the ergodic approximation, valid when the orbits in the trap are fast compared to the mean Q-He collision time. The simulations also assume the Q density is small compared to the He density, so Q-He collisions need to be considered but Q-Q collisions can be mostly ignored. In all collisions, the small effect of the trap light on the buffer gas atoms is ignored. We also approximate all collisions as equally likely to occur at all times so that the ergodic theorem is relevant for drawing initial conditions, even though, as in Section V A, collisions are technically more likely to occur at times of Q's motion where its speed is higher.
In our simulations, we approximate the trap potential of a single lattice site j as a truncated harmonic potential, in order to efficiently sample initial conditions for collisions from the energy hypersurfaces.
Here, the ω i represent the trap angular frequencies ω x = ω y = 4η(j)k B T /m Q w(j) 2 = ω r and ω z = 8π 2 η(j)k B T /m Q λ 2 . We also ignore the constraint that orbits in the trap conserve the z component of angular momentum. We stress that we are still making the ergodic approximation, we are simply taking initial conditions for collisions from a simpler physical system than the fully nonlinear trapping potential. The simulation volume for lattice site j, centered on the antinode at z j , is defined as the region where the approximate trap potential in Eq. (C1) is nonzero, and consists of an ellipsoid with volume V e (z j ) = w 2 (z j )λ/3, assuming the lattice site is near the trap center where wavefront curvature can be ignored.
In action-angle coordinates, sampling from the energy hypersurface amounts to drawing three random angle variables and three random actions J i under the constraint that E + ηk B T = 3 i=1 ω i J i . This is efficiently done by drawing two random numbers between 0 and 1 and using them as a partition of the interval into three randomly drawn pieces. These are translated back to regular phase-space coordinates to draw initial conditions for bound Q molecules in collisions.
There are two forms of the ergodic loading simulations, namely, that which computes the loading fractions f i , and that which computes the mean number of Q-He collisions needed to eject a trapped molecule, k Q-He , which are discussed in Section V B 2. The results of both are required to construct the loading model in Section V B 3. For the simulation of the loading fractions f i , an untrapped Q molecule is first spawned randomly in the simulation volume. Its velocity is drawn from a Maxwell distribution at the buffer gas temperature T = 1.5 K and its position x is drawn uniformly within the simulation volume, before its kinetic energy is increased due to the local trap potential U (x). Its speed is computed as v = 2(E − U )/m Q , and its direction is randomized. The initial condition for a first billiard-ball collision is now set by drawing a buffer gas atom from a Maxwell distribution at T . If the molecule is bound by the collision, only its energy needs to be retained, since before every subsequent collision new initial conditions are drawn by randomly sampling the energy hypersurface within the ergodic approximation. There is no need to integrate the Q trajectory. The simulation continues until Q's energy either becomes positive and the molecule escapes, or it stays negative and falls below the energy of a molecule drawn from a Boltzmann distribution in the trap at T .
For the simulation of the number of Q-He collisions needed to eject a trapped molecule, k Q-He , a molecule is initially drawn with an energy from a Boltzmann distribution within the trap, and with phase-space coordinates drawn ergodically from the corresponding energy hypersurface. Billiard-ball collisions are repeatedly enacted with buffer gas atoms until the molecule is eventually ejected from the trap (E > 0), giving k Q-He .

D. Trap heating caused by laser noise
Laser noise can cause trap heating in two ways [174,175]. Firstly, laser intensity noise (RIN) at twice the trap angular frequencies ω r = 2π × 2 MHz and ω z = 2π × 68 MHz lead to parametric driving of trapped molecules, resulting in a heating rate Γ RIN at which energy grows exponentially in the trap: where S (ω) is the one-sided power spectrum (in units of dBc/Hz) of the relative intensity noise ∆I/I at an angular frequency ω [174].
We have measured the RIN of our laser oscillator (NKT Koheras Adjustik Y10) to be the measured RIN values outside the cavity. We note that the value of Γ RIN z is within a factor of two of the shot-noise of the circulating power inside the cavity. We also note that imperfect spatial mode matching will reduce the cavity's suppression. In particular, 2ω r = 2π × 4 MHz is close to the difference in resonance frequencies of ∆f 10 = 3.6 MHz between the TEM 01 /TEM 10 modes and the fundamental TEM 00 mode. Finally, we do not expect that the high-power fiber amplifier to be used in the experiment adds substantial RIN at the frequencies of interest.
Secondly, movement of the trap center by an amount x(t) results in a linear growth in average energy, which now could lead to significant heating due to the quartic dependence on the trap frequency. Here, S x is the one-sided power spectrum of position fluctuations in the trap center [174]. However, we expect acoustic noise at MHz frequencies to be sufficiently small.
E. Two-body loss models

Universal loss
In optical trapping experiments on bialkali molecules, two-body loss is a result of so-called "universal loss", in which a large fraction of collisions between trapped molecules lead to loss.
One argument suggests that our trap will not experience universal loss. In this argument, the loss mechanism involves strong perturbation of molecular energy levels by a partner, which leads to an excited state of the two-body complex located one photon energy of the trap light above the ground state [176]. For long enough lifetimes of the complex, the probability to reach this excited state, causing subsequent loss of the molecules, then approaches unity.
However, the interactions between SCS molecules are often of order 20 meV [177,178], which is two orders of magnitude too small to perturb the excited states by several eV and unlock a single-photon excitation of the complex. Moreover, the complex lifetimes of SCS molecules (based on a semiclassical model presented in [153]) are only of order picoseconds, so the notion of "sticky molecular collisions" seems inapplicable to our trap. Hence, one expects not to see universal loss for Q.
However, recent literature has cast some doubt over this explanation [129]. Therefore, we consider a worst-case scenario in which universal loss occurs for all Q-Q collisions in the trap. This worst-case scenario provides an upper bound for the two-body loss rates expected from general two-body loss processes that are difficult to estimate in the cold, high-intensity conditions of our trap. We do not, however, treat Q-He collisions as lossy, as the interactions are only O(5 meV) [126].
In universal loss, the functional form of the two-body loss coefficient β Q-Q depends on whether or not an s-wave two-body interaction is allowed, which in turn depends on whether the molecules can be regarded as indistinguishable and Bosonic. Since the s-wave loss rates are higher than the p-wave loss rates, we must assume in a worst-case estimate that we will be limited by s-wave scattering, in which [179] where the Van der Waals length a is here defined by the C 6 coefficient through a = 0.4778(2µC 6 / 2 ) 1/4 [180], and µ is the reduced mass of the two bodies.
For many small molecules, a typical value of C 6 is a few tens to hundreds of atomic units (E h a 6 0 ), and the C 1/4 6 dependence of β Q-Q means a good estimate is not critical. Using C 6 coefficients from an approximate formula in [181], we obtain, including an extra factor of two as noted in [179], values of β Q-Q slightly less than 1 × 10 −10 cm 3 /s for many molecules, and assign β Q-Q = 10 −10 cm 3 /s for Q, which is consistent with typical values in [179]. In the main text, we keep in mind the strong likelihood that universal loss does not occur for Q, in which case two-body losses may not limit the trapped density. By experimentally observing the two-body loss rate for molecules in our trap, we will be able to provide useful insight into the nature of cold collisions and the mechanism of universal loss.

Collision-induced absorption
In high-density gases, optical transitions that are forbidden for single molecules can become allowed through collision-induced absorption (CIA) [182]. To our knowledge, of the molecules in Table I, only O 2 exhibits CIA at 1064 nm. From Fig 1 (b) of [183], the absorption of this feature is α CIA = 10 −6 cm −1 am −2 = 1.39 × 10 −45 cm 5 /molecule 2 . Assuming every CIA event leads to loss, the resulting two-body loss coefficient is β = α CIA I/ ω = 2.2 × 10 −15 cm 3 /s, which is small compared to other two-body loss coefficients considered in Fig. 6. It therefore appears that even in the worst case of CIA features centered at 1064 nm, the densities required to induce absorption are high compared to those we expect to see in our experiment. CIA effects can therefore be ignored in our experiment.

F. Trap heating and loss from rotational Raman scattering
Here, we compute the rotational Raman scattering (RSS) loss term Γ RR of Eq. (7). RRS from the J = 0 rotational ground state to the J = 2 excited state increases a molecule's internal energy by ∆E RR = 6B, which is ∼15 K, and hence greater than T trap , for many molecules in Table I. Rotationally inelastic collisions typically occur once per ∼ 10 elastic collisions [123], which is 10 4 s −1 , much faster than RRS rates in Table I. We therefore assume molecules are only rotationally excited once before being immediately rotationally cooled to the ground state, whereby their rotational energy is converted to translational kinetic energy.
In Q(J = 2)-Q(J = 0) inelastic collisions, which occur at a rate ν Q-Q,i , on average ∆E RR is split evenly between the molecules, so most likely a molecule will be ejected from the trap.
On the other hand, in Q(J = 2)-He inelastic collisions, which occur at a rate ν Q-He,i , only µ/m Q = 1/8.5 of ∆E RR goes to the molecule's translational energy, the rest being carried away by the much lighter He atom. Therefore, while the fraction ν Q-Q,i /(ν Q-He,i + ν Q-Q,i ) of RRS events simply lead to loss, the remaining fraction ν Q-He,i /(ν Q-He,i + ν Q-Q,i ) do not. They instead lead to heating of the sample of N trapped molecules, increasing the total kinetic and potential energy of the trapped sample at a rateĖ h = µ m Q ∆E RR W RR N . This heating will lead to some molecule loss, which we can bound above by assuming, conservatively, that the buffer-gas cooling does not compensate forĖ h at all, and the heat is instead balanced exclusively by molecule loss. Every time a thermalized molecule is lost from the trapped sample, it carries away on average (η + κ − 3)k B T of energy (κ = (1 − Γ(5, η))/(ηΓ(3, η) − 4Γ(4, η)), where Γ is the incomplete Gamma function [134]). BalancingĖ h with this loss and combining with the Q(J = 2)-Q(J = 0) case, an upper bound for the overall loss rate due to RRS is This loss term has a nonlinear dependence on N as ν Q-Q,i depends on the trapped density.
To evaluate Eq. (F1), inelastic collision cross sections are required. For N 2 (Fig. 6 (d)), we approximate He-N 2 inelastic collision cross sections using the elastic cross section from [125] and the inelastic-to-elastic collision ratio from [184], and we take the N 2 -N 2 inelastic collision cross sections from [185,186]. We estimate the inelastic N 2 -N 2 and N 2 -He cross sections at 1.5 K to be 100Å 2 and 40Å 2 , respectively.

G. Laser-induced breakdown
Laser-induced breakdown occurs when an electron created in the laser beam, for instance by multiphoton ionization (MPI), is heated by inverse bremsstrahlung in collisions with neighboring molecules to above the ionization energy of a molecule, so it can free new electrons by impact ionization. This can lead to exponential growth in the electron number.
The process is well understood in high-intensity, pulsed lasers and at high gas densities [104,[187][188][189][190][191]. However, the treatment of breakdown in continuous-wave (CW) laser beams is often left as an afterthought, because CW intensities are rarely high enough to induce breakdown.
Fortunately, the theory of breakdown in CW laser beams is straightforward. In this case, breakdown can be treated as an electron diffusion problem: breakdown occurs when an electron spawns in the beam and causes more than one ionization event before it leaves.
Once outside the beam, the electrons will cool down rather than be collisionally heated, and eventually recombine with ions, and so will no longer contribute to further ionization events.
We do not consider the ion motion as contributing to breakdown, because the collisional heating rate of the ions is suppressed by their much higher mass.
At the relatively low buffer gas densities we are interested in, macroscopic treatments of breakdown [189] are not required, and it is instead numerically tractable to simulate entire electron trajectories through the laser beam. Electrons are spawned at the center of the beam in the radial direction, and within a Rayleigh range of the focus in the axial direction.
Based on [187], their initial energies are drawn from a Boltzmann distribution with mean energy I 0 /3, I 0 being the He ionization energy. For self-consistency, we have confirmed in our simulations that electrons formed from impact ionization events at intensities of order 100 GW/cm 2 have roughly this energy. The electron motion is modeled step-wise, with step sizes drawn from an exponential distribution of mean size l e , where l e is the electron mean free path, computed as 1 l e = gas i process j where i indexes between Q and He, n i represent the gas number densities, and σ ij is the cross section for electron collisions with gas i that lead to process j. The gas densities are taken to be consistent with Eq. (7) during the loading phase, and we assume there is no ion density, as electron multiplication begins before there is a large number of ions in the beam. The processes we consider are elastic collisions, electronically exciting collisions, and ionizing collisions. Exciting collisions are counted as ionizing collisions as the trap light would likely ionize any electronically excited molecule. In each case, we treat the gas particles as stationary targets due to their much lower velocity, and model the cross sections as spherically symmetric for simplicity.
In each collision, a number of photons may be absorbed due to inverse bremsstrahlung, which modifies the cross sections. For an exciting or elastic process with a laser-free cross section σ (0) E 1 →E 2 (E 1 ), the cross section is modified to [192][193][194][195] σ (0) where n represents the number of photons absorbed in the process, e the electron charge, E 0 the electric field amplitude, m e the electron mass, ω the optical angular frequency, J n the n th -order Bessel function of the first kind, and ∆k x the change (along the axis of optical polarization, here the x axis) in the electron's k vector in the collision. We do not modify cross sections in which a third body is created in this way.
In our simulations, we use a particular target, CH 4 , to represent Q, because electronimpact cross section data is difficult to choose for Q. The energy-dependent, laser-free cross sections for helium are compiled from [196] (elastic), [197] (exciting), and [198] (ionizing).
CH 4 cross sections are taken from Song [199], except for dissociation cross sections (which we take to be inelastic but not ionizing) which are found in [200] and multiplied by 2 to fit with Song's sparsely populated but likely more accurate data. The binary encounter dipole model ( [201] for helium and [202] for methane) is used to draw electron velocities after impact ionization events to continue their trajectories. Other molecules, like N 2 and O 2 , have similar orders of magnitude for their electron-impact cross sections [203,204], so we believe that our conclusions for breakdown in a mostly helium He-CH 4 mixture will generalize to other species we would like to trap.
One challenge with using these cross sections is that l e depends on the local light intensity and gas density. One would need to update l e depending on where the electron moves in one step, which is itself determined by l e . We overcome this using a "feedforward" method, in which l e is determined by the local intensity and densities at the previous step. We expect this to only cause small errors, and this approximation is needed to make the problem tractable.
Armed with the relevant cross sections, at each step of mean distance l e and random direction, we stochastically decide which process occurs and how many photons are absorbed with relative weights n i σ (n) ij . From the chosen process, we update the electron energy, and count whether an ionization event has occurred. We also update the electron energy by the ponderomotive potential, which causes minimal change to the simulation results. We stop an electron's simulation when it has moved three waists from the beam's central axis. We also track the total time taken for the electron to diffuse, and we can add electrons spawned from multiphoton ionization (MPI) during this time to compute the breakdown condition when MPI is included. We count how many electrons are created from each electron we spawned, and if this number is greater than one, we decide that breakdown has occurred.
The results are shown in Fig. 7. It is known that in the diffusive breakdown regime, the intensity threshold I th depends on the gas pressure p as [190] I th ∼ p −2/m , where m is a constant near 1. For breakdown dominated by MPI instead of electron diffusion, this constant m becomes quite large, and the pressure dependence becomes small.
We can see both of these regimes in Fig. 7. When MPI electrons are ignored, there is a strong dependence on pressure, but when they are added, I th becomes a slow-varying function of pressure at high enough intensities. Breakdown appears to occur at densities and intensities Simulations determine when one electron spawned in the beam creates more than one electron in the laser beam before it diffuses out. The solid curve includes electrons created independently through multiphoton ionization (MPI) before the electron diffusion out of the beam. The dashed curve ignores these electrons and purely counts electrons formed by impact ionization. The red-filled region is that bounded below by the solid curve, indicating the region where breakdown occurs. Inset are the CH 4 densities used in the simulation, which are taken from the steady state of Eq. (7), despite the fact that the model is not strictly valid at the comparatively high densities used in this simulation. Here, we fix n CH 4 ,LR /n He = 1/100 and assume universal loss in Eq. (7). Our conclusions on breakdown are not sensitive to these choices. much higher than we will require, so it will not limit our experiment.