H₂ is the most common molecule in the interstellar medium. Its formation is assumed to occur on grains as the gas phase formation is too slow to explain its abundance. In diffuse clouds the grains are bare and they are usually classified as silicates or carbonaceous (Wakelam et al., 2017). In dense interstellar clouds on the other hand the grains are often covered by a mantle of ice, which may contain for instance H₂O, H₂CO, N₂, O₂, CO, CO₂, H₂O₂, NH₃, CH₄ and CH₃OH (Tielens and Hagen 1982; Gibb et al., 2000). The ice is expected to primarily be in the form of amorphous solid water (ASW) or polycrystalline ice (PCI).

The process of forming H₂ on an icy grain surface is believed to involve hydrogen atom physisorption followed by diffusion such that another hydrogen atom is encountered whereby the two atoms can react and form a hydrogen molecule. Large amounts of energy is released as H₂ forms, which together with the weak interaction of the hydrogen molecule with the surface makes it plausible that it quickly desorbs and ends up in the gas phase. The process just described is essentially an example of how heterogeneous catalysis works. The hydrogen atoms on the grain surface do however not have to come from the gas phase. It is also possible that hydrogen atoms form in the bulk of the ice due to photodissociation. It should also be mentioned that if there is a hydrogen atom on the surface it could be hit directly by a hydrogen atom approaching from the gas phase, which is the Eley-Rideal mechanism. This mechanism is not so likely for a grain with low coverage of H atoms.

Much experimental effort has gone into investigating the diffusion of hydrogen and deuterium atoms on ASW and PCI at temperatures at and around 10 K (Manicò et al., 2001; Hornekær et al., 2003; Perets et al., 2005; Amiaud et al., 2007; Matar et al., 2008; Watanabe et al., 2010; Hama et al., 2012; Kuwahata et al., 2015). The diffusion is found to be much faster on PCI than on ASW. Recent work of Kuwahata et al. (2015) studies the diffusion of D and H atoms on both PCI and ASW. The authors find that the kinetic isotope effect (KIE) for diffusion of H atoms compared to D atoms is larger on PCI than on ASW. They also find that the KIE is larger for their larger fluxes of atoms towards the surfaces than for lower fluxes. Kuwahata et al. (2015) interpret their results as being due to a larger contribution of tunnelling to the diffusion rate on PCI than on ASW. For ASW they find that thermal hopping dominates although tunnelling diffusion may also partly contribute.

There has also been a substantial amount of theoretically based work on absorption diffusion on grains in the context of astrochemistry (Gould and Salpeter, 1963; Williams, 1968; Augason, 1970; Hollenbach and Salpeter, 1970; Hollenbach and Salpeter, 1971; Lee, 1972; Watson and Salpeter, 1972; Jura, 1975; Goodman, 1978; Smoluchowski, 1979, 1981, 1983; Pirronello et al., 1997a, 1997b, 1999; Katz et al., 1999; Karssemeijer et al., 2014; Senevirathne et al., 2017; Fredon and Cuppen 2018). Important for H₂ formation on a surface is that the residence time of the adsorbed species is long enough that the adsorbed hydrogen atoms have time to diffuse on the surface before desorbing. In early model calculations, using two types of adsorption sites (i.e. weak and strong adsorption sites), by Hollenbach and Salpeter (1970,1971) the tunneling was efficient enough to allow the adsorbed hydrogen atoms to diffuse over most of the grain surface before desorbing. Pirronello et al. (1997a, 1997b, 1999) measured hydrogen formation on olivine and carbonaceous substrates and found much lower hydrogen atom mobility and H₂ formation rate than Hollenbach and Salpeter and pointed to the temperature dependence of this rate. The experimental results were also modelled by Katz et al. (1999). Smoluchowski (1981, 1983) noticed based on quantum mechanical calculations that adsorbed hydrogen atoms would on an amorphous quickly get trapped in the deep sites, which lowered the H₂ production efficiency by orders of magnitude compared to Hollenbach and Salpeter’s.

In the present work we theoretically study the kinetic isotope effect between deuterium and protium (i.e. the common hydrogen isotope) atom diffusion by estimating the transition rates of H and D atoms hopping between adjacent local minima on the grain surface. This is done by employing simple models that are intended to represent the hexagonal ice and the amorphous solid water. In this way we obtain kinetic isotope effects for the transition rates and thereby also for the diffusion rates. We aim to offer some insight as to why Kuwahata et al. (2015) found larger values on hexagonal ice (I _h ) than on ASW.

