Quantum dynamics of CO-H$_2$ in full dimensionality

Accurate rate coefficients for molecular vibrational transitions due to collisions with H$_2$, critical for interpreting infrared astronomical observations, are lacking for most molecules. Quantum calculations are the primary source of such data, but reliable values that consider all internal degrees of freedom of the collision complex have only been reported for H$_2$-H$_2$ due to the difficulty of the computations. Here we present essentially exact full-dimensional dynamics computations for rovibrational quenching of CO due to H$_2$ impact. Using a high-level six-dimensional potential surface, time-independent scattering calculations, within a full angular-momentum-coupling formulation, were performed for the deexcitation of vibrationally excited CO. Agreement with experimentally-determined results confirms the accuracy of the potential and scattering computations, representing the largest of such calculations performed to date. This investigation advances computational quantum dynamics studies representing initial steps toward obtaining CO-H$_2$ rovibrational quenching data needed for astrophysical modeling.

of the CO fundamental band probe the physical conditions in the inner disk region, ∼10-20 astronomical units (AU), the site of planet formation. Detailed modeling of such environments requires state-to-state inelastic rovibrational excitation rate coefficients for CO due to H 2 collisions, but current simulations are limited to approximate scaling methods due to the dearth of explicit data [18,23].
In this Article, we address the two issues outlined above by advancing the state-of-theart for inelastic quantum dynamics with a full-dimensional investigation using an accurate potential energy surface relevant to this scattering process, with a particular emphasis on the important region of the van der Waals complex. This was made possible through the accurate computation and precise fitting of a 6D CO-H 2 potential energy surface (PES) in the relevant region of the formaldehyde tetra-atomic system and the further development of the inelastic diatom-diatom time-independent scattering code, TwoBC [24], which performs full angularmomentum coupling, the CC formalism [1], including vibrational degrees-of-freedom. We first briefly describe the new CO-H 2 PES and its testing through comparison of rotational excitation calculations using the 6D PES and 4D PESs (neglecting vibrational motion) and available experimental results. The full-dimensional (6D), essentially exact, computations for rovibrational quenching of CO(v 1 = 1) due to H 2 collisions are presented resolving a two-orders-of-magnitude discrepancy between earlier 4D calculations which adopted various approximations [25][26][27][28]. Finally, the current results are consistent with the rovibrational quenching measurements for the CO-H 2 system, performed at the Oxford Physical Chemistry Laboratory from 1976-1993, which no prior calculation has yet been able to adequately explain [29][30][31].

