Rotationally inelastic collisions of SiO with H₂

ABSTRACT

The silicon monoxide (SiO) molecule is a key species for the study of the interstellar medium as it is used to trace warm shocked gas. A large number of transitions, including high rotational levels, are observed, and the modelling of these emission lines can provide valuable information on the chemical and physical conditions of the observed regions. In these environments, where the local thermodynamical equilibrium approximation is not valid, an accurate modelling requires collisional rate coefficients with the most abundant species. We focus on the calculation of rate coefficients of SiO in its ground vibrational state in collision with para- and ortho-H₂ using a new high accurate 4D potential energy surface. Dynamical calculations of pure rotational (de)excitation of SiO were performed for the lowest 21 rotational levels using the close-coupling (CC) approach, while the coupled-state (CS) approximation was used to derive rate coefficients among the first 30 rotational levels. State-to-state rate coefficients were obtained for temperatures ranging from 5 to 300 K in the CC calculations and for temperatures up to 1000 K in the CS approximation. Propensity rules show that rate coefficients for Δj₁ = 1 transitions are dominant for both para- and ortho-H₂ colliders. The rotational rate coefficients are compared with recent results obtained for j₁ ≤ 5 levels in a full dimensionality approach. These new data will help to model emission lines in warm environments such as shocked layers of molecular outflows in star-forming regions.

molecular data, molecular processes, scattering, ISM: abundances

1 INTRODUCTION

The silicon monoxide (SiO) molecule was first detected by Wilson et al. (1971) in the interstellar medium. Then, SiO was observed in a large variety of objects and it was shown that it is mostly present in regions associated with warm, dense, and shocked gas. Rotational lines of SiO are thus used as tracers of molecular outflows surrounding star-forming regions (Gusdorf et al. 2008). Maser lines involving several vibrational levels of SiO were observed in circumstellar environments of asymptotic giant branch (AGB) stars (Tercero et al. 2011; Desmurs et al. 2014; Wong et al. 2016) and high vibrational (⁠|$v$| = 0 − 3) SiO mega-masers were detected near the centre of Seyfert 2 galaxy NGC 1068 (Wang et al. 2014). Agúndez et al. (2012) observed in the inner layers of IRC+10216 the j₁ = 2 − 1 to j₁ = 8 − 7 rotational transitions of the |$v$| = 0 state and three transitions in the |$v$| = 1 state, as well as several transitions of ²⁹SiO, Si¹⁸O, and Si¹⁷O. Using the Heterodyne Instrument for the Far Infrared (HIFI) instrument aboard the Herschel Space Observatory, Justtanont et al. (2012) have identified in O-rich AGB stars several high-j₁ rotational transitions of SiO, up to the 16–15 line in the |$v$| = 0 state and the 23–22 line in |$v$| = 1. In order to study the chemical processes at different positions in the circumstellar envelope of the carbon star IRS+10216, Fonfría et al. (2014) performed interferometric observations of several molecules among which SiO (⁠|$v$| = 0, 1, j₁ = 6 − 5). Observations of a large number of rotational transitions covering a wide range of excitation energies coupled to radiative transfer calculations put constraints on chemical and physical conditions. However, as pointed out by Fonfría et al. (2014), a proper modelling of the emission lines remains challenging due to the lack of accurate collisional rate coefficients, in particular with the most abundant species, that is generally H₂.

Dayou & Balança (2006) calculated SiO rotationally inelastic rate coefficients for collisions with He for SiO rotational levels up to j₁ = 26 and temperature up to T = 300 K in a full close-coupling (CC) approach with a highly accurate rigid rotor potential energy surface (PES). An estimation of rate coefficients for collisions with H₂ were also derived from the same rigid-rotor PES but using the SiO-H₂ reduced mass in the dynamics calculations (Dayou & Balança 2006). More recently, using the vibrational close-coupling rotational infinite order sudden (VCC-IOS) approach and an accurate 3D-PES, Balança & Dayou (2017) have calculated ro-vibrational excitation of SiO by He for the first six SiO vibrational levels. Finally, very recent results for the ro-vibrational excitation of SiO by collision with H₂ were reported by Yang et al. (2018). In this work, a highly accurate full dimensional (6D) PES was computed and full quantum CC calculations of rotational and ro-vibrational excitation of SiO (⁠|$v$| = 0, 1) were performed for the six lowest SiO rotational levels. These calculations show that vibrational excitation has little effect on pure rotational excitation, the pure rotational cross-sections being almost identical in |$v$| = 0 and |$v$| = 1. Moreover, as previously reported in the VCC-IOS study of SiO-He ro-vibrational excitation (Balança & Dayou 2017), the rate coefficients corresponding to the vibrationally inelastic processes are several orders of magnitude lower than the pure rotational coefficients, what justifies to study pure rotational excitation at relatively low energies.

The aim of this paper is to present rate coefficients for the rotational (de)excitation of SiO by collision with both para- and ortho-H₂ species. Fully quantum scattering calculations were carried out using a new, 4D, ab initio PES. The CC method was used to compute rate coefficients for the lowest 21 rotational levels of SiO and for the temperature range 5–300 K, while the coupled states (CS) approximation (McGuire & Kouri 1974) extended the calculations to the first 30 rotational levels and temperatures up to 1000 K. The paper is organized as follows: Sections 2 and 3 describe the theoretical methods used for the PES and scattering calculations. The results are presented and discussed in Section 4. Conclusions are given in Section 5.

2 SiO-H₂ potential energy surface

