Quantum tunneling during interstellar surface-catalyzed formation of water : the reaction H + H 2 O 2 − → H 2 O +

The final step of the water formation network on interstellar grain surfaces starting from the H+O2 route is the reaction between H and H2O2. This reaction is known to have a high activation energy and therefore at low temperatures it can only proceed via tunneling. To date, however, no rate constants are available at temperatures below 200 K. In this work, we use instanton theory to compute rate constants for the title reaction with and without isotopic substitutions down to temperatures of 50 K. The calculations are based on density functional theory, with additional benchmarks for the activation energy using unrestricted single-reference and multireference coupled-cluster singlepoint energies. Gas-phase bimolecular rate constants are calculated and compared with available experimental data not only for H+H2O2 −−→ H2O+OH, but also for H+H2O2 −−→ H2+HO2. We find a branching ratio where the title reaction is favored by at least two orders of magnitude at 114 K. In the interstellar medium this reaction predominantly occurs on water surfaces, which increases the probability that the two reactants meet. To mimic this one, two, or three spectator H2O molecules are added to the system. Eley-Rideal bimolecular and Langmuir-Hinshelwood unimolecular rate constants are presented here. The kinetic isotope effects for the various cases are compared to experimental data as well as to expressions commonly used in astrochemical models. Both the rectangular barrier and the Eckart approximations lead to errors of about an order of magnitude. Finally, fits of the rate constants are provided as input for astrochemical models.


Introduction
In the dense and cold regions of the Interstellar Medium (ISM), water is known to be formed on the surface of dust grains via sequential hydrogenation of O, O 2 , or O 3 . The full water surface reaction network consists of ∼ 15 reactions and depending on density, temperature and H, H 2 , and O abundance of the interstellar region, different reaction pathways towards the formation of water are important. 1,2 In regions with the highest density, the absolute amount of oxygen becomes sufficiently high for the following reaction pathway to contribute significantly to the formation of H 2 O and OH: The final step in this reaction route is known to have a large activation barrier in the gas phase, see Baulch et al. 3 and references therein. It has also been studied experimentally in the solid phase at low temperatures and a kinetic isotope effect was found. 4 This Only the first reaction produces H 2 O, but the total reaction rate constant, k tot. = k 1 + k 2 , determines the destruction of H 2 O 2 on the surface. Therefore, both reactions need to be considered for an accurate description of the surface process. Furthermore, since tunneling is involved, the kinetic isotope effect should also be studied explicitly.
Here, we present a theoretical study of the reactions R1 and R2 within the concept of studying reactions on a surface that proceed via tunneling. Calculations are performed using a DFT functional and basis set combination that is benchmarked to single-reference and multireference coupled-cluster single-point energies as outlined in Sections 2.1 and 3.1. Rate constants are calculated with instanton theory which is briefly described in Section 2.2. We give activation barriers and rate constants for three different cases: the gas-phase reaction (Section 3.2), the reaction with several spectator H 2 O molecules (small clusters) to mimic a surface related to a bimolecular reaction (Section 3.3), and the reaction with the same small clusters related to a unimolecular reaction (Section 3.4). Previous studies have focused only on the pure gasphase reaction, did not take into account the differences between hydrogenation and deuteration, have not benchmarked their DFT functional, neglected tunneling and/or calculated rate constants only down to 200 K. [9][10][11] Finally we explain how the calculated reaction rate constants can be implemented in astrochemical models (Section 4) and give more general conclusions (Section 5).