RESULTS
The CO-H 2 potential energy surface. The CO-H 2 interaction has been of considerable interest to the chemical physics community for many decades with one of the first 4D surfaces for the electronic ground state constructed by Schinke et al. [32], which was later extended by Baĉić et al. [25]. The group of Jankowski, Szalewicz, and coworkers performed accurate 5D and 6D electronic energy calculations, but averaged over monomer vibrational modes, and performed several fits to obtain a series of 4D rigid-rotor surfaces, referred to as the V98 [33], V04 [34], and V12 [35] PESs. A 6D PES for formaldehyde was constructed earlier by Zhang et al. [36], but it was developed for reactive scattering applications and consequently, limited attention was given to the long-range CO-H 2 van der Waals configuration. Therefore as a prerequisite to 6D inelastic dynamics studies, we carried out an unprecedented potential energy calculation including over 459,756 energy points (see Methods for details). The potential energy data were then fit using the invariant polynomial method with Morse-type variables in terms of bond-distances [37,38]. The resulting 6D PES, referred to as V6D, is shown in Fig. 1 for one sample configuration. Some features of V6D are also illustrated in Supplementary Figs. S2 -S4.
Cross sections and rate coefficients. Time-independent quantum scattering calculations were performed using the CC formulation of Arthurs and Dalgarno [1] as implemented for diatom-diatom collisions in the 4D rigid-rotor approximation in MOLSCAT [39] and extended to full-dimensional dynamics as described by Quéméner, Balakrishnan, and coworkers [12,13] in TwoBC. In the first set of scattering calculations, the new 6D PES was tested for pure rotational excitation from CO(v 1 = 0, j 1 = 0, 1), where v 1 and j 1 are the vibrational and rotational quantum numbers, respectively. The crossed molecular beam experiment of Antonova et al. [40], who obtained relative state-to-state rotational inelastic cross sections, is used as a benchmark. The experimental cross sections were determined at three centerof-mass kinetic energies (795, 860, and 991 cm −1 ), but with an initial state distribution of CO estimated to be 75 ± 5% for j 1 = 0 and 25 ± 5% for j 1 = 1. Antonova et al. normalized the relative cross sections to rigid-rotor calculations done with MOLSCAT using the V04 PES of Jankowski and Szalewicz [33]. Comparison of the experiment with new 4D rigidrotor calculations on V12 and full-dimensional calculations on V6D are shown in Fig. 2. We find no difference in the excitation cross sections when using V04 or V12, while the RMS cross section errors between the normalized experiment and the V12 and V6D calculations range from 0.56-0.89×10 −16 cm 2 and 0.55-0.95×10 −16 cm 2 , respectively (See Supplementary   Table S1 and Supplementary Note 1). Clearly, a 4D rigid-rotor treatment of the dynamics is sufficient for describing rotational excitation at these relatively high energies.
The importance of full dimensionality for rotational excitation becomes more evident as the collision energy is reduced (see also Supplementary Figs. S5 and S6). Low-energy excitation cross sections for the process (1) or (0000)→(0100), using the notation defined in Methods, were measured by Chefdeville et al. [41] in a crossed-beam experiment. They obtained the excitation cross section for center of mass kinetic energies from 3.3 to 22.5 cm −1 . Though their energy resolution was limited, three broad features were detected near 6, 13, and 16 cm −1 attributed to orbiting resonances. Computed cross sections for process (1) using the 4D V12 PES and the fulldimensional V6D PES are presented in Fig. 3a. While both computations reveal numerous resonances, the resonances are generally shifted by 2-3 cm −1 between calculations. The energy and magnitude of the resonances are very sensitive to the details of the PESs, but differences may also be due to the relaxation of the rigid-rotor approximation with the use of V6D in TwoBC. In Fig. 3b, the two calculations are convolved over the experimental energy resolution and compared to the measured relative cross sections. Except for the peak near 8 cm −1 , the current 6D calculation appears to reproduce the main features of the experiment. RMS errors are found to be 0.355 and 0.228×10 −16 cm 2 for the V12 and V6D PESs, respectively. Further details on the rotational excitation calculations can be found in Supplementary Note 1. Now that the improved performance of the V6D potential for pure rotational excitation is apparent, we turn to rovibrational transitions. As far as we are aware, prior experimental [29-31, 42, 43] and theoretical [25][26][27][28] studies are limited to the total quenching of CO(v 1 = 1). In the scattering calculations of both Baĉić et al. [25,26] and Reid et al. [27], the 4D potential of Baĉić et al. [25] was adopted, in which two coordinates were fixed (r 2 = 1.4 a 0 , φ = 0), and various combinations of angular-momentum decoupling approximations for the dynamics were utilized (e.g., the infinite order sudden, IOS, and coupled-states, CS, approximations; see Supplementary Note 2). More recently, Flower [28] performed CC calculations on a parameterization of the 4D PES of Kobaysashi et al. [44]. These four sets of 4D calculations for quenching due to para-H 2 (j 2 = 0) are compared in Fig. 4 to the current 6D/CC calculations on the V6D surface for the case of j 2 = j 2 = 0. Corresponding state-to-state and total cross sections for collisions with ortho-H 2 are given in Supplementary   Figs. S7-S9 with further details in Supplementary Note 2 A cursory glance at Fig. 4 reveals a more than two orders of magnitude discrepancy among the various calculations. The large dispersion for the previous calculations is due to a combination of reduced dimensionality and decoupled angular momentum which makes it difficult to assess the reliability of each approximation. The current results, however, remove these uncertainties by utilizing i) a full-dimensional (6D) PES, ii) full-dimensional (6D) dynamics, and iii) full angular-momentum coupling. The sharp peaks in the cross sections over the 1-10 cm −1 range in the 6D/CC results are due to resonances [45,46] [29][30][31] reported by Reid et al. [27]. The comparison is not straight-forward because i) the measurements correspond to an initial thermal population of H 2 rotational states, ii) the initial rotational population of CO(v 1 = 1, j 1 ) was unknown, iii) the experimental rate coefficients for ortho-H 2 are estimated from para-and normal-H 2 measurements, and iv) the contribution from a quasi-resonant channel, dominates the para-H 2 case for T > ∼ 50 K. Fig. 5a displays the CO(v 1 = 1) rovibrational deexcitation rate coefficients from the current 6D/CC calculations for collisions of ortho-H 2 for j 2 =1 and 3, separately. These rate coefficients are summed over all j 1 and j 2 =1, 3, and 5. Assuming a Boltzmann average at the kinetic temperature T of the H 2 rotational levels j 2 = 1 and 3, presumed to be representative of the experimental conditions, rate coefficients are computed and found to be in good agreement with the measurements above 200 K. From   Fig. 5a, the contribution from j 2 =3 is seen to be only important above 150 K. The remaining difference with experiment at low temperatures may be due to the fact mentioned above that the ortho-H 2 rate coefficients are not directly measured, but deduced from normal-H 2 and para-H 2 experiments. In particular, Reid et al. assume that the ortho/para ratio in the normal-H 2 measurements is 3:1, i.e. statistical, and independent of temperature. Further, as stated above, the experimental CO rotational population distribution in v 1 = 1 is unknown.
Nevertheless, the current 6D/CC computations are a significant advance over the 4D results of Reid et al. which also correspond to a Boltzmann average of rate coefficients for j 2 =1 and 3.
As indicated above, the situation for para-H 2 collisions is more complicated due to the quasi-resonant contribution (2), a mechanism not important for ortho-H 2 . Boltzmannaveraged rate coefficients are presented in Fig. 5b including j 2 =0 and 2 summed over j 2 = 0, 2, and 4, with and without the quasi-resonant contribution, j 2 = 2 → j 2 = 6. While the current 6D/CC results and the 4D calculations of Reid et al. are in agreement that the quasi-resonant contribution becomes important for T > ∼ 50 K, the relative magnitude compared to the non-resonant transitions from the 6D/CC calculation is somewhat less than obtained previously with the 4D potential. This is partly related to the fact that the 6D/CC rate coefficients for j 2 =0 are significantly larger than those of Reid et al. (see also the corresponding cross sections in Fig. 4 and in the Supplementary Fig. S9). Compared to the experiment, we obtain excellent agreement for T < ∼ 150 K, but are somewhat smaller at higher temperatures. This small discrepancy may be related to the unknown CO(v 1 = 1, j 1 ) rotational population in the measurement. Nevertheless, it is only the 6D/CC computations, i.e., dynamics in full dimensionality with full angular momentum coupling, which can reproduce the measurements for both ortho-and para-H 2 . Further, the computations for the quasi-resonant process (2) were the most challenging reported here due to the requirement of a very large basis set (see Supplementary Table S2) which resulted in long computation times, a large number of channels, and usage of significant disk space (∼0.5 TB per partial wave). In total, the cross sections given here consumed >40,000 CPU hours.