2.1 Electronic structure calculations

In the present work, we focus on inelastic processes for rotational levels of SiO lower than the threshold of the first excited vibrational mode (1230 cm⁻¹, i.e. 1769 K for SiO |$v$| = 1). In a recent discussion on the importance of the vibrational motion on the rotational excitation process, Faure et al. (2016a) conclude that, in the CO-H₂ test case, a full 6D treatment is not needed for rotational excitation at low energies, as an equally good agreement with the experimental results (Chefdeville et al. 2015) is achieved using an approximate 4D treatment where the vibrational motion of both molecules is not taken into account. However, as previously suggested by Mas & Szalewicz (1996), the most adequate choice of rigid-body geometry is the vibrationally averaged geometry and not the equilibrium geometry. This was also found by Valiron et al. (2008) in a study of the H₂O-H₂ system. In the present study, the two collisional partners were considered as rigid species, with bond lengths set equal to their average value in the ground vibrational state. For H₂, we used |$r_{\mathrm{H}_{2}} = 1.448736$| bohr (Le Roy & Hutson 1987). For SiO we built a Rydberg–Klein–Rees (RKR) diatomic potential using the procedure of Le Roy (Le Roy 2004) and the spectroscopic constants of Lovas, Maki & Olson (1981). The vibrationally averaged value r_SiO = 2.860462 bohr was obtained from vibrational wavefunctions computed by the Fourier Grid Hamiltonian method (Marston & Balint-Kurti 1989).

The four-dimensional (4D) PES was built using the coordinate system shown in Fig. 1. The geometry of the SiO-H₂ complex is described by three angles (θ₁, θ₂, φ₂), and the distance R between the centres of masses of SiO and H₂. The polar angles of SiO and H₂ with respect to the space-fixed Z-axis chosen along |${\boldsymbol R}$| are denoted by θ₁ and θ₂, respectively, while φ₂ is the dihedral angle between the half-planes containing the SiO and H₂ bonds. With the chosen convention θ₁ = θ₂ = 0° corresponds to the linear O-Si-H-H configuration.

Figure 1.

4D coordinate system (R, θ₁, θ₂, φ₂) employed to describe the interaction between the SiO and H₂ molecules.

The PES was computed for the ground electronic state of SiO-H₂ by means of ab initio quantum chemistry calculations, using the coupled-cluster method at the CCSD(T) level of theory (Hampel, Peterson & Werner 1992; Watts, Gauss & Bartlett 1993). The calculations were performed using standard aug-cc-pVXZ basis sets (Woon & Dunning 1994) with double (X = D), triple (X = T), and quadruple (X = Q) zeta quality. These three basis sets are designated below as aVDZ, aVTZ, and aVQZ, respectively. Each basis set was employed together with the [3s2p1d] bond functions of Williams et al. (1995), located at the position |${\boldsymbol R}_{b}$| given by Akin-Ojo, Bukowski & Szalewicz (2003):

\begin{eqnarray*} {\boldsymbol R}_{b} = \frac{1}{2} \frac{ \sum _{ij}w_{ij}\left({\boldsymbol R}_{i}+{\boldsymbol R}_{j}\right)}{\sum _{ij} w_{ij}}, \end{eqnarray*}

(1)

where |$w_{ij} = r_{ij}^{-6}$|⁠, and |$r_{ij} = |{\boldsymbol R}_{j}-{\boldsymbol R}_{i}|$| is the distance between atoms i and j with positions |${\boldsymbol R}_{i}$| and |${\boldsymbol R}_{j}$|⁠, respectively. At each geometry (R, Ω) ≡ (R, θ₁, θ₂, φ₂), the total energy |$E_{{\rm SiO-H_2},X}$| calculated for a given basis set aVXZ was corrected from the basis set superposition error (BSSE) following the counterpoise procedure of Boys and Bernardi (Boys & Bernardi 1970):

\begin{eqnarray*} V_X({R, \Omega })& = &E_{{\rm SiO-H_2},X}({R, \Omega }) \nonumber \\ & & -E_{{\rm SiO},X}({R, \Omega })-E_{{\rm H_2},X}({R, \Omega }), \end{eqnarray*}

(2)

where the SiO and H₂ species are computed with the SiO-H₂ basis set, and V_X is the SiO-H₂ interaction potential corrected from BSSE. All calculations were carried out using the molpro suite of programs (Werner et al. 2012). The grid for the radial coordinate R included 44 values, with R varying from 2.5 to 9.5 bohr by steps of 0.25 bohr, from 10 to 14 bohr by steps of 0.5 bohr, and then 15, 16, 18, 20, 25, and 30 bohr. The calculations were performed at 1000 angular geometries for each R value for both the aVDZ and aVTZ basis sets, while a smaller subset of 310 angular geometries was considered for the aVQZ basis set. The values of the angular coordinates Ω = (θ₁, θ₂, φ₂) were selected following a random sampling procedure (Valiron et al. 2008).

Following the strategy developed for the H₂O-H₂ (Valiron et al. 2008) and the SO₂-H₂ (Spielfiedel et al. 2009) systems, a complete basis set (CBS) extrapolation scheme was employed to determine two distinct PESs, V_DT and V_TQ. Each PES corresponds to the CBS limit of the interaction potential obtained by extrapolation of the correlation energies computed for given aVXZ basis sets, the aVDZ and aVTZ basis sets for V_DT, and aVTZ and aVQZ basis sets for V_TQ (see Appendix  A for details). The V_DT and V_TQ CBS PESs (obtained from equation A1 in Appendix  A) are finally employed to build the SiO-H₂ PES as follows:

