Quantum mechanical study of the high-temperature $\mathrm{H}^+ + \mathrm{HD} \to \mathrm{D}^+ + \mathrm{H}_2$ reaction for the primordial universe chemistry

TL;DR

This study employs quantum mechanics to analyze the high-temperature H+ + HD reaction, providing detailed cross sections and rate coefficients relevant for primordial universe chemistry, highlighting the reaction's dependence on initial states and vibrational levels.

Contribution

It presents the first comprehensive quantum-mechanical calculations of state-to-state reaction rates for H+ + HD at high temperatures, considering various vibrational and rotational states.

Findings

Rate coefficients vary significantly with initial rotational states below 4000 K.

Reaction probability decreases with higher final vibrational levels.

Reaction dynamics depend on whether the process is endothermic or exothermic.

Abstract

We use the time-independent quantum-mechanical formulation of reactive collisions in order to investigate the state-to-state $\mathrm{H}^+ + \mathrm{HD} \to \mathrm{D}^+ + \mathrm{H}_2$ chemical reaction. We compute cross sections for collision energies up to 1.8 electron-volts and rate coefficients for temperatures up to 10000 kelvin. We consider HD in the lowest vibrational level $v=0$ and rotational levels $j=0$ to 4, and H$_2$ in vibrational levels $v'=0$ to 3 and rotational levels $j'=0$ to 9. For temperatures below 4000 kelvin, the rate coefficients strongly vary with the initial rotational level $j$, depending on whether the reaction is endothermic ($j\le 2$) or exothermic ($j\ge 3$). The reaction is also found less and less probable as the final vibrational quantum number $v'$ increases. Our results illustrate the importance of studying state-to-state reactions, in the context of the chemistry of the primordial universe.

Tables

Table 1. Table 1: Parameters of the intervals into which the energy grid is split: minimum E min subscript 𝐸 min E_{\mathrm{min}} , maximum E max subscript 𝐸 max E_{\mathrm{max}} and increment Δ ​ E Δ 𝐸 \Delta E in eV, as well as the maximum total angular momentum J max subscript 𝐽 max J_{\mathrm{max}} .

$(v,j)$	$E_{\min }$	$E_{\max }$	$\Delta E$	$J_{\max }$
$(0,0)$	0.031	0.1	0.003	28
$(0,0)$	0.11	0.3	0.01	42
$(0,0)$	0.31	1	0.03	62
$(0,0)$	1.1	1.8	0.1	75
$\left(0,j\ge 2\right)$	0.0034	0.01	$6\times {10}^{-4}$	24
$\left(0,j\ge 2\right)$	0.011	0.03	0.001	28

Equations

$${\sigma }_{vj,v^{\prime }{j}^{\prime }}(E)=\frac{\pi {\hslash }^2}{2\mu (E-E_{vj})(2j+1)}\sum _J(2J+1){\left|S_{vj,v^{\prime }{j}^{\prime }}^J(E)\right|}^2$$

$$k_{vj,v^{\prime }{j}^{\prime }}(T)=\sqrt{\frac{8}{\pi \mu {\left(k_{B}T\right)}^3}}{\int }_0^{+\infty }d{E}_c{\sigma }_{vj,v^{\prime }{j}^{\prime }}(E_c)E_c{e}^{-E_c/k_{B}T},$$

$$k(T)=\frac{1}{Z_{\mathrm{rot}}}\sum _{j=0}^4(2j+1)e^{-E_{vj}/k_{B}T}\sum _{j^{\prime }=0}^9{k}_{0j,0j^{\prime }}(T)$$

Full text

††pubyear: 2019††pagerange: Quantum mechanical study of the high-temperature $\mathrm{H}^{+}+\mathrm{HD}\to\mathrm{D}^{+}+\mathrm{H}_{2}$ reaction for the primordial universe chemistry–Quantum mechanical study of the high-temperature $\mathrm{H}^{+}+\mathrm{HD}\to\mathrm{D}^{+}+\mathrm{H}_{2}$ reaction for the primordial universe chemistry

