Roaming pathways and survival probability in real-time collisional dynamics of cold and controlled bialkali molecules

Perfectly controlled molecules are at the forefront of the quest to explore chemical reactivity at ultra low temperatures. Here, we investigate for the first time the formation of the long-lived intermediates in the time-dependent scattering of cold bialkali \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{23}\hbox {Na}^{87}$$\end{document}23Na87Rb molecules with and without the presence of infrared trapping light. During the nearly 50 nanoseconds mean collision time of the intermediate complex, we observe unconventional roaming when for a few tens of picoseconds either NaRb or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Na}_2$$\end{document}Na2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Rb}_2$$\end{document}Rb2 molecules with large relative separation are formed before returning to the four-atom complex. We also determine the likelihood of molecular loss when the trapping laser is present during the collision. We find that at a wavelength of 1064 nm the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Na}_2\hbox {Rb}_2$$\end{document}Na2Rb2 complex is quickly destroyed and thus that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{23}\hbox {Na}^{87}$$\end{document}23Na87Rb molecules are rapidly lost.


Results
Potential energy surface and reaction paths. We have performed non-relativistic coupled-cluster (CC) calculations with single, double, and non-iterative triple (CCSD(T)) excitations using the Karlsruhe def2-TZVPP basis set 58 in order to find stationary points and reaction paths on the ground electronic state of Na 2 Rb 2 . We correlated all inner and valence electrons in the CC calculations. This state has a zero total electron spin. The CC calculations show that the stationary points predominantly occur when the four atoms are located at C s symmetries with all four atoms in a single plane. Figure 1 shows the potential energies at stationary points with C s symmetry as well as the corresponding planar locations of the four atoms. The reactants NaRb+NaRb are shown on the left-hand side of the figure, while the products Na 2 +Rb 2 are shown on the right-hand side. The product energy is only hc × 37 cm −1 higher Scientific Reports | (2021) 11:10598 | https://doi.org/10.1038/s41598-021-90004-0 www.nature.com/scientificreports/ than that of the reactants and compares favorably with the spectroscopically determined value of hc × 47 cm −1 derived in Ref. 35 . For the energetics of this system it is also worth realizing that the quantum-mechanical zeropoint energy of the individual Na 2 , NaRb, and Rb 2 dimers are hc × 80 cm −1 , hc × 53 cm −1 , and hc × 29 cm −1 in their ground X 1 + g , X 1 + , and X 1 + g state, respectively. This corresponds to a nearly identical combined zero-point energy of the reactants and products that is significantly larger than the hc × 37 cm −1 endothermicity of our CC calculations. On the other hand, the difference between the combined zero-point energies is much smaller than the endothermicity. Figure 1 further shows four minima. At the global minimum of the potential the atoms have D 2h symmetry as well. D 2h symmetry corresponds to planar geometries where the center of masses of the two Na atoms and of the two Rb atoms coincide and the Na-Na and Rb-Rb interatomic axes make a 90 • angle. The extrema have multiple equivalent atom geometries or configurations with the same potential energy reached by reflections through perpendicular planes or interchange of the Na atoms (or Rb atoms.) Only for the local minima at hc × −1840 cm −1 relative to the reactant state, do we show equivalent states as those will be helpful in describing reaction paths below. Our potential energies at the minima agree well with the predictions in Ref. 59 .
The five saddle points or transition states in Fig. 1 separate out into three "conventional" planar transition states (TS) and two roaming transition states (r-TS). For both types of transition states the curvature of the potential along one or more normal-mode directions is negative, i.e. has a so-called imaginary frequency. Conventional or tight transition states have imaginary frequencies that are on the order of typical zero-point energies, here ∼ hc × 100 cm −1 from the dimer zero-point energies. Roaming transition states have at least one imaginary frequency that is much smaller than these typical values. Only one of the "conventional" transition states has D 2h symmetry.
The two r-TS in our system have potential energies close to the reactant and product states and have geometries that closely resemble those of the reactant and product states. They are thus in the threshold regions of the PES. From Refs. 50,60 we learned that such states can initiate dynamically distinct collisions. For Na 2 Rb 2 , these transition states will play an essential role in roaming collisional dynamics as discussed below. The only outof-plane stationary point (except for equivalent geometries reached by identical atom interchange, reflections through the center of mass, etc.) is a second-order saddle point at an energy of hc × − 1070 cm −1 relative to the reactant state and has two directions along which the curvature is negative. It is not shown in Fig. 1.
Finally, Fig. 1 shows planar reaction pathways. These reaction paths are one-dimensional curves with C s symmetry connecting the reactants with the products along points of lowest energy through intermediate saddle points and local and global minima, i.e. minimum energy paths (MEPs). Two pathways connecting transition states and minima exist. The first pathway passes through the global minimum with its D 2h symmetry and crosses over the tight TS2 transition state at hc × −1468 cm −1 to reach the local minimum at hc × −2736 cm −1 and the product state. The TS2 saddle point can also direct the first pathway back to the global minimum and reactant state. The second pathway has a flatter energy landscape associated with the looser saddle points labelled D 2h TS and TS1 as well as the two equivalent local minima at hc × −1840 cm −1 . Both pathways encounter our two roaming transition states.
We end this subsection by noting that when not only the positions of the atoms lie in a single plane but also their momenta are directed in this plane then there exist no forces from the interatomic potentials that can break 2 Rb 2 ground-state potential based on coupled-cluster calculations with all four atoms located in the same plane, i.e. have C s symmetry. All energies are given in units of hc × 1 cm −1 . The energy is zero for two 23 Na 87 Rb molecules at their equilibrium separation of the X 1 + state. This state is shown on the left-hand side of the figure. The right-hand side shows the endothermic di-atomic 23 Na 2 + 87 Rb 2 ground-state limit at the respective equilibrium separations. Other stationary states, minima and saddle points, in the PES (red horizontal lines) are shown in between. Saddle points are labeled by the abbreviations TS and r-TS corresponding to transition states and roaming transition states, respectively. For all stationary states a diagram shows the location of the four atoms. Na and Rb atoms are indicated by small cyan and large magenta spheres, respectively. Dashed black and blue lines connect stationary states and define reaction paths. Figure was prepared in Adobe Illustrator CS 6 with molecular structure pictures prepared using GaussView 5.09 for Mac. www.nature.com/scientificreports/ C s symmetry. Similarly, if the positions and momenta of the four atoms lie on a single line, then they will always remain or move along this line.