\begin{eqnarray*} V({R, \Omega }) = V_{\rm DT}({R, \Omega })+V_c({R, \Omega }), \end{eqnarray*}

(3)

where V_c = (V_TQ − V_DT) is an energy difference that serves as a correction to the CBS PES V_DT obtained with the smaller basis sets. The V_DT PES is known for a larger set of angular geometries (1000 geometries) than the PES correction V_c (310 geometries). However, the anisotropy of the PES correction is weak, which allows a rapid convergence of the angular expansion employed to represent this correction term.

2.2 Analytical representation

In order to solve the quantum CC equations for the scattering process, an analytical representation of the SiO-H₂ PES of equation (3) was determined by expanding separately the CBS PES V_DT and the correction term V_c over a set of suitable angular functions for intermolecular distance R (Green 1975):

\begin{eqnarray*} V\left(R,\theta _1,\theta _2,\varphi _2\right) = \sum _{l_{1}l_{2}l}v_{l_{1}l_{2}l}(R) A_{l_1 l_2 l} \left(\theta _1, \theta _2, \varphi _2\right), \end{eqnarray*}

(4)

where the expansion coefficients |$\nu _{l_{1}l_{2}l}(R)$| are radial functions, and the basis functions |$A_{l_1 l_2 l} (\theta _1, \theta _2, \varphi _2)$| are products of associated Legendre polynomials P_lm:

\begin{eqnarray*} A_{l_1 l_2 l} \left(\theta _1, \theta _2, \varphi _2\right) &=& \sqrt{\frac{2l_1+1}{4\pi }} \Big \lbrace \langle l_10l_20\mid l0\rangle P_{l_10}(\theta _1)P_{l_20}(\theta _2) \nonumber \\ &&+\,\, 2 \sum _{m}(-1)^{m} \langle l_1ml_2-m\mid l0\rangle \\ &&\times \, P_{l_1m}(\theta _1)P_{l_2m}(\theta _2) \cos (m \varphi _2)\Big \rbrace, \end{eqnarray*}

(5)

where 〈....|....〉 is a Clebsch–Gordan coefficient. Here l₁, l₂ are associated with the rotational motion of SiO and H₂, respectively, the index l₂ being restricted to even values due to the homonuclear symmetry of the hydrogen molecule.

Owing to the large anisotropy of the SiO-H₂ interaction at short intermolecular distances that leads to steric hindrance in the angular expansion, we adopted the regularization procedure proposed by Wernli (2006) and Wernli et al. (2007) using the following formula:

\begin{eqnarray*} \begin{array}{l l l l}V_{\rm reg} & = & V & (V\le V_{a})\\ & = & V_{a}+(V_{b}-V_{a})\displaystyle {f\left(\frac{V-V_{a}}{V_{b}-V_{a}}\right)} & (V_{a} \le V \le V_{b}) \\ & = & V_{a} +\displaystyle {\left(\frac{2}{\pi }\right)^2}(V_{b}-V_{a}) & (V \ge V_{b}) \end{array}, \end{eqnarray*}

(6)

where the scaling function is given by

\begin{eqnarray*} f(u) = \left(\frac{2}{\pi }\right)^{2}\sin \left[\frac{\pi }{2} \sin \left(\frac{\pi u}{2}\right)\right]. \end{eqnarray*}

(7)

We have selected as threshold values |$V_{a} = 10\, 000$| cm⁻¹ and |$V_{b} = 50\, 000$| cm⁻¹. The above regularization procedure was applied separately to the V_DT and V_TQ CBS PESs.

At each intermolecular distance, the regularized PES V_DT and the correction term V_c were then developed over the angular expansion (equation 4) using a standard linear least-square procedure. We included all anisotropies up to l₁ = 12, l₂ = 6, and l = 18, resulting in 174 expansion functions. Only significant terms were selected using a Monte Carlo error estimator (defined in Rist & Faure (2012)), resulting in a final set of 68 expansion functions for V_DT, with anisotropies up to l₁ = 12, l₂ = 6, and l = 13. For the correction V_c, 21 expansion functions were retained, with anisotropies up to l₁ = 9, l₂ = 4, and l = 9. A cubic spline radial interpolation of the expansion coefficients was finally employed over the whole intermolecular distance range (R = 2.5 − 30 bohr) and it was smoothly extrapolated (using exponential and power laws at short- and long-range, respectively) as in Valiron et al. (2008). The root mean square (rms) residual was found to be lower than 1 cm⁻¹ from the long-range part of the interaction potential down to the regions of the global and secondary minima (R ≥ 5 bohr).

The relatively large number of l₁, l₂ terms required to obtain a good accuracy of the fit is due to the large anisotropy of the interaction as illustrated in Fig. 2 that shows the R-dependence of the PES for different sets of (θ₁, θ₂, φ₂) angles. The global minimum of the fitted PES is equal to −300.472 cm⁻¹ and has been found for the collinear configuration Si-O-H-H (θ₁ = 180°, θ₂ = φ₂ = 0°) with R = 7.431 bohr. The large anisotropy is also pointed out in Fig. 3 that shows a 2D contour plot of the PES as a function of θ₁ and θ₂ for fixed value of φ₂ = 0° and R = 7.431 bohr. Comparison with the 6D PES of Yang et al. (2018) is discussed in Section 4.3 below.

Figure 2.

