Rate Constants for Fine-structure Excitations in O–H Collisions with Error Bars Obtained by Machine Learning

Of the elements heavier than helium, oxygen is the most abundant in the interstellar medium. Interstellar oxygen is predominantly in the form O i as it is shielded from ionizing radiation by the more abundant hydrogen atoms. The fine-structure transitions in collisions of O i with atomic hydrogen,

$$\begin{eqnarray}&&{\rm{O}}{(}^{3}{P}_{j})+{\rm{H}}{(}^{2}S)\to {\rm{O}}{(}^{3}{P}_{j^{\prime} })+{\rm{H}}{(}^{2}S),\end{eqnarray} \tag{ 1 }$$

are hence an important cooling process in the interstellar medium (Dalgarno & McCray 1972; Cowie & Songaila 1986; Pequignot 1990; Shaw et al. 2006; Glover & Clark 2013) and may be a dominant mechanism at low temperatures (Dalgarno & McCray 1972).

The relative intensities of the O i spectral lines serve as an important diagnostic probe of the atomic environment (Cowie & Songaila 1986; Jenkins & Tripp 2011). The absorption lines can be used as a source of information on the column densities for the ${}^{3}{P}_{j}$ levels, and the strong emission lines in diffuse regions can be used to study a diverse range of nebulae (Pequignot 1990). Accurate rate constants for fine-structure transitions (1) are needed for understanding cooling in early galaxies (Greif et al. 2010), atomic-line diagnostics and gas modeling in protoplanetary disks (Meijerink et al. 2008; Tilling et al. 2012), and the spectroscopy of stellar atmospheres (Fabbian et al. 2009). The rates of the reaction (1) have been used in models of turbulence in the interstellar medium (Lesaffre et al. 2013; Guillard et al. 2015) and studies of both the ionized and neutral parts of the warm interstellar medium (Dong & Draine 2011; Jenkins 2013).

Most of these models rely on the theoretical calculations of collision-induced fine-structure transition rates by Launay & Roueff (1977) and, more recently, by Abrahamsson et al. (2007). The interaction of O(³P) with H is governed by four molecular potentials. Abrahamsson et al. (2007) used the interaction potentials computed by Parlant & Yarkony (1999). However, the potentials used by Abrahamsson et al. (2007) were originally constructed for the scattering calculations of the collision rates for the electronic excitations O(³P) $\to$ O(¹D). As we will show below, adapting the potentials to the collision problem (1) resulted in an error, rendering the calculations of Abrahamsson et al. (2007) inaccurate. In addition, the interaction potentials of Parlant & Yarkony (1999), as with any intermolecular potentials, are subject to uncertainties associated with errors of the ab initio calculations. It is important to evaluate the effect of these uncertainties on the rate constants. However, because the collision problem is parametrized by four potentials and it is not known which of these potentials play a dominant role in determining the fine-structure transitions, estimating the errors of the rate constants stemming from uncertainties in the potentials presents a challenge.

In the present paper, we develop an approach using a combination of coupled channel scattering calculations with a machine-learning technique based on Gaussian Process regression to determine the sensitivity of the rate constants to variations of the four interaction potentials. Using this approach, we repeat the calculations of Abrahamsson et al. (2007) to report accurate values of the rate constants based on the potentials of Parlant & Yarkony (1999) and compute the error bars of the rate constants corresponding to reasonable variation of the ab initio potentials. We show that this method can be used to determine the relative importance of the individual potentials in determining the outcome of different fine-structure-changing collisions.

The collisions of O(³P) with H(²S) are determined by four interaction potentials of the following symmetries: X²Π, ²Σ, ⁴Π, and ⁴Σ. These potentials have been computed most recently by Parlant & Yarkony (1999). Figure 1 shows the calculations of Parlant & Yarkony (1999) as symbols. Instead of using the fixed interaction potentials for the quantum scattering calculations as was done by all previous authors, in the present work we ask the following question: if each of the potentials is uncertain to within the gray band shown in Figure 1, what is the corresponding uncertainty in the collision rates?