Electronic structure
The core system is small (19 electrons) and could be well treated by high-accuracy methods, but the actual focus is a reaction on a surface of water molecules and hence the method of choice to describe the electronic structure is density functional theory (DFT). A suitable functional and basis set needs to describe the interaction between hydrogen and hydrogen peroxide, as well as between water molecules. We perform a benchmark study for the activation and reaction energies of reactions R1 and R2 with respect to two coupled cluster methods. The first benchmark method is spin-orbital based (unrestricted) coupledcluster theory with singles, doubles, and perturbative triples clusters (UCCSD(T)) 12-14 and explicitly-correlated geminal functions (UCCSD(T)-F12) 15,16 employing a restricted Hartree-Fock (RHF) reference function and the cc-pVTZ 17 or cc-pVTZ-F12 18 basis set. To test for multireference character, the internally contracted multireference coupled-cluster method, again with singles, doubles, and perturbative triples clusters (icMRCCSD(T)) was used. 19,20 These computations were carried out with a cc-pVTZ basis set. 17 A complete active space self-consistent field (CASSCF) reference was used for this method. The CAS used for H 2 O 2 , the transition state, H 2 O, OH, and HO 2 are (6e,6o), (7e, 7o), (4e,4o), (3e,3o), and (5e,5o). The active space was chosen to describe the unpaired electron and all bonding and antibonding sigma orbitals of the system. These benchmark computations were carried out for single geometries, as obtained from DFT computations optimized on MPW1B95/MG3S 21,22 and M05-2X/MG3S 22,23 levels for reactions R1 and R2, respectively. This choice of functional is based on finding the best match to coupled cluster energetics. Data are given in Section 3.1.
All geometry optimizations are performed using DL-find 35 within the Chemshell 36,37 framework and NWChem versions 6.3 and 6.6 38 . The single energy points are calculated with Molpro 39 for UCCSD(T) and UCCSD(T)-F12 and GeCCo 19 for icMR-CCSD(T). VMD version 1.9.2 40 and wxMacMolPlt version 7.7 41 are used for visualization.

Reaction rate constants
The hydrogenation reactions are initially modeled in the gas phase. To investigate the influence of the surface more specifically, we added one, two, or three water molecules to the previously optimized structures and re-optimized the resulting configuration. Reaction rate constants are calculated using instanton theory [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60] which has been shown to provide accurate tunneling rates down to very low temperature and is increasingly used to predict rate constants. 56,58, Instanton theory treats the quantum effects of atomic movements by Feynman path integrals. The main tunneling path, the instanton, is described by a closed Feynman path, which connects the reactant and product valleys of the potential energy surface. The instanton represents the tunneling path with the highest statistical weight at a given temperature. It is located by a Newton-Raphson optimization scheme. 56,57 A semiclassical approximation results in the rate constants. More details on our implementation of instanton theory are given elsewhere. 56,57 Rotational and translational partition functions were approximated by their classical analogues (J-shifting approximation), which is generally accepted as a good approximation at the temperature scale considered here. In bimolecular cases, the product of the partition functions of the separated reactants was used, in unimolecular cases, the partition function of the encounter complex. The Feynman paths were discretized to 60 images. Instanton theory is applicable below the crossover temperature T c , which is defined as where ω b is the absolute value of the imaginary frequency at the transition state. The instanton represents the tunneling path with the highest statistical weight at a given temperature. Instantons were optimized to a residual gradient (derivative of the effective energy of the instanton with respect to the mass-weighted atomic coordinates) below 10 −8 atomic units (hartree bohr −1 m −1/2 e ). This and other parameters were chosen equivalently to previous Table 1 DFT functional/basis set combination benchmark with respect to UCCSD(T)-F12/cc-pVTZ-F12, UCCSD(T)/cc-pVTZ and frozen-core icMRCCSD(T)/cc-pVTZ single-point energies for reactions R1 and R2, respectively. Reaction energies computed from the HEAT protocol are given, too. Values are given in kJ/mol excluding zero-point energies. First, we consider the pure gas-phase reaction. This may not be very relevant in terms of astrochemistry, but it is important to understand the simplest case first. In the gas phase, a hydrogen atom can approach the molecule on both oxygen atoms, since they are equivalent. In other words, the rate constants need to be multiplied with a rotational symmetry factor. The symmetry factor used here is 2, resulting from the pointgroup C 2 for the hydrogen peroxide molecule. 84 Concerning the surface reaction, there are typically two reaction mechanisms taken into consideration: the Eley-Rideal (ER) and the Langmuir-Hinshelwood (LH) processes. The ER mechanism describes one species to be adsorbed on the surface and the other approaching from the gas phase, i.e., an overall bimolecular reaction. For the LH mechanism both species are adsorbed on the surface, approach each other via diffusion and form an encounter complex of H and H 2 O 2 on the surface. This encounter, or pre-reactive, complex can then decay to yield the reaction products in a unimolecular process. Moreover, to extend the results to the solid phase it is key to realize that rotational motion on the surface is restricted. Therefore, rate constants calculated for both ER and LH mechanisms need to keep the rotational partition function constant between the reactant and transition state. Moreover, the surface structure breaks the gas-phase symmetry, hence no symmetry factor is required.
For astrochemical modelers to be able to easily implement the calculated rate constant, we fitted these to the rate expression 85 The parameters α, β , γ, and T 0 are all fitting parameters, where α has the units of the rate constant, β regulates the lowtemperature behavior, and γ and T 0 can be related to the activation energy of the reaction. Instanton rate calculations were used for the fits at low temperature, below T c , while rate constants obtained from transition state theory including quantized vibrations and a symmetric Eckart model for the barrier were used above T c .