Interaction potential for the SiO-H₂ system as a function of the intermolecular distance R for selected orientations (θ₁, θ₂, φ₂) of the SiO and H₂ molecules (see Fig. 1).

Figure 3.

Contour plot for the SiO-H₂ system as a function of θ₁, θ₂ with φ₂ = 0° and R = 7.431 bohr. Contours are drawn by step of 50 cm⁻¹.

It is interesting to analyse the contribution of (l₁, l₂, l) terms to the angular expansion of the PES, as the anisotropic expansion coefficients |$v_{l_1l_2l}(R)$| are the coupling terms responsible for the rotationally inelastic SiO(j₁) →SiO(⁠|$j^{\prime }_1$|⁠) transitions (Green 1975). Taking the expansion at R = 7.5 bohr (near the global minimum), one may observe in Table 1 that, except the isotropic term (000), the dominant terms are the (101), (123) (202), and (224), and at R = 20 bohr, the dipole–quadrupole (123) and the quadrupole–quadrupole (224) terms are the dominant ones. The contribution of these two terms to |$j_{1}\rightarrow j^{\prime }_{1}$| transitions is missing when the H₂ rotational basis set is restricted to j₂ = 0. It thus seems necessary to include j₂ = 2 in the H₂ basis set to properly describe the rotational excitation of SiO by para-H₂. The results of Table 1 can also explain some differences between the excitation cross-sections (and rate coefficients) corresponding to para-H₂ and ortho-H₂ colliders (see Section 4.2).

Table 1.

Dominant expansion terms (in cm⁻¹) of the SiO-H₂ PES, at R = 7.5 and R = 20 bohr.

l1	l2	l	|$v_{l_{1}l_{2}l}(R = 7.5)$|	|$v_{l_{1}l_{2}l}(R = 20.0)$|
0	0	0	−2329.097	−10.633
1	0	1	773.244	−0.107
1	2	3	704.797	13.153
2	0	2	391.267	−0.534
2	2	4	−342.839	−2.526
3	0	3	251.484	0.078

3 SCATTERING CALCULATIONS

In this paper, we focus on the following process:

\begin{eqnarray*} {\rm SiO}(j_{1})+{\rm H_{2}}(j_{2}) \rightarrow {\rm SiO}(j^{\prime }_{1})+{\rm H_{2}}(j_{2}), \end{eqnarray*}

(8)

where j₁ and j₂ are the rotational angular momenta of SiO and H₂, respectively. We only consider |$j_{1}\rightarrow j^{\prime }_{1}$| transitions where j₂ remains fixed. However, the calculations include the possibility of H₂ excitation during the collision with a low probability at the considered energies. Considering the important contribution of the H₂ quadrupole terms to the angular expansion of the PES (see Table 1), the H₂ basis set was taken equal to j₂ = 0, 2 for para-H₂ and j₂ = 1, 3 for ortho-H₂. CC calculations using a smaller H₂ basis set (j₂ = 0 and j₂ = 1 for para and ortho-H₂, respectively) were performed for comparison. The results show differences up to 20 per cent for collisions with para-H₂, whereas the cross-sections differ by less than 5 per cent for collisions with ortho-H₂. This shows that the contribution of H₂ quadrupole terms to |$j_{1}\rightarrow j^{\prime }_{1}$| transitions is important and needs to be taken into account. This quadrupole contribution is included for collisions with ortho-H₂j₂ = 1 but not for collisions with para-H₂j₂ = 0. As a consequence, all CC calculations for para-H₂ were performed with the j₂ = 0, 2 basis set while, to avoid prohibitively long calculations, CC cross-sections for ortho-H₂ at energies larger than 800 cm⁻¹ were obtained using the j₂ = 1 basis set.

The energies of the rotational levels of SiO were obtained from the rotational and distortion constants of Lovas et al. (1981) (B_e = 0.7267 cm⁻¹, α_e = 5.0379 × 10⁻³ cm⁻¹, and D_e = 9.9232 × 10⁻⁷ cm⁻¹). For H₂, we used B_e = 60.853 cm⁻¹, α_e = 3.062 cm⁻¹, and D_e = 0.0471 cm⁻¹. Full quantum CC calculations include all (de)excitation processes among the 21 lowest SiO rotational levels. For CC calculations, the basis set for SiO includes, at each collisional energy, all open channels plus several closed channels. For the largest energies, the basis set included up to j₁ = 30 SiO rotational levels to insure convergence of the cross-sections. A basis set of 40 SiO rotational levels were considered in the CS calculations to obtain converged cross-sections and rate coefficients up to j₁ = 30. The energies of the lowest 41 rotational levels of SiO are presented and compared with CDMS experimental data (Endres et al. 2016) in Table 2.

Table 2.

Theoretical and CDMS (Endres et al. 2016) SiO rotational energy levels |$\epsilon _{j_{1}}$|⁠, in cm⁻¹.

j1	Theory	CDMS	j1	Theory	CDMS
0	0.000	0.000	21	334.384	334.383
1	1.449	1.449	22	366.208	366.207
2	4.345	4.345	23	399.474	399.473
3	8.691	8.691	24	434.183	434.182
4	14.484	14.484	25	470.332	470.331
5	21.726	21.726	26	507.926	507.921
6	30.416	30.416	27	546.953	546.951
7	40.554	40.554	28	587.423	587.421
8	52.140	52.140	29	629.332	629.329
9	65.173	65.173	30	672.679	672.676
10	79.654	79.653	31	717.463	717.460
11	95.581	95.581	32	763.684	763.680
12	112.956	112.956	33	811.340	811.336
13	131.778	131.778	34	860.432	860.427
14	152.045	152.045	35	910.958	910.953
15	173.759	173.759	36	962.918	962.912
16	196.918	196.918	37	1016.310	1016.303
17	221.522	221.522	38	1071.134	1071.126
18	247.572	247.571	39	1127.389	1127.380
19	275.065	275.065	40	1185.073	1185.064
20	304.003	304.003

