Rotational relaxation of HCO+ and DCO+ by collision with H2

The HCO + and DCO + molecules are commonly used as tracers in the interstellar medium. Therefore, accurate rotational rate coefﬁcients of these systems with He and H 2 are crucial in non-local thermal equilibrium models. We determine in this work the rotational de-excitation rate coefﬁcients of HCO + in collision with both para-and ortho-H 2 , and also analyse the isotopic effects by studying the case of DCO + . A new four-dimensional potential energy surface from ab initio calculations was developed for the HCO + –H 2 system, and adapted to the DCO + –H 2 case. These surfaces are then employed in close-coupling calculations to determine the rotational de-excitation cross-sections and rate coefﬁcients for the lower rotational states of HCO + and DCO + . The new rate coefﬁcients for HCO + + para-H 2 were compared with the available data, and a set of rate coefﬁcients for HCO + + ortho-H 2 is also reported. The difference between the collision rates with ortho-and para-H 2 is found to be small. These calculations conﬁrm that the use of the rate coefﬁcients for HCO + + para-H 2 for estimating those for HCO + + ortho-H 2 as well as for DCO + + para-H 2 is a good approximation.


I N T RO D U C T I O N
Although molecular cloud temperatures can range from 5 K up to a few tens of K, and densities from 10 2 up to 10 6 cm −3 , the temperature and densities of most molecular clouds are low, 10-15 K and 10 2 -10 4 cm −3 , respectively, and therefore, the collisions are rare. As a result, the rotational population of molecules in these regions cannot be described by a Boltzmann distribution. The physico-chemical conditions in the molecular clouds should be determined then by nonlocal thermal equilibrium (non-LTE) models. Such models require not only the knowledge of the Einstein coefficients, usually known, but also the state-to-state rate coefficients of the observed molecule with the most common colliders in the interstellar medium (ISM), named H 2 , He, e, and H. HCO + , discovered by Buhl & Snyder (1970), was the first cation identified in the cold ISM (Kraemer & Diercksen 1976). It has been detected almost in all the regions of the ISM, including dense cores, molecular outflows, and even diffuse clouds (Lucas & Liszt 1996). Its deuterated species, DCO + , is also considered as one of the primary tracers of the deuterium chemistry in the ISM (Huang et al. 2017). HCO + and DCO + are often used as tracers of ionization in dense cores (Caselli et al. 1998) and DCO + is easily observed in the molecular layer of protoplanetary discs (e.g. Dutrey et al. 2014). Accurate rate coefficients of these systems with He and H 2 are thus important for a correct description of the physico-chemical conditions in these astrophysical media. E-mail: otoniel.denis@uautonoma.cl The rotational rate coefficients of HCO + with He were reported by Buffa, Dore & Meuwly (2009). For the collision of HCO + with H 2 , the first set of rate coefficients was computed by Monteiro (1985) from a potential energy surface (PES) determined at the selfconsistent field (SCF) level of theory and limited to low temperatures (T ≤ 30 K) and the lower rotational states. The data that can be found in the LAMDA (Schöier et al. 2005) and BASECOL ) data bases are those reported by Flower (1999), who extended the calculations of Monteiro (1985). However, such rate coefficients were computed using the same PES. More recently, Yazidi, Ben Abdallah & Lique (2014) calculated a new set of rate coefficients for the lower 21 transitions of HCO + from a PES averaged over the H 2 orientation. Such data showed significant differences with the rate coefficients of Flower (1999). The use of the averaged PES in the determination of the rate coefficients was validated by the calculations of Massó & Wiesenfeld (2014). These authors employed a four-dimensional PES for the HCO + -H 2 complex at a high level of theory and computed the cross-section and rate coefficients for the lower states of HCO + . However, all available sets of rate coefficients, large enough for being included in non-LTE calculations, are restricted to the collision with para-H 2 . The only study that considered the collision with ortho-H 2 was limited to the three lower rotational states of HCO + (Massó & Wiesenfeld 2014). Therefore, a large set of rate coefficients for the collision of HCO + with ortho-H 2 deserves to be computed.
The case of DCO + is quite different. Its collision with He was also studied (Buffa 2012), and differences of about 20 per cent were found between the rate coefficients of HCO + + He and DCO + + He. However, there is no set of rate coefficients for its collision with Internal coordinates used to describe the HCO + + H 2 system. The azimuthal angle ϕ is undefined when θ 1 or θ 2 is equal to 0 • or 180 • . H 2 . Pagani, Bourgoin & Lique (2012) estimated the hyperfine rate coefficients for DCO + + H 2 but using the rate coefficients for HCO + + H 2 and ignoring the isotopic substitution of HCO + . In the main data bases for collisional data, namely LAMDA (Schöier et al. 2005) and BASECOL , the rate coefficients for DCO + + H 2 are assumed to be the same as for HCO + + H 2 , and astronomers usually employ these rate coefficients in the non-LTE models used to analyse the observations of DCO + . However, in the case of the isoelectronic molecule HCN in collision with He, the ratio between the rate coefficients for both isotopologues, namely HCN and DCN, varies from 0.4 to 3.9 (Denis-Alpizar, Stoecklin & Halvick 2015). This approximation used for DCO + thus deserves to be checked.
The main goal of this work is to determine the first set of rotational de-excitation rate coefficients of HCO + in collision with both paraand ortho-H 2 that can be used in non-LTE models from a fourdimensional PES, also developed in this work. Furthermore, this PES will be modified to describe the DCO + + H 2 system, and the dynamics of DCO + + H 2 will also be studied. This paper is organized as follows: In the next section, the methods employed are described, while the results are analysed in Section 3. Finally, the summary of this work is provided in Section 4.