Results and Discussion
Here we first present the results of the DFT benchmark (Section 3.1) and subsequently give the results for the three cases of the reaction that we study: the gas-phase reaction (Section 3.2), the bimolecular Eley-Rideal reaction (Section 3.3), and the unimolecular Langmuir-Hinshelwood reaction (Section 3.4). All values for the rate constants are given in the Electronic Supplementary Material. Table 1 gives an overview of the activation and reaction energies without zero-point energy (ZPE) corrections for reactions R1 and R2. All DFT calculations comprise full geometry optimizations and all stationary points were verified by their appropriate number of imaginary frequencies, i.e., zero for the reactants/products and one for the transition states. The coupled cluster values are single-point energies on MPW1B95/MG3S and M05-2X/MG3S levels for reactions R1 and R2. Finally, reaction energies have also been calculated with the use of atomization energies computed with the HEAT-456QP protocol (as tabulated in Harding et al. 24 ). CCSD(T) results are commonly used in computational chemistry as a gold standard for activation and reaction energies. This is, however, valid only for species where a single reference wavefuction is a good approximation. We found values of the T1 and D1 diagnostics of 0.022 and 0.062 (R1) and 0.031 and 0.109 (R2) in our CCSD(T)-F12 calculations. Thefore, additional tests with the icMRCCSD(T) method seemed important. The icMR-CCSD(T) method translates the accuracy of CCSD(T) to multireference cases and has been successfully applied to predicting barriers of reactions. 86 An estimate of the multireference effects can be obtained from the results of icMRCCSD(T) and UCCSD(T) calculations using the cc-pVTZ basis. The computations indicate, that multireference effects only slightly lower the activation energy (−2.8 kJ/mol for R1 and −1.3 kJ/mol for R2). We hence conclude that the multireference character of the transition state is not pronounced. This can also be seen from the contributions Comparing the coupled cluster values to those obtained with DFT, it is clear that the barrier of reaction R1 is best described by the combination MPW1B95/MG3S and reaction R2 by M05-2X/MG3S, which is in full agreement with previous findings. 10 The slight overestimation of the barrier can result in an underestimation of the calculated reaction rate constants, but we will show that the spread in the activation energies resulting from the interaction with spectator H 2 O molecules is much larger. The need for using two different functionals for the O-O bond breaking and H-abstraction reactions is somewhat disatisfactory from a purist's point of view. At present, however, our choice is dictated by the absence of any practial functional that is good as describing both reaction paths. 10 The issue is alleviated to some degree by the fact that the two processes take place on different parts of the potential energy surface.