The reduced mass of SiO-H₂ was taken equal to 1.927303 amu. CC calculations were performed for total energies up to 2000 cm⁻¹ for collisions with para- and ortho-H₂. Convergence of the cross-sections required total angular momentum channels up to J = 128 to be included at the largest collision energies. We used a rather fine energy grid of 0.05 cm⁻¹ for E ≤ 9.0cm⁻¹, 0.1 cm⁻¹ for 9.0 ≤ E ≤ 250 cm⁻¹, 0.2 cm⁻¹for 250 ≤ E ≤ 480 cm⁻¹, and 10 cm⁻¹ for 480 ≤ E ≤ 550 cm⁻¹, plus 15 energy values up to 2000 cm⁻¹. Several points were added near the SiO rotational energy thresholds. The same energy grid was used for ortho-H₂ and CS calculations for E ≤ 2000 cm⁻¹. CS calculations were extended up to 10 000 cm⁻¹ to obtain rate coefficients up to T = 1000 K. All the scattering calculations were performed with the quantum scattering code molscat (Hutson & Green 1994). The modified log-derivative Airy propagator of Alexander & Manolopoulos (1987) with a variable step size was used to solve the scattering equations from R = 4 to R = 50 − 100 bohr depending on the collision energy domain.

From the rotationally inelastic cross-sections |$\sigma _{j_1 \rightarrow j_1^{\prime }}(E)$|⁠, one can obtain the corresponding thermal rate coefficients at temperature T by an average over the collision kinetic energy E_k:

\begin{eqnarray*} k_{j_1 \rightarrow j^{\prime }_1}(T) &=& \left(\frac{8k_{{\rm B}}T}{\pi \mu }\right)^{\frac{1}{2}}\left(\frac{1}{k_{{\rm B}}T}\right)^{2}\\ &&\times \,\int _{0}^{\infty }E_{k} \sigma _{j_1 \rightarrow j^{\prime }_1} (E_{\rm k}) e^{\frac{-E_{\rm k}}{k_{{\rm B}}T}}dE_k \end{eqnarray*}

(9)

where k_B is the Boltzmann constant. Fully converged CC rate coefficients were obtained for the 21 lowest SiO rotational levels and temperatures up to 300 K, while CS calculations give converged rate coefficients for the 30 lowest SiO levels and temperatures up to 1000 K. While we report in the following de-excitation cross-sections and rate coefficients, the results for the reverse transitions are obtained from the detailed balance relation:

\begin{eqnarray*} k_{j^{\prime }_1 \rightarrow j_1}(T) = \frac{2j_1+1}{2j^{\prime }_1+1} \exp \left(\frac{\epsilon _{j^{\prime }_1}-\epsilon _{j_1}}{k_{{\rm B}}T}\right) k_{j_1\rightarrow j^{\prime }_1}(T), \end{eqnarray*}

(10)

where |$\epsilon _{j_1}$| denotes the energy of the rotational level j₁.

4 RESULTS

4.1 Cross-sections

Cross-sections for the rotational de-excitation of SiO by para-H₂ and ortho-H₂ are presented in Fig. 4 as a function of kinetic energy for |$\Delta j_1 = j_1-j^{\prime }_1 = 1$| transitions. The para- and ortho-H₂ cross-sections are similar in magnitude and exhibit many resonances at kinetic energies below 100 cm⁻¹. This relates to the presence of quasi-bound states formed into the potential well of the PES before the collisional complex dissociates. For SiO in collision with ortho-H₂, the cross-sections appear to present a smoother energy dependence than the para-H₂ ones. This is due to the fact that the number of coupling terms |$v_{l_1l_2l}(R)$| involved in a given transition |$j_{1}\rightarrow j^{\prime }_{1}$| is larger for ortho-H₂ than for para-H₂, what gives rise to many more overlapping resonances in the cross-sections in the case of ortho-H₂. Cross-sections for other Δj₁ transitions exhibit a similar behaviour.

Figure 4.

State-to-state CC cross-sections for the rotational de-excitation of SiO by para-H₂ (upper panel) and ortho-H₂ (lower panel) for |$\Delta j_1 = j_1-j^{\prime }_1 = 1$| as function of the kinetic energy. The cross-sections are labelled according to the related transition |$j_{1}\rightarrow j^{\prime }_{1}$| between the rotational states of SiO.

4.2 Rate coefficients

Rate coefficients for the rotational de-excitation of SiO by para-H₂ and ortho-H₂ are shown as a function of temperature in Figs 5 and 6 for Δj₁ = 1 and Δj₁ = 2 transitions, respectively. For a given value of Δj₁, the rate coefficients display a similar temperature dependence, different in the para- and ortho-H₂ cases. Their magnitude increases with the initial j₁ state, and more rapidly for the low j₁ states than for the higher ones. For a given initial state j₁, the rate coefficients decrease with increasing Δj₁. These general features can be understood by considering the coupling terms |$v_{l_1l_2l}(R)$| involved in each |$j_{1}\rightarrow j^{\prime }_{1}$| transition (Green 1975). Indeed, for a given Δj₁, the main contribution is due to coupling terms |$v_{l_1l_2l}(R)$| for which l₁ ≥ Δj₁, with even/odd l₁ values for even/odd Δj₁. Furthermore, the number of coupling terms involved increases with j₁, since |$l_1\le j_1+j^{\prime }_1$|⁠, but their magnitude decreases with increasing l₁ (see Table 1). The case ortho-H₂ differs from para-H₂ with the contribution of coupling terms associated with l₂ ≥ 2 for each |$(j_1,j^{\prime }_1)$| pair of states, whereas they contribute through indirect coupling with other rotational states in the case para-H₂.