Zoom In Zoom Out Reset image size

This question could be answered by a quantum scattering calculation on a grid of parameters determining the variation of each of the four potentials. However, the number of such scattering calculations would have to be unfeasibly large. For example, if the variation of each potential was determined by n parameters, the rate constants would be $4n$-dimensional functions, requiring a total of $\sim {10}^{4n}$ scattering calculations of collision rates at each temperature in order to determine the global response of rates to the variation of the potentials. Instead, we propose a machine-learning approach based on training a statistical learning method with the results of quantum scattering calculations, reducing the required number of computations to a linear function of the number of parameters.

In order to account for the variation of the potentials, we first produce the analytical fits of the numerical data of Parlant & Yarkony (1999), shown in Figure 1 as solid curves. We then modulate each of the analytical fits of the repulsive potentials for the excited molecular states ²Σ, ⁴Π, and ⁴Σ by a scaling factor varied from 0.8 to 1.2. The interaction potential of the ground state X²Π is modulated by means of two parameters: the overall scaling factor inducing a 20% change in the magnitude of the potential and a factor translating the interatomic distance by ±10%, leading to shifts in the position of the potential minimum, a change in the slope of the repulsive part of the interaction potential, and a change of the classical turning point in the limit of zero collision energy. We thus parametrize the collision problem by five parameters ${\boldsymbol{x}}=\{{x}_{1},...,{x}_{5}\}$.

We use the Latin Hypercube Sampling (LHS) scheme to generate Q sets of these parameters at random (Stein 1987; McKay et al. 1979). LHS ensures that these sets sample the five-dimensional space evenly. For each set of the five parameters ${\boldsymbol{x}}$, we compute the scattering rates for the fine-structure transitions (1). The computations are performed using the rigorous coupled channel method described in multiple previous papers, including the work by Launay & Roueff (1977) and Launay (1977). The rate constants are obtained by integrating the energy dependence of the scattering cross sections with the Maxwell–Boltzmann distribution. The cross sections for a fixed collision energy are computed by integrating a system of coupled channel equations presented by Launay (1977). As recommended by NIST,¹ we adapt the energies of 158 cm⁻¹ for O(${}^{3}{P}_{j=1}$) and 227 cm⁻¹ for O(${}^{3}{P}_{j=0}$) relative to the energy of O(${}^{3}{P}_{j=2}$).

The coupled channel calculations produce Q values of the rate constants corresponding to Q sets of the parameters ${\boldsymbol{x}}$ at each temperature. We treat each fine-structure transition separately and compute only the rate constants for the collision-induced excitations: $j=2\to j^{\prime} =1$, $j=2\to j^{\prime} =0$ and $j=1\to j^{\prime} =0$. These rate constants are used as input for Gaussian Process (GP) regression.

GP regression is an efficient, statistical learning technique for interpolation in multidimensional spaces (Cressie 1993; Neal 1998, p. 6; Williams 1998, p. 599; Papritz & Stein 2002; Stein 1999; Quinonero-Candela & Rasmussen 2005; Rasmussen & Williams 2006). The application of GP regressions to molecular physics problems has been described in our recent work (Cui & Krems 2015, 2016; Cui et al. 2015), where more details on the implementation of the method can be found. In brief, the interpolation starts by estimating the statistical correlations between rate constants corresponding to different points ${\boldsymbol{x}}$ in the five-dimensional space considered here. We approximate the correlation function as

$$\begin{eqnarray}R({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} ) & = & \left\{\displaystyle \prod _{i=1}^{q=5}\left(1+\displaystyle \frac{\sqrt{5}| {x}_{i}-{x}_{i}^{\prime }| }{{\omega }_{i}}+\displaystyle \frac{5{({x}_{i}-{x}_{i}^{\prime })}^{2}}{3{\omega }_{i}^{2}}\right)\right.\\ & & \left.\times \exp \left(-\displaystyle \frac{\sqrt{5}| {x}_{i}-{x}_{i}^{\prime }| }{{\omega }_{i}}\right)\right\},\end{eqnarray} \tag{ 2 }$$