In calculating the KIE’s, we will begin by using harmonic transition state theory. First we apply it in completely classical form and thereafter with quantized partition functions. We will also consider a simple Wigner tunnelling correction to the transition state theory (TST) rate constant. As far as ASW is concerned, TST turns out to suffice to reproduce the experimentally derived KIE. For I _h , however, the TST derived KIE is not compatible with the experimental result. For I _h we will therefore pay more attention to tunnelling as it may be enhanced due to its more regular structure than in the ASW case. For this purpose we will treat the physisorbed hydrogen atom as moving in a double well model. First we will treat the double well as isolated and thereafter as weakly interacting with the ice lattice. Finally we will employ a band model by Kehr (1997), where we assume that the I _h is perfectly periodic over large distances.

The paper is structured such that this Introduction is followed by a brief description in Section 2 of the experiments of Kuwahata et al. (2015). Then in Section 3 the ice model that is used for the present calculations is introduced. Thereafter calculations and results are presented in Section 4 (ASW) and Section 5 (I _h ). Section 6 contains a discussion and conclusions.

Our calculations are inspired by the experimental work of Kuwahata et al. (2015) on H and D atom diffusion on PCI and ASW. Their work was performed at (or close to) 10 K. In the experimental work H (or D) atoms are deposited on a PCI (or ASW) surface, at varying fluxes. During deposition, atoms are photodesorbed by weak laser radiation and then detected by time-delayed resonance enhanced multi-photon ionization (REMPI). The time-delay between the laser pulses was varied and the H/D REMPI intensity ratios were summed over the delay-time. In this way KIE’s for the diffusion rate could be determined. This was done by assuming that the REMPI signal due to H (D) atoms depend on how many of the H (D) atoms were desorbed from the surface by the weak desorption laser, which in turn would depend on how many H (D) atoms had already disappeared by forming H₂ (D₂) due to diffusion on the surface.

The experiments were repeated for a range of deposition fluxes. The idea behind this is that at low flux, the surface coverage of atoms will be low and the atoms will be far apart. As they diffuse they are likely to encounter irregularities or defects in the ice surfaces. Therefore Kuwahata et al. (2015) suggest that thermal diffusion will be rate limiting. For high flux on the other hand the atoms will encounter each other after shorter travel. For the I _h this could mean that a more periodic surface is encountered, supposedly allowing enhanced tunnelling. In the Discussion we will return to this.

By reading from Figure 2 in Kuwahata et al. (2015) the following results are found. Their measurements on PCI give at high flux a KIE of 120 and accounting for their error bars the KIE is in the range 65–175. Their corresponding results for ASW give a KIE of 16 and accounting for error bars the KIE is in the range 10–25. The KIE is taken to be the ratio of the rate for the lighter isotope to that for the heavier one. At lower fluxes the KIE’s were smaller than at large fluxes, both for ASW and for PCI. Kuwahata et al. (2015) conclude that for PCI the much faster surface diffusion of the H atoms than the D atoms at high flux cannot be explained by classical thermal hopping. They also note that their work is the first experimental evidence of quantum-tunnelling diffusion on a PCI surface. For the ASW surface they conclude that tunnelling diffusion may partly contribute, but the diffusion is predominantly thermal hopping given the smaller KIE than that observed on PCI.

We end this section by noticing that the experimentally determined KIE represents the KIE of the diffusion constants for H and D, which may be related to hopping rates between potential minima as follows. The diffusion constant D on a surface can be obtained from the relationship

where < r ² > is the average distance that an atom/molecule diffuses on the surface in time t _diff , which is the sum of the individual hopping times, t _i , between potential minima on the surface (that occur in the total time t _diff ). Since t _i = 1/k _i where k _i is the rate constant for a jump between two minima we can write D = < r 2 > 4 Σ i 1 k i . (2)

If all k _i are increased by the same factor, we see that in that case D increases by the same factor. This motivates that we below assume that a certain KIE in the hopping rates between minima will be reflected in the same KIE in the diffusion constants. ¹

For the purpose of the present paper we will need thermal rate constants for H and D atoms hopping between adjacent wells on ice surfaces. Here we will treat the underlying ice as frozen and thus the hydrogen atom is the only moving species, thereby making the problem three-dimensional. This is the same model as was used by Senevirathne et al. (2017). In deriving parameters for the models to be described below we use the interaction potential of Andersson et al. (2006), noticing the minor corrections described in Senevirathne et al. (2017). We note that since the ice is treated as frozen, self-trapping will not occur.