Figure 5.

State-to-state CC rate coefficients for the rotational de-excitation of SiO by para-H₂ (upper panel) and ortho-H₂ (lower panel) for Δj₁ = 1 as function of temperature. The rate coefficients are labelled according to the related transition |$j_{1}\rightarrow j^{\prime }_{1}$| between the rotational states of SiO.

Figure 6.

Same as Fig. 5 for Δj₁ = 2.

We display in Fig. 7 the de-excitation rate coefficients for rotational transitions out of the initial state j₁ = 20 of SiO at two selected temperatures, T = 10 K and 100 K. As expected the rate coefficients decrease with increasing Δj₁, since the main contribution is due to coupling terms of decreasing magnitude (l₁ ≥ Δj₁) for both para- and ortho-H₂ species. In the case of para-H₂ and low Δj₁ values, there is a slight propensity rule favouring odd Δj₁ transitions over even Δj₁ ones. This can be related to the marked asymmetry of the PES regarding the orientation of the SiO molecule (variation of the PES along the coordinate θ₁, see Fig. 3), which gives rise to an important magnitude of coupling terms associated with odd l₁ values. Such a propensity rule is not observed for ortho-H₂ due to the contribution of additional coupling terms with l₂ ≥ 2. At low temperature the ortho-H₂ rate coefficients are generally larger than the para-H₂ ones due to a dominant contribution of l₂ = 2 quadrupole terms in the long-range part of the interaction. At higher temperature the para- and ortho-H₂ rate coefficients are in closer agreement. Such a similarity was also found for collisions of H₂ with HC₃N (Faure, Lique & Wiesenfeld 2016b), C₆H⁻ (Walker et al. 2017) and CN⁻ (Kłos & Lique 2011), and can be considered as due to the presence in these systems of a large well depth where the orientation of the H₂ fragment is averaged during the collision (Walker et al. 2017).

Figure 7.

State-to-state CC rate coefficients for the rotational de-excitation of SiO(j₁ = 20) by para-H₂ and ortho-H₂ as a function of |$\Delta j_1 = j_1-j^{\prime }_1$| at temperatures T = 10 K (upper panel) and T = 100 K (lower panel).

In order to emphasize more globally the differences between the para- and ortho-H₂ rate coefficients, we compare in Fig. 8 the two sets of coefficients for all the de-excitation transitions with initial states up to j₁ = 20 at T = 10 K and 100 K. The horizontal axis represents the rate coefficients with para-H₂ and the vertical axis those with ortho-H₂. At 100 K, it appears clearly that the two sets of data agree generally within a factor of 2, the largest rate coefficients being, for both sets of data, those for Δj₁ = 1, 2. At 10 K, most of ortho-H₂ rate coefficients are larger than the para-H₂ ones, reflecting the important contribution of l₂ = 2 long-range quadrupole terms at low energies.

Figure 8.

Comparison between para-H₂ and ortho-H₂ CC rate coefficients for the rotational de-excitation of SiO(j₁ ≤ 20) at temperatures T = 10 K (upper panel) and T = 100 K (lower panel). The horizontal axis corresponds to para-H₂ rate coefficients and the vertical axis to ortho-H₂ ones. The two dashed lines define the region where the rate coefficients differ by less than a factor of 2.

Full quantum CC calculations are computationally very expensive at high kinetic energies when the number of open channels increases that requires large basis sets. The CS approximation (McGuire & Kouri 1974) should provide a reasonable estimate for calculations at high energies far enough from the transition thresholds. This approximation assumes that off-diagonal Coriolis couplings are negligible. The accuracy of this approximation was studied in different collisional situations, rotational excitation of HCl (Heil, Kouri & Green 1978), H₂O (Dubernet et al. 2009), SO₂ (Balança, Spielfiedel & Feautrier 2016), C₆H⁻ (Walker et al. 2016) in collision with H₂, and recently ro-vibrational excitation of H₂ (Bohr et al. 2014) and CO (Forrey et al. 2015) in collision with H₂. Differences between the CS and CC cross-sections and rate coefficients are highly dependent on the systems, the transitions and the range of collision energy considered.

Fig. 9 displays the comparison between the CC and CS rate coefficients at 100 K for all the de-excitation transitions corresponding with states j₁ ≤ 20. One can observe a relatively good agreement between the two sets of rate coefficients for collisions with para-H₂. For ortho-H₂, discrepances up to 50 per cent are found for the rate coefficients of large magnitude (those corresponding to transitions with small Δj₁ values), while they can be as large as a factor of 2 or 3 for the rate coefficients of very small magnitude (corresponding to large Δj₁ values). This can be explained by the largest number of angular momentum couplings involved in collisions with ortho-H₂. From this comparison, we estimate that, for SiO in collision with para- and ortho-H₂, the CS results are reliable to within 50 per cent for the main transitions at high temperature and give the correct trends of the rate coefficients for all transitions. The complete sets of CC and CS para- and ortho-H₂ de-excitation rate coefficients will be made available through the LAMDA (Schöier et al. 2005) and BASECOL (Dubernet et al. 2013) databases.

