An Efficient Statistical Method to Compute Molecular Collisional Rate Coefficients

A general property of astronomical molecular spectra is that the populations of the energy levels are rarely at local thermodynamical equilibrium (LTE). Indeed, in space, the density is usually such that the frequency of collisions is neither negligible nor large enough to maintain LTE. Deviations from LTE, including strong deviations like population inversions, are thus the rule rather than the exception. In such conditions, interpreting spectra requires one to solve simultaneously the radiative transfer equation and a set of statistical equilibrium equations for the molecular energy levels. Solving the statistical equilibrium in turn necessitates the availability of the energy levels, the state-to-state rates for spontaneous radiative decay, and the state-to-state rate coefficients for collisional (de)excitation. The collisional (inelastic) rate coefficients are extremely difficult to measure in the laboratory and radiative transfer models rely almost exclusively on theoretical estimates.

Nowadays theoretical data can be obtained for neutral—closed shell—diatomic and polyatomic molecules colliding with atoms or light diatomic perturbers like H₂ at the full quantum time-independent close-coupling (CC) level of theory. In addition, accuracy can be inferred by comparing calculations directly with state-to-state experimental data thanks to recent developments in double resonance (Mertens et al. 2017) and crossed-beam techniques (Chefdeville et al. 2015; Vogels et al. 2015). Such comparisons have demonstrated the high accuracy of modern computations based on state-of-the-art potential energy surfaces. For (polyatomic) reactive radicals and ions or polyatomic projectiles, however, the numerical solution of the coupled equations becomes impractical due to the memory and CPU requirements.⁴ The problem is particularly acute for systems with deep potential wells because an excessively large number of closed channels is needed to achieve convergence. As a result, for collisions where the intermediate complex is a stable molecule or ion, CC calculations have been restricted to the very simplest atom–diatom systems; see, e.g., Stoecklin et al. (2015) and Bulut et al. (2015) for OH⁺+H.

Various approximative methods, however, exist to treat the most demanding systems. A standard approach to reduce the number of coupled channels is to use angular momentum decoupling approximations such as the coupled-states and the infinite-order-sudden methods (Kouri 1979). Such approximations, however, become reliable only at high collision energies with respect to the potential well and/or to the rotational spacings, which in practice restricts their applicability to temperatures above, typically, 300 K. Quantum time-dependent (wavepacket) methods suffer similar limitations, although recent progress has been made with the Multi Configuration Time Dependent Hartree method (Ndengué et al. 2017). Semi-classical and quasi-classical methods offer another alternative, their common feature being that at least one degree of freedom in the collisional system is treated classically. The good agreement observed between such methods and full quantum calculations is encouraging (Faure et al. 2016; Semenov & Babikov 2017), but there are still intrinsic limitations, especially when resonances or interferences are important.

In the present Letter, we explore another class of approximation to treat purely inelastic collisions: the statistical method. The statistical approach was initially developed for chemical reactions (see Pechukas et al. 1966; González-Lezana 2007; Park & Light 2007; Dagdigian & Alexander 2018 and references therein). Statistical models are expected to be accurate for reactive or inelastic collisions proceeding through formation and decay of a strongly bound collision complex. Although the basic hypotheses of statistical theories are similar, the different models introduce various dynamical constraints. Recent examples include the treatment of the reactive collision ${{\rm{H}}}_{3}^{+}+{{\rm{H}}}_{2}$ by Park & Light (2007), later extended by Gómez-Carrasco et al. (2012), or the study of rotationally inelastic collisions of CH with H₂ (Dagdigian 2016) and H (Dagdigian 2017) using the quantum statistical method of Rackham et al. (2001). Here, we employ a simpler approach inspired by the statistical adiabatic channel theory of Quack & Troe (1974, 1975). The method is checked against full quantum time-independent calculations for a series of benchmark systems, including ionic and neutral targets, and for a large range of temperatures. In Section 2, we describe the statistical approach and the determination of collisional rate coefficients. In Section 3, we compare the results obtained with this method to exact quantum calculations. A discussion on the possible use of our model concludes this Letter.

The CC method is based on an expansion of the total wavefunction Ψ into an angular basis set $\{| \alpha \rangle \}$:

$$\begin{eqnarray}&&{\rm{\Psi }}={R}^{-1}\displaystyle \sum _{\alpha }| \alpha \rangle {\chi }_{\alpha }(R),\end{eqnarray} \tag{ 1 }$$