1 Introduction

It is now well established that the molecules composing the primordial universe, that contained neutral or singly-ionised hydrogen (H), deuterium (D), helium (4He) and lithium (7Li), were out of thermodynamic equilibrium (Coppola et al., 2011, 2013, 2016). In other words, the populations in their rovibrational levels did not follow a Maxwell-Boltzman distribution. Therefore, in order to understand the physical and chemical evolution of the primordial gas, which led in particular to the formation of the first stars, it is necessary to have quantitative information about the chemical reactions involving different rovibrational levels of the reactants and the products, so-called state-to-state reactions (Galli & Palla, 2013; Bovino & Galli, 2018). In particular, rate coefficients are needed on a wide range of temperatures, up to a few thousand kelvin (K). Theoretical calculations of rate coefficients are welcome, since they allow for covering such a wide range, in a more accessible way than experimental measurements.

In this context, the molecules H2 and HD play a central role, as they act as coolants of the primordial gas, through the mechanism of collisional excitation to higher rovibrational levels followed by the spontaneous emission of a photon. This cooling process has a strong influence on the gravitational collapse leading to the first structures. Below 500 K, HD is the main coolant, due to its electric dipole moment (about $10^{-3}$ debye) and its smaller rovibrational spacings compared to H2 (see for instance Galli & Palla (2002); Ripamonti (2007); Kreckel et al. (2010)). Therefore, it is crucial to investigate the chemical reactions involving HD, especially with neutral and ionised atoms, in order to characterise the cooling dynamics. If the reaction with the most abundant species H has been widely studied (see for instance (Flower & Roueff, 1999; Ely et al., 2016; Desrousseaux et al., 2018)), very little is known about its ionic counterpart $\mathrm{H}^{+}+\mathrm{HD}$, which is the subject of the present article.

In the literature, the scarce results given on the $\mathrm{H}^{+}+\mathrm{HD}\to\mathrm{D}^{+}+\mathrm{H}_{2}$ reaction (Henchman et al., 1981; Smith et al., 1982; Millar et al., 1989; Gerlich & Schlemmer, 2002; Jambrina et al., 2012) generally belong to articles focused on the inverse reaction $\mathrm{D}^{+}+\mathrm{H}_{2}\to\mathrm{H}^{+}+\mathrm{HD}$ (Fehsenfeld et al., 1974; Gerlich, 1992; Jambrina et al., 2009; Honvault & Scribano, 2013; González-Lezana et al., 2013; Sahoo et al., 2014, 2015; Lara et al., 2015; Bhowmick et al., 2018). In Refs (Henchman et al., 1981; Smith et al., 1982), thermal rates of $1.1\pm 0.2$ and $1.7\pm 0.2\times 10^{-10}$ cm3.s*-1* are measured at 205 and 295 K respectively using the selected ion flow tube (SIFT) technique. Millar et al. (1989) derived an Arrhenius-type formula from those measurements. Later, Gerlich & Schlemmer (2002) also gave an Arrhenius-type formula fitted from most dynamically biased (MDB) statistical calculations between 30 and 130 K. Finally, Jambrina et al. (2012) give thermal rate coefficients calculated both with a time-independent quantum-mechanical (TIQM) method and with variants of the quasi-classical trajectory (QCT) method on the so-called ARTSM potential-energy surface (PES) by Aguado et al. (2000). The agreement of the TIQM results with the above-mentioned experiments is good, even if the latter are better reproduced by the statistical QCT calculations at 295 K. It is worthwhile noting that the reactive collision between Rydberg hydrogen atoms and HD has been investigated experimentally (Yu et al., 2014).