Figure 9.

Comparison between CC and CS rate coefficients for the rotational de-excitation of SiO(j₁ ≤ 20) by para-H₂ (upper panel) and ortho-H₂ (lower panel) at the temperature T = 100 K.

4.3 Comparison with previous results

A full-dimensional study of the SiO-H₂ collisional system has been published very recently by Yang et al. (2018). From a 6D interaction potential, these authors calculated pure rotational (j₁ ≤ 5) as well as rovibrational (⁠|$v$| = 0, 1, j₁ ≤ 5) SiO-H₂ cross-sections and rate coefficients for both para- and ortho-H₂ species. The effect of vibrational excitation is shown to have little impact on pure rotational transitions. This allows us to compare the present 4D and the 6D (Yang et al. 2018) results for rotational transitions in |$v$| = 0.

The first ingredient for the treatment of the SiO-H₂ scattering is the interaction potential. The two PESs were both calculated using high-level quantum chemistry methods. The global minimum reported for the 6D PES is equal to −284.2 cm⁻¹ (ab initio) and −279.5 cm⁻¹ (fit) at R = 7.4 bohr for the linear Si-O-H-H arrangement and monomer equilibrium values (Yang et al. 2018). Using their fitted PES at the vibrationally averaged SiO and H₂ geometries considered in our 4D ab initio PES, we found a minimum of −293.65 cm⁻¹ to be compared to our value −300.47 cm⁻¹. The two values agree within 2 cm⁻¹ if we consider that their fitted PES is∼ 5 cm⁻¹ higher than the ab initio value in the region of the global minimum (in our case the difference is ∼ 0.4 cm⁻¹ at the global minimum). The R-dependence of the SiO-H₂ interaction at selected (θ₁, θ₂, φ₂) orientations, as well as the angular dependence around the global minimum, appear very similar for the two PESs (compare present Figs 2 and 3 with the ones of Yang et al. (2018)). Note the different convention used for the angular coordinates: θ₁(4D-PES) |$= \pi +\theta _1$| (6D-PES).

We compare in Fig. 10 our rate coefficients with those reported by (Yang et al. 2018) for the rotational de-excitation of SiO by para-H₂ as a function of temperature for three selected initial j₁ states. As expected from the similarity of the two PESs, we observe a good global agreement for both the magnitude and the temperature dependence of the rate coefficients. This is also a confirmation that the vibrational motion of SiO, not included in the present study, has a quite small influence on pure rotational energy transfers Yang et al. (2018); Balança & Dayou (2017). The slightly larger differences observed at low temperatures could be attributed to small differences in the description of long-range interactions or the treatment of cross-sections just above thresholds. We cannot perform such a comparison for the ortho-H₂ case as the corresponding rate coefficients values were not reported by Yang et al. (2018).

Figure 10.

State-to-state CC rate coefficients for the rotational de-excitation of SiO(j₁) by para-H₂ as a function of temperature for the initial states j₁ = 1 (upper panel), j₁ = 3 (middle panel), and j₁ = 5 (lower panel).

It is also of interest to compare the present rate coefficients obtained for SiO in collision with para-H₂ with the data obtained for collisions with He (Dayou & Balança 2006). Fig. 11 presents a comparison between the SiO-para-H₂ de-excitation rate coefficients and the SiO-He ones, multiplied by a factor of 1.4. This approximation, often used to get an estimate of rate coefficients with para-H₂ (j₂ = 0) from those with He, considers identical cross-sections values for the two colliding systems and applies the 1.4 scaling factor to account for the associated different reduced masses. As can be seen in Fig. 11, important discrepances are obtained between the para-H₂ and scaled He rate coefficients. There is no systematic trend observed, the discrepances being highly dependent on the initial state and the transition considered. The propensity rule in favour of odd Δj₁ transitions is well predicted, but the discrepances in magnitude can grow up to a factor of 3 even for the highest rate coefficients values corresponding to Δj₁ = 1, 2 transitions. Similar discrepances are obtained whatever the temperature range investigated. Such differences are mainly explained by the distinct interaction potentials as the cross-sections are very sensitive to the shape and depth of the PES well. As stated previously, differences up to 20 per cent also exist between the cross-sections for collisions with para-H₂(j₂ = 0, 2) and para-H₂(j₂ = 0) due to the contribution of H₂ quadrupole terms in the PES expansion. The use of the mass-scaling factor is thus expected to provide only the correct order of magnitude of rate coefficients for SiO in collision with para-H₂. Accordingly, the new set of SiO-H₂ rate coefficients presented in this work could have a significant influence on the analysis of astrophysical environments based on the modelling of SiO emission lines.

Figure 11.

Comparison of the SiO-para-H₂ CC rotational de-excitation rate coefficients for SiO(j₁ ≤ 20) at the temperature T = 100 K with previous CC results (Dayou & Balança 2006) obtained for SiO-He and scaled by a factor of 1.4. The two dashed lines define the region where the rate coefficients differ by less than a factor of 2.

5 CONCLUSION