Our main interest is in studying the difference in hopping rates between H and D on ASW and I _h . ASW can be expected to have a larger spread in barrier heights than I _h (Senevirathne et al., 2017). On the other hand, typical pathways for diffusion likely avoid the highest barriers. This would make the barriers that are actually passed differ less in height than would otherwise be the case (Senevirathne et al., 2017).

Here we will pick out a single well and transition state from an I _h surface and base all of our calculations on that. The I _h surface was obtained by Andersson et al. (2006) by taking snapshots from molecular dynamics simulations. We choose a well that can reasonably well characterize H and D atom hopping for both our ASW and I _h calculations. Thus, possible differences that we find in hopping rates between ASW and I _h jumps should thus be related to how we treat effects of symmetry and periodicity of the I _h lattice, which are absent in ASW. Some relevant data regarding the potential well and transition state can be found in Table1.

In this section we will investigate what transition state theory (TST) tells us about the H to D diffusion rate ratios on ASW and I _h . We assume that the diffusion rates are proportional to the rate constants for hopping between adjacent minima on the ice surfaces, see motivation at the end of Section 2. Our interest will thus be in finding the thermal rate constants for such hops by applying TST. In general we assume the underlying ice surfaces to be frozen such that only the physisorbed H (D) atoms move.

In TST the rate constant may be written:

where k _B is the Boltzmann constant, T is the temperature in Kelvin, h is the Planck constant, Q ^‡ is the partition function at the transition state (with the imaginary frequency vibration excluded), Q ^R is the reactant partition function and E _b is the barrier height on the potential energy surface.

If we evaluate the partition functions classically and assume uncoupled harmonic oscillators we get a contribution k _B T/ℏω _i from each vibrational mode i. Considering other approximations involved we may as well set the deuterium mass to be twice that of protium. In a fully classical TST calculation of the KIE, most factors will then cancel out, resulting in a KIE of 2 ≈ 1.4 . This value is independent of temperature, as long as the temperature is high enough to motivate a classical treatment. This is however often not the case, and certainly not in the present case. Our result for the KIE is also independent of the number of harmonic oscillators that are included. Thus, for this case we do not need to assume the atoms in the ice to be frozen (ignoring self-trapping).

Let us now take TST a step further by evaluating the partition functions quantum mechanically, again assuming uncoupled harmonic oscillators and let us set the reference level for each individual oscillator to be its zero point energy (ZPE) level rather than the bottom of the potential well. Accounting for this the TST rate constant may be expressed as

where the partition functions are evaluated with the ZPE level as reference and ΔE ^VAG is the vibrationally adiabatic ground state barrier height Δ E V A G = E b + ∑ i = 2 3 ℏ ω i ‡ / 2 − ∑ i = 1 3 ℏ ω i / 2 , (5)

where ℏ is the reduced Planck constant.

The smallest real vibrational angular frequency involved in the evaluation of k(T) is 12.3 ps⁻¹ (Table1), which is one of the deuterium frequencies at the transition state. This individual oscillator gives a contribution to the partition function of ≈1.0001 (with the ZPE level as reference) when evaluated at 10 K, which will be the most interesting temperature as there we have experiments to compare with. All other individual oscillator contributions to the overall partition functions will be smaller (closer to unity), and we can thus approximate them all as essentially being unity. From Table 1 we find that ΔE ^VAG = 3.04 meV for H and 5.56 meV for D. This gives a difference in VAG barrier height of 2.52 meV between the isotopes, which yields (using Eq. (4)) a KIE of 18.6 at 10 K.

For illustrative purposes we now perform the TST calculation for a one-dimensional case. Then Q Z P E ‡ = 1 and in addition Q Z P E R = 1 at low temperature. We also obtain

and

where ω _1,H = 21.8 ps⁻¹ is the relevant angular frequency for the hydrogen atom and other notation should be intuitive. Again taking the deuterium mass to be twice that of protium we find Δ E D V A G − Δ E H V A G = ℏ ω 1 , H 2 1 − 1 / 2 . (8)

Using this to evalute the KIE [using Eq. (4)] at 10 K we find that it becomes 11.4. This result is independent of the choice of potential beyond the values used for ω _1,H and ω _1,D .