The weak interest to the $\mathrm{H}^{+}+\mathrm{HD}\to\mathrm{D}^{+}+\mathrm{H}_{2}$ reaction is probably due to its endothermiciity (for HD and H2 in their ground rovibrational levels $(v,j)=(0,0)$ and $(v^{\prime},j^{\prime})=(0,0)$). This endoothermicity, equal to 39.5 meV, mainly comes from the difference in zero-point energies between H2 and HD (35.8 meV), but also from the difference in ionisation potentials of H and D (3.7 meV). But firstly, this activation energy is widely overcome at the temperatures that we consider here; and secondly, the title reaction becomes exothermic for HD in rovibrational levels higher than $(v=0,j=3)$. On the other hand, for electronic energies larger than 1.83 eV above the dissociation limit $\mathrm{H}^{+}+\mathrm{H}_{2}(X^{1}\Sigma_{g}^{+})$, charge transfer becomes possible towards the channel $\mathrm{H}(^{2}S)+\mathrm{H}_{2}^{+}(X^{2}\Sigma_{g}^{+})$. However, time-dependent wave-packet (TDWP) calculations for $\mathrm{D}^{+}+\mathrm{H}_{2}$ (Ghosh et al., 2015) and $\mathrm{H}^{+}+\mathrm{H}_{2}$ (Ghosh et al., 2017) on the three lowest diabatic PESs of H${}_{3}^{+}$ (Viegas et al., 2007) have shown that the charge transfer processes are much less probable than the reactive one without charge transfer.

In this article, we compute the cross sections and the rate coefficients characterizing the reaction $\mathrm{H}^{+}+\mathrm{HD}(v,j)\to\mathrm{D}^{+}+\mathrm{H}_{2}(v^{\prime},j^{\prime})$, with HD in the lowest vibrational level $v=0$ and rotational levels $j=0$ to 4, and H2 in vibrational levels $v^{\prime}=0$ to 3 and rotational levels $j^{\prime}=0$ to 9. Our calculations are performed with the TIQM method for reactive collisions, based on hyperspherical coordinates, which take into account the indistinguishability of the two H nuclei. Within the framework of the Born-Oppenheimer approximation, we characterize the motion of the two H and the D nuclei on the so-called VLABP ground state global PES of H${}_{3}^{+}$ calculated by Velilla et al. (2008), an improved version of the ARTSM PES (Aguado et al., 2000) that cautiously takes into account long-range interactions. Because we ignore the hyperfine interactions, we do not account for the difference in ionisation potentials between H and D, which introduces an uncertainty of 3 meV on collisional energy. But in this study, we are more interested in large collision energies, up to 1.8 eV above the lowest channel $\mathrm{H}^{+}+\mathrm{HD}(0,0)$, for which we can ignore charge transfer. Because our cross sections are very low at 1.8 eV, we can give converged rate coefficients for temperatures up to 10000 K for the lowest rovibrational levels of H2. We also calculate thermal rate coefficients, assuming HD rovibrational levels in thermodynamic equilibrium, in order to compare our results with literature, and find larger rates. We possibly attribute the discrepancies with previous TIQM results to differences in the asymptotic region of the underlying PESs. Regarding previous experimental results, rate coefficients are given at only two temperatures, and so we think that additional measurements would be particularly relevant.

This article is organised as follows. Section 2 describes our TIQM method, giving in particular the expressions of the cross section and rate coefficient. Section 3 presents our results, dealing with HD in the rovibrational ground level (Sec. 3.1), in rotationally-excited levels (Sec. 3.2), and in thermodynamic equilibrium (Sec. 3.3), to allow comparison with literature results. Section 4 contains concluding remarks.

2 Method

In this article, we focus on two quantities characterizing the reactive collision $\mathrm{H}^{+}+\mathrm{HD}(v,j)\to\mathrm{D}^{+}+\mathrm{H}_{2}(v^{\prime},j^{\prime})$: the state-to-state cross section $\sigma_{vj,v^{\prime}j^{\prime}}$ given as a function of the total energy $E$,

[TABLE]

and the state-to-state rate coefficient $k_{vj,v^{\prime}j^{\prime}}$ given as a function of temperature $T$,

[TABLE]

where $\hbar$ is the reduced Planck’s constant, $k_{B}$ is Boltzman’s constant, $\mu$ is the reduced mass of the reactants, $J$ is the total angular momentum of the H2D+ system, $S_{vj,v^{\prime}j^{\prime}}^{J}(E)$ are the elements of the scattering matrix, and $E_{c}=E-E_{vj}$ is the collision energy in the entrance channel $\mathrm{H}^{+}+\mathrm{HD}(v,j)$.