DFT benchmark
As we will outline in Section 3.2 the main focus of the paper is on reaction R1 and therefore we focus our attention concerning the description of dispersive interactions on this reaction only. The functional MPW1B95 has been shown to have a good performance for hydrogen bonding and weak interaction calculations. 21 Therefore it is not obvious if an empirical dispersion correction should be applied additionally. We have tested the per-  and such a surface is important for our astrochemical application, we will not use an additional dispersion correction. Figure 1 shows the bimolecular reaction rate constants calculated for reactions R1 and R2 for the pure gas-phase situation. Note that the black curve has been calculated with the MPW1B95/MG3S combination and the red curve with M05-2X/MG3S following earlier statements. Instanton theory is not applicable above the crossover temperature, which are 275 and 463 K for R1 and R2 respectively. This is why the two curves are cut off at the higher temperature range. From the figure the rate constants calculated at the highest temperatures can be compared to the expressions summarized by Baulch et al. 3 based on experimental and modeling work. [87][88][89][90] These expressions are recommended down to 280 or 300 K and a reasonably good correspondence to our calculated values can be seen.

Gas-phase reaction
Already at a temperature of 114 K the difference between the rate constants of R1 and R2 is more than two orders of magnitude. This can be explained by the much higher activation energy, despite the fact that tunneling might be expected to dominate H-abstraction at low temperatures more than the breaking of the O-O bond. Therefore, we exclude it from further investigation, i.e., the branching ratio R1:R2 at lower temperature is at least 100:1. Figure 3 shows the instanton paths at 114 K for both reactions. The path essentially shows the delocalization of the  atoms involved in the reaction. This deviates from the classical picture of overcoming a barrier and visualizes tunneling through a barrier. The lower the temperature, the more delocalization is usually visible, since tunneling then plays a larger role.
Several possible isotopic substitutions can be made: O vs. 18 O and H vs. D. In Fig. 2 the resulting rate constants for such substitutions are depicted. Firstly, changing the oxygen atom to a heavier isotope results in a lowering of the rate constant by a factor 1.4 at 50 K, because the O-O bond needs to be broken. In Fig. 3 it is visible why: the oxygen atoms are somewhat delocalized as well, meaning that they take part in the tunneling process. Exchanging a protium for a deuterium on the peroxide molecule has a similar small effect, i.e., decreasing the rate constant by a factor 1.3 at 50 K. However, for the hydrogen that approaches the peroxide and is to be added to the oxygen there is a strong kinetic isotope effect of a factor 229 at 50 K, see Table 3, again visualized by strong delocalization.

Eley-Rideal surface reaction mechanism
The Eley-Rideal reaction mechanism in surface chemistry corresponds to a reaction between one reactant that is adsorbed on a surface, and a second reactant that approaches directly from the gas-phase, hence to a bimolecular reaction between the Table 2 Activation energies, E a for the hydrogenation and deuteration versions of reaction R1 with respect to two separated reactants (bimolecular ER) and to an encounter complex (unimolecular LH) in kJ/mol including zero-point energy. The cross-over temperature indicating when tunneling dominates the reaction rate in is given in K.  114  36  26  50  24  30  64  97  58  53  87  39  56  120  84  91  87  146  63  84  215  74  117  128  204  85  130  326  66  169  148  266  127  176  461  59  169  176  269  133  -510  54  199  202  306  163  340  -50  229  202  522  197  -1194 adsorbate-surface system and the incoming atom. In our case, the adsorbate is the H 2 O 2 molecule, the incoming atom is the H atom, and the surface is represented by spectator H 2 O molecules. The opposite case where H is pre-adsorbed and peroxide comes in from the gas phase is much less likely given the low H 2 O 2 gasphase abundance in the ISM. Adding one (A or B), two (A and B), or three (A, B, and C) spectator water molecules, as depicted in Figure 4, helps to elucidate the effect of the interaction between the water surface and and the reactants. In Fig. 4, a transition structure with three added water molecules is shown, including their hydrogen bonded structure and respective labeling. Table 2 summarizes the activation energies for the reaction between hydrogen peroxide and protium or deuterium. Adding water molecules to the reactive site indeed has an influence on the activation energies, which consequently span a range between 25.1 and 32.5 kJ/mol. The activation energy seems to correlate with the O-O and O-H (incoming) distances in the transition state: the higher the activation energy, the larger the O-O and the smaller the O-H bond lengths, i.e., the later the transition state is.
The corresponding reaction rate constants are presented in Fig. 5. Here, rate constants are calculated keeping the rotational partition function constant and without a symmetry factor. Therefore, the black curves without water have been recalculated and  differ from those in Figs. 1 and 2. Comparing the five curves to their respective activation energies we can see that the higher the barrier, the lower the rate constant, as can be expected. The spread at the lowest temperature, 50 K, is about 1.5 orders of magnitude. However, all curves do seem to follow the same trend, which indicates that the way the reaction proceeds it not altered by the presence of water molecules. This means that although hydrogen bonds can and are formed, the influence thereof lies only in the height of the activation energy. The rate constants seem to level off around a value of 10 −19 cm 3 s −1 at 50 K. Again, the same correspondence between the activation energy and rate constant is found for the deuterated reactions. Furthermore, a similar spread in the rate constants is observed at 50 K. The ratios between the reaction rate constants for hydrogenation and deuteration of H 2 O 2 for the ER and LH rate constants are summarized in Table 3. This ratio is defined here as the kinetic isotope effect, φ KIE ,