The difference between the KIE of 11.4 and the classical value of 1.4 may be viewed as resulting from the quantum mechanical zero point energy. We would typically see this as the effective barrier height to reaction becoming lower when the ZPE is accounted for and this effect is obviously larger for H than for D and thereof the increase in KIE compared to the classical case. The further change to find a KIE of 18.6 when treating three vibrational motions for the reactant, rather than a single one, is related to the additional changes in zero point energies.

The simple evaluation that we have just performed accounts harmonically for quantum effects in all degrees of freedom. For the reaction coordinate we are however lacking dynamic tunnelling through the barrier, which results from the wave function never decreasing to zero within the barrier. Since the cross over temperature (Hänggi et al., 1990) T _c = 7.2 and 5.1 K for H and D respectively, tunnelling is however not expected to be dominant at 10 K. A simple Wigner tunnelling correction (Wigner, 1932) would however give an increase in rate constant with a factor of 1.86 for H and 1.43 for D at 10 K (using imaginary angular frequencies from Table 1). This would increase the KIE to 24. We should also point out that reflection from the barrier (or recrossing classically) is not accounted for in TST, which would act to reduce the rate constants for H and D and more for H than D, so as to oppose the effect of tunnelling on the KIE.

The simple calculations just described give KIE’s that agree with the experimental results for ASW where the KIE is in the range 10–25 for high flux. While this may be a coincidence, it still suggests that there is some realism to the calculations, particularly as we have noticed that the calculated KIEs do not explicitly depend on the barrier height. The KIE results obtained by TST are summarized in Table 2. We consider these results final for ASW but not for I _h where we need to consider the effect of tunnelling in more detail. In the next section we focus on the tunnelling aspects for hexagonal ice.

The discussion above regarding TST estimation of the KIE for ASW applies just as well to hexagonal ice and calculations would give exactly the same result for the same potential. Kuwahata et al. (2015) however observed a substantial increase in KIE for I _h compared to ASW. This difference must result from the differences between potential energy landscapes. Here we will investigate some effects of symmetry and periodicity on the potential energy surface for the hexagonal ice, but otherwise let the single well be the same as used for ASW.

Three simple models from the literature will be briefly described in the context of H and D atom transitions between adjacent minima on hexagonal ice. The reader is reminded that in all cases the underlying ice is treated as frozen and thus the hydrogen atom is the only moving species, thereby making the problem three-dimensional.

As mentioned the present work has been inspired by the experimental work of Kuwahata et al. (2015) at 10 K, where they for instance state that tunnelling diffusion may occur for a short distance within each single crystal, whereas it should be highly suppressed for long-distance diffusion. This statement is fully supported by our results obtained using the Kehr model (Kehr 1997), described in Sec. 5.3. We find that the model is not applicable at 10 K, which supports that coherent tunnelling over several wells is not likely to happen.

We next discuss tunnelling over short distances on crystalline ice. The isolated symmetric double well model described in Sec. 5.1 represents an idealized case where we can numerically exactly calculate the time for a transition from one well to the other. For protium this time is 3.3 × 10^–11 s (for deuterium the corresponding time is 6.9, ×, 10^–10 s), assuming that we begin in a superposition of the two lowest eigenstates. If we start from a superposition of the next pair of vibrational eigenstates this time decreases markedly to become 2.8 × 10^–12 s (1.2 × 10^–10 s for D).

For a diffusion process to occur, the hydrogen atom must spend much longer time in a well waiting to jump to another, than the time that the actual jump takes. Otherwise the jumps may be correlated and not correspond to a diffusion process. The transition times mentioned in the previous paragraph correspond to the time the transition takes for a symmetric double well potential. In reality self-trapping and coupling to the thermal motion of the lattice would cause the transition time to be longer as we would have to wait for the lattice to thermally end up in a geometry that would correspond to a symmetric double well potential, thereby allowing fast tunnelling. In this way a diffusion process can occur.