The scattering matrix is calculated for a given total energy $E$ and total angular momentum $J$ using a fully Coriolis-coupled TIQM method based on the body-fixed democratic hyperspherical coordinates, and described in details in (Honvault & Launay, 2004). At each hyperradius $\rho$, the scattering wave function is expanded on a set of appropriate hyperangular basis functions. The $\rho$-dependent coefficients are solutions of a set of coupled second-order differential equations, which are solved using the Johson-Manolopoulos log-derivative propagator (Manolopoulos, 1986). The scattering wave function is computed up to the hyperradius $\rho_{\mathrm{max}}=17.5\,a_{0}$, where the $S_{vj,v^{\prime}j^{\prime}}^{J}(E)$ matrix elements are extracted for many rovibrational levels of the reactant HD and the product H2. This method has been successfully applied to the isotopic variants $\mathrm{H}^{+}+\mathrm{H}_{2}$ (Honvault et al., 2011b, a; Rao et al., 2014; González-Lezana et al., 2014; González-Lezana & Honvault, 2017), $\mathrm{D}^{+}+\mathrm{H}_{2}$ (Honvault & Scribano, 2013; González-Lezana et al., 2013, 2014; Lara et al., 2015) and $\mathrm{H}^{+}+\mathrm{D}_{2}$ (González-Lezana et al., 2009). Here we use the VLABP PES calculated by Velilla et al. (2008), and which accurately describes the long-range interactions between H+ and H2.

Because the charge-transfer channel $\mathrm{H}+\mathrm{H}_{2}^{+}$ is located at an electronic energy of 1.83 eV above $\mathrm{H}^{+}+\mathrm{H}_{2}$, we perform our scattering calculations up to 1.8 eV above the lowest entrance channel $\mathrm{H}^{+}+\mathrm{HD}(v=0,j=0)$. Our energy grid, which contains about 80 points, is denser for low collision energies. The energy grid is split into intervals inside which we take the same maximum total angular momentum $J_{\mathrm{max}}$ (see Table 1). Regarding the ground rovibrational level of HD, the lowest collision energy that we take is 0.031 eV, for which the reaction is impossible, but we see the opening of the reactive channel at 0.037 eV. For each excited rotational level $J\geq 2$, we also consider collision energies down to 0.003 eV, a value below which our matching distance $\rho_{\mathrm{max}}=17.5\,a_{0}$ does not allow for a satisfactory convergence. On the other hand, to get a good convergence for an energy of 1.8 eV above the $\mathrm{H}^{+}+\mathrm{HD}(v=0,j=0)$, we need to include 252 and 98 (open and closed) rovibrational levels for HD and H2. In the case of HD, $j$ ranges from 0 to 29 for $v=0$, and $v$ ranges from 0 to 11 for $j=0$. As a consequence, we get 350, 2875 and 3892 coupled channels for $J=0$, 10 and 20, respectively.

3 Results

In this section, we present cross sections and rate coefficients for selected initial (HD) and final (H2) rovibrational levels. The collision energies are expressed in electron-volts (eV), cross sections in units of squared Bohr radius $a_{0}^{2}$, temperature in kelvin (K), and rate coefficients in cm3.s*-1*.

3.1 HD in $(v=0,j=0)$

Figure 1 shows the cross sections of the reactions $(v=0,j=0)\to(v^{\prime}=0,j^{\prime})$ for $j^{\prime}=0$ to 9, as functions of the collision energy. For all those endothermic reactions, the logarithmic scale enables us to see the opening of the reactive channel associated with each rotational level $j^{\prime}$ of H2. For energies slightly above that threshold, the cross sections abruptly increase, while they slowly decrease at high energies, down to values smaller than 2 $a_{0}^{2}$ at 1.8 eV. In addition, the curves present resonant peaks for energies lower than 0.5 eV. Those peaks are higher and more numerous for the low rotational levels of H2: for instance, the highest peak is observed for $j^{\prime}=1$ (54 $a_{0}^{2}$ at 0.085 eV). The curve with $j^{\prime}=0$ also presents a narrow and high resonance (49 $a_{0}^{2}$ at 0.094 eV); but that curve is generally below the $j^{\prime}=1$ one.