Langmuir-Hinshelwood surface reaction mechanism
The Langmuir-Hinshelwood reaction mechanism assumes two already adsorbed species on the surface that find each other via diffusion, sit next to each other, i.e., form a pre-reactive complex, and then attempt a reaction. This is thus related to a unimolecular reaction starting out from the encounter complex; note the different units. Table 2 gives the activation energies for the reactions with respect to this encounter complex. It is clear from the small differences in the barriers between the first and the second column of the table that the interaction between the hydrogen atom and the H 2 O 2 (+(H 2 O) n ) system is very weak. In fact the interaction energies are slightly positive, 1-2 kJ/mol. This could be caused by the poor description of the Van der Waals interaction by DFT,  however, even an additional D3 correction leaves the interaction energies positive. The result is that at lower temperatures the rate constants become more noisy. This weak interaction of the H atoms shows that it is not trivial to define an encounter complex, since most likely one should take an ensemble of possible configurations into account. This will be dealt with in a future study.
The reaction rate constants are presented in Fig. 6. Again the rotational partition function has been kept constant. Furthermore, a similar trend as for the ER bimolecular case can be observed, with a spread of about an order of magnitude. The rate constants seem to level off around a value of 2 × 10 4 s −1 at 50 K.
Experimentally, under ultra-high vacuum conditions, it is assumed that the LH mechanism plays the dominant role. Early work by Miyauchi et al. 92 , Ioppolo et al. 93 , Cuppen et al. 94 already specifically mentions the role of the reaction H + H 2 O 2 − − → H 2 O + OH within the water network, but it was only by Oba et al. 4 that the reaction and the role of tunneling were explored to a greater detail. In their paper, they show a clear kinetic isotope effect between hydrogenation and deuteration of H 2 O 2 and also of D 2 O 2 . Since experimentally it is not possible to determine the amount of H/D atoms residing on the surface, the reported φ KIE 's are based on effective rate constants and have a value of approximately 50. A more exact treatment would have to include the exact H:D flux ratio as well as diffusion and recombination rates. Our low-temperature φ KIE is at least a factor of four higher, which is caused by the fact that the theoretical φ KIE is based solely on the rate constant for the reactions, excluding any additional effects. Since the longer lifetime of D atoms on the surface with respect to that of H atoms on the surface indeed is expected to lower the kinetic isotope effect, the agreement to experiments is actually reasonable. Note that in the unimolecular case, the addition of spectator molecules strongly influences the kinetic isotope effect, the reason for this is at present unclear. Eley-Rideal mechanism: Comparison between rate constants calculated with transition state theory including quantized vibrations, the Eckart approximation, and instanton theory. The fit to Eqn 2 is also shown.