DISCUSSION
The current investigation of the CO-H 2 inelastic collision system has been performed with the intent of minimal approximation through the computation of a high-level potential energy surface, robust surface fitting, and full-dimensional inelastic dynamics with full angular-momentum coupling. That is, within this paradigm for studying inelastic dynamics, we have advanced the state-of-the-art for diatom-diatom collisions through this unprecedented series of computations. The approach has been benchmarked against experiment for pure rotational and rovibrational transitions giving the most accurate results to date within the experimental uncertainties and unknowns. The accuracy and long-range behavior of the 6D potential energy surface is found to be comparable to previous, lower-dimensional surfaces. The agreement of the current computation for the CO(v 1 = 1) rovibrational quenching with measurement resolves a long-standing (more than two decades) discrepancy and justifies the requirement of a full-dimensional approach. This methodology can now, though with significant computational cost, be applied to a large range of initial rotational levels for v 1 = 1 and for higher vibrational excitation to compute detailed state-to-state cross sections unobtainable via experiment.
This advance in computational inelastic scattering is particularly timely as ground-based (e.g., the Very Large Telescope (VLT)) observations have focused on CO rovibrational emission/absorption in a variety of astrophysical objects, while related observations are in the planning stages for future space-based telescopes (e.g., NASA's James Webb Space Telescope). In particular, we are in an exciting era of investigation into the properties of protoplanetary disks (PPDs) around young stellar objects [47]. PPDs provide the material for newly-forming stars and fledgling planets. The CO fundamental band (|∆v 1 | = 1) is a tracer of warm gas in the inner regions of PPDs and, with appropriate modeling, gives insight into disk-gas kinematics and disk evolution in that zone where planets are expected to be forming.
A recent survey of 69 PPDs with the VLT [19] detected CO vibrational bands in 77% of the Remarkably, rotational excitation as high as j 1 = 32 was observed. However, the modeling of PPDs, and other astrophysical sources with CO vibrational excitation, is hindered by the lack of rate coefficients due to H 2 collisions. We are now in an excellent position to provide full-dimensional, state-to-state CO-H 2 collisional data which will not only have a profound impact on models characterizing these intriguing environments that give birth to planets, but also aid in critiquing current theories used to describe their evolution.

METHODS
Potential energy computations. The potential energy computations were performed using the explicitly correlated coupled-cluster (CCSD(T)-F12B) method [48,49], as implemented in MOLPRO2010.1 [50]. The cc-pcvqz-f12 orbital basis sets [51] that have been specifically optimized for use with explicitly correlated F12 methods and for core-valence correlation effects have been adopted. Density fitting approximations [49] were used in all explicitly correlated calculations with the AUG-CC-PVQZ/JKFIT and AUG-CC-PWCVQZ-MP auxiliary basis sets [52,53]. The diagonal, fixed amplitude 3C(FIX) ansatz was used, which is orbital invariant, size consistent, and free of geminal basis set superposition error (BSSE) [54,55]. The default CCSD-F12 correlation factor was employed in all calculations and all coupled cluster calculations assume a frozen core (C:1s and O:1s). The counter-poise (CP) correction [56] was employed to reduce BSSE. Even though the explicitly correlated calculations recover a large fraction of the correlation energy, the CP correction is still necessary, mainly to reduce the BSSE of the Hartree-Fock contribution. Benchmark calculations at the CCSD(T)-F12b/cc-pcvqz-f12 level were carried out on selected molecular configurations and results were compared with those from the conventional CCSD(T) method using aug-cc-pV5Z and aug-cc-pV6Z basis sets. Results showed that the CP corrected interaction energy agrees closely with those derived from CCSD(T)/aug-cc-pV6Z.
To construct the potential energy surface (PES), the computations were performed on a six-dimensional (6D) grid using Jacobi coordinates as shown in Supplementary Fig. S1. R is the distance between the center-of-mass of CO and H 2 . r 1 and r 2 are the bond lengths of CO and H 2 , respectively. θ 1 is the angle between r 1 and R, θ 2 the angle between r 2 and R, The PES fit. The CO-H 2 interaction PES has been fitted in 6D using an invariant polynomial method [37,38]. The PES was expanded in the form, c n 1 ···n 6 y n 1 1 y n 6 6 [y n 2 2 y n 3 3 y n 4 4 y n 5 5 + y n 5 2 y n 4 3 y n 3 4 y n 2 5 ], where y i = e −0.5d i is a Morse-type variable. The internuclear distances d i between two atoms are defined as The total power of the polynomial, N = n 1 + n 2 + n 3 + n 4 + n 5 + n 6 , was restricted to 6. Expansion coefficients c n 1 ···n 6 were obtained using weighted least-squares fitting for potential energies up to 10,000 cm Scattering theory and computational details. The quantum scattering theory for a collision of an S-state atom with a rigid-rotor was developed [2,[57][58][59] based on the close-coupling (CC) formulation of Arthurs and Dalgarno [1]. Details about its extension to diatom-diatom collisions with full vibrational motion can be found in Refs. [12,13]. In this approach, the interaction potential V (R, r 1 , r 2 , θ 1 , θ 2 , φ) is expanded as, with the bi-spherical harmonic function expressed as, where 0 ≤ λ 1 ≤ 10, 0 ≤ λ 2 ≤ 6 was used in the scattering calculations. Due to the symmetry of H 2 , only even values of λ 2 contribute.
For convenience, the combined molecular state (CMS) notation is applied to describe a combination of rovibrational states for the two diatoms. A CMS represents a unique quantum state of the diatom -diatom system before or after a collision. The CMS will be denoted as (v 1 j 1 v 2 j 2 ). v and j are the vibrational and rotational quantum numbers.
The rovibrational state-to-state cross section as a function of collision energy E is given by, where (v 1 j 1 v 2 j 2 ) and (v 1 j 1 v 2 j 2 ) are, respectively, the initial and final CMSs of CO-H 2 , the wave vector k 2 = 2µE/h 2 , and S is the scattering matrix. l is the orbital angular momentum and J the total collision system angular momentum, where J = l + j 12 and j 12 = j 1 + j 2 .
Thorough convergence testing was performed in the scattering calculations by varying all relevant parameters. The CC equations were propagated for each value of R from 4 to 18.0 a 0 using the log-derivative matrix propagation method of Johnson [60] and Manolopoulos [61], which was found to converge for a radial step-size of ∆R = 0.05 a 0 . The convergence tests of the v 1 = 1 → 0 vibrational quenching cross section of CO with respect to the number of v 1 = 1 rotational channels found that at least 13-15 channels have to be included in the v 1 = 1 basis set, especially for low-energy scattering. Based on convergence tests with respect to the adopted maximum R for the long range part of the PES, we found that the cross sections are converged down to the lowest collision energy of 0.1 cm −1 . This value also guarantees that the rate coefficients are converged for temperatures greater than 1 K. The number of discrete variable representation points N r 1 and N r 2 ; the number of points in θ 1 and θ 2 for Gauss-Legendre quadrature, N θ 1 and N θ 2 ; and the number of points in φ for Gauss-Hermite quadrature, N φ , which were applied to project out the potential expansion coefficients were all tested for convergence with the final adopted values given in Supplementary Table S2. The basis sets and the maximum number of coupled channels are also presented in Supplementary Table S2.
The resulting integral cross sections were thermally averaged over a Maxwellian kinetic energy distribution to yield state-to-state rate coefficients as function of temperature T , where m is the reduced mass of the CO-H 2 complex, β = (k B T ) −1 , and k B is Boltzmann's constant.
where σ is the cross section and N the number of data points. The RMS errors for the 4D version of V6D with EE and VV bond lengths, the 4D V12 PES, and the V6D PES are shown in Supplementary Table 1 with the cross sections for the latter two shown in Fig. 2 of the main text. The current V6D results for the final CO j 1 distributions (Fig. 1) are seen to be in excellent agreement with experiment, validating the accuracy of the fitted potential as well as the present scattering calculations.
To compare with the measurements, the computed cross sections were convolved with the experimental kinetic energy distribution according to, where δ is given by [S7], The detailed cross sections and discussion are presented in the main text. [S10] used a mixed IOS-CS formalism in which the rotational motion of CO was neglected, but that of H 2 retained, allowing for significant reduction of basis sets. The cross sections from Flower [S11] are likely too large due to an insufficient CO basis, an effect studied by Baĉić et al. [S8]. In fact, using Flower's basis set in CC calculations on V6D, we obtained significantly larger cross sections than our converged results, verifying the primary cause for the discrepancy. However, in most cases differences are related to the adopted PESs and dimensionality of the scattering approaches.

Quenching rate coefficients
In the calculations of Baĉić et al. [S8, S9] and Reid et al. [S10], the total quenching cross section from state (v 1 = 1, j 2 ) was obtained by summing the state-to-state quenching cross sections over the final rotational state of H 2 , j 2 , allowing for H 2 rotational inelastic transitions (with v 2 = v 2 = 0 and no information on j 1 or j 1 as the IOS approximation was adopted for CO in both previous calculations.) By thermally averaging the cross sections over the collision energy, the quenching rate coefficients from v 1 = 1, j 2 were obtained.
The total thermal rate coefficients for CO v 1 = 1 → v 1 = 0 vibrational quenching were obtained by summing k(v 1 = 1, j 2 → v 1 = 0)(T ) over all initial H 2 rotational states j 2 weighted by their populations assuming a Boltzmann distribution at T . As observed by Millikan and Osburg [S12], para-H 2 is a more efficient collision partner than ortho-H 2 for vibrational quenching of CO due to the quasi-resonant transition: CO(v 1 = 1) + H 2 (v 2 = 0, j 2 = 2) → CO(v 1 = 0) + H 2 (v 2 = 0, j 2 = 6). Guided by the measurements, this quasiresonant process was also theoretically investigated by Baĉić Supplementary Fig. S11 gives a similar comparison for ortho-H 2 collisions though a quasiresonant process is not operable due to much larger asymptotic energy differences. Otherwise, the trends are very similar to those noted for para-H 2 collisions, though the current V6D/CC results fall somewhat between the earlier two 4D calculations. Note that we did not consider the H 2 inelastic channel j 2 = 1 → 5 as Bacîć et al. [S9] find the cross sections to be more than an order of magnitude smaller than the j 2 = 1 → 3 transition. This is in contrast to Reid et al. who find the j 2 = 1 → 5 transition to dominate the quenching of (1001). As a consequence, the thermally-averaged rate coefficients for the quenching by to draw conclusions from one-to-one comparisons of the various calculations. However, the main deficiency in the IOS approach is that it does not resolve the final CO(j 1 ) channels