Figure 2 shows the rate constants as functions of the temperature between 100 and 10000 K, for the same reactions as Fig. 1. All the curves evolve in a similar way with temperature, showing a fast increase, a maximum, and then a slow decrease. The low-temperature increase is faster with smaller $j^{\prime}$, due to the lower energy thresholds observed on Figure 1. Moreover, as $j^{\prime}$ increases, the high-temperature decrease in $k(T)$ is less pronounced; for example for $j^{\prime}=9$, the bump is hardly visible since the maximum rate coefficient is equal to $6.06\times 10^{-11}$ cm3.s*-1* (at 8700 K), while it is equal to $5.98\times 10^{-11}$ cm3.s*-1* at 10000 K.

For levels $j^{\prime}=5$ to 9, the rate coefficients are smaller on the whole range of temperatures of Fig. 2(b) when $j^{\prime}$ increases. But for lower values of $j^{\prime}$, the hierarchy is less clear. Except for $T<300$ K, the reaction towards $j^{\prime}=0$ is never the dominant one, as expected with cross sections (see Fig. 1). The highest rate is obtained for $j^{\prime}=1$ ($3.33\times 10^{-10}$ cm3.s*-1*at 1450 K), but for higher temperatures, this rate becomes smaller than those for $j^{\prime}=2$ to 5. At the temperature of 6000 K for instance, the dominant reactions are toward $j^{\prime}=3$, 2 and 4.

It is important to stress that the high-temperature decrease of Fig.  2 is not due to a bad convergence of the integral in Eq. (2). To check it, we have removed the last two cross sections (at 1.7 and 1.8 eV) in the calculation of rate coefficients. At 10000 K, we found a relative difference of 0.7, 1.2 and 8.7 % for $j^{\prime}=0$, 3 and 9, respectively. The convergence is less good for $j^{\prime}=9$ because the cross section at 1.8 eV, equal to 1.46 $a_{0}^{2}$ represents one third of its maximum value (4.30 $a_{0}^{2}$ at 0.91 eV); its removal has thus a sizeable influence on the integral (2).