where ${\omega }_{i}$ are the unknown parameters representing the characteristic length scales of the correlations. Note that the specific choice of the correlation function in Equation (2) is not unique. The mathematical form of the correlation function determines the efficiency of the GP regression, i.e., the number of rate constant calculations needed to obtain the global five-dimensional dependence on the parameters. However, given any physically reasonable mathematical form of R, GP regression should become accurate, when based on a large enough number of rate constant calculations.

For a set of Q points in the five-dimensional space, we construct the correlation matrix ${\boldsymbol{A}}$ whose elements are given by $R({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} )$, with ${\boldsymbol{x}}$ and ${\boldsymbol{x}}^{\prime}$ covering all pairs of rate constants in the computed set. To find the estimates of the parameters ${\boldsymbol{\omega }}={({\omega }_{1},{\omega }_{2},\cdots ,{\omega }_{5})}^{\top }$, we maximize the log-likelihood function

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }({\boldsymbol{\omega }}| {\boldsymbol{Y}})=-\displaystyle \frac{1}{2}[Q\mathrm{log}{\hat{\sigma }}^{2}+\mathrm{log}(\det ({\boldsymbol{A}}))+Q]\end{eqnarray} \tag{ 3 }$$

by iteratively computing the determinant and the inverse of the matrix ${\boldsymbol{A}}$. In Equation (3),