Collision time of the intermediate complex.
In this subsection we determine the mean duration or mean collision time of the intermediate four-atom Na 2 Rb 2 complex. We do so by computing an ensemble of classical trajectories of two colliding 23 Na 87 Rb molecules using the Mercury/Venus96 code 44,61 . For these classical simulations we use an analytical ground-state PES based on the dimer-in-molecule (DIM) approach described in Methods. Essentially, this four-atom potential is constructed from the spectroscopically-accurate singlet X 1 + and triplet a 3 + pair-wise potentials for each of the dimers in the complex such that the total electron spin of the four doublet alkali-metal atoms is zero. Coupled-cluster calculations, described in the previous subsection, are too computationally demanding within the Mercury/Venus96 code. The DIM potential can be quickly computed for any geometry of the four atoms, while still giving a reasonable description of the four-body potential. We set the total energy in the classical simulations equal to that of two 23 Na 87 Rb molecules, each in their v = 0, J = 0 ro-vibrational state plus a fixed relative kinetic energy of k × 1 mK, where k is the Boltzmann constant. We use micro-canonical sampling for the initial relative orientation of the molecules and initial relative momenta within each molecule given the constraints from the total energy and zero point energies of the 23 Na 87 Rb dimers. We mainly present data for simulations with an impact parameter equal to zero. At cold collision energies collisional lifetimes (and other "global" properties) are nearly independent of impact parameter up to impact parameters about five times larger than the equilibrium separation of ground state 23 Na 87 Rb. We postpone further discussion of the role of the impact parameter until the end of the last subsection of Results. Finally, the calculations are started at a molecule-molecule separation of R s = 57a 0 , where the potential energy between the molecules is approximately hc × −1 cm −1 . The fixed propagation time step is 0.1 fs. Our calculations show that the total energy is conserved to hc×0.1 cm −1 on the order of the rotational constant of hc×0.07cm −1 . This change in total energy is sufficiently small for the purposes of finding collision times and excitation rates.
For the purpose of defining the duration or lifetime of a collision, a collision starts when for the first time the total kinetic energy is larger than two times the initial total kinetic energy. It ends when for the last time the total kinetic energy is less than two times this initial total kinetic energy. The time difference is the lifetime of a trajectory τ . The numerical simulations halt when the separation between products reaches R f = R s − 0.5a 0 (The choice of halting separation slightly less than R s is an artifact of the fact that for cold collision and our choice of R s products can not separate to infinity. We have verified that the conclusions in this section do not change when we use initial relative kinetic energies larger than hc × 1 cm −1 .) We then analyze the geometry of the four-atom system. Mainly, the products are simply those of the initial state or those where the two sodium atoms (or the two rubidium atoms) have interchanged. In rare cases the classical simulations lead to the homonuclear Na 2 +Rb 2 products. The endothermicity of the system, shown in Fig. 1, is smaller than the zero point energy of both NaRb+NaRb and Na 2 +Rb 2 . Thus with micro-canonical sampling the product Na 2 +Rb 2 can sometimes appear. In any quantum mechanical simulation this product cannot form. They do not significantly affect our value for the collision lifetime. Figure 2(a) shows the cumulative distribution function (CDF) of lifetimes τ as determined from approximately 1200 zero-impact-parameter trajectories on a double logarithmic scale. The cumulative distribution function at argument τ is the probability that the collision time is less than or equal to τ . We observe that there exist three distinct classes of trajectories. Class I includes trajectories with τ < 0.07 ns. This group corresponds to approximately 4% of the trajectories. Class II includes the trajectories with τ between 0.07 ns and 1 ns and has fewer members. Class III consists of long-lived trajectories with τ larger than 1 ns with a majority of cases around 50 ns.
The inset of Fig. 2(a) depicts the cumulative distribution function for paths with τ > 1 ns. This restricted CDF is well represented by the function 1 − exp(−τ/�τ �) corresponding to an exponential or Poissonian probability function of lifetimes. In fact, the mean lifetime τ is 50 ns from a least-squares fit. Figure 2(b) shows examples of the total kinetic energy as functions of time for one trajectory out of each of the three groups. The behavior of the kinetic energy for short-lived, class-I trajectories is intriguing. An in-depth investigation showed that for most of these trajectories the two NaRb molecules approach each other in near collinear orientation and collide head-to-tail, as sketched in Fig. 2(a), and their motion is quasi-one-dimensional with small excursions along perpendicular directions weakly breaking collinear symmetry. Moreover, the atoms come close only once, corresponding to the near zero kinetic energy around time 0.1936 ns. The existence of these trajectories follows from the observation that when the positions and velocities of the atoms are in a single common direction the forces that would move the atoms out of alignment are zero. The relatively-large 4% probability of such collisions is due to the small positive curvature of the potentials along all transverse directions. The smaller this curvature or frequency, the larger the acceptance angle for quasi-1D collisions. The rare class-II trajectories are initially head-to-head collisions with larger rotational energy leading to their multiple re-collisions with various orientations.
Finally, in long-lived, class III collisions the atoms undergo chaotic motion with all orientations, geometries, and velocities appearing at random times. Short intervals of regular motion, however, do appear once in a while and will be discussed in the next section.