where ${\chi }_{\alpha }(R)$ are radial functions describing the nuclear motion and the form of the angular functions depend on the type of collision (Flower 2007). Integrating the Schrödinger equation over the angular variables leads to the set of coupled differential radial equations,

$$\begin{eqnarray}&&-\displaystyle \frac{{{\hslash }}^{2}}{2\mu }\displaystyle \frac{{d}^{2}{\chi }_{\alpha }}{{{dR}}^{2}}+\displaystyle \sum _{{\alpha }^{{\prime} }}\left\langle {\alpha }^{{\prime} }\left|{H}_{\mathrm{int}}+V+\displaystyle \frac{{{\boldsymbol{L}}}^{2}}{2\mu {R}^{2}}\right|\alpha \right\rangle {\chi }_{\alpha }=E{\chi }_{\alpha },\end{eqnarray} \tag{ 2 }$$

to be solved for a given collision energy E and for each value of the total angular momentum J. In this equation, ${H}_{\mathrm{int}}$ is the Hamiltonian describing the internal ro-vibrational motion of the molecule(s), V denotes the potential energy surface of the system, and ${\boldsymbol{L}}$ is the angular momentum describing the relative motion of the two colliding partners. The total angular momentum ${\boldsymbol{J}}={{\boldsymbol{j}}}_{1}+{{\boldsymbol{j}}}_{2}+{\boldsymbol{L}}$, where ${{\boldsymbol{j}}}_{1}$ and ${{\boldsymbol{j}}}_{2}$ are the angular momenta of the colliding molecules (${{\boldsymbol{j}}}_{2}=0$ for atom–molecule collisions), can be used to conveniently reformulate the coupled equations. The total angular momentum is conserved during a collision, and the coupled equations become block-diagonal in J. Solving these equations with appropriate boundary conditions gives access to the collision matrix $S(E,J)$ and the inelastic cross sections are obtained by summing the contributions of all values of J until convergence is reached.

In our approach, we diagonalize the Hamiltonian excluding the nuclear kinetic term (i.e., only the second term of Equation (2)) for each value of J, which results in a set of adiabatic curves. It is important to note that the present approach requires an accurate ab initio potential energy surface for the collisional complex. Each adiabatic curve can be associated asymptotically with an internal state of the colliding molecules, and the number of adiabatic curves is determined by the values that the quantum number L can take according to the standard angular momentum coupling rules. Examples of adiabatic channels are represented in Figure 1. The basic assumption is that if the kinetic energy is larger than the centrifugal barrier in the entrance channel, an intermediate molecular complex can be formed and an inelastic collision can take place. The probability of inelastic transition is determined solely by the number $N(E,J)$ of open exit channels based on an analysis of the adiabatic curves. According to statistical theory, all open channels are assigned equal weight, and the collision ${\boldsymbol{S}}$-matrix elements are given by $| {S}_{{if}}(E,J){| }^{2}\,=1/N(E,J)$ for open channels, and zero otherwise. We note that in the case of reactive (exothermic) systems, the present approach can be easily extended to include the reactive channels in the counting of open channels. The cross section is calculated as

$$\begin{eqnarray}&&{\sigma }_{{if}}(E)=\displaystyle \frac{\pi }{(2{j}_{1}+1)(2{j}_{2}+1){k}^{2}}\displaystyle \sum _{J=0}^{\infty }(2J+1)| {S}_{{if}}(E,J){| }^{2}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

The rate coefficients are then computed by integrating the cross section weighted by a Maxwell–Boltzmann distribution of energies,

$$\begin{eqnarray}&&{k}_{{if}}(T)={\left(\displaystyle \frac{2}{{k}_{{\rm{B}}}T}\right)}^{3/2}\displaystyle \frac{1}{\sqrt{\pi \mu }}{\int }_{0}^{\infty }{{Ee}}^{-E/{k}_{{\rm{B}}}T}{\sigma }_{{if}}(E)\ {dE}.\end{eqnarray} \tag{ 4 }$$

Note that the present approach allows to save an enormous amount of both CPU time and memory compared to standard CC calculations as the adiabatic curves are independent of the collision energy. Moreover, the convergence of the adiabatic curves with respect to the ro-vibrational basis is much faster than the convergence of CC cross sections so that it is not necessary to include closed channels in the basis, even for strongly bound collisional systems. We also note that the present approach can be used for open-shell (polyatomic) molecules showing a complex rotational structure.