Implementation in astrochemical models
The rate constants that are presented in this work can be used as input for astrochemical models. To that purpose, we have fitted rate constants to equation 2 with fitting parameters α, β , γ, and T 0 . The values of these parameters are given in Table 4. Note that in all cases, we used the pure gas-phase geometries, i.e. without spectator molecules, but calculated the ER and LH rate constants keeping the rotation partition function constant. In Fig. 5 we showed that the rate constants without spectator molecules lie nicely within the range spanned by the various geometries and can thus be interpreted as an average value. The full spread is about 1-1.5 orders of magnitude and therefore the exact values of the rate constant fits should be considered to have an uncertainty of approximately a factor 5. We stress here that in astrochemical grain-surface reaction modeling this accuracy is in fact quite sufficient. Uncertainties of that order can even be found in gas-phase reaction networks although these are typically better constained. 95 As outlined above, this reaction can, in principal, take place both in the gas phase and on a (water) surface. As a result of the low interstellar abundances of H 2 O 2 in the gas phase, the gas-phase route is unlikely. A gas-phase hydrogen atom could directly strike a reaction partner that is adsorbed to a water or grain surface, as in a bimolecular process. However, the timescales involved in dark cloud chemistry are large enough for a hydrogen atom to scan the surface and meet a H 2 O 2 molecule. Therefore, unimolecular rate constants have a direct correspondence to the process as it occurs in the interstellar medium. For the sake of completeness we provide fit parameters for all three possibilities.
The rate constants presented here are not directly comparable to the expressions commonly used in astrochemical modelling. Usually in such models the tunneling probability, P tunn. , is de-  scribed by the rectangular barrier approximation where a is interpreted as the barrier width, µ as the effective mass and E a the activation energy of the reaction including ZPE. This is a very convenient expression, because is it implementable in rate equation models. Most of these take diffusion into account in the calculation of the LH reaction rate: Here, the rate of diffusion, R diffusion , includes the diffusion of both species, with k diff the unimolecular diffusion rate constant, N sites the number of surface sites and n X the concentration of species X. The probability for reaction is composed of the competition between reaction, diffusion out of the site, and desorption rate constants, In models, at low temperature, k react. is often approximated by a trial frequency, ν, multiplied by P tunn. , whereas what we calculated corresponds directly to k react. = k tunn. . Therefore, the best way to compare the instanton rate constants to the rectangular barrier approximation is to look at the kinetic isotope effect, see Table 5. Furthermore, Taquet et al. 11 improved on the rectangular barrier approximation by the use of an Eckart correction to describe tunneling. In the same Table 5

Conclusion
With the study of the reaction between hydrogen and hydrogen peroxide the final step of the water formation sequence starting from molecular oxygen is further quantified. In particular, attention is paid to the low-temperature behavior, the kinetic isotope effect, and the influence of spectator water molecules mimicking an icy grain surface. Specifically, • prior to calculating rate constants the method of choice (DFT) has been benchmarked to single-reference and multireference coupled cluster single-point energies for the stationary points on the potential energy surface, • a branching ratio between the rate constants for O-O bond breaking, R1, and H-abstraction, R2, of at least 100:1 was established at a temperature of 114 K, • rate constants that apply to the gas-phase, surface Eley-Rideal, and surface Langmuir-Hinshelwood mechanisms are now available down to 50 K, • the 50 K results, or an extrapolation down to at least 30 K via the fitted expression Eqn 2 (Table 4), can be used as a reasonable guess for even lower-temperature surface processes thanks to the rate constants leveling off with decreasing temperature, • quantitative agreement with experimental gas-phase data and qualitative agreement with the experimental surface kinetic isotope effects is found, • the addition of spectator molecules indeed influences on the reaction rate constant and kinetic isotope effect, mainly by influencing the transition state structure, which leads to a change in the activation energy of the reaction, • general trends such as the asymptotic behavior of the curves, don't seem to be strongly affected by the addition of spectator molecules, but it is important to note that the surface aids in bringing the two reactants together and allows for heat dissipation of the exothermicity, • a comparison between the rectangular barrier approximation, the Eckart approximation, and instanton rate constants shows that both approximations leads to large errors (more than an order of magnitude).