Quantum CC calculations of cross-sections and rate coefficients for the rotational (de)excitation of SiO by collision with para-H₂ and ortho-H₂ are presented in this work. The scattering calculations were carried out on a new highly accurate 4D ab initio PES, with the SiO and H₂ bond lengths held fixed at vibrationally averaged values. The CC method was employed for scattering calculations for the first 21 levels of SiO and temperatures in the 5–300 K range, while the CS approximation was used to treat the rotational (de)excitation for the first 30 levels of SiO and temperatures up to 1000 K. CC and CS rate coefficients agree relatively well for collisions with para-H₂, whereas differences by factors up to 50 per cent are found in the ortho-H₂ case for transitions associated with large rate coefficients values (transitions with small Δj₁). Propensity for dominant |$\Delta j_{1} = j_1-j^{\prime }_1 = 1$| rate coefficients is obtained for both para- and ortho-H₂ species, with a rapid decrease of the rate coefficients as Δj₁ increases. A propensity rule in favour of odd Δj₁ transitions is also observed for para-H₂, but not in the ortho-H₂ case. The para-H₂ and ortho-H₂ rate coefficients can differ significantly at low temperatures, with ortho-H₂ rates generally larger than the para-H₂ ones, but are of comparable magnitude at high temperatures. State-to-state rate coefficients were compared with recent data obtained with a full-dimensional 6D approach as well as with previous approximate results obtained using a SiO-He potential and a mass-scaling factor. The present set of SiO-H₂ rate coefficients will contribute to help to accurately model astrophysical observations of SiO emission lines.

ACKNOWLEDGEMENTS

This work was supported by the CNRS programmes PCMI (Physique et Chimie du Milieu Interstellaire) and PNPS (Programme National de Physique Stellaire). The PES calculations were performed on work stations at the Centre Informatique of Paris Observatory and all the dynamics calculations were performed using HPC resources from GENCI/IDRIS (grant N°2015047344). LW and CB thank COST action CM1401 ‘Our Astrochemical History’ for some travel support. The fit of the SiO-H2 PES was performed on the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER0713 CIRA: http://www.ci-ra.org).

APPENDIX A: CBS EXTRAPOLATION

The SiO-H₂ PESs V_X (equation 2) computed using the three aVXZ basis sets with X = (D,T,Q) were employed to determine two distinct PESs, V_DT and V_TQ, by applying a CBS extrapolation scheme. Each PES corresponds to the CBS limit of the interaction potential obtained by extrapolation of the correlation energies E_{corr, X} = E_{CCSD(T), X} − E_{HF, X} computed for given aVXZ basis sets, the aVDZ and aVTZ basis sets for V_DT, and aVTZ and aVQZ basis sets for V_TQ. The correlation energies were extrapolated according to (Schwartz (1962) and Helgaker et al. (1997)) E_{corr, X} = E_{corr, CBS} + Ax⁻³, where E_{corr, CBS} is the CBS limit of the correlation energy, and x = (2, 3, 4) for X = (D, T, Q), respectively. Since the Hartree–Fock energy converges much faster than the correlation energy with the basis-set size, the Hartree–Fock energies were not extrapolated and we used the results obtained for the larger of the two basis sets. In this case, at each geometry (R, Ω) ≡ (R, θ₁, θ₂, φ₂), the CBS limit of the interaction potential is given by

\begin{eqnarray*} V_{\rm CBS}(R,\Omega) & = & V_{\rm HF,X_2}(R,\Omega)\nonumber\\ & & +\displaystyle {\frac{x_2^3}{x_2^3-x_1^3}\left[V_{\rm X_2}(R,\Omega)-V_{\rm HF,X_2}(R,\Omega)\right]}\nonumber\\ & & -\displaystyle {\frac{x_1^3}{x_2^3-x_1^3}\left[V_{\rm X_1}(R,\Omega)-V_{\rm HF,X_1}(R,\Omega)\right],} \end{eqnarray*}

(A1)

where x₂ > x₁, V_{HF, X} is the Hartree–Fock interaction potential, and V_X is the CCSD(T) interaction potential. Both interaction potentials are corrected from BSSE following equation (2). The V_DT and V_TQ CBS PESs employed to build the composite PES of equation (3) were both obtained from equation (A1).

In order to illustrate the above procedure, we display in Fig. A1 as a function of the intermolecular distance R the interaction potentials obtained at several levels of calculation for the linear configuration Si-O-H-H corresponding to the global minimum of the SiO-H₂ system. As can be seen, the ab initio interaction energies V_X are continuously lowered on increasing the size of the aVXZ basis set, while a continuous decrease of the equilibrium distance is observed. The V_DT and V_TQ CBS PESs follow the correct trend regarding the variation of the interaction energy and equilibrium distance with the basis-set size. In this sense, the V_TQ PES is expected to approach the ab initio result we would have obtained by increasing further the basis-set size. The analytical PES, which was built by fitting separately the V_DT PES and the correction term V_c (equation 3), is found in close agreement with the V_TQ result. Additional calculations were performed for several others orientations of SiO and H₂ and a similar level of agreement is obtained (note that these geometries and those of Fig. A1 were not included in the fitting procedure). This indicates that the anisotropy of the correction term V_c is sufficiently weak to be properly described from a limited number of angular geometries.

Figure A1.

Interaction potential for the SiO-H₂ system as a function of the intermolecular distance R at the angular geometry (θ₁ = 180°, θ₂ = φ₂ = 0°) corresponding to the global minimum. Comparison between the ab initio interaction potentials V_X computed from three distinct aVXZ basis sets (equation 2), the V_DT and V_TQ CBS PESs obtained from the basis set extrapolation procedure (equation A1), and the analytical potential obtained by fitting the composite PES of equation (3).