Ab initio calculations
The HCO + + H 2 complex was described employing the coordinates shown in Fig. 1, where R connects the centres of mass of the H 2 and HCO + molecules and θ 1 , θ 2 , and ϕ describe the relative angular orientations. The HCO + and H 2 molecules were considered as rigid rotors. The use of this approximation has allowed an excellent prediction of experimental results for several molecule-diatom collisions, e.g. O 2 -H 2 ) and CO-D 2 (Stoecklin et al. 2017). Furthermore, recent close-coupling calculations considering the bending motions for the collision with a triatomic system have shown that the rigid rotor approximation is valid for collisional energies lower than the bending frequency (Stoecklin et al. 2013;Denis-Alpizar, Stoecklin & Halvick 2014;Denis-Alpizar et al. 2015). The bond length of H 2 was taken as the vibrationally averaged value in the rovibrational ground state r H-H = 0.767 Å (Jankowski & Szalewicz 1998), while for the linear HCO + molecule, the CO and CH distances were set to the experimental equilibrium values r CO = 1.105 Å and r CH = 1.097 Å (Woods 1988).
The ab initio calculations were performed at the explicitly correlated coupled-cluster level, including single, double, and perturbative triple excitations [CCSD(T)-F12a] using the MOLPRO package Figure 2. Ab initio energies computed at the CCSD(T)-F12a/aug-cc-pVTZ, CCSD(T)/aug-cc-pVQZ, and CCSD(T)/aug-cc-pV5Z level of theory and those extrapolated to the CBS limit, at θ 1 = 90 • and θ 2 = 0 • . (Werner et al. 2012). In such method, the direct inclusion of F12 terms in the triples is not available, and an improvement of the energy of the triples was obtained by scaling the triples' energy contribution (Knizia, Adler & Werner 2009) as implemented in MOLPRO (Werner et al. 2012). This correction gives energies close to the completed basis set (CBS) limit for the augmented correlationconsistent polarized triple zeta basis set (aug-cc-pVTZ) of Dunning (1989) (Knizia et al. 2009). Therefore, the aug-cc-pVTZ basis set was then employed in this work. The non-size consistency of the CCSD(T)-F12a method was corrected by shifting up the ab initio interaction energies to make it vanish at R = 350 Å. The basis set superposition error was corrected using the counterpoise procedure of Boys & Bernardi (1970). Fig. 2 shows the energies computed in the region of the minimum of the complex, and those employing the standard CCSD(T) method with the aug-cc-pVQZ and aug-cc-pV5Z basis sets. The two-point extrapolation formula where X and Y are the cardinality of the basis set (X = 5 for aug-cc-pV5Z and Y = 4 for aug-cc-pVQZ; Halkier et al. 1999), was used for obtaining the energies at the CBS limit. The agreement of the energies computed with the CCSD(T)-F12a/aug-cc-pVTZ method and those at the CBS limit is excellent.
In total, 49 038 ab initio energies were computed. The radial grid includes 29 values of R from 1.8 to 20.0 Å, while θ 1 and θ 2 vary from 0 • to 180 • in steps of 15 • . The azimuthal angle ϕ varies in steps of 20 • in the [0, 180] • interval.

Analytical PES
The ab initio energies were fitted to the analytical function with the form where the angular part is represented by a product of normalized associated Legendre functions and a cosine function. Only even values of l 1 are included in this function due to the symmetry of H 2 . For each value of R, the ab initio energies were fitted using a least-square method. The coefficients f l 1 ,l 2 m were fitted using the Reproducing Kernel Hilbert Space procedure (Ho & Rabitz 1996)  1. Rotational de-excitation cross-sections (in a 2 0 ) of HCO + by collision with para-H 2 using one (j H 2 = 0) and two (j H 2 = 0, 2) rotational states of H 2 in the calculations at the collisional energies of 10 and 200 cm −1 for selected transitions.
where N is the number of radial points of the grid R k . The q 2, 2 (R, R k ) is the one-dimensional kernels defined by Ho & Rabitz (1996), which ensures the long-range R −3 behaviour (Soldán & Hutson 2000). The α l 1 ,l 2 ,m k coefficients were found by solving the system of linear equations q(R i , R j )α = f l 1 ,l 2 m (R i ), where i and j label the different radial geometrical configuration of the grid.