In our supposedly most realistic model of the ones employed, a symmetric double well (DW) interacts weakly with the lattice. In this model we have approximately accounted for the interaction with the lattice vibrations, but not for self-trapping (or electron drag, should that matter). For this case, from the transition rates, we obtain lifetimes for how long time it takes for half of the population to move from one well to the other, i.e. the half-life. For protium the half-lives are 6.2 × 10^–10 s, 4.8 × 10^–10 s and 2.3 × 10^–11 s at 0 K, 10 and 25 K respectively. For TST the half-lives are 1.1 × 10^–10 and 9.7 × 10^–12 at 10 and 25 K respectively, while at 0 K the half-life goes to infinity. These half-lives, together with the tunnelling times for the isolated DW, suggest that at these temperatures the H atoms predominantly move by diffusion on the surface. As the temperature is further increased, however, a more fluid-like motion rather than a diffusive motion becomes more and more important (as does desorption).

The DW model accounts for tunnelling through the barrier while TST, without tunnelling correction, accounts for transitions above the barrier. Summing these two contributions up should thus roughly give us the total rate constant for transition from one well to the next. In Figure 1 this is illustrated for protium. We notice that k (T) approaches a constant value as T decreases. This can be seen to be due to the rate becoming limited by tunnelling from the lowest pair of eigenstates in the double well. This happens just above 7 K, which agrees with the cross over temperature being 7.2 K. ² As the temperature increases above 7 K, TST makes the dominant contribution to the thermal rate constant, indicating that over-barrier jumps dominate.

In Figure 2, k (T) for deuterium is shown. As can be expected, in this case TST contributes even more than what tunnelling does in comparison to the H atom case. From the total rate constants for H and D transitions, we can obtain values for the KIE for various temperatures as illustrated in Figure 3. We notice that the KIE obtained is extremely sensitive to the temperature. At 10 K the calculated KIE is much lower (about a factor of six) than the experimentally estimated value. At 7 K, however, the calculated value has increased to coincide with the experimental value.

In Figure 3 we also show the KIE that is obtained from the TST analyses only, i.e. over the barrier transitions. Since the vibrational partition functions are very close to unity (ZPE level as reference), this KIE is determined by the difference in VAG barrier height. Even if the classical potential barrier height is in error, this does not affect the KIE calculated by TST, beyond changes in vibrational frequencies. The TST contribution is thus likely less sensitive to the potential than what the tunneling contribution is. For the ASW, we expect much smaller tunnelling contribution than for the crystalline ice. We note that for the ASW, the TST KIE agrees with the experimentally derived KIE.

One explanation for the fact that we cannot reproduce the experimentally derived KIE on I _h could be that the height and/or shape of the employed potential is inaccurate. The calculated KIE for ASW hinges on the difference in TST rate constants for H and D, which in turn hinges on the difference in ΔE ^VAG , which hinges on the vibrational frequencies and only indirectly on the barrier height itself. The tunnelling contribution to the KIE, which is more important for I _h , depends strongly on the barrier height and width. In addition, in the present double well calculations the KIE depends on at which temperature contributions from the excited pair of eigenstates kick in noticeably. If the energy level of the second pair of eigenstates is raised, its contribution, which lowers the KIE, decreases at a given temperature and thus the KIE will increase at a given temperature. In order to reproduce the experimentally deduced KIEs, the tunnelling contribution would need to be larger. This suggests that an increase in barrier height and/or width, which would increase the KIE for I _h much more than for ASW, would lead to better agreement with the experimental results.

The work on model potentials presented here has been aimed at giving basic insights into, and qualitative understanding of, the kinetic isotope effects observed by Kuwahata et al. (2015). Clearly, for quantitative results, other calculations using accurately determined potentials rather than the model potentials employed here would be useful. It should also be kept in mind that even the fairly smooth I _h surface does in reality show disorder even at short range, so that a single well with a single barrier height and no self-trapping will not capture all aspects of the adsorbate-ice interaction.

In conclusion we have found that transition state theory explains the experimentally found kinetic isotope effect on ASW. This supports the conclusion of Kuwahata et al. (2015) that for the ASW surface diffusion is predominantly thermal hopping rather than tunnelling. On the other hand, transition state theory does not explain the larger kinetic isotope effect on PCI. Our calculations indicate that a substantial tunneling contribution is required to explain the KIE on Ih. This agrees with the conclusion of Kuwahata et al. (2015). Further, Kuwahata et al. (2015) indicate that diffusion by tunnelling is temperature independent. From Figure 1 we see that this agrees with the present calculations up to say 10 K for H atom tunnelling, but for higher temperatures the tunnelling rate becomes more and more dependent on temperature. For the D atom, the transition from temperature independent to temperature dependent tunnelling occurs at a slightly lower temperature than for the H atom.