State-to-State Quantum Mechanical Calculations of Rate Coefficients

Mar 1, 2013 - The dynamics of the D+ + H2 → HD + H+ reaction on a recent ab initio potential energy surface (Velilla, L.; Lepetit, B.; Aguado, A.; B...
0 downloads 0 Views 998KB Size
Download PDF
Article pubs.acs.org/JPCA

State-to-State Quantum Mechanical Calculations of Rate Coefficients for the D+ + H2 → HD + H+ Reaction at Low Temperature P. Honvault*,† and Y. Scribano‡ †

Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR CNRS 6303, Université de Bourgogne, 21078 Dijon Cedex, and UFR Sciences et Techniques, Université de Franche-Comté, 25030 Besançon cedex, France ‡ Laboratoire Univers et Particules de Montpellier, Université de Montpellier, LUPM-UMR CNRS 5299, 34095 Montpellier cedex 05, France ABSTRACT: The dynamics of the D+ + H2 → HD + H+ reaction on a recent ab initio potential energy surface (Velilla, L.; Lepetit, B.; Aguado, A.; Beswick, J. A.; Paniagua, M. J. Chem. Phys. 2008, 129, 084307) has been investigated by means of a time-independent quantum mechanical approach. Cross-sections and rate coefficients are calculated, respectively, for collision energies below 0.1 eV and temperatures up to 100 K for astrophysical application. An excellent accord is found for collision energy above 5 meV, while a disagreement between theory and experiment is observed below this energy. We show that the rate coefficients reveal a slightly temperature-dependent behavior in the upper part of the temperature range considered here. This is in agreement with the experimental data above 80 K, which give a temperature independent value. However, a significant decrease is found at temperatures below 20 K. This decrease can be related to quantum effects and the decay back to the reactant channel, which are not considered by simple statistical approaches, such as the Langevin model. Our results have been fitted to appropriate analytical expressions in order to be used in astrochemical and cosmological models. association, H + D → HD + hν,1 even if stimulated effects due to the background radiation field are considered.8 So the title reaction plays an important role in the chemistry of the early universe, especially in the deuterium fractionation. It is worthwhile to mention that it is also one of the main sources of HD in interstellar clouds.9−11 Direct efficient formation of HD on dust grains does not proceed because of the absence of dust grains at the considered epoch, which were only created later after the first stars appeared. In addition, at the appearance of first grains, the D+ + H2 → HD + H+ process is still the main mechanism of HD formation helped by the fact that the H2 abundance is strongly boosted by grain surface reactions for all considered metallicities.12 The reverse reaction HD + H+ → D+ + H2 is endothermic in contrast with the title reaction, which is exothermic by 39.37 meV (456.9 K) and thus has no energy threshold. So, at low temperatures ( 0 for D+ + H230,31 for collision energies from 1.7 to 2.5 eV on new diabatic PES.32 Classical dynamical computations employing the trajectory surface hopping method were also carried out31,33 for very high collision energies (>2.5 eV). Given its cosmological importance, we have carried out a theoretical study of dynamics of the D+ + H2 (v = 0, j = 0,1) → HD + H+ reaction by means of a TIQM approach on an updated version34 of the ab initio PES of Aguado et al.16 Total ICSs have been calculated as a function of the collision energy up to 0.1 eV. At this energy and below, only the rovibrational states (v = 0, j = 0,1) of H2 and (v′ = 0, j′ = 0, 1, 2, 3, 4) of HD can be populated. State-to-state ICSs have also been reported for D+ + H2 (v = 0, j = 0) → HD (v′ = 0, j′) + H+ and D+ + H2 (v = 0, j = 1) → HD (v′ = 0, j′) + H+. Rate coefficients are then calculated up to 100 K using an integration of ICSs over a Maxwell−Boltzmann distribution of collision energy. Here, we only consider collision energies well below 1.8 eV, so D+ + H2 → HD + H+ is the only reactive channel, and we do not study any nonadiabatic effects. Though the charge transfer channel leading to H2+ + H(D) is not accessible at the energies involved in the study (much lower than 1.8 eV where it becomes open), nonadiabatic interaction with the charge transfer channel will modify the energy level spectrum of H3+ (H2D+), and it may influence the dynamics. The present work follows a recent study by our group of the ortho−para conversion of H2 by proton exchange,35−38 and to the best of our knowledge, it gives the first quantum mechanical total and state-to-state rate coefficients at low temperature for D+ + H2.

cooling requires a complete set of the state-to-state rate coefficients as well as the total rate coefficient. Besides its importance in astrophysics, the D+ + H2 reaction has been investigated by several theoretical and experimental groups because it appears as a prototype for dynamical studies of ion−molecule reactions. In 2005, Gonzalez-Lezana et al.15 performed time-dependent wave packet (TDWP) and statistical quantum-mechanical (SQM) calculations on an ab initio potential energy surface (PES) published in 2000 for the H3+ system16 to study the reactive noncharge transfer ion− molecule collisions D+ + H2. Reaction probabilities have been computed up to 1.4 eV. They found that the SQM results reproduce very well the reaction probability at a zero total angular momentum J, while for higher J some discrepancies are found. Integral cross-sections (ICS) and differential crosssections have also been obtained using different approximations, such as the centrifugal sudden approach. In 2010, Jambrina et al.17 reported time-independent quantum-mechanical (TIQM) calculations using the close-coupled hyperspherical method implemented in the ABC code18 for D+ + H2 on the same PES reported by Aguado et al.16 They compared the TIQM reaction probabilities at J = 0 and J > 0 and ICS with those obtained from two approximate methods: quasi-classical trajectory (QCT) and statistical quasi-classical trajectory (SQCT). Very recently, Jambrina et al.19 compared the results obtained from the three methods mentioned above with available experimental measurements. For the D+ + H2 reaction, the rate coefficients were computed in the 100−500 K temperature range and compared with the measurements done at 205 and 295 K with normal H2 using the temperatureselected ion flow tube (SIFT)20,21 or a flowing afterglow22 technique. The theoretical rate coefficients are in good agreement with experiment. They reproduce the decrease with increasing temperature between 200 and 300 K. The cross-sections multiplied by the relative velocity were also presented between 1 meV and 1.2 eV for H2 (v = 0, j = 0,1,2,3). They have been averaged over the thermal H2 rotational state distribution at Trot = 300 K, and a good agreement was obtained with the experimental results.23,24 About the experimental results, Gerlich measured23 the rate coefficients for D+ + (normal)H2 based on the merged beams technique for temperatures between 180 and 350 K. He found a rate that is nearly constant in this temperature range with a value of 1.6 × 10−9 cm3 molecule−1 s−1. An analytical expression for the rate coefficient between 30 and 130 K have been derived by Gerlich25,26 from the most dynamically biased (MDB) statistical theory.25 In addition, Gerlich23 also presents experimental cross-sections multiplied by the relative velocity for collision energy down to 1 meV. McCarroll,27 assuming again that the reaction proceeds via the formation of a long-lived complex, employed a statistical mixing model including the nuclear symmetry constraint to compute the rate coefficient between 10 and 400 K for ortho, para, and normal H2. The results are in excellent agreement with the measurements at 205 and 295 K,20−22 and the dependence of the rate coefficients on the temperature is well reproduced. Nonadiabatic dynamical calculations to investigate the charge transfer mechanism were undertaken above the crossing with the first excited PES by means of exact time-independent closecoupling formalism with hyperspherical coordinates at J = 028,29 on different PESs built from the diatomics-in-molecule

II. TIME-INDEPENDENT QUANTUM METHOD Although the calculations for J > 0 are very expensive in CPU time and memory due to the presence of a deep well (4.60 eV) in the PES, we have carried out three-dimensional reactive scattering calculations using a fully Coriolis-coupled TIQM method based on the body-frame hyperspherical democratic coordinates approach. An important advantage of this method is that it yields the entire S matrix at a given energy, and so state-to-state scattering attributes, such as cross-sections and rate coefficients, can be easily extracted. TIQM methods are also best suited for the calculations at low collision energies, in contrast with the TDWP methods because of the large de Broglie wavelengths involved in this energy regime. In addition, TIQM methods are appropriate for an accurate study of complex-forming reactions. On the contrary, TDWP calculations are more difficult and even unfeasible in some cases for dynamics involving long-lived resonances because the propagation time can become very long. All these aspects are presented in detail in a recent review on the complex forming reactions.39 The hyperspherical method was presented in detail in ref 40, and only a brief description is given here. It has already been used for the isotopic variant H+ + H2 as mentioned above. This method has also previously proved successful in describing complex-forming reactions, such as the OH + O,41 H + O2,42 and OH + N43,44 reactions and ultracold alkali-dialkali collisions.45 At each hyper-radius, the scattering wave function is expanded on a set of hyperspherical adiabatic states of a reference Hamiltonian H = T + V, which incorporates the kinetic energy T arising from deformation at fixed hyper-radius and the potential energy V. The expansion coefficients are the B

dx.doi.org/10.1021/jp3124549 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

solution of a set of coupled second-order differential equations, which are solved using the Johnson−Manolopoulos logderivative propagator.46 In this work, for J = 0, 200 states dissociate at large hyper-radius into the H2(20,18,16,14,12,10,6) rovibrational set (this notation indicates the largest rotational level j for each vibrational manifold v = 0,1,...,6) and the HD(25,23,20,18,16,14,11,8,1) rovibrational set for v′ = 0,1,...,8. So a large number of closed channels needed for convergence has been included in the present study. Propagation goes from 0.5a0 up to the asymptotic matching distance of 30a0, where the S matrix is extracted. When computing J > 0 partial waves, we have considered the Ω components from 0 to 25 in the close-coupling expansion to obtain accurate ICS. Thus, the number of coupled equations increases from 200 for J = 0 to 1938 for J > 25. In addition to the inclusion of many closed channels at 0.1 eV (the highest collision energy considered in the present work), we checked convergence using larger basis sets and larger asymptotic distances. The main errors might come from the PES used in the reactive scattering calculations. As already mentioned, we used a recent high quality ab initio PES of the ground electronic state of H3+.34 It is invariant under all permutations of the nuclei and includes the long-range electrostatic interactions analytically. New refinements of a previous PES16 have been added in the present PES, such as more ab initio points and the inclusion of a functional form of the long-range electrostatic interaction in the analytical representation. So, the new PES should properly include the long-range potential in the H(D)+ + H2 channel, which is crucial for quantum dynamics calculations at low collision energies. The present PES is therefore the one assumed as more accurate and should give the best possible description of the existing potential interactions.

Figure 1. Total rate coefficients as a function of the translational energy Et for D+ + H2 (v = 0, j = 0) → HD + H+ (black line) and D+ + H2 (v = 0, j = 1) → HD + H+ (red line). The measurements found by Gerlich23 are also shown (blue circles). The lower panel presents the same results with a logarithmic scale for Et.

III. RESULTS AND DISCUSSION We have computed TIQM ICS for the reaction D+ + H2 (v = 0, j) → HD (v′ = 0) + H+ for the two initial rotational states j = 0,1 and integrated over all final rotational states of HD. Figure 1 (upper panel) shows the energy dependence of the rate coefficient kj(Et) obtained by multiplying the TIQM crosssections for D+ + H2 (v = 0, j = 0) or H2 (v = 0, j = 1) with the relative velocity. Et denotes the translational energy. We can see that the j = 0 and j = 1 initial rotational states have a similar contribution to the total rate coefficient. The TIQM rate coefficients reveal a highly oscillating behavior with Et, which indicates a large number of resonances. We can associate the observed narrow peaks to the formation of long-lived intermediate complex, which can be formed in the deep well (4.60 eV) of the H2D+ PES. The TIQM rate coefficient is therefore very sensitive to the translational energy Et. However, in average it keeps a constant value around 1.7 × 10−9 cm3 molecule−1 s−1 down to about 0.01 eV, and it decreases below this energy as Et decreases. It is worthwhile to recall that a simple capture model gives a dependence of the cross-section in Et−1/2 for an ion−neutral molecule reaction, yielding a constant rate coefficient. However, at very low energy, the capture model can lead to wrong results as already found in the alkalin−dialkalin collisions (see for instance Figure 4 in ref 47). Gerlich measured rate coefficients23,26 for Et between 1 meV and 1 eV using a merged beam technique with normal H2. As can be seen in the upper panel of Figure 1, for Et > 0.006 eV, the measured rate coefficients are in excellent agreement for both the magnitude and the energy dependence with the

present theoretically determined values, even if the resonance structure is not seen in experiment. However, below 0.006 eV, the experimental and theoretical results predict totally different dependences on Et . The experimental rate is still independent of Et, and in the same range, it is above the theoretical results. At 1 meV (the lowest Et in experiment), the TIQM rate coefficient is a factor of 3 smaller than that predicted by experiment. This result needs to be confirmed by new experiments especially at very low energies. In addition, this discrepancy below 0.006 eV may originate from a possible shortcoming of the PES used in the present study. The computed rate coefficient, which is too small at very low energy, might be the consequence of a long-range potential, which is not sufficiently attractive, although the analytical representation of the PES included an explicit model for the long-range interaction. Dynamical effects could also exist and other parts of the PES may be involved. Moreover, as mentioned above, the D+ + H2 → HD + H+ reaction is exothermic by 39.37 meV (456.9 K) and by 36.58 meV if we include the zero-point energies of H2 and HD only. As we employed a PES built for the H3+ system, this PES does not take into account the difference in the ionization potentials of H and D. Figure 1 (lower panel) presents the same data but using a log-scale for Et. This figure clearly shows that, after a decrease C

dx.doi.org/10.1021/jp3124549 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

down to 10−5 eV, the TIQM rate coefficient becomes independent of Et with a constant value. This behavior is expected from the Wigner threshold laws.48,49 For D+ + H2 below 10−5 eV, only the s-wave (l = 0) contributes to the reactivity and the rate coefficient, which varies, as Elt must be constant for l = 0. So the TIQM result has the correct dependence on Et at very low energy. Above 10−5 eV, other partial waves contribute yielding another dependence, and for Et > 0.006 eV, the large number of partial waves give the expected behavior for an ion−molecule reaction, i.e, a rate coefficient that is a constant. TIQM state-to-state ICS for the reaction D+ + H2 (v = 0, j) → HD (v′ = 0, j′) + H+ are reported in Figures 2 and 3 for,

Figure 4. State-to-state rate coefficients as a function of the temperature for D+ + H2 (v = 0, j = 0,1) → HD (v′ = 0, j′) + H+ for j′ = 0,1,2,3.

highest rate coefficient is obtained for the transition group j = 0,1 → j′ = 0,1,2. All the reaction rate coefficients slightly depend on the temperature and are located between 0.9 × 10−10 cm3 molecule−1 s−1 and 0.8 × 10−9 cm3 molecule−1 s−1. The decrease with decreasing temperature is more pronounced for temperatures below 20 K. The values of rate coefficients are highly below the Langevin rate coefficient, 2.1 × 10−9 cm3 molecule−1 s−1. The TIQM rate coefficients of the reaction j = 0,1 → j′ = 3 are significantly smaller than the rates cited above for temperatures below 60 K. Above this temperature, the rate coefficient of the reaction j = 1 → j′ = 3 increases and then reaches a similar value than the reaction rate coefficient obtained for j = 1 → j′ = 0. The rate coefficient for the reaction D+ + H2 (v = 0, j) → HD + H+ is reported in Figure 5 using the TIQM ICS. The specific rate coefficients, kj(T), for j = 0,1 have a similar amplitude, and

Figure 2. State-to-state integral cross-sections as a function of the collision energy for D+ + H2 (v = 0, j = 0) → HD (v′ = 0, j′) + H+ for j′ = 0,1,2,3,4.

Figure 3. State-to-state integral cross-sections as a function of the collision energy for D+ + H2 (v = 0, j = 1) → HD (v′ = 0, j′) + H+ for j′ = 0,1,2,3,4.

respectively, j = 0 and j = 1. For the reactions j = 0 → j′ = 0,1,2 and j = 1 → j′ = 0,1,2, on average, ICS decrease smoothly with the collision energy. This is consistent with a barrierless entrance channel. In contrast, the reactions j = 0 → j′ = 3,4 and j = 1 → j′ = 3,4 show an energy threshold corresponding to the energy difference between the rotational levels. The TIQM state-to-state rate coefficients using the previous ICS are displayed in Figure 4 for temperature up to 100 K. The

Figure 5. Theoretical rate coefficients computed from different theoretical methods. TIQM for D+ + H2 (v = 0, j = 0) → HD + H+ (black solid line, this work), TIQM for D+ + H2 (v = 0, j = 1) → HD + H+ (red solid line, this work), TIQM for thermal rate (green solid line), statistical thermal rate from Gerlich25,26 (dashed blue line), TIQM thermal rate computed19 using the ABC code18 (blue square symbol), and experimental thermal rate (white circles with error bars).51 The Langevin value, 2.1 × 10−9 cm3 molecule−1 s−1, is also shown (black dotted line). D

dx.doi.org/10.1021/jp3124549 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

the thermal rate coefficient is close to the rate computed from j = 0. We have also included data coming from other studies. The single point computed by Jambrina et al.19 with the ABC code18 at T = 100 K is in good agreement with our computed value. The statistical method used by Gerlich et al.25 to compute the thermal rate coefficient is in good agreement with our TIQM thermal result for temperatures around 100 K. However, a slight disagreement appears as the temperature decreases. Our quantum result gives a rate coefficient smaller than the one given by the pure statistical approach. We can also see in Figure 5 that the constant Langevin value is much higher than the TIQM rate computed in this work. Another model based on the Langevin approach50 was recently proposed by McCarroll,27 which is adapted to take into account the statistical molecule properties and the rotation of H2. His model gives rate coefficients for j = 0 and j = 1, which are, respectively, larger and smaller than our predictions by a factor between 2 and 4 depending on the temperature. These results yield a thermal rate coefficient that is relatively close to the Langevin value and higher than the TIQM prediction for temperature below 100 K. The good agreement of McCarroll results at higher temperature with the experimental values (at 205 and 295 K) seems to illustrate again that, even if the nuclear spin statistical properties are well incorporated in this modified Langevin model, it does not take into account the fact that collision complex can decay back to the reactant channel, especially at low temperature. It should be also mentioned that the Gerlich measurements (1.6 × 10−9 cm3 molecule−1 s−1) between 180 and 350 K are consistent with older measurements of Feshenfeld et al.,51 yielding a value of 1.0(−0.25; +0.5) × 10−9 cm3 molecule−1 s−1, which was estimated to be independent of the temperature for the 80−278 K temperature range. Figure 5 shows that the TIQM thermal rate coefficient is in good accord with the latter experimental value, and it is consistent with the measurements performed by Gerlich at higher temperature. Nevertheless, the disagreement between experiment and theory observed at very low collision energy (Figure 1) suggests that the TIQM rate coefficient is converged above 50 K, but maybe not fully converged at lower temperatures due to some possible deficiencies in the PES. For practical purposes, fits of the TIQM rate coefficients to the analytical formula, k(T) = α(T/300)β exp(−γ/T), have been performed. These new rate coefficients will be soon added in the chemical database KIDA52 in order to be used by the astrochemical community. Values for the parameters α, β, and γ are shown in Table 1. The fitting expressions using this formula provide a fairly good description over the entire temperature range, with errors below 9% and an average value of 2.5%. These values are consistent with respect to the disagreement

observed between theory and experiment at very low collision energy.

IV. CONCLUSIONS We have performed TIQM calculations of the total and stateto-state integral cross-sections and rate coefficients for the D+ + H2 → HD + H+ reaction at low collision energies and temperatures (