Scattering
The dynamics was performed with the DIDIMAT code (Guillon et al. 2008) employing the explicitly correlated four-dimensional PES developed in this work. In this code, the close-coupling equations were solved in the space-fixed frame using the logderivative propagator (Manolopoulos 1988) for collisional energies from 10 −2 up to 1000 cm −1 . The rotational constants of H 2 were fixed to B H 2 = 60.853 cm −1 and for HCO + this value was B HCO + = 1.4875 cm −1 . The minimum propagation distance used was 50 a 0 , and the convergence was checked as a function of the total angular momentum and maximum propagation distance for each collision energy.
The basis set used in the close-coupling calculations included 32 rotational states of HCO + . Table 1 shows the influence of the number of rotational states of H 2 included in the basis set. At 10 cm −1 , and for the largest j transitions, differences can seem to be relatively large. They are, however, expected in this region where the resonances have the strongest effect on the magnitude of the cross-sections. In contrast, the agreement of the cross-sections at 200 cm −1 , which is still below the well depth, is much better as the contribution of the resonances in the magnitude of the cross-section is decreasing when the collision energy increases. Consequently, the effect of the limitation of the basis set of H 2 to its lowest rotational state only marginally affects the state-to-state rate coefficients as also found by Massó & Wiesenfeld (2014). Therefore, only one rotational state of H 2 was included in the basis set.