Roaming dynamics and conical intersections.
Analysis of the long-lived, class III trajectories discussed in the previous section has led us to the observation of unconventional roaming pathways, previously observed in unimolecular and atom-molecule reactions with light atoms 48,50,62 . A main characteristic of roaming pathways is that they go through near-threshold roaming transition states (r-TS). They exhibit bond stretching with large amplitude relative motion of molecules. www.nature.com/scientificreports/ We observe multiple appearances of roaming pathways in the complex-forming Na 2 Rb 2 trajectories. Motion in the complex is once in a while interrupted by dimer molecules moving away from the complex with large amplitude angular motion. This lasts for a few tens of ps after which they return to again form a four-atom intermediate. This process continues until products or reactants have enough kinetic energy to fly away from each other. An example of such event is shown in Fig. 3. Panels (a) and (b) of this figure show the total kinetic energy, E k , of a long-lived trajectory with lifetime τ = 200 ns as functions of time. Panel (a) shows a time slice just at the beginning of the association of two NaRb into a tetramer. Panel (b) displays a roaming event in the trajectory. Panels (c) and (d) show the relevant pair separations. The bond between Na and Rb atoms is temporarily weakened, their separations are stretched out to ≈ 30a 0 , while at the same time Na 2 and Rb 2 molecules are formed near the roaming-TS on the right-hand side of Fig. 1. These homonuclear dimers slowly roam about each other in the van der Waals region of the potential. After about seven vibrational periods of the homonuclear dimers the roaming molecules return back to reform the tetramer. In other trajectories with long lifetimes, we have also observed roaming to NaRb+NaRb over similarly brief time periods corresponding to the roaming-TS on the left-hand side of Fig. 1.
Previous studies of roaming with hydrogen-containing molecules 63,64 emphasized the role that seams of conical intersections (CIs) between two PESs 45 play. Here, CIs are the geometries where the potential energies of molecular states with the same properties under symmetry operations are equal. The existence of CIs in alkali-metal tetramers was previously predicted by Refs. 30,65 . For a better understanding of roaming in Na 2 Rb 2 , we studied the potentials of the singlet ground state and the first singlet excited state within the dimer-in-atom model (See Methods for the calculation of the excited-state potential.) The dimer-in-atom model puts all CIs at D 2h geometries. Figure 4 shows that there exist three continuous CI seams in D 2h geometries, uniquely described by the homonuclear pair separations R NaNa and R RbRb . Two of these seams occur when either the R NaNa or R RbRb is significantly larger than their dimer ground-state equilibrium separation. The CI geometry where the potential energy is lowest on any of the three seams is also shown in the figure. This minimum-energy conical intersection (MECI) is about hc × 500 cm −1 higher in energy than our entrance-channel energy and, thus, does not overlap the classically allowed region for the quasi-classical trajectories. We must conclude that roaming does not require the presence of CIs. The complex landscape of the six-dimensional potential energy surface with its chaotic motion is sufficient for the appearance of roaming.
Light-induced transition to excited states. In this subsection, we describe the time-dependent quasiclassical dressed-state scattering of two cold 23 Na 87 Rb molecules in the presence of trapping laser light with a wavelength of 1064 nm. We also estimate the probability that the laser excites the four-atom complex during the 50 ns mean collisional lifetime found in the previous subsections.
For these goals we do not only need the electronic ground state but also excited states that have significant electric-dipole coupling with the ground state, i.e. those that are optically-active and can be resonantly excited. The dimer-in-molecule model cannot be used to determine all relevant optically-active excited-state potentials www.nature.com/scientificreports/ and electric dipole moments. Instead, we use density functional theory (DFT) to obtain the on-the-fly groundstate potential, as the four atoms classically evolve under the forces due to this potential, and time-dependent density functional theory (TD-DFT) to compute splittings between ground-and excited-state potentials as well as transition dipole moments. A description of these electronic-structure calculations is given in Methods.   www.nature.com/scientificreports/ Figure 5 shows ground-and spin-singlet excited-state potential energies as functions of time up to 8 ps for four of our trajectories as well as their initial state sampled from the v = 0, J = 0 X 1 + state of NaRb. Their initial relative kinetic energy is 1 mK at an initial dimer separation of 38a 0 . The trajectories are calculated using a variable time step size no larger than 4 fs. For the trajectories the collision did not finish before we stopped the simulations. In the figure we also show the energy of the dressed ground state corresponding to the groundstate potential plus the energy of a 1064 nm photon. The four-atom complex can be optically-excited when the potential energy of an excited state equals that of the dressed ground state. In fact, for the 1064 nm laser the potential of only the two energetically lowest excited states "cross" that of the dressed ground state. The figure shows multiple crossings within our ps time windows.
The computation of on-the-fly DFT ground and excited potentials is numerically demanding and we had to decrease the initial dimer separation from that used in the previous sections and limit the time evolution to about 10 ps. Error bands around excited-state potential energies in Fig. 5 are our estimates of the standard-deviation uncertainty of splittings between potentials based on additional EOM-CCSD calculations for the three energetically lowest potential surfaces.
The crossings shown in Fig. 5 immediately suggest a means to determine the probability of non-adiabatic transitions between the potentials. The physics reduces to a time-dependent two-state system, considered by Landau 56 and Zener 57 , with diabatic energies that vary linearly in time and a time-independent coupling due to the molecule-laser interaction. In the spatial domain this corresponds to a transition from one potential to another that is "vertical", i.e. where the atomic positions and velocities are unchanged in the transition. Landau-Zener theory then states that the transition probability of excitation from the ground state at the i th crossing, occurring at time t i and i = 1, 2, 3, · · · , is where energy is the coupling matrix element that depends on transition electric dipole moment d(t i ) to the relevant excited state and laser intensity I. Here, ǫ 0 is the vacuum permittivity. The rate of change of the energies at time t i is α(t i ) = v i β i , where v i is the velocity along the classical path R cl (t) of the atomic coordinates at t = t i and (1) www.nature.com/scientificreports/ is the spatial derivative of the potential energy difference between the ground-and excited-state potential, V gr (R) and V exc (R) , respectively, along the classical path at t = t i . Here, direction n cl = R cl /|R cl |.
Once excited the atoms accelerate and move according to the forces generated by the electronically excited state potential. The excited state can decay back to the ground state by spontaneous emission or undergo other Landau-Zener transitions, either back to the ground state or up to doubly-excited states whenever dressed-state potentials are resonant. In either case the system's molecular total energy has changed sufficiently from that of the initial collision state so that the molecules are no longer trapped by the laser and lost.
For our system typical values of |β(t i )| are ∼ hc × 10 −2 cm −1 /a 0 , while those of |d(t i )| range from 10 −3 ea 0 to 1 ea 0 , where e is the elementary charge. Moreover, for a laser wavelength of 1064 nm the mean time interval between crossings δt = �t i+1 − t i � is 0.23 ps based on averaging within a trajectory and our sample of trajectories. In fact, we find that the intervals t i+1 − t i have a Poisson distribution. For a laser at this wavelength and an intensity of 10 kW/cm 2 , the mean excitation probability for a crossing p = �p i � is 2.1 × 10 −5 , orders of magnitude smaller than one. From Eq. (1) it follows that p is proportional to laser intensity. The probability distributions for t i+1 − t i and p i , which we find to be uncorrelated, are not expected to change for propagation times significantly larger than 10 ps.
As we are unable to run the on-the-fly calculations out to collision times of order of tens of nanoseconds, we must derive a coarse grained model for the excitation process based on the inequality δt ≪ τ . First, we realize that after N crossings for each trajectory the total or cumulative likelihood to be excited, and thus lost, is after the first excitation. Next, we assume that t ≫ δt (and thus N ≫ 1 ) and sample p i and t i+1 − t i from their respective probability distributions. The averaged cumulative excitation probability at time t is then well approximated by by replacing p i and t i+1 − t i by their mean values, using the geometric series to evaluate the finite sum, and using that p ≪ 1 . On the right hand side of the equation, we defined the physically relevant laser-induced survival rate κ =p/δt of the collisional complex, which for typical laser intensities is proportional to the laser intensity.
For the laser with a wavelength of 1064 nm and an intensity of 10 kW/cm 2 the decay time 1/κ = 12 ns, much larger than δt , justifying our derivation. On the other hand the 12-ns decay time is on the order of the mean lifetime τ of the collisional complex. These two characteristic times of the collisional complex can be combined to give the likelihood that a NaRb molecule does not survive the collision. This likelihood is P( τ ) , the averaged cumulative excitation probability at time τ . For a 1064 nm laser at I = 10 kW/cm 2 its value is 1 − 0.016 , close to one. A NaRb molecule is not likely to survive a collision in the presence of 1064 nm light.
Finally, we make some observations regarding the relevance of the impact parameter on our conclusions so far. All reported data have been for zero impact parameter. From additional trajectory simulations (not shown) we found that for impact parameters as large as 40a 0 , approximately five times the di-atomic X-state equilibrium separations, the duration of the collision and mean excitation probabilities fluctuate by no more than 20%. (Impact parameters larger than 40a 0 are inconsistent with our initial molecular separation and have not been considered.) The independence of the relevant quantities with respect to impact parameter is consistent with the theory behind spiraling collisions, Langevin or orbiting cross sections, as well as the chaotic motion in the sixdimensional ground-state potential. Firstly, because for low collision energies and the attractive long-range vander-Waals interaction between two NaRb molecules, the maximum expected contributing impact parameter 66 is several times larger than 40a 0 . Moreover, chaotic motion with its inherent sampling of all allowed phase space implies that once the short-ranged four-body complex is formed at the beginning of the collision the typical or mean collision time is independent of impact parameter. Thus, we can conclude that the light-induced inelastic cross section is approximately P( τ ) times the Langevin cross section.