Now we consider the product H2 in its first excited vibrational level $v^{\prime}=1$. The cross sections are plotted in Fig. 3, in an energy range between 0.5 and 1.8 eV and in linear scale. In that range, the energy resolution does not allow for distinguishing the thresholds of the reactions $(0,0)\to(1,0)$ and $(0,0)\to(1,1)$; but it is sufficient in order to understand the overall evolution of the cross section, which is relevant for the calculation of rate coefficients. Globally, the cross sections for $v^{\prime}=1$ are much smaller than for $v^{\prime}=0$ (see Fig. 1. This trend is also visible for the rate coefficients, plotted in Fig. 4. For example at 1000 K, the rate coefficients for $v^{\prime}=1$ are on the order of $10^{-13}$ cm3.s*-1*or below, that is three orders of magnitude lower than those for $v^{\prime}=0$. This gap decreases with temperature, and at 10000 K, there is approximately a one-order-of-magnitude difference between $v^{\prime}=0$ and $v^{\prime}=1$. As for the $j^{\prime}$-dependence of the rates, at 6000 K, the dominant reactions involve intermediate $j^{\prime}$-values (3, 4 then 5), while $j^{\prime}=0$ has the second smallest rate.

We have also checked the convergence of the rate coefficients given in fig. 4 by removing the last two cross sections of Fig. 3. This reduces the rate coefficients by 5.9 % and 42 % for $j^{\prime}=0$ and 9 respectively. The latter value is thus less accurate, but one should keep in mind that corresponding reaction is by far not the dominant one driving to the destruction of HD($v=0,j=0)$.

When $v^{\prime}=2$ and 3, the rate coefficients are even smaller, as shown in fig. 5: at 6000 K, the rates lose one order of magnitude when $v^{\prime}$ increases by one unity. Regarding the convergence of equation (2), for $v^{\prime}=3$, there are only five non-zero cross sections, from 1.4 to 1.8 eV; removing the last value decreases the rate coefficients by 38 %.

3.2 HD in excited rotational levels $(v=0,j)$

For the first rotationally excited level $j=1$, the results, very similar to the previous ones, are given in the supplementary material. In particular, at 6000 K, the largest reactions rate involve $j^{\prime}=3$, 4 then 2. The case $j=2$ is a little particular, since the reaction $\mathrm{H}^{+}+\mathrm{HD}(v=0,j=2)\to\mathrm{D}^{+}+\mathrm{H}_{2}(v^{\prime}=0,j^{\prime}=0)$ is endothermic by only 3.4 meV (neglecting the difference in ionisation energies between H and D). So for a given $j^{\prime}$-value, the energy thresholds for $(0,2)\to(0,j^{\prime})$, see Fig. 6, are lower than for $(0,0)\to(0,j^{\prime})$, see Fig. 1, and the maximum of the curve has a smaller energy and a larger cross section. The consequence on the rate coefficients can be seen in Fig. 7: they are all larger for the reaction $(0,2)\to(0,j^{\prime})$ at low temperatures; but as $T$ increases, the ratio $k_{02,0j^{\prime}}(T)/k_{00,0j^{\prime}}(T)$ decreases, and becomes smaller than unity for the lowest $j^{\prime}$-values (see insets of Fig. 7).

The case $j=3$ is particularly interesting, since the reactions $(0,3)\to(0,0)$ and $(0,3)\to(0,1)$ are exothermic by 0.029 and 0.015 eV. Figure 8 shows the corresponding cross sections down to an energy of 3 meV, a value below which our value of $\rho_{\mathrm{max}}=17.5~{}a_{0}$ does not allow for a full convergence of the cross section. The reaction $(0,3)\to(0,1)$ strongly dominates over $(0,3)\to(0,0)$, for instance by a factor of 3 at 3 meV. This predominance is also visible on rate coefficients, see Fig. 9. On the whole range of temperature, $k_{03,01}$ is larger than $k_{03,00}$ by a factor of $\simeq 3$; but is exceeded by $k_{03,02}$ for $T\geq 400$ K and then by other rates. At large temperatures, the results are similar to other $j$-values: the rates slowly decrease, and at 6000 K, the dominant reactions are towards $j^{\prime}=4$, 3 and 5.

In order to check the convergence at low temperature, we have compared the previous rates with those calculated by removing the first two cross sections (at 3 et 3.4 meV) in equation (2). At 100 K, the rates $k_{03,00}$ and $k_{03,01}$ decrease by 6.5 % and 6.9 % respectively. At large temperature, the convergence is very similar to the one of $k_{00,00}$.

In the case $j=4$, there are three exothermic reactions, for $j^{\prime}=0$, 1 and 2. The cross sections and rate coefficients, which are given in the supplementary material, look like Figs. 8 and 9. At low energies, the cross section $\sigma_{04,02}$ is the largest ($169~{}a_{0}^{2}$ at 3 meV), followed by $\sigma_{04,01}$ ($128~{}a_{0}^{2}$) and $\sigma_{04,00}$ ($36~{}a_{0}^{2}$). This trend shows up on low-temperature rates, while the high-temperature rates are similar to the other $j$-values: at 6000 K, the largest rates are towards $j^{\prime}=4$, 3 and 5.

3.3 HD in thermodynamic equilibrium

In order to compare our results with literature, we have computed the thermal rate coefficient $k(T)$ of the title reaction, assuming that the rovibrational levels of HD are populated according to a Maxwell-Boltzman distribution. Since we deal with rather low temperatures in this paragraph, we can assume that only the five lowest rotational levels $(v=0,j\in[0;4])$ are significantly populated. For instance at 295 K, the population in $(0,0)$ is equal to 20.1 %, that of $(0,1)$ is 39.1 %, while those of $(0,4)$ and $(0,5)$ are 2.4 % and 0.4 % respectively. Moreover, as shown in Fig. 7, the reactions producing vibrationally excited H2 can safely be ignored. Therefore, our thermal rate coefficient is given by

[TABLE]

where $Z_{\mathrm{rot}}=\sum_{j=0}^{4}(2j+1)e^{-E_{vj}/k_{B}T}$ is the rotational partition function restricted to the five lowest rotational levels of HD.

By doing so, we find a thermal rate coefficient that can be very well fitted with an Arrhenius formula $A\times\exp(-B/T)$, between 100 and 400 K, with $A=1.84\times 10^{-9}$ cm3.s*-1* and $B=416$ K. More specifically, we obtain $k(T=205~{}K)=2.4$ and $k(T=295~{}K)=4.5\times 10^{-10}$ cm3.s*-1*, which are approximately twice as large as the measurements of Refs. (Henchman et al., 1981; Smith et al., 1982), and also significantly larger the rates calculated by Jambrina et al. (2012). In view of these discrepancies, we have checked the quality of our scattering matrices, by using them to calculate cross sections and rate coefficients of the reverse reaction $\mathrm{D}^{+}+\mathrm{H}_{2}\to\mathrm{H}^{+}+\mathrm{HD}$: our results are in very good agreement with former calculations (Honvault & Scribano, 2013), themselves in very good agreement with experiments for collision energies larger than 0.01 eV (Gerlich, 1992).

The discrepancies obtained with the TIQM calculations of Jambrina et al. (2012) may come from the PESs employed. Indeed, Rao et al. (2014) have shown that the cross sections and rate coefficients of the $\mathrm{H}^{+}+\mathrm{H}_{2}$ reaction are always significantly larger with the VLABP PES than with the KBNNPES by Kamisaka et al. (2002). But as shown in (González-Lezana et al., 2009) on the $\mathrm{H}^{+}+\mathrm{D}_{2}$ reaction, there is not such a heavy trend between the KBNN and ARTSM PESs. Extending those results to the title reaction thus suggests that the ARTSM and VLABP PESs are likely to give significantly different cross sections and rate coefficients. Because the calculations of Jambrina et al. (2012) were performed with the ARTSM PES, this may explain why our rate coefficients overcome those of Jambrina et al. (2012). In any case, additional measurements, with more collision energies or temperatures, would be necessary.

4 Conclusion

We have computed cross sections and rate coefficients of the state-to-state $\mathrm{H}^{+}+\mathrm{HD}(v,j)\to\mathrm{D}^{+}+\mathrm{H}_{2}(v^{\prime},j^{\prime})$ reaction, where HD is in the vibrational ground level. The rate coefficients are computed for temperatures from 100 to 10000 K. They are significantly smaller for vibrationally excited H2 molecules, approximately losing one order of magnitude when $v^{\prime}$ increases by one. For temperature below 4000 K, rate coefficients also strongly depend on the initial rotational level, as a function of which the reaction is endothermic or exothermic. Investing state-to-state $\mathrm{H}^{+}+\mathrm{HD}\to\mathrm{D}^{+}+\mathrm{H}_{2}$ reactions is therefore crucial to model the chemical evolution of the primordial gas. In a future work, we will consider HD in excited vibrational levels, for which the reaction is exothermic for all rotational levels.

In order to compare our results with the few literature ones, we have also calculated thermal rate coefficients, assuming the HD rovibrational levels in thermodynamic equilibrium. We find rate coefficients roughly twice as large as published values. Regarding previous theoretical results, we possibly attribute those discrepancies to differences in the asymptotic region of the potential-energy surfaces used for the calculations. In any case, more experimental results are necessary, since thermal rate coefficients have only been measured for two temperatures.

Acknowledgements

M.L. acknowledges the financial support of “Région Bourgogne Franche Comté” under the projet 2018Y.07063 “ThéCUP”. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.