Explanation of efficient quenching of molecular ion vibrational motion by ultracold atoms

Buffer gas cooling of molecules to cold and ultracold temperatures is a promising technique for realizing a host of scientific and technological opportunities. Unfortunately, experiments using cryogenic buffer gases have found that although the molecular motion and rotation are quickly cooled, the molecular vibration relaxes at impractically long timescales. Here, we theoretically explain the recently observed exception to this rule: efficient vibrational cooling of BaCl+ by a laser-cooled Ca buffer gas. We perform intense close-coupling calculations that agree with the experimental result, and use both quantum defect theory and a statistical capture model to provide an intuitive understanding of the system. This result establishes that, in contrast to the commonly held opinion, there exists a large class of systems that exhibit efficient vibrational cooling and therefore supports a new route to realize the long-sought opportunities offered by molecular structure.

T he internal structure of molecules offers a host of scientific and technological opportunities 1 , including the manipulation of quantum information, critical insight into quantum chemistry and improved tests of the Standard Model. To utilize this potential of molecules typically requires the preparation of molecular samples at very low temperatures, where only a single quantum state is occupied. Unfortunately, experiments attempting to reach these temperatures by buffer gas cooling have found that although the molecular motion and rotation are quickly cooled to the cryogenic temperature 2,3 , the molecular vibration relaxes at impractically long timescales 4 . However, in a recent study 5 , sympathetic cooling through collisional interaction with laser-cooled atoms was demonstrated to be an alternative and efficient approach to quench molecular ion vibrational motion. Although this result was predicted based on a semi-classical argument 6 , it had not been verified by detailed quantum mechanical calculations. Further, it was not known how widely this technique could be applied, nor how to estimate its efficiency for other systems.
Here, we perform a detailed theoretical study of the Ca-BaCl þ system and compare to the recent experimental results, as well determine efficiency criteria to predict vibrational quenching rates in similar systems. Specifically, we first build an analytical model of the potential energy surface (PES) of the Ca-BaCl þ collision using a large grid of ab initio points, taking special care to accurately describe its long-range behaviour, which has an important role at very low collision energy. We then perform close-coupling calculations of the vibrational quenching of BaCl þ by collisions with Ca atoms. Although such calculations are nowadays straightforward and relatively fast, it is highly computationally demanding in the case of Ca-BaCl þ owing to the large mass and bonding energy. We therefore give a brief account of the method used to make the calculations feasible in a reasonable amount of computer time and compare the results with experiment. In addition, we use a scattering model based on quantum defect theory (QDT) with generalized short-range boundary conditions to gain an insight into the vibrational cooling of the molecular ions by the ultracold atoms. Finally, we compare the close-coupling vibrational quenching results with those obtained for four other similar systems: He-N 2 þ , He-NO þ , He-CH þ and Ar-NO þ and propose a very simple statistical capture model, which reproduces the close-coupling results and provides a simple means to estimate the efficiency of vibrational quenching for a given system. This result establishes that, in contrast to the commonly held opinion, there exists a large class of systems that exhibit efficient vibrational cooling and therefore supports a new route to realize the long-sought opportunities offered by molecular structure.