The method is expected to yield the best results for collisions in which an intermediate molecular complex can be formed. As an example, we consider rotationally inelastic OH⁺–H collisions, for which accurate CC results are available in the literature (Bulut et al. 2015; Stoecklin et al. 2015). In this application, we only consider rigid-rotor collisions on the quartet potential energy surface, thus neglecting the (low probability) reactive channels (Bulut et al. 2015). The intermediate complex, H₂O⁺, has a large depth (about 0.5 eV or 4000 cm⁻¹) that should favor a statistical approach. Using the methodology described above, we have determined the OH⁺–H rotational rate coefficients in the temperature range 10–1000 K. Results are presented in Figure 2 for four sample temperatures between 50 and 1000 K and they are compared to the (rigid-rotor) CC results of Bulut et al. (2015). At low temperatures (below 300 K), we observe that the present statistical method predicts the rate coefficients with a very good accuracy over many orders of magnitude. Up to 300 K, most rate coefficients are accurate to within 50%, which is satisfactory for most astrophysical applications. We note, however, that larger discrepancies appear as the temperature increases for the dominant transitions. This is expected since at high energies, the collision time becomes too short to assume formation of a long-lived complex, so that the hypothesis of an equiprobable distribution of exit channels is not guaranteed. However, differences remain below a factor of ∼3 up to 1000 K, which is still quite reasonable. In Figure 3, we compare the rate coefficients for the initial state ${j}_{1}=6$ at a temperature of 200 K. We observe that the adiabatic model slightly underestimates the transitions with ${\rm{\Delta }}{j}_{1}=\pm 1$ compared to the CC results, but all other transitions and propensities are nicely reproduced. The same behavior is observed at other temperatures. We note that strong resonance and/or interference effects would necessarily be missed in the statistical approach.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the discrepancies between the CC results and the statistical results, we compute the weighted mean error factor (WMEF) for each temperature, defined as

$$\begin{eqnarray}&&\mathrm{WMEF}=\displaystyle \frac{{\displaystyle \sum }_{i}{k}_{i}^{\mathrm{CC}}{r}_{i}}{{\displaystyle \sum }_{i}{k}_{i}^{\mathrm{CC}}},\end{eqnarray} \tag{ 5 }$$

where ${k}_{i}^{\mathrm{CC}}$ is the rate coefficient for transition i obtained with the CC method and r_i is the ratio of the rate coefficients obtained with the two methods, defined as ${r}_{i}=\max ({k}_{i}^{\mathrm{CC}}/{k}_{i}^{\mathrm{adia}},{k}_{i}^{\mathrm{adia}}/{k}_{i}^{\mathrm{CC}})$ so that ${r}_{i}\geqslant 1$. We use a weighted average to underline the fact that for radiative transfer applications, the largest rate coefficients are the most important. It is clear that the method is very efficient to describe collisions in strongly bound systems such as OH⁺–H with a WMEF of only 1.18 at 50 K. At such low temperature, the approximate wavepacket and semi- or quasi-classical approaches would fail by much larger factors. The WMEF increased to 2.20 at 1000 K, which is still reasonable. These good results confirm that the statistical method is well adapted to treat systems with "covalent" interactions, as expected.

In many astrophysical environments, such as the cold interstellar medium, the most abundant colliding partners are not hydrogen atoms but H₂ and He. In this case, the corresponding intermediate complexes may have much shallower well depth, and one might wonder whether the statistical method is still applicable. We first consider collisions of molecular ions with He or H₂. These collisions provide an intermediate case as a non-covalent but strong van der Waals interaction occurs due to the charge-induced dipole potential. Figure 4 illustrates the results for a cationic system, OH⁺–He, and an anionic system, CN⁻–H₂. The well depths of the potential energy surfaces are very similar, 730 cm⁻¹ and 750 cm⁻¹, respectively. The rotational rate coefficients have been calculated with the CC method for OH⁺–He (Gómez-Carrasco et al. 2014) and the coupled-states method for CN⁻–H₂ (Kłos & Lique 2011). For this latter system, CC calculations have shown that the coupled-states approach is accurate to better than 50%. In the case of OH⁺–He, we observe a reasonable agreement (better than a factor 3 for most transitions) for the three temperatures considered, although at 300 K the dominant rate coefficients are not well reproduced. For CN⁻–H₂, the agreement is not that good, which might reflect the fact that H₂ has a large rotational constant and does not behave statistically. This is particularly the case for transitions with small rate coefficients ($\lt {10}^{-11}$ cm³ s⁻¹). On the other hand, the dominant transitions are rather accurately reproduced at all temperatures. As for OH⁺–H collisions, we note that the agreement between the statistical and the coupled-channel results decreases with increasing temperature. We conclude that the statistical approach can be useful also for (non-covalent) strongly bound van der Waals systems, provided that low temperatures are considered.

Zoom In Zoom Out Reset image size

Finally, we consider collisions in which the intermediate complex is bound by weak van der Waals forces. Due to its large abundance in many astrophysical environments, we used CO in collision with He and H₂, for which CC rate coefficients are available (Cecchi-Pestellini et al. 2002; Yang et al. 2010). In both cases, the depth of the potential energy surface is very small, 24 and 94 cm⁻¹, respectively. Quite unexpectedly, the statistical method still performs relatively well, with most dominant transitions (i.e., rate coefficients above 10⁻¹¹ cm³ s⁻¹) being reproduced with an error less than a factor of 3 and WMEF in the range 1.5–3, as illustrated in Figure 5. On the other hand, the scatter of points is much larger than for the previous systems, as expected from the much smaller well depths. In addition, for CO–H₂ the WMEF decreases with increasing temperature, which might suggest the method is out of its applicability domain.

Zoom In Zoom Out Reset image size

In summary, the statistical method is found to provide average accuracies better than 50% for strongly bound (covalent) collisional systems up to room temperature. At higher temperature, or for systems with weaker interactions, larger inaccuracies are observed, as expected from the basic assumption of the method.

In this Letter, we have presented a new statistical method to compute inelastic rate coefficients for interstellar molecules. The method has been benchmarked versus "exact" coupled-channel calculations for five collisional systems (OH⁺–H, OH⁺–He, CN⁻–H₂, CO–He, and CO–H₂). These five systems include a strongly bound molecular system, two strongly bound van der Waals systems, and two weakly bound systems. For the strongly bound case, we found the statistical results to be in excellent agreement with the coupled-channel results, especially at low temperatures where the lifetime of the collision complex is the longest. When the temperature increases above room temperature, the agreement decreases as expected. For the two strongly bound van der Waals systems, the agreement is still (perhaps surprisingly) rather satisfactory. For the two weakly bound systems, the accuracy decreases significantly, as could be anticipated.

Whereas accurate collisional rate coefficients are now available for many non-reactive molecules (Roueff & Lique 2013), very few such reliable data exist for inelastic collisions involving reactive radicals such as NH or reactive ions such as ${{\rm{H}}}_{3}^{+}$, H₂O⁺, and H₃O⁺, for which reactive and inelastic channels are competing. In order to provide the astrophysical community with collisional data for these systems, the helium atom is generally employed as a substitute for H or H₂ for simplicity. However, it is now well established that significant differences exist between He, H, and H₂ collisional rate coefficients (Roueff & Lique 2013). Collisions with He atoms can thus provide guidance at the order-of-magnitude level but they can certainly not reproduce H or H₂ to within 50% accuracy. The present approach leads to much more accurate results for inelastic collisions and should be considered as a useful alternative to prohibitive CC calculations. Future applications will include exothermic reactive collisions such as the triatomic CH⁺+H system for which full-dimensional CC calculations are available (Werfelli et al. 2015).

Finally, we note that in the coma of comets and in protoplanetary disks, new partners enter the collision game. Thus, in comets, molecules are excited by H₂O, which is the most abundant neutral species in the coma. The computation of rates for the excitation of molecules due to H₂O collisions is not yet numerically accessible through quantum coupled-channel calculations due to the large potential wells (Sánchez & Dubernet 2017; Kalugina et al. 2018). The present approach has no such CPU limit and as soon as a potential energy surface is available, collisional data can be easily predicted. As a result, the present statistical method will allow us to study the collisional excitation of cometary molecules such as CO, H₂O, or HCN by H₂O and to provide the first collisional rate coefficients for H₂O–molecule systems.

We thank T. Stoecklin, M. Alexander, and P. Dagdigian for useful discussions. This research was supported by the French National Research Agency (ANR) through a grant to the Anion Cos Chem project (ANR-14-CE33-0013). We also acknowledge the financial support of the Region Normandie through the Actions Thématiques Stratégiques project managed by Normandie University.