$$\begin{eqnarray}&&{\hat{\sigma }}^{2}({\boldsymbol{\omega }})=\displaystyle \frac{1}{Q}{({\boldsymbol{Y}}-{\boldsymbol{\beta }})}^{\top }{{\boldsymbol{A}}}^{-1}({\boldsymbol{Y}}-{\boldsymbol{\beta }}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\hat{\beta }({\boldsymbol{\omega }})={({{\boldsymbol{I}}}^{\top }{{\boldsymbol{A}}}^{-1}{\boldsymbol{I}})}^{-1}{{\boldsymbol{I}}}^{\top }{{\boldsymbol{A}}}^{-1}\,{\boldsymbol{Y}},\end{eqnarray} \tag{ 5 }$$

Y is the vector of all computed rate constants for a given fine-structure transition, ${\boldsymbol{\beta }}=\beta {\boldsymbol{I}}$, ${\boldsymbol{I}}$ is the identity vector, and the hat over the symbol denotes the maximum likelihood estimators (i.e., expressions giving the maximum likelihood estimates).

Given the best values of ${\boldsymbol{\omega }}={({\omega }_{1},{\omega }_{2},\cdots ,{\omega }_{q})}^{\top }$ thus found, we can construct the final correlation matrix ${\boldsymbol{A}}$. The prediction of the rate constant at any point x₀ of the five-dimensional space ${\boldsymbol{x}}$ can then be expressed in terms of this matrix and a vector

$$\begin{eqnarray}&&{{\boldsymbol{A}}}_{0}=\left[\begin{array}{c}R({{\boldsymbol{x}}}_{0},{{\boldsymbol{x}}}_{1})\\ R({{\boldsymbol{x}}}_{0},{{\boldsymbol{x}}}_{2})\\ \vdots \\ R({{\boldsymbol{x}}}_{0},{{\boldsymbol{x}}}_{Q})\end{array}\right]\end{eqnarray} \tag{ 6 }$$

as follows:

$$\begin{eqnarray}&&k({{\boldsymbol{x}}}_{0})=\beta +{{\boldsymbol{A}}}_{0}^{\top }{{\boldsymbol{A}}}^{-1}({\boldsymbol{Y}}-{\boldsymbol{\beta }}).\end{eqnarray} \tag{ 7 }$$

The vector A₀ describes the statistical relationship between the rate constant at the point of interest x₀ and all other computed rate constants assembled in the vector ${\boldsymbol{Y}}$. Equation (7) gives the five-dimensional dependence of the rates on the parameters modulating the OH interaction potentials.

Given Equation (7), we can compute the rates corresponding to the original potentials of Parlant & Yarkony (1999), calculate the error bars arising from the variation of the interaction potentials, and examine the sensitivity of the rates to each of the five different parameters. To do that, we apply a variance-based sensitivity analysis (Mcrae et al. 1982; Saltelli & Bolado 1998; Saltelli et al. 1999; Helton & Davis 2002; Cacuci & Ionescu-Bujor 2004). As illustrated in the following section, this analysis can be a powerful tool to explore the mechanisms of inelastic scattering by providing information on which parts of the underlying interaction potentials are predominantly responsible for the collisional energy transfer.

We begin by illustrating the accuracy and efficiency of GP regression for the prediction of rate constants for reaction (1). Figure 2 presents a sample of 200 rate constants for the $j=1\to j^{\prime} =0$ transition at T = 300 K randomly selected in the five-dimensional space of the potential parameters. For each point, we compute the rate constant with the coupled channel method and with Equation (7) based on three different numbers of training points: Q = 10, Q = 100, and Q = 1000. The results illustrate that 1000 scattering calculations produce an accurate GP model giving the five-dimensional dependence of the rate constants on the potential parameters. Note that the efficiency of GP regression is known to improve with increasing dimensionality of the problem, i.e., the total number of points divided by the number of dimensions should decrease with increasing dimensionality (Deisenroth 2010). The results of Figure 2 thus indicate that it is feasible to construct a GP model of rate constants parametrized by as many as 50 parameters. Constructing such a model would require an iterative computation of the inverse of the correlation matrix of size less than 10,000 × 10,000. For the following calculations, we use the GP model trained by Q = 2000 scattering calculations. This number should be contrasted with 10⁵, which would be the number of scattering calculations to construct the global dependence of the rate constants on the five potential parameters, if the calculations were performed on a grid with 10 points per parameter.

Zoom In Zoom Out Reset image size

Figure 3 presents the temperature dependence of the rate constants for fine-structure excitations with the error bars corresponding to the variations of the interaction potentials described in the previous sections and depicted in Figure 1. These temperature dependencies are obtained from the GP model. As an additional check of the accuracy of these results, we computed the rate constants with the coupled channel approach using the potentials of Parlant & Yarkony (1999). The results, shown by blue triangles in Figure 3, are graphically indistinguishable from the means of the uncertainty intervals represented by the red solid lines.

Zoom In Zoom Out Reset image size

Figure 3 also includes the data of Abrahamsson et al. (2007) represented by the line with green squares. As can be seen, their data underestimate or overestimate the current calculation by almost a factor of 2 at low temperatures and ∼10%–30% at elevated temperatures. To illustrate the discrepancy between the previous calculations and the current results, we plot the difference in the rate constants versus temperature in Figure 4. We have traced the origin of the discrepancy to the incorrect long-range behaviour of the interaction potentials used by Abrahamsson et al. (2007) resulting from an extrapolation of the data of Parlant & Yarkony (1999). The ²Σ, ⁴Π, and ⁴Σ potentials of Abrahamsson et al. (2007) become degenerate at interatomic separations $R\gt 6$ Bohr, as a result of matching the potentials to the long-range analytical form with the same dispersion coefficient. As illustrated in the inset of Figure 1, the three degenerate excited-state potentials used by Abrahamsson et al. (2007) have an abnormally large van der Waals minimum at long range. Recently, Dagdigian et al. (2016) demonstrated that these molecular potentials in fact have a very shallow long-range minimum, not exceeding 27 cm⁻¹.

Zoom In Zoom Out Reset image size

For future reference, we present the numerical values of the rate constants for the fine-structure excitation transitions in Tables 1–3. The tables include the data obtained by Launay & Roueff (1977) and by Abrahamsson et al. (2007), as well as the uncertainty intervals shown in Figure 3.

Since GP regression produces a global dependence of rates on the potential energy parameters, it can be used to explore the sensitivity of the rate constants to different interaction potentials. To do this, we perform the variance-based sensitivity analysis as described by multiple authors (Mcrae et al. 1982; Saltelli & Bolado 1998; Saltelli et al. 1999; Helton & Davis 2002; Cacuci & Ionescu-Bujor 2004). This analysis produces five coefficients that quantify the contribution of each underlying parameter to the variance of the rate constants. Figure 5 presents the temperature dependence of these coefficients for the three fine-structure excitations.

Zoom In Zoom Out Reset image size

Surprisingly, we find that the $j=1\to j^{\prime} =0$ transition is predominantly determined by a single potential of the ²Σ electronic states. For the other transitions, we find that the dynamics of fine-structure excitations at low temperatures is determined by a combination of all four potentials, each with a varying degree of importance. However, at temperatures above 500 K, the rate constants for the $j=2\to j^{\prime} =1$ and $j=2\to j^{\prime} =0$ transitions are much more sensitive to variations of the ⁴Π potential than to variations of any of the other potentials. We have verified these results by direct coupled channel calculations of the rates as functions of different potential parameters.

The sensitivity analysis is useful for practical purposes. It can be used to guide further refinement of the rate constant calculations. For example, if the rate constants for the $j=1\to j^{\prime} =0$ transition need to be computed with high precision, our results show that the ab initio calculations must focus on determining the interaction potential for the ²Σ state. In order to refine the rate constants for the $j=2\to j^{\prime} =1$ and $j=2\to j^{\prime} =0$ transitions at high temperatures, the ab initio calculations must refine the ⁴Π potential. On the other hand, the attractive ${}^{2}{\rm{\Pi }}$ potential is important for the rate constants at low temperatures. Note that moving the position of the minimum or the repulsive wall of this potential does not appear to have a significant effect on the rate constants at any temperature.

In this work, we have computed the rate constants for the fine-structure transitions in collisions of O(${}^{3}{P}_{j}$) with H(²S), which is of importance in astrophysics, using the most accurate potentials available to date. We have developed an approach to analyze the variation of the rate constants upon modulation of the interaction potentials for collision systems with complex potentials. The fine-structure transitions in ${\rm{O}}({}^{3}{P}_{j})+{\rm{H}}$ collisions represent an example of non-adiabatic dynamics determined by multiple adiabatic interaction potentials. The method described here can be used to obtain information on the sensitivity of the rate constants to each of the underlying potentials and to compute error bars of the rate constants corresponding to a certain variation of each potential. In this work, we have assumed that the potentials are known to better than 20% and computed the distributions of rate constants corresponding to the 20% variation of all four potentials determining the collision dynamics. The approach presented here can be easily extended to systems with more potentials.

We note that the exact uncertainty of the molecular potentials is difficult to assess. In a recent work, Faure et al. (2016) computed the interaction potential for the H₂ + CO molecular system with extremely high precision. The authors estimate the uncertainty of the ab initio calculations to be within 1%. To the best of our knowledge, the calculation of Faure et al. (2016) is one of the most accurate quantum chemistry calculations to date. To compute the interaction potentials for the OH molecular systems, Parlant & Yarkony (1999) used a state-of-the-art multireference configuration interaction method. However, the open-shell electronic structure of the OH molecule makes these calculations difficult. The calculations of Parlant & Yarkony (1999) should thus be expected to be less accurate than those of Faure et al. (2016). However, given that the ab initio calculations can produce results accurate to close to 1%, we believe that the error of 20% adopted here for the OH molecular potentials is an overestimate. Our results show that the 20% variation of the potentials affects the rate constants by less than 50% at low temperatures and less than 20% at elevated temperatures. These errors are generally acceptable for models used in astrophysics. Further refinement of the O(³P) + H(²S) scattering rates, if necessary, must focus on assessing and improving the accuracy of the molecular interaction potentials, guided by the results of Figure 5.

We thank Professor Paul Dagdigian for the discussions that have led us to discover the issue with the potentials used by Abrahamsson et al. (2007). The work was supported by NSERC of Canada.