Results
Ca-BaCl þ PES calculation. To calculate the Ca-BaCl þ vibrational quenching rate, we first built an analytical model of the PES of the Ca-BaCl þ collision using a large grid of ab initio points. The electronic ground state of the CaBaCl þ complex, a 1 A 0 state, was calculated with the MOLPRO programme suite 7 (see the Methods for details) on a three-dimensional (3D) grid of points in the Jacobi coordinates space r, R and y. Here, r represents the BaCl þ bond length, R the distance between Ca and the centre of mass of BaCl þ , and y the angle between r and R, with the linear structure Ba-Cl-Ca corresponding to y ¼ 0°. The functional form V(r,R,y) of this PES is defined as the sum of the interaction energy V I between Ca and BaCl þ and the potential of the isolated diatomic BaCl þ , V BaCl : These terms are found by interpolation over the ab initio calculations (see the Methods for details). The results of these calculations are shown in Fig. 1 along two dimensions in the Jacobi space. Figure 1a shows the existence of a relatively deep potential well in good agreement with the charge-transfer nature of the bonding within this ionic complex. Figure 1b reveals the existence of two minimal structures and two saddle points connecting these equilibrium structures. It also shows that the potential is strongly anisotropic. Although the Ba-Cl bond length is only slightly extended by the interaction with the calcium, we observe in Fig. 1a that the vibrational potential of BaCl þ is significantly modified by the latter interaction. This indicates there is a significant coupling between the vibration of BaCl þ and the other modes of motion. This coupling is expected to promote vibrationally inelastic collisions.
Close-coupling calculation of the vibrational quenching rate. Calculation of the vibrational quenching rate on such a complex PES is normally done through the use of approximation Contour energies are regularly spaced by 1,000 cm À 1 . (b) Contour plot of the three-body interaction energy for r ¼ 5 a 0 . Below À 500 cm À 1 , contour energies are regularly spaced by 500 cm À 1 . Red contours correspond to positive energies, and blue to zero and negative energies.
techniques (see the Methods for details), which are not applicable at ultracold temperatures. Therefore, the use of the close-coupling method 8 is compulsory, despite the fact that several other features of this system make the calculation tremendously intense-in fact, to our knowledge, no similar calculation has ever been performed before. Owing to the strong long-range and anisotropic potential, the small value of the rotational constant of BaCl þ , and the large relative mass, the calculation must be performed to extraordinarily large distances (up to 2,000 a o ) with matrices for a given value of J and parity that are as large as 10 4 Â 10 4 . Therefore, several theoretical and numerical improvements were required to make the calculations possible (see the Methods for details), including the development of a new scattering code using asynchronous task parallelization to calculate the vibrational quenching rates. Figure 2b shows the calculated vibrational quenching rate coefficients for a selected set of rovibrational states of BaCl þ . A thick horizontal line represents the experimental measurement 4 of the population averaged vibrational relaxation rate for v ¼ 1 and v ¼ 2, where the length of the line is representative of the energy range in the experiment. The ion-induced-dipole Langevin law 9 is shown as a dashed line. The calculated values compare very well with the experiment, whereas the Langevin law is roughly double the experimental value. The close-coupling rate follows the Langevin law in the temperature domain of the experiment and departs from it at lower temperature-as discussed later, this departure is not due to quantum suppression 10 as one might expect. In Figures 2a, 2c, 2d, the vibrational and rotational quenching are compared for several initial rotational levels belonging, respectively, to the vibrational levels n ¼ 1, 2 and 3. The vibrational quenching is always larger than the rotational quenching. This very unusual result is due to the low value of the vibrational frequency of BaCl þ and the deep potential well, together yielding a strong coupling between many vibrational levels. This is in contrast with previously studied atom-diatom van der Waal neutral or ionic complexes [11][12][13] , where the potential well is usually not deep enough to couple even two different vibrational levels of the diatom. The only other possibility to obtain vibrational quenching comparable to rotational quenching is when the bond length of the complex is smaller than expected with a pure Van der Waals interaction, indicating the rise of chemical bonding induced by electron sharing between monomers. This is, for example, the case of the He-CH þ complex 14 .
QDT calculation of the vibrational quenching rate. Despite its obvious utility, the close-coupling calculation does not lend itself to an intuitive understanding of the collision physics. Therefore, we have also performed a QDT calculation [15][16][17] of BaCl þ -Ca vibrational quenching, where the radial Schrodinger equation is solved for the long-range, isotropic R -4 induction term from large separation and matched with boundary conditions at short range describing the amplitude and phase of flux returning from the chemical bonding region Z cm (E) and d cm (E), respectively (see the Methods for details).   (a) Comparison between the vibrational and rotational quenching rate coefficients for BaCl þ in the initial states (n ¼ 1, j ¼ 0,1,2,3,4). The first and second numbers designate, respectively, the initial vibrational and rotational quantum number of BaCl þ . (b) Comparison between the vibrational quenching rate coefficients of several excited rovibrational levels (n, j) of BaCl þ resulting from collisions with Ca with the experimental results and with the Langevin law. The first and second numbers designate, respectively, the initial vibrational and rotational quantum number of BaCl þ . (c) Comparison between the vibrational and rotational quenching rate coefficients for BaCl þ in the initial states (n ¼ 2, j ¼ 0,1,2,5). The first and second numbers designate, respectively, the initial vibrational and rotational quantum number of BaCl þ . (d) Comparison between the vibrational and rotational quenching rate coefficients for BaCl þ in the initial state (n ¼ 3, j ¼ 2). The label of each curve designates the final vibrational level. Figure 3a shows the quenching rate coefficient in the so-called universal limit, where all collisions reaching short range lead to quenching, that is, Z cm (E) ¼ 0. This rate agrees reasonably well with experiment and theory at the experimentally relevant energies, but dramatically overestimates the quenching rate at low energies. We thus conclude that the reduction in quenching rate at low energy is not due to quantum suppression as would be expected in systems with shorter ranged potentials 8 . Therefore, we match the QDT result to the close-coupling calculation by parameterizing the short-range boundary condition as d cm (E) when 0rZ min rZ cm (E)r1 and Z cm (E) ¼ Z min or 1 otherwise, where Z 0 ,Z min , Z 1 and d 0 are determined by least-squares fitting (see the Methods for details). Figure 3b shows the good agreement between this fit and the close-coupling calculation and indicates that the suppression of the rate coefficient at low energies comes from lower quenching probabilities for small partial waves, as might be expected since the vibrational quenching is driven by asymmetry in the interaction potential.
Statistical capture model for the vibrational quenching rate. In addition to these calculations, it is desirable to develop a model for vibrational quenching, which provides useful rules-of-thumb to both better understand the collision physics and guide future experiments. Unfortunately, the available models of vibrational quenching 18 are not appropriate in this regime. The Landau-Teller model 19 , which is the most common model, is restricted to high collision energy. The only model available at low collision energy, owing to Dashevskaia and Nikitin 20 , cannot be applied to ionic systems as it makes the hypothesis that the long-range interaction potential is weak. Therefore, in order to rationalize our results, we introduce a statistical capture model following the recent work of Lara et al. on the ultra-low-temperature reactivity of D þ þ H 2 reaction 21 . This model divides the quenching process into two steps: the formation of the collision complex through the long-range interaction potential followed by its fragmentation to produce a vibrationally quenched diatomic cation. This scheme relies on the statistical ansatz, which is valid for collisions involving deep intermediate wells and/or ultra low temperature. The capture process is described by the Langevin rate, whereas the probability of vibrational quenching is proportional to the number of accessible vibrational channels, which is roughly measured by the ratio D e /o e , where D e is the triatomic molecular ion dissociation energy and o e the diatomic vibrational frequency. Therefore, we define the statistical capture rate constant k SC 10 , independent of the temperature, as: The validity of this model can be tested by comparing to known atom-ion quenching rates. Figure 4 shows a comparison of the statistical capture model to the quenching rate calculated from the Wigner threshold law for Ca-BaCl þ as well as several other species, which have been analysed by some of the authors [9][10][11][12] . Given the strong correlation with statistical capture model, we can create rules-of-thumb to aid future experiments realize efficient vibrational cooling. Namely, the vibrational quenching efficiency depends both on the strength of the long-range interaction potential, given by a, and on the density of states of the complex, given by D e /o e . The increase of the state density increases the lifetime of the complex and facilitates vibrational quenching. Loss rate coefficients based on the QDT as functions of collision energy. The black solid line corresponds to the loss rate coefficient assuming universal or completely absorbing short-range boundary conditions. Coloured lines correspond to loss rates contributions from individual partial waves and its projections with lr3. Finally, the dashed line is the loss rate coefficient found by the Langevin capture theory. (b) Loss rate coefficients K loss as functions of collision energy or temperature for a partial-wave dependent QDT optimized to agree with the coupled-channels BaCl þ (v ¼ 0, j ¼ 1) þ Ca vibrational quenching rate coefficient. The blue curve shows K loss (E) as a function of collision energy, whereas the black curve shows the corresponding thermally averaged rate coefficient. The red curve corresponds to the close-coupling results 'Vib 1 0' shown in Fig. 2. Shape resonances in K loss (E) are assigned by their partial wave l. The inset shows the short-range amplitude with parameters Comparison between the close-coupling vibrational quenching rate coefficients k CC 10 T ð Þ with the statistical capture rate k SC 10 T ð Þ for five different colliding systems with the diatomic cation in the initial state Discussion Among the possible explanations of the very large rate coefficients for the vibrational quenching, which is even larger than those for rotational quenching, one could be tempted to think about quasi resonant vibration-rotation energy transfer. This type of transfer is well documented for neutral systems like the collisions of Li 2 and HF with neutral noble gases at high temperatures 22 and has also been predicted theoretically to occur for the neutral He-H 2 collision at very low temperatures 23 . Quasi resonant vibrationrotation energy transfer, however, is expected for molecules in highly excited initial rotational states and takes place when the rapidly rotating diatom is stretched to its outer turning point and collinear with the atom 24 . In the present case, the calculated vibrational quenching rate coefficients are already unusually large for the lowest j ¼ 0-5 rotational quantum numbers of BaCl þ . This mechanism then cannot be at play.
The calculations presented here confirm the observation of efficient vibrational cooling in the Ca-BaCl þ system and provide a simple means to predict new systems that will exhibit efficient cooling. As all laser coolable atoms exhibit large polarizabilities, the only requirement on future experiments is to choose molecular ions that are strongly bound. This requirement is not particularly restrictive as other requirements of the effort, for example, non-reactivity with the ultracold atom and the existence of a large dipole moment, are typically only satisfied for strongly bound molecular ions. Therefore, sympathetic cooling of molecular ions with ultracold atoms appears to be poised to provide a generic and robust route to harness the potential of molecular structure for science and technology.
Methods PES calculation. The electronic ground state of the CaBaCl þ complex is a 1 A 0 state. In preliminary computations, we found that the wavefunction is dominantly described by a single electronic configuration. Therefore, the interaction energy between Ca and BaCl þ entities was computed using the coupled cluster method based on single and double electronic excitations and a perturbative treatment of triple excitations, CCSD(T). The counterpoise procedure was used to correct the interaction energy for the basis set superposition error. All electronic calculations were performed with the MOLPRO suite of programmes 7 . We used the def2-QZVPPD basis set 25 , which was augmented by diffuse functions. For chlorine, this results on (21s,15p,5d,2f,1g) primitives, which were contracted into (10s,7p,5d,2f,1g) and for calcium (24s,18p,6d,3f) contracted into (11s,6p,4d,3f) basis sets. For the barium atom, only the valence electrons were described by the (8s,8p,5d,3f) primitives contracted into (7s,5p,3d,3f), whereas we used also the quasi-relativistic 10-valence-electron pseudopotential ECP46MWB 26 for the 46 inner electrons.
In this manner, the interaction energy was calculated on a 3D grid of points in the Jacobi coordinates space. The 3D grid is a direct product of the three 1D grids spanning the Jacobi coordinates r, R and y. Here, r represents the BaCl þ bond length, R the distance between Ca and the centre of mass of BaCl þ , and y the angle between r and R, with the linear structure Ba-Cl-Ca corresponding to y ¼ 0°. For the r bond length, we selected 15 values ranging from 3.85 to 6.05 a 0 . The distance R took 50 values varying from 4 to 50 a 0 , with a variable step increasing from 0.2 to 10 a 0 . The y angle was uniformly distributed from 0°to 180°by step of 10°.
The evolution of the potential energy curve (PEC) and the dipole moment of the isolated BaCl þ diatomic along the r coordinate were calculated with the Davidson corrected multi-reference configuration interaction (MRCI þ Q) method 27 using the molecular orbitals calculated by the complete active space self consistent field (CASSCF) technique 28 . State-averaged CASSCF calculations (over 8 states with equal weight) were performed, along with 12 active orbitals and 6 active electrons. The basis set and the core pseudopotential defined above were used. A grid of 26 points from 3.4 to 20 a 0 was calculated.
The parallel and perpendicular static polarizabilities of BaCl þ were calculated using the finite-field procedure. For that purpose, we computed the ground-state CASSCF/MRCI þ Q energy, the CASSCF calculations being performed with the ground-state orbitals only. The polarizabilities were calculated on a grid of 17 points from 3.7 to 8 a 0 . The ionization energy of BaCl þ was obtained as the energy difference between the energy of BaCl þ and of BaCl þ þ , both taken in their respective electronic ground state. The potential of BaCl þ þ was calculated at the CASSCF/MRCI þ Q level where the CASSCF was state-averaged on two states and using 5 active electrons distributed in 12 active orbitals.
The functional form V(r,R,y) of the PES is defined as the sum of the interaction energy V I between Ca and BaCl þ and the potential of the isolated diatomic BaCl þ , V BaCl : V r; R; y ð Þ¼V I r; R; y ð ÞþV BaCl r ð Þ: ð3Þ The potential V BaCl is obtained by cubic splines interpolation of the ab initio data. The functional form of the interaction energy V I (r,R,y) is defined as the sum of the short-range V S and long-range V L contributions, combined with the switching function S(R): with the short-range function defined by where P l (cosy) are normalized Legendre polynomials. The coefficients C kl (R) have been first calculated at the points of the R-grid by the linear least squares method, using k max ¼ 10, l max ¼ 10 and r 0 ¼ 5 a 0 . Then, a cubic splines interpolation is used to obtain the coefficients for any value of R inside the limits of the R-grid. The switching function, which ensures a smooth connection between the long-range and the short-range functions, is defined by with R 0 ¼ 20 a 0 , and A 0 ¼ 0.14 a 0 À 1 . The long-range part is defined as a the sum of the leading terms of induction and dispersion energies 29 where the dipole moment m, the parallel a || and perpendicular a > polarizabilities of BaCl þ were obtained by polynomial interpolation of the ab initio data and where the coefficients a k , b k and c k have been determined by the linear least square method. The polarizability of calcium 30 is a Ca ¼ 168.71 a 0 3 , the ionization energy of calcium 28 I Ca ¼ 6.113 eV, and the ionization energy of BaCl þ I BaCl ¼ 13.9 eV. Let us note that, although using the long-range interaction model without any dependence in the internal coordinates is a customary approximation, we have chosen here to take in account the vibrational coordinate because this work focuses on vibrational quenching.
The spectroscopic properties of the fitted diatomic potential for BaCl þ (X 1 S þ ) are r e ¼ 4.899 a 0 , D e ¼ 39,278 cm À 1 , o e ¼ 333.8 cm À 1 and o e x e ¼ 0.837 cm À 1 . These values are close to the theoretical data previously reported 31 . The dipole moment calculated at the equilibrium bond length is m e ¼ 8.927 Debye. This calculation was done with the origin of coordinates fixed at the centre of mass, using the average atomic masses.
We display in Fig. 1a the contour plot of the total potential energy along the r and R coordinates, and in Fig. 1b the contour plot of the interaction energy along the R and y coordinates. The plot in Fig. 1a shows the existence of a relatively deep potential well in good agreement with the charge-transfer nature of the bonding within this ionic complex. Figure 1b reveals the existence of two minimal structures and two saddle points connecting these equilibrium structures. It also shows that the potential is strongly anisotropic. The existence of these stationary points on the ground potential of the Ca-BaCl þ system was checked by ab initio geometry optimizations where all coordinates were kept free. The Jacobi coordinates of the global minimum of V(r,R,y) are r ¼ 5.08 a 0 , R ¼ 6.54 a 0 and y ¼ 54.7°. This means that the calcium atom is bonded to the chlorine end of BaCl þ , with bond length r CaCl ¼ 5.34 a 0 and angle y CaClBa ¼ 87.1°. The electronic dissociation energy is D e T ¼ 6,263.8 cm À 1 . If we consider the interaction potential alone, the coordinates of the global minimum are then r ¼ 5.41 a 0 , R ¼ 6.49 a 0 and y ¼ 52.8°, and the corresponding dissociation energy is 7,441.9 cm À 1 .
Although the Ba-Cl bond length is only slightly extended by the interaction with the calcium, we observe in Fig. 1a that the vibrational potential of BaCl þ is significantly modified by the latter interaction. This indicates there is a significant coupling between the vibration of BaCl þ and the other modes of motion. This coupling is expected to promote vibrationally inelastic collisions.
At long range, the induction term that scales as R À 5 is the leading anisotropic interaction. The angular dependence (see equation 6) of this long-range interaction potential tends to align the three atoms Ca-Ba-Cl. Thus, at low collision energies, we expect a propensity for the alignment of the reactants, resulting in pendular states of BaCl þ around the linear structure. It is only for Ro10 a 0 that the propensity for the linear structure begins to disappear gradually with decreasing R. The minimum of the potential energy of the linear structure is found at r ¼ 4.94 a 0 NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11234 ARTICLE and R ¼ 9.50 a 0 , with a dissociation energy of 5,387.6 cm À 1 . This minimum is a saddle point, as bending the Ca-Ba-Cl structure leads to the global minimum.
Close-coupling calculations. The use of the long established close-coupling method 6 to study the dynamics of atom-diatom rovibrational inelastic collisions seems at first to be a simple task. However, the deep potential well of the system and both the very small vibrational and rotational quanta 29 of BaCl þ make the size of such calculations tremendous and one could be tempted to consider the use of more approximate alternative approaches like the infinite order sudden approximation 32 or the coupled states approximation 33,34 methods. Unfortunately, a prerequisite for use of the infinite order sudden approximation is that the rotational spacing of the diatom has to be negligible compared with the collision energy, which would become true above 10 cm À 1 , but is not satisfied in the energy domain of the ultracold experiment. Further, this approximation requires that the well depth be small compared with the collision energy, such that couplings with the closed channels remain small. With a well depth of 7,442 cm À 1 , the use of this approximation is then restricted to very high collision energies for this system. Similarly, the coupled states approximation method is expected to be valid for rotor states whose relative kinetic energy is large compared with the well depth; it is therefore not applicable to Ca-BaCl þ at very low collision energy.
The use of the close-coupling method then appears to be compulsory for this collision, despite the fact that several other features of this system that make such calculations especially intense. First, the strong angular anisotropy of the potential and the small value of the rotational constant of BaCl þ result in the coupling of many rotational levels. Second, the strong long-range potential requires propagation of the calculation to very long intermolecular distances. Third, the large relative mass necessitates many values of the total angular momentum quantum number J to reach convergence. The size of the matrices that have to be propagated for a given value of J and parity can then be as big as 10 4 Â 10 4 (above a collision energy of 1 cm À 1 ), which make the calculations unusually heavy, even at very low collision energy.
Therefore, several theoretical and numerical improvements were required to make the calculations possible. First, the vibrational coupling matrices hj vj |(r-r 0 ) k |j v 0 j 0 i, where j nj denotes the asymptotic rovibrational eigenfunctions, were calculated once and stored before the close coupling calculation. Thus, the large vibrational quadrature (about 100 points) necessary to compute these coupling terms was reduced to a summation over 11 terms (k runs from 0 up to 10) in the close-coupling calculation. Second, the calculation of the cross-sections for a single collision energy, even the lowest one, converged as a function of the total angular momentum, takes several hundred of hours of CPU time and reaches quickly (around a collision energy of 1 cm À 1 ) a thousand of hours. Therefore, we developed an MPI version of our Newmat code 35 using asynchronous task parallelization. The elementary task is the propagation of the wavefunction at one particular collision energy. The MPI code distributes N tasks over M processors. Because the tasks are independent, this parallelization scheme requires no overhead and feeds efficiently all the processors.
We included 25 rotational levels in each of the 19 vibrational levels used to perform the calculations. The propagation step size was taken to be 0.015 a 0 and the maximum distance of propagation was 400 a 0 . The calculations were performed in the collision energy interval (10 À 6 , 10 0 ) cm À 1 . The relative convergence criterion of the inelastic cross-sections as a function of the total angular momentum was taken to be 0.1% for the lowest energies, 1% for intermediate energies and 5% for the highest energies around 1 cm À 1 .
Quantum defect theory. Despite its obvious utility, the close-coupling calculation does not readily lend itself to an intuitive understanding of the collision physics. Therefore, we have also developed a QDT for scattering between BaCl þ and Ca based on solutions of a single radial Schrödinger equation where the potential between the molecular ion and neutral atom is dominated by its longest-ranged and isotropic induction potential, V(R) ¼ À C 4 /R 4 . At R ¼ R short , where R short is a characteristic short-range separation defined more precisely below, these solutions are uniquely specified by a boundary condition with a limited number of parameters that summarize the reflections (elastic collisions) and absorption (quenching collisions) in the complex, anisotropic and possibly chaotic evolution at short separations. By carefully choosing the collision-energy and partial-wave dependence of the boundary condition, we reproduce the temperature dependence of the rate coefficients for rovibrational relaxation as obtained by close-coupling calculations. The behaviour of these boundary conditions gives a simple, intuitive picture of the collision processes.
The long-range induction potential is much longer ranged than the van der Waals potential between two neutral atoms. As a result, many partial waves contribute to the rates; in fact only for temperatures of 100 nK are collisions s-wave dominated. Nevertheless, as we will show below, even for mK collision energies quantum mechanical effects profoundly affect the scattering of BaCl þ and Ca and play an important role in the description of the interplay between inelastic and elastic collisional processes.
In our QDT, we numerically solve the single-channel radial Schrödinger equation for collision energy E, partial wave l and projection m from separation R ¼ R short to . The functions exp(±i[y À p/4]) can be recognized as Wentzel-Kramers-Brillouin solutions of a À C 4 /R 4 potential at zero collision energy and partial wave. The short-range amplitude Z lm ðEÞe 2idlmðEÞ with real-valued functions Z cm (E) and d cm (E) determines the flux returning from the chemical bonding region where all three atoms are separated by no more than R short ooR 4 ; we use R short E30a 0 , just outside the separations where the electron wave functions have significant overlap. The short-range amplitude, in principle, are determined by the close coupling simulations. Flux conservation requires that 0rZ cm (E)r1, where Z cm (E) ¼ 0 corresponds to the case where no flux is returned from short range and Z cm (E) ¼ 1 corresponds to the case where everything is reflected back. The phase d cm (E) describes the relative phase shift of the returning flux. In the limit R-N, the wavefunction approaches c lm R ð Þ ! e À iðkR À lp=2Þ À S ii E; lm ð Þe iðkR À lp=2Þ , where E ¼ ' 2 k 2 =ð2mÞ, k is the relative collision wavevector, and S ii (E,lm) is the diagonal S-matrix element from which elastic and total inelastic rate coefficients can be determined. In fact, the total inelastic rate coefficient is K loss ðEÞ ¼ P lm K lm loss ðEÞ with where v r is the relative velocity. Here, we have used flux conservation or the unitarity of the S-matrix to rewrite the loss rate coefficient solely in terms of the diagonal S-matrix element, S ii (E,lm). At ultracold temperatures, only a few partial waves l contribute as for higher l the centrifugal barrier prevents the atom and molecule from approaching each other and K lm loss ðEÞrapidly goes to zero with increasing l. The À C 4 /R 4 potential is the largest energy scale at R ¼ R short and, consequently, we initially assume that Z cm (E) and d cm (E) are independent of collision energy E, partial wave l and projection m. In fact, ultracold reactions between neutral KRb and K, with a long-range van der Waals dispersion potential, forming K 2 have been successfully modelled 36 with Z cm (E)oo1. It is then natural to first study the universal limit of Z cm (E) ¼ 0, that is, all short-range collisions lead to quenching, for our induction potential. Figure 3a shows the universal loss rate coefficient as a function of E from E/k ¼ 1 nK to 0.01 K. At very low energy, the rate is dominated by s-wave scattering and is a decreasing function of E. At higher energy, other ARTICLE partial waves contribute creating a rate coefficient that weakly oscillates around the rate coefficient predicted by Langevin capture theory. Although the QDT universal approximation roughly agrees with both the experimental determination and the close-coupling calculation for energies above 1 mK, the universal result dramatically over estimates the quenching rate for temperatures below 1 mK. Therefore, in order to reproduce the temperature dependence of the coupled-channel calculation with QDT, we must use a collision-energy-and partial-wave-dependent short-range amplitude Z lm ðEÞe 2idlmðEÞ . For practical purposes, the number of parameters that describe Z cm (E) and d cm (E) must be limited and, here, we use the fact that QDT aims to only describe cross-sections over a small range of collision energies and limited number of partial waves. In our case, we need to represent collision energies below 0.1 K and partial waves up to 35. We then observe that at R ¼ R short the rotational energy ' 2 lðl þ 1Þ=ð2mR 2 short Þ for l ¼ 35, although much smaller than the À C 4 =R 4 short potential energy, is much larger than the relevant range of E. Hence, we can assume that the short-range parameters only depend on l. In fact, we further limit the parameterization to d cm (E) ¼ d 0 and Z cm (E) ¼ Z 0 þ Z 1 l(l þ 1) when 0rZ min rZ cm (E)r1 and Z cm (E) ¼ Z min or 1 otherwise. There are four free variables Z 0 , Z min , Z 1 and d 0 . Figure 3b shows the best fit of QDT to the BaCl þ (v ¼ 1,j ¼ 0) þ Ca vibrational quenching rate coefficient. The figure shows both its energy dependence and the corresponding thermally averaged value. The partial wave dependence of Z cm (E) is shown in the inset. The energy-dependent rate coefficient is less than 10 À 10 cm 3 s À 1 for E/kE1 mK. For these energies only l ¼ 0 and 1 contribute, Z cm (E) is close to one, and losses are suppressed. The rate coefficient then increases, interrupted by multiple resonances, to a rate coefficient of about one half of the Langevin rate coefficient near E/kC0.1 K, where nearly 35 partial waves contribute and Z cm (E) ¼ Z min . In the thermally averaged rate coefficient, the resonances, except for those below T ¼ 100 mK, have washed out. Rate coefficients for other rovibrational levels (v,j) can be found in a similar manner and lead to slightly different parameters Z 0 , Z min , Z 1 and d 0 .
The location and origin of the series of resonances explains much of our successful fit. They are shape resonances behind the even partial-wave centrifugal barriers. Analysis of the analytic E ¼ 0 solution of the À C 4 /R 4 potential with short-range amplitude Z cm (E ¼ 0) ¼ 1 and l-independent phase d cm (E ¼ 0) has shown that if an E ¼ 0 bound state occurs for partial wave l then such bound state also exist for the ..., l-4, l-2, l+2, l+4, ... partial waves. This condition is only satisfied for two values of the short-range phase d 0 , one corresponding to all even partial waves, one to the odd waves. For different values of d 0 as well as Z cm (E)o1, these bound states become shape resonances where the inelastic loss is resonantly enhanced. For BaCl þ þ Ca collisions with its Z cm (E)EZ 0 E1 for small partial waves, we find that near the optimal d 0 resonances occur for odd partial waves and, in particular, near E/k ¼ 10 mK the rate coefficient is sensitive to the location of the l ¼ 3, f-wave resonance. For E/k40.5 mK, where lZ8 wave collisions become prominent, the amplitude Z cm (E) deviates sufficiently from one, such that resonances from even partial wave collisions are observed. Finally, for the optimal conditions, the absolute value of the elastic scattering length a for s-wave BaCl þ þ Ca collisions is larger than the natural length scale of the C 4 potential. In this case, the finite Wigner-threshold prediction for the loss rate coefficient is only reached for collision energies below E/k ¼ 1 nK not shown in Fig. 3b. Figure 5 predicts the ratio of the QDT elastic and loss rate coefficient as a function of temperature. Ratios much larger than one indicate that the kinetic energy in the centre-of-mass motion of the molecular ion is more efficiently relaxed than its vibrational energy. For example, for BaCl þ ions at T ¼ 0.1 K, the ratio is close to one. We anticipate two cooling stages for BaCl þ ions in collisions with ultracold and highly polarizable Ca atoms. In the first stage cooling or relaxation of the rovibrational states occurs, as the small ratio K elastic /K loss implies that translational cooling in this region will be inefficient. Only when the molecules are in the lowest rovibrational state, the second stage, will elastic collisions cool the external motion of the molecules.
Traditional cooling schemes often involve collisions of ionic molecules with He gas at cryogenic temperatures. Here, we analyse the ratio K elastic /K loss for two ionic systems He-N 2 þ and He-CH þ to demonstrate a different cooling mechanism than for cooling with Ca. First, we fit the He-N 2 þ and He-CH þ vibrational quenching rate coefficients obtained in the coupled-channel calculations 10,12 to our QDT theory with partial wave-dependent short-range parameters. The parameters of our best fit are Z 0 ¼ 1-5.7*10 À 6 , Z 1 ¼ À 2*10 À 6 , d cm (E) ¼ 0.76p for He-N 2 þ and Z 0 ¼ 1-7*10 À 3 , Z 1 ¼ 0 and d cm (E) ¼ 0.06p for He-CH þ . Figure 6 shows the temperature-dependent ratio K elastic /K loss for these two systems. It is evident that for both systems the elastic rate coefficient is much larger than that of the inelastic processes. Hence, for a molecular ion in a given excited rovibrational level its translational motion will be cooled first. Occasionally, an inelastic process relaxes this internal state at the cost of rapid increase of the translational temperature. Elastic collisions will then start the cooling all over again.
Statistical capture model. To rationalize our results, we compare in Table 1 five collisions involving a diatomic cation and a rare gas. We report in this table the dissociation energy of the complex, the vibrational frequency of the diatomic cation, the relative mass of the system, the polarizability of the rare gas atom and the imaginary part of the scattering length b 10 associated with the collisional vibrational quenching of the lowest rotational level of the first excited vibrational level of the diatomic cation. In the Wigner regime, the quenching rate coefficient is independent of the temperature and directly proportional to b 10 : The values reported in the table were obtained from our close-coupling calculations for the collision energy of 10 À 6 cm À 1 , where only the s-wave contributes to the collision. The propagation at such low collision energy has to be performed up to very large values of the intermolecular separation coordinate. For Ca-BaCl þ , this maximum distance had to be extended to 2,000 a 0 . Because of the long-range ion induced dipole potential, which furthermore requires applying modified effective range theory 37 , the use of the scattering length approximation is limited to even lower collision energies. The usual formulas linking the real and imaginary parts of the scattering length with the scattering S matrix are in this case valid only for collision energies lower that 10 -10 cm À 1 . The properties of the Ca-BaCl þ system, namely a very low vibrational frequency, a deep potential well, a large value of the three quantities-equilibrium diatomic distance, relative mass and atomic rare gas polarizability-are seen to differ strongly from those of the other systems. But what makes BaCl þ especially singular is the unusually small value of its vibrational frequency (330 cm À 1 )-compared with B2,000 cm À 1 for the other diatoms considered. As seen in Table 1, the vibrational quenching scattering length b 10 monotonously increases as a function of the triatomic ion well depth D e and also as a function of D e divided by the diatom vibrational frequency o e . This correlation was noticed long ago by many authors 18  b 10 is the value of the Close Coupling imaginary part of the scattering length for the collisional vibrational quenching rate coefficient of the (n ¼ 1, j ¼ 0) state of the cationic diatom, computed at the collision energy of 10 -6 cm À 1 . For each system, the relative mass m, the equilibrium intermonomers distance R e , the well depth D e of the triatomic ABC þ complex, the diatom vibrational frequency o e , the diatom rotational constant B rot and the polarizability of the impinging atom are reported. ARTICLE Landau-Lifshitz probability for vibrational quenching of neutral molecules used by Dashevskaya and Nikitin 18 , for example, is mainly dependent on both parameters o e and D e . At very low collision energy, the long-range part of the interaction potential plays a major role, whereas the ratio D e /o e offers a crude measurement of the efficiency of the vibrational energy redistribution. Therefore, we define the statistical capture rate constant k SC 10 , independent of the temperature, by simply multiplying the Langevin rate coefficient by this ratio D e /o e : This relation can be formally introduced through a capture statistical approach, which was presented in the body of this manuscript. This simple relation makes the vibrational quenching proportional not only to the D e /o e ratio but also to the polarizability of the ultra cold atom and to the relative mass of the system. These two last results were expected as they express the dependence of the quenching efficiency on the strength of the long-range interaction potential and on the density of states of the complex. The increase of the state density increases the lifetime of the complex and facilitates vibrational quenching. We notice, however, that the 4 He-CH þ system seems to behave differently than the other systems. This is effectively what can be seen on Fig. 7 where the Close coupling rate coefficients were reported for these five systems in the [10 À 7 , 1] Kelvin interval. The 4 He-CH þ rate coefficient is the only one for which the rate coefficient decreases instead of increasing above the Wigner regime. This behaviour was shown in our paper dedicated to this system 12 to be due to virtual state scattering 38 . This is indeed the only one of these five systems for which the real part of the scattering length is negative. Virtual state scattering increases the close coupling value of the quenching rate coefficient and our simple model then allows predicting a lower bound for the zero temperature limit for all the studied systems.