R E S U LT S
The large number of ab initio energies, together with the use of the least-square and kernel procedure, allows a fit of high quality. For negative energies, the root mean square error deviation (RMSD) is 1.99 × 10 −2 cm −1 , while for energies in the interval 0 < E < 1000 cm −1 , the RMSD is 7.78 × 10 −2 cm −1 . Fig. 3 shows the contour plots for the HCO + + H 2 complex at several configurations. The global minimum of the analytical PES, named −1509.35 cm −1 , was found at R = 3.44 Å, θ 1 = 90 • , and θ 2 = 0 • . This well is slightly deeper (1.46 per cent) than the one reported by Massó & Wiesenfeld (2014) (−1487.5 cm −1 ), which was found at the same angular configuration and R = 3.49 Å. Such a small different value in the minimum of the PES arises from the fact that we used the explicitly correlated CCSD(T)-F12a method instead of the standard CCSD(T) used by Massó & Wiesenfeld (2014). This value is also in good agreement with previous studies for the HCO + + H 2 complex. Bieske et al. (1995) found a value of −1443 cm −1 from the quadratic configuration interaction method, including single and double excitations with perturbative treatment of triples [QCISD(T)], while Li et al. (2008) reported a well of −1482 cm −1 computed with the CCSD(T) method. However, the PES from Monteiro (1985), computed only at the SCF level, and especially from Yazidi et al. (2014), that used a CCSD(T) PES averaged over the H 2 orientation exhibits much larger differences (−168 and −913 cm −1 , respectively). If our surface is averaged over the H 2 orientation, the global minimum (−971 cm −1 ) is also deeper than the one of Yazidi et al. (2014), but shows similar anisotropy. This difference is because the basis set used by Yazidi et al. (2014) produces results comparable to the energies computed at the CCSD(T)/aug-cc-pVQZ level, while in our case the energies are close to the CBS limit, as can be seen from Fig. 2.
The PES previously presented was employed for determining the cross-sections for the lowest rotational states of HCO + (up to j = 21) in collision with both para-and ortho-H 2 . The state-to-state rate coefficients were determined by the average of the corresponding cross-sections over a Maxwell-Boltzmann distribution of the collision energy.
The rotational rate coefficients computed here for HCO + + H 2 and those computed by Monteiro (1985), Flower (1999), Yazidi et al. (2014), andMassó &Wiesenfeld (2014) are shown in Table 2. The agreement with the calculations of Yazidi et al. (2014) is excellent. This is a new confirmation that the use of an averaged PES is a good approximation for studying the collision of a linear molecule with para-H 2 (j = 0). It is worth noting that the averaged PES approximation employed by Yazidi et al. (2014) transforms a fourdimensional PES to a two-dimensional surface by an average over the angular orientation of H 2 , and the collision with para-H 2 (j = 0) is then treated as the scattering of a linear rigid rotor with an atom. This is justified because in the dynamics of the collision with para-H 2 (j = 0), the coupling matrix elements are non-zero only if l 1 = 0 (Cabrera-González, Páez-Hernández & Denis-Alpizar 2020). The main limitation of such approximation is that it is valid only for H 2 (j = 0). Also, from Table 2, it can be seen that the rate coefficients for the collision with para-and ortho-H 2 are very close in magnitude. Fig. 4 shows the maximum percentage difference between both sets of rate coefficients at all | j| considered in this work. The behaviour of the per cent difference is almost independent of the temperature, and even if it increases after j = 13, the values remain very close between them. The ratio between the rate coefficients for HCO + + ortho-H 2 and HCO + + para-H 2 varies from 0.75 up to 2.14. This similarity in the rate coefficients or cross-sections for collisions with para-and ortho-H 2 has also been observed for other molecular ions, e.g.   Table 2. Rotational de-excitation rate coefficients (×10 −10 cm 3 molecule −1 s −1 ) of HCO + by collision with para-H 2 (p-H 2 ) and ortho-H 2 (o-H 2 ) at several temperatures. Rate coefficients for HCO + + para-H 2 from (Monteiro 1985) (M85), (Flower 1999 Walker et al. (2017) pointed out that the similarities come from the fact that the effects in the long-range interaction outweigh those from the short range. However, in the study of C 3 N − + H 2 , Lara-Moreno et al. (2019) found that the similarities decrease at low temperatures, where the long range dominates the dynamics. These authors showed that the short-range interaction appears to be dominated by the attractive isotropic A 0 00 term that gives non-zero contributions in the matrix element of the potential for collisions involving both para-H 2 (j = 0) and ortho-H 2 (j = 1) (Lara-Moreno  et al. 2019). The case of HCO + + H 2 is particular because the similarities are observed at both the lower and larger temperatures analysed, which can be explained by a combination of the arguments mentioned above, one acting at low temperatures and the other at high temperatures.
A small modification in the PES allows studying the collision of DCO + + H 2 . The coordinates for DCO + + H 2 differ from those for HCO + + H 2 only by small displacement of the centre of mass in the DCO + molecules. The old coordinates were expressed as a function of the new one from geometrical considerations (Sultanov, Guster & Adhikari 2012), and a new grid of energies was generated and fitted using the same expression 1. The influence of this correction in the PES is small (see Fig. 5), and it is not expected to play an essential role in the dynamics.
The modified surface was employed in close-coupling calculations. The same parameters used in the dynamics for HCO + + H 2 were employed. The only different data used were the mass of deuterium and the rotational constant of DCO + named B DCO + = 1.201 cm −1 (Caselli & Dore 2005). Fig. 6 shows the rotational de-excitation rate coefficients for the collision of DCO + with para-and ortho-H 2 from j = 5. The difference between both sets of rate coefficients is also very small, Figure 6. De-excitation rate coefficients of DCO + with para-H 2 (solid lines) and ortho-H 2 (dashed lines) from j = 5. The rate coefficients for HCO + + para-H 2 (dash-dotted line) are also included. like in the case of HCO + , as expected due to the very small differences between both PESs. The rate coefficients for HCO + + para-H 2 are also included in Fig. 6, and the differences with DCO + are very small. For all the rate coefficients computed in this work, the ratio between the DCO + + para-H 2 and HCO + + para-H 2 rate coefficients varies from 0.44 up to 1.98. Using the HCO + + para-H 2 rate coefficients for reproducing those of DCO + + H 2 then appears to be approximately valid. The DCO + + para-H 2 rate coefficients that can be found in the supplementary materials should, however, be preferred for analysing the DCO + observations as these differences are not completely negligible.

S U M M A RY
A new PES for the HCO + + H 2 complex at the CCSD(T)-F12a/augcc-pVTZ level of theory was developed in this work. Close-coupling calculations were performed using this PES. A new set of rate coefficients for HCO + + para-H 2 was reported and compared with previous calculations. Furthermore, the first set of rate coefficients for a large number of rotational states of HCO + in collision with ortho-H 2 was determined. Similarities between the rate coefficients for collisions with para-and ortho-H 2 were found. The collision of DCO + with H 2 was also investigated, and a set of rate coefficients for DCO + + para-H 2 was also reported. These rate coefficients are very close to those for HCO + + H 2 , which justifies the use of the latter to analyse the observations of the DCO + .

AC K N OW L E D G E M E N T S
Computer time for this study was provided by the Mésocentre de Calcul Intensif Aquitain, which is the computing facility of Université de Bordeaux et Université de Pau et des Pays de l'Adour. We acknowledge the support from the ECOS-SUD-CONICYT project number C17E06 (Programa de Cooperación Científica ECOS-CONICYT ECOS170039).

DATA AVA I L A B I L I T Y
The data underlying this article are available in the article and in its online supplementary material.