Discussion
We have performed a detailed examination of classical collision trajectories in 23 Na 87 Rb+ 23 Na 87 Rb interactions both with and without the presence of laser light trapping the cold molecules. We computed the mean collision duration of the four-body complex without laser excitations and showed evidence of a long-lived "collision complex" that includes roaming motion to briefly form homonuclear 23 Na 2 + 87 Rb 2 . We also determined the survival or collision time of the four-body complex in the presence of 1064-nm laser light, which can induce transitions to two electronically excited states, and found that it is of the same order of magnitude as the duration of the four-body complex for a typical experimental laser intensity of 10 kW/cm 2 . Under these circumstances 23 Na 87 Rb molecules will be rapidly removed from the optical trap.
Our investigation, which includes roaming pathways within the collisional complex, leads to different, muchshorter lifetimes than those from the RRKM theory, which does not have a concept of roaming pathways. We believe that the existence of roaming events during collisions can significantly modify the collisional lifetime and change the reaction outcome. One example of such effect is given in Ref. 67 , which reported the several orders of magnitude decrease of the rate constant for H-atom formation in photodissociation of C 2 H 5 molecule due to the roaming dissociation channel.
Recent experiments 37,38 , using chopped ODTs with a time duration of laser pulses of order of 100 µs and/or box-shaped potentials, attempted to decrease the molecular loss rate. No change was observed. Our calculations www.nature.com/scientificreports/ giving short collision times of 50 ns would predict that colliding molecules are excited by the trapping light and mostly lost from the trap. Other loss mechanisms must be present when the laser light is not present. Several open questions remain regarding the validity of the approximations used in the paper and must be answered with future investigations. These concerns include the need for full quantum simulations of the collision on the ground-state potential, treatment of the CIs, and more sophisticated models for laser excitations. We also envision a study of a wavelength dependence of the optical excitation. Longer wavelengths should lead to fewer excitations as, for example, excitations to the second excited electronic state become energetically forbidden. Those to the first excited state with its CIs with the ground state potential are always possible. Nevertheless, we feel that the general conclusions will not be invalidated.

Methods
Dimer-in-molecule potentials. For the determination of the mean duration of 23 Na 87 Rb+ 23 Na 87 Rb collisions and the tetramer conical intersections we require an easy to compute but at the same time reasonably accurate electron-spin-singlet ground-state PES. We use the six-dimensional non-relativistic potential based on the dimer-in-molecule theory of Ref. 43 . In this theory the electron wavefunctions of the tetramer are approximated by superpositions of products of four spin-1/2 ground-state alkali-metal-atom electron wavefunctions, such that the total electron spin is zero. Two (geometry-or position-independent) singlet electron wavefunctions can be constructed with straight-forward application of angular momentum algebra 68 . The corresponding 2 × 2 tetramer potential matrix is determined from matrix elements for the sum of six pairwise potential operators. Each pairwise potential operator is given by 1 S ij =0 P S ij V S ij (R ij ) , where S ij is the electron spin of atom pair (i, j), operator P S ij is the projection operator on states with pair spin S ij with P 0 + P 1 equal to the identity operator, and V S ij (R ij ) is the isotropic pair-wise X 1 + (or X 1 + g ) Born-Oppenheimer potential at pair separation R ij for S ij = 0 and the triplet a 3 + (or a 3 + u ) potential for S ij = 1 . The three matrix elements of the 2 × 2 matrix are sums of V S ij (R ij ) with weights given by matrix elements of a unitary transformation reflecting the different ways to couple four spin-1/2 atoms into pairs. These weights are proportional to a nine-j symbol and phase factors 68 .
We use the spectroscopically-accurate dimer potentials from Refs. [69][70][71] with adjustments to their repulsive walls at small pair separations, where the potential energy is larger than that at the dissociation limit. These adjustments removed unphysical trends going from the light Na 2 , to NaRb, to the heavier Rb 2 dimers. The extrapolation of spectroscopic data to small separations can not reliably specify the shape of the inner walls of potentials.
The energetically-lowest eigenvalue of the 2 × 2 potential matrix corresponds to the ground-state tetramer potential and is used in classical trajectory calculations. This PES is not separable but, nevertheless, does not include so-called non-additive contributions. We assume that the effects of such contributions on the mean lifetime are small. Here, we rely on the results of time-independent quantum-mechanical studies of reactive atom-dimer collisions with alkali-metal atoms and molecules 72,73 . They found that total reaction rate coefficients are not affected when the non-additive part of the potential is or is not included. Similarly, we assume that the lifetimes are not affected by the precise nature of the adjustments to the inner walls of the pair potentials. In addition, long-range four-body electric dipole-dipole and quadrupole-quadrupole interactions between the polar NaRb molecules are absent in our PES. We do note that our cold entrance channel NaRb+NaRb conditions only include the non-rotating J=0 ground state only. Then to first order these four-body interactions are averaged out and should not affect our estimated value of the collision time.
Conical intersections in 2 × 2 dimer-in-molecule model occur when the two eigenvalues of the matrix are equal. A necessary condition for such occurrence is that the off-diagonal matrix element is zero. Inspection of the corresponding weighted sum of V S ij (R ij ) for Na 2 Rb 2 shows that this matrix element is zero for all geometries with D 2h symmetry. Thus, for such geometries conical intersections occur when the diagonal elements of the 2 × 2 matrix are equal. Figure 4 shows the locations of the conical intersections in the dimer-in-molecule model.
The on-the-fly optical excitation model. The on-the-fly optical excitation model for 23 Na 87 Rb+ 23 Na 87 Rb collisions in the presence of a laser field has two steps. First, we performed classical trajectory calculations using the on-the-fly ground-state electronic potential energies and gradients from a density functional theory (DFT) method. Second, at each of the geometries along a trajectory we applied time-dependent density functional theory (TD-DFT) calculations to obtain energy splittings between ground and electronically excited states as well as the corresponding electric dipole moment. The molecular dynamics equations have only been propagated to about 10 ps.
The DFT calculations have been performed for total electron spin S = 0 using the hybrid functional wb97xd 74 . This functional is a long-range corrected functional that allows for a good description of the long-range van-der-Waals interaction between atoms and molecules. We use a correlation-consistent triple zeta basis (cc-pvtz) for the Na atom and the Stuttgart/Cologne-group effective-core potential ECP28MDF 75 for the Rb atom with Hill and Peterson's augmented correlation-consistent quadruple-zeta polarization potential basis (aug-cc-pvqz-pp) 76 . This choice of basis allows for reasonable computational times of the on-the-fly dynamics.
The TD-DFT calculations, which found about ten excited states, are based on the spin-unrestricted Coulombattenuated B3LYP functional (uCAM-B3LYP). We compare the lowest two TD-DFT excitation energies with those from coupled-cluster equation-of-motion with singles and doubles (EOM-CCSD) calculations using the same basis set and find a nearly 5% difference as shown in Fig. 5. This small difference is sufficient for our purposes.