Long-Time Scale Simulations of Tunneling ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Long-Timescale Simulations of Tunneling-Assisted Diffusion of Hydrogen on Ice Surfaces at Low Temperature Vilhjálmur Ásgeirsson, Hannes Jonsson, and Kjartan Thor Wikfeldt J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b10636 • Publication Date (Web): 01 Jan 2017 Downloaded from http://pubs.acs.org on January 3, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Long-Timescale Simulations of Tunneling-Assisted Diffusion of Hydrogen on Ice Surfaces at Low Temperature † ´ V. Asgeirsson, H. Jónsson,¶ and K. T. Wikfeldt∗,†

†Science Institute, University of Iceland, 107 Reykjav´ık, Iceland ‡Mulliken Center for Theoretical Chemistry, University of Bonn, Beringstr. 4, 53115, Bonn, Germany ¶Faculty of Physical Sciences and Science Institute, University of Iceland, 107 Reykjav´ık, Iceland E-mail: [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Hydrogen (H) atom diffusion on dust grain surfaces is the rate-limiting step in many hydrogenation reactions taking place in interstellar clouds. In cold (10–30 K) molecular clouds, the dust grains are coated by amorphous water ice. Therefore, H adatom mobility on ice surfaces is of fundamental importance in this context. We have calculated H atom adsorption and diffusion on both crystalline and amorphous ice surfaces using an analytic interaction potential for H2 O-H. Tunneling rates for H atom hops between adsorption sites are explicitly calculated, the kinetic Monte Carlo method is used to simulate long timescale evolution and the diffusion coefficient, D, is evaluated for the temperature range 5-120 K. For ice Ih we find D = 1.6 · 10−7 cm2 /s at 10 K and below that temperature tunneling becomes the dominant diffusion mechanism. On the amorphous ice surface, the mobility of H is much slower than for ice Ih, D = 5.8 · 10−11 cm2 /s at 25 K. Below 25 K the H adatom becomes trapped in the deepest adsorption site which is located in a small surface pore. Furthermore, by blocking this site the diffusion increases by several orders of magnitude. H2 formation is thus likely to take place in deep adsorption sites at the low coverage and low temperature characteristic of molecular clouds. Deep adsorption sites can also explain experimental observations indicating that tunneling does not significantly contribute to the diffusivity, even though H adatoms can tunnel between shallow adsorption sites, the rate-determining transitions out of deep sites require thermal activation. Furthermore, in the context of coarse-grained astrochemical models, we find the ratio of the activation energy of diffusion and the adsorption energy of H to be 0.64 on amorphous water ice.

Introduction The adsorption and diffusion of hydrogen atoms on the surfaces of solids is a fundamental process that contributes to the chemistry of a diverse range of technological, environmental and astrophysical systems. However, while the electronic structure of H is significantly 2


Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simpler than that of heavier atoms, the opposite is true when it comes to nuclear dynamics. This is because of the small mass of the H atom resulting in a notably larger impact of nuclear quantum effects (NQEs) e.g., quantum tunneling and zero-point motion, which complicates the theoretical treatment of H atom dynamics. NQEs of H atoms have been found to play an important role not only for surface diffusion, but also associative desorption of H2 1 and in a diverse range of systems including enzymatic reactions, 2 high-pressure ice 3 and other H-bonded crystals and ferroelectrics. 4,5 Surface diffusion on metals has been extensively investigated both experimentally and in theoretical simulations 6–9 providing valuable insights into the tunneling behavior of H atoms. However, while H adsorption and diffusion on metals is an important process in various technological applications, quantum tunneling typically becomes dominant only at a temperature below ∼100 K, which is far below the working conditions of these technologies. In contrast, tunneling is believed to be the dominant mechanism in various chemical reactions taking place on dust grain surfaces, in interstellar molecular clouds, where the temperature is typically around 10 to 30 K and thermally activated processes are therefore strongly suppressed. In the case of H atoms, however, there has been a long–standing debate on the exact nature of H atom mobility on grain surfaces in molecular clouds, 10 particularly since gas-phase formation routes of H2 are too inefficient to explain the observed abundance. 11 Hence, surface diffusion of H atoms will be the ratelimiting step in the formation of H2 . 12,13 The interaction between physisorbed H and the surface is weak, so diffusion barriers are low and thermal diffusion cannot be excluded. In addition, a “hot-atom” mechanism has been proposed based on the idea that thermalization of impinging H atoms may be sufficiently slow for the H atom to rapidly scan the surface and find reaction partners before dissipating its excess kinetic energy. Understanding the underlying mechanisms of H atom mobility on dust grain surfaces is crucial not only for H2 formation, the most common and fundamental of all astrochemical reactions, but also for various hydrogenation reactions. These sub-micron sized dust grains contain a carbonaceous or a silicate core covered with an icy mantle predominantly

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

composed of H2 O molecules in an amorphous solid water (ASW) phase. Many previous studies of H2 formation have therefore focused on H diffusion on ASW surfaces, although crystalline ice surfaces have sometimes been investigated in tandem, for comparison. Moreover, the morphology of interstellar ASW is not known, and both porous and non-porous ASW surfaces have been investigated. An early picture proposed by Smoluchowski 14,15 suggested that H atom tunneling should not be effective on ASW due to the poor ground state energy level matching, resulting from the wide distribution of well depths. While on crystalline ice surfaces tunneling should proceed efficiently. Experiments based on temperature programmed desorption (TPD) both support and contradict this picture. 10 However, it has been pointed out that TPD may not be able to clearly characterize the H diffusion mechanism on ice surfaces. 16 In recent work, 16–18 Watanabe and coworkers have developed a new experimental procedure combining photostimulated desorption (PSD) with resonanceenhanced multiphoton ionization (REMPI) to overcome the limitations of TPD and measure the diffusion-limited recombination of H and D atoms on both ASW and polycrystalline ice surfaces. Even at a temperature as low as 8 K, only a weak isotope effect was observed on ASW surfaces, 17 suggesting that thermal hops dominate diffusion on ASW. In their most recent work, on the other hand, the first clear experimental evidence was presented for tunneling-dominated diffusion on polycrystalline ice surfaces at 10 K. 16 Molecular simulations can complement experiments by providing detailed atomistic information on the diffusion of H atoms on ice surfaces, and a number of previous studies have applied molecular dynamics (MD) simulations to investigate the sticking coefficient, adsorption and diffusion of H on ASW. 19–21 MD is useful to determine the temperaturedependent sticking coefficient as a function of temperature, but two key aspects of H atom diffusion on ice are beyond the capabilities of standard MD methods. First, MD is based on a classical treatment of the atomic nuclei and therefore does not include NQEs. Second, MD is limited to short timescales and has therefore not been used to investigate dynamics much beyond the timescales of the hot-atom diffusion, which is on the order of picoseconds.

4


Page 4 of 31

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


NQEs in the adsorption and diffusion of H on ASW have been investigated theoretically, starting with the work of Buch and coworkers, who calculated the nuclear wave functions of adsorbed H atoms. 22 More recently, Senevirathne et al. 23 have performed harmonic quantum rate theory (’instanton’ theory) calculations of H on ASW and ice Ih surfaces. They report tunneling rates for typical H diffusion barriers and find the average crossover temperature to tunneling-dominated diffusion to be around 10 K for both ASW and ice Ih, suggesting that thermal hops dominate over tunneling at T>10 K. The motivation behind the present work is to go beyond these previous studies by applying a novel simulation approach that explicitly accounts for both NQEs and long-timescale dynamics in the adsorption and diffusion of H atoms on ice surfaces. We report distributions of adsorption energy, residence time, diffusion barriers and hop lengths. These data are then used as input to kinetic Monte Carlo (KMC) simulations of long-timescale diffusion, where the hopping rates between neighboring minima have been obtained by explicitly considering quantum tunneling contributions. In agreement with recent experimental results, 16 we find that tunneling contributes significantly to the diffusion on crystalline ice, only at temperature below approximately 10 K. However, the negligible tunneling contributions for the diffusion on ASW is not the absence of tunneling per se, but rather the fact that H becomes trapped in deep adsorption sites on the ASW surfaces. This article is organized as follows. In Section 2 we describe the methodology used in this work, including descriptions of the crystalline and amorphous model ice surfaces, the analytic interaction potential employed, classical and quantum transition (or hop) rates of H atoms and finally of the KMC algorithm, which is used for calculating the diffusion coefficients. In Section 3 we present our results on H atom adsorption sites, the minimum energy paths connecting these sites, rate constants and residence times, and classical and quantum diffusion constants. Finally, concluding remarks are given in Section 4.

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Methodology Two structure models for ice surfaces are used in this work, the basal plane of crystalline hexagonal ice (ice Ih) and a surface of low-density amorphous solid water (ASW). The two models consist of 96 water molecules with a unit cell of identical size, or approximately 13.5 ˚ A and 15.6 ˚ A in x- and y- directions, respectively. The ice Ih model is adapted from the work of Hayward and Reimers. 24 Among the possible ice bilayers that can be used to generate a basal plane surface, we select one which has a near-linear pattern of dangling hydrogen atoms, because this configuration has been found to be particularly stable. 25,26 15 ˚ A of vacuum are added normal to the basal plane and the resulting slab is structurally optimized using PBE-D3 27,28 as implemented in the CP2K program. 29,30 Furthermore, the calculation employs a plane-wave basis set, with a cutoff of 280 Ry and a molecularly-optimized DZVP (mDZVP) Gaussian basis set. 31 During the optimization procedure, the lower half of the molecules are kept frozen in the bulk geometry. The ASW model is generated using the classical SPC/E 32 force field by carrying out NPT molecular dynamic simulations of liquid water at ambient temperature and pressure (1 bar) and rapidly cooling the system until it glassifies as an amorphous solid. As in the case of crystalline ice, 15 ˚ A of vacuum are added normal to the ASW surface, to create a slab. The molecular structure of the slab is then refined by minimizing the energy obtained using the same procedure as discussed above, with the lower half of the water molecules fixed. The reason for using rather small model systems is to make it possible to carry out periodic DFT calculations of H atom diffusion and H2 formation reactions in the future.

Computational approach We employ an analytic interaction potential to compute adsorption energy, activation energy and rate constants. Specifically, we use the H2 O-H interaction potential developed by Andersson et al. 33 from fits to the ab initio YZCL2 potential energy surface 34 for the H2 O+H↔OH+H2 reaction. A correction 23 for a slight error in the expression for the inter6


Page 6 of 31

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


action potential in the original publication 33 was used in this work. The H2 O-H potential includes dispersion and repulsion terms between the single H atom and all atoms of the water molecule. In addition, a Morse potential is employed between the H atom and the oxygen atoms. Using this interaction potential, all adsorption sites are located by relaxing the H atom position from a set of starting points defined by periodic grids on the model crystalline and amorphous ice surfaces. Nudged elastic band (NEB) 35–37 calculations are then performed between all adjacent energy minima and the highest energy image on the minimum energy path (MEP) is used to create a good initial guess for a saddle point search using the minimum-mode following algorithm. 38,39 Additionally, vibrational analysis is performed at all stationary points (minima and saddle points) to yield the vibrational frequencies.

Thermal and quantum-corrected rate constants To calculate thermal rate constants for H adatom hopping between energy minima sites on the surfaces, we rely on harmonic transition state theory (HTST). In HTST the rate constant is expressed as

k

HTST

Q3N m 1 −Ea i=1 ωi = Q −1 SP exp 2π 3N kB T j=1 ωj

(1)

where Ea is the activation energy of of the transition, ωSP and ωm are the vibrational frequencies obtained from the eigenvalues of the Hessian matrix at the relevant stationary points, T is the temperature and kB is the Boltzmann constant. Generally, rate constants obtained by HTST can be considered to be an accurate, upper-bound estimate of the exact rate constant. However, we are interested in a temperature range extending down to 5 K, where quantum tunneling competes with thermal hops. Therefore, the rate constants obtained with HTST have to be corrected. An established approach to describe tunneling is the semi-classical Wentzel-Kramer-Brillouin (WKB) method, which estimates the quantum transmission probability through a one dimensional potential barrier. A quantum correction

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 31

factor is calculated as the ratio of the thermally averaged quantum transmission probability and classical transmission probability over the same potential barrier. The WKB rate is then expressed as the product of the correction factor and the HTST rate, or by

kWKB = R ∞ 0

R∞

P (E) exp(−Ea /kB T )dE kHTST . H(E − VSP ) exp(−Ea /kB T )dE 0

(2)

H is the Heaviside function and is included because the classical transmission probability is zero for potential energy lower than the saddle point energy. P (E) denotes the quantum transmission probability and is given by the relation

P (E) =

1 1 + exp (2θ(E))

(3)

where θ(E) is referred to as the barrier penetration integral and expressed as 1 θ(E) = ~

Z

s2 s1

p 2µ(s)(V (s) − E)ds.

(4)

The integral extends over the classically forbidden region V (s1 ) = V (s2 ) = E and µ(s) represents the effective mass along the reaction coordinate. In this work only a single H atom is considered to participate in the tunneling process (the substrate remains fixed) and therefore the effective mass is replaced with the atomic mass of H. From the equations above, it is apparent that the WKB approach accounts for the dependence of tunneling on the mass and energy of the tunneling atom and the height and width of the potential barrier. A second more elaborate approach to describe quantum tunneling is harmonic quantum transition state theory (HQTST), which is also referred to as instanton theory, 40,41 and is based on statistical Feynman path-integrals. HQTST can first of all be used to estimate the crossover temperature, Tc , below which tunneling becomes the dominant mechanism. Below Tc the ring-polymer spreads over the barrier, describing the tunneling path of a quantum particle. Tc can be calculated from the curvature along the minimum energy path using the

8


Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


expression 42 Tc =

~ω 2πkB

(5)

where ω is the magnitude of the imaginary frequency at the saddle point. As the temperature is lowered below Tc , the Feynman paths open up and tend to slide down towards the energy minima, reducing the effective tunneling barrier. Furthermore, the tunneling path is often found to be shorter than the minimum energy path, an effect known as corner-cutting. The computational effort of estimating the transition rate using HQTST is substantially larger than the WKB calculation. While we present estimates of Tc for H atom diffusion below, the transition rates are evaluated using WKB in the present simulations.

Average residence time If the rate constants are known for all transitions that transport the H adatom away from a particular adsorption site, it becomes possible to calculate the average residence time at the site, by the expression 1 tres = PN

i=1

ki

(6)

where N is the total number of possible transitions. We calculate the average residence times of the H adatom using both classical HTST and quantum WKB-corrected rates.

Long-timescale dynamics and diffusion coefficients Given that the rate constants are known for all transitions out of every adsorption site on the surface, it becomes possible to simulate the time evolution of the system using a kinetic Monte Carlo (KMC) algorithm. We use KMC as implemented in the EON software 43 to simulate H diffusion via thermal hops and tunneling on both ice Ih and ASW. The application of KMC in our case is particularly simple since diffusion of only a single H atom is simulated, corresponding to a realistic description of conditions in molecular clouds where the flux of gasphase species onto the surface is extremely low. 10,53 For each system, multiple simulations

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

are performed at a temperature ranging from 5 K to 120 K, and all simulations are continued for 107 MC steps. The timescale reached by each simulation is governed by the input rate constants and hence the simulation temperature. As an example, the simulation time for ice Ih is 10−6 seconds at 120 K and 102 seconds at 5 K. From a KMC simulation, the temperature dependent diffusion coefficient, D, is calculated from the mean squared displacement (MSD) as

D=

hr(t) − r(0)i2 2dt

(7)

where the dimensionality is d = 2, appropriate for surface diffusion. KMC simulations have non-equidistant time steps, which makes the MSD evaluation more complicated than from regular MD trajectories. In a first step we generate a complete trajectory of the H adatom, where hops outside of the simulation box are allowed. The resulting trajectory is partitioned to ten equidistant parts, where each part is sampled by selecting all successive points along the (partitioned) trajectory and calculating the relative distances to all succeeding points, the distances are counted and averaged in a predefined time bin. The vertical error bars are estimated from the standard deviation of the 10 partitions.

Results and discussion The two different ice surfaces used in this study, i.e. the (0001) basal plane of hexagonal ice (ice Ih) and the amorphous (ASW) surface, are both visualized in Figs. 1(a)–(b). Oxygen atoms in ice Ih are arranged in a crystalline structure, but the orientation of the water molecules is random (apart from obeying the ice rules 44 ) and the hydrogen atoms are therefore disordered. Half of the upper half-bilayer surface molecules have one of their H atoms pointing out of the surface, these atoms are referred to as dangling H atoms and their ordering has been shown to affect the stability of ice Ih surfaces. 25,45 A perfect Fletcher-phase is a linear ordering of dangling H atoms and has been found 10


Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to be the lowest energy structure of the basal plane surface. This structure minimizes the dangling-H order parameter COH 25 and is assigned a value of 2.0, while a completely random (and less stable) arrangement gives COH ≈ 3.0. The basal plane surface shown in Fig.1(a) has COH ≈ 2.2 and should therefore represent a realistic hexagonal ice surface. Because of the proton disorder, even the ice Ih crystal gives rise to an interaction energy surface for the diffusing H adatom that resembles that of an amorphous solid, as has been noted previously for the case of H2 O and CO diffusion. 46–48 The ASW surface is, on the other hand, completely disordered and contains two shallow surface pores. Several previous studies 19,51 have investigated the sticking coefficient of H on ASW and ice Ih. A typical dependence of the sticking coefficient on the incident energy was found, where the sticking probability increases with decreasing incidence energy of the H atom. The sticking coefficient is close to unity at low incidence energy (∼ 10 K). All results reported in the following are based on a physical picture in which an H atom lands on the ice surface and becomes thermalized in a stable adsorption site.

Adsorption sites The adsorption sites found on the ice Ih surface are visualized in Fig.1(c). A total of 24 adsorption sites are found but only the deep adsorption sites are visualized. Figure 2(a) shows the distribution of the adsorption energy. It is important to note that the adsorption energies have not been corrected to account for zero point energy and are thus not directly comparable to experiments. The adsorption sites can be divided into two groups based on the spatial position of the H adatom, and Fig. 2(a) shows a distinct difference in the average adsorption energy. The first group consists of 12 deep minima which are located above the center of each hexamer, shown by green spheres in Fig. 1(c). The average adsorption energy for these sites is 41 meV. The results agree qualitatively with Al-Halabi et al., 19 who found the energy minima to be located at the center of each hexamer on ice Ih with an average adsorption energy of 35±1 meV using a different interaction potential developed by Buch and 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 12 of 31

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Zhang. 51 The second group consists of more shallow energy minima located above oxygen atoms at the corner of each hexamer, where bifurcations occur in MEPs emanating from the deep adsorption sites. On the ASW surface, a total of 12 adsorption sites are found with an average adsorption energy of 46 meV (see Fig. 1(d) and Fig. 2(b)). The deepest adsorption site is located in a surface pore on the ASW surface where the H atom interacts with many neighboring water molecules. The distribution in adsorption energy on ASW is significantly wider than for the ice Ih surface, reflecting the more heterogeneous nature of the ASW surface.

(a)

Figure 2: surfaces.

(b)

Distribution of adsorption energy for a H adatom on (a) ice Ih and (b) ASW

The results outlined above are in agreement with the work of Senevirathne et al. (2014), 23 who employed the same analytic interaction potential on larger structural models (360 H2 O molecules) of ice Ih and ASW. They also found two types of energy minima on ice Ih with an adsorption energy ranging from 42-47 meV and 34-38 meV. Their calculations on the ASW surface yields adsorption energy in the range of 39 to 89 meV, close to what is found here.

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Minimum energy paths Activation energy for diffusion hops, Ea , are key to the overall surface mobility of H and we now turn to an analysis of the minimum energy paths obtained for the ice Ih and ASW surfaces. Fig. 1(c)-(d) show the connectivity between all adsorption sites on both surfaces. The corresponding distribution in Ea is shown in Fig. 3. In the case of ice Ih, all energy barriers are seen to be low, ranging from 6 to 12 meV. In the case of ASW, the calculated activation energy of diffusive hops varies from (close to) 0 to 56 meV. The reported energy barriers have not been corrected to account for zero point energy. It is clear, that a much wider distribution in activation energy is found for ASW compared to that of ice Ih. Again, reflecting the disordered structure of ASW. Watanabe et al. 18 concluded from their experiments that two groups of adsorption sites were present for H on ASW, where one group has an average activation energy of diffusion of 18±2 meV and the other an activation energy of 50 meV. Hama et al. 17 further refined the work of Watanabe et al. and suggested there were three distinct types of H adsorption sites on ASW: shallow, middle and deep sites, with an average activation energy of diffusion of 18±2, 22±1 and exceeding 30 meV, respectively. These experimental findings, as well as the results of other experimental work, 52 are within the range of the results presented here. However, we find no evidence for different well-defined subclasses of H adsorption sites on ASW surfaces. Senevirathne et al. 23 report activation energy of diffusion of up to 14 meV for ice Ih and up to 69 meV for ASW. The difference can be attributed to the larger surface size employed by Senevirathne et al. and possibly to differences in surface morphology. The activation energy determines the hopping rates of H (along with the vibrational partition functions of the initial minima and saddle points) within the framework of HTST, but when quantum tunneling starts to contribute significantly at low temperature, the shape and width of the potential energy barrier also becomes important, as outlined in Section 2. Distributions of MEP lengths (i.e. hop distances) for both ASW and ice Ih are shown in Fig. 3(c). The total length of the reaction coordinate, referred to as a hop distance, on ice 14


Page 14 of 31

Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ih is on average 5.5 ˚ Awith a very narrow spread, while the hop distance on ASW can range from approximately 2.9 to 9.4 ˚ A.

(a)

(b)

(c)

Figure 3: Distribution of activation energy for diffusion between two adjacent sites on (a) ice Ih and (b) ASW. The distribution of H atom hop distances (total length of the MEPs) for the two surfaces is visualized in (c). Having computed the adsorption energy at the various sites along the ice Ih and ASW surfaces and the activation energy for all transitions (hops) out of each adsorption site, it is interesting to consider how these two quantities correlate. The correlation may find use in coarse-grained astrochemical models, where the energy barrier for diffusion are often roughly approximated from readily available desorption barriers, i.e. EA,diff = cEdes , where c is a

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

fixed constant, usually assigned values between 0.3 and 0.8. To obtain a first measure of c for both surfaces, the relationship between the adsorption energy and the lowest diffusion barrier are investigated and a linear fit is made, see Fig.4. The coefficients of determination (R-squared) are found to be 0.71 and 0.74 and the respective slopes of the fit (or the values of c) are found to be 0.68 and 0.64 for ice Ih and ASW, respectively. Interestingly, the values differ slightly from a previously reported value, c = 0.78, by Katz et al., 49 where they fitted TPD data to rate-equations for a H adatom on an amorphous carbon surface. More recently, Karssemeijer and Cuppen 50 predict c to be 0.34 and 0.42 for CO and CO2 on a proton disordered ice Ih surface, where the water molecules have been frozen. Their prediction is significantly lower than our corresponding values for H on a disordered ice Ih surface.

(a)

(b)

Figure 4: Correlation between adsorption energy and the lowest activation energy for a diffusion hop from each adsorption site, shown for the (a) ice Ih and (b) ASW surfaces. The lines represent linear fits with a slope of 0.68 for ice Ih and 0.64 for ASW.

Classical and quantum rates In the low temperature range of cold molecular clouds (10–30 K), most thermally activated processes are inactive and quantum tunneling becomes the dominant reaction mechanism. However, for for the transitions with the lowest activation energy, thermal hops are likely 16


Page 17 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to compete with tunneling. Recent experiments 16 report the co-existence of tunneling and thermal diffusion of H atoms on ASW at 10 K, where thermal diffusion is proposed to be the dominant mechanism. Fig. 5 shows the rate constants obtained using HTST and WKB methods, outlined in Section 2, as a function of temperature for a selected typical diffusion hop on the ASW surface, both the forward and backward rates are depicted. The WKB calculation clearly captures the typical behavior of quantum tunneling with a crossover to tunneling-dominated hops occurring at around 10 K. At high temperature the WKB-corrected and HTST rate constants are approximately equal, as expected, but as the temperature is lowered the WKB becomes significantly higher than the classical HTST rate. The transition where the H atom hops from the weaker to the stronger binding site has a temperature-independent hop rate, for low temperature. In this deep tunneling regime the tunneling rate is orders of magnitude larger than the rate of the exponentially decreasing thermal hops. In the opposite direction, as H hops from the deeper to the more shallow minimum, an activation energy of roughly 15 meV is required even in the deep tunneling regime, because the final state is higher in energy than the initial state. The temperature where kWKB and kHTST start to diverge corresponds to the crossover temperature Tc of a barrier. This crossover temperature can also be estimated independently from the Feynman path integral formulation using Eq. 5, as stated in Section 2. The estimated Tc clearly agrees well with the change in slope of the WKB calculated rates. By applying this estimate to all the diffusion hops on the ASW and ice Ih surfaces, we find that Tc ranges from (close to) 0 K up to 9 K for ASW, and up to a maximum Tc of 11 K for ice Ih. This suggests that although tunneling is likely to contribute to some extent above 10 K, thermal activation is still the dominant diffusion mechanism in the typical temperature range of dense molecular clouds.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5: The HTST rate constant (dashed lines) and WKB corrected rate constants (solid lines) are plotted for a single selected reaction path along the ASW surface as a function of temperature. Both the forward and backward rates are shown individually, and denoted by numbers 1 and 2 respectively. The dashed line shows the estimated crossover temperature between classical and quantum tunneling mechanisms, as estimated from the Feynman path integral formulation, Eq. 5.

Residence times The average residence time of H at the adsorption sites contains valuable information on the mobility of H on the two surfaces and allows us to analyze further the impact of tunneling at low temperature. Fig. 6 show the average residence time of H in each adsorption site as a function of temperature. The visualized residence times are calculated separately using WKB-corrected and HTST rate constants. The average residence times on ice Ih have a narrow spread among the different adsorption sites, which are also all very short, e.g., the average residence time of H in all of the adsorption sites, at a temperature of 5 K, is below 1 ms, regardless of whether tunneling is accounted for or not. Residence times on ASW follow a different trend and range from roughly 10−10 s to 1025 s at 5 K (for HTST), showing that H atoms can become trapped at a low temperature in deep adsorption sites. Such deep adsorption sites, not found on ice Ih surfaces, are a consequence of the heterogeneous nature of the ASW surface and the resulting wide range

18


Page 18 of 31

Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of adsorption energies, as discussed above. Furthermore, Fig. 6 allows us to further inspect the importance of tunneling for each system by comparing residence times obtained from the WKB-corrected and HTST rates. As expected, tunneling lowers the average residence times and thus increases the H mobility on both surfaces. However, in the case of ice Ih the impact of tunneling on the average residence time of H is rather small and only becomes noticeable below 10 K. In contrast, for ASW the residence times calculated from WKB-corrected and HTST rate constants can differ by many orders of magnitude at low temperature. These results reveal that tunneling can have a large effect on individual H atom hops on ASW. This, however, is unrelated to long-distance surface diffusion, as will be discussed next.

(a)

(b)

Figure 6: Average residence time, tR , of H in each adsorption site as a function of temperature for (a) ice Ih and (b) ASW. The same logarithmic time scale is used for both surfaces, for clarity. Red curves are calculated using the WKB-corrected rate constants and the blue curves using the HTST rate constants. At and below 10 K, the residence time can be significantly reduced by the inclusion of tunneling.

KMC and diffusion coefficients Long time-scale simulations of both thermal (classical) and tunneling-assisted H diffusion were carried out using KMC, for both ice Ih and ASW. The results are shown in Fig.7. The diffusion coefficient for a H adatom on ice Ih is estimated to be D = 5.5·10−3 cm2 /s at 120 K. 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The explicit inclusion of tunneling does not affect the calculated diffusion coefficient down to a temperature of roughly 10 K, where D = 1.6 · 10−7 cm2 /s. At this point, the HTST and WKB predicted diffusion coefficients begin to deviate from one another, as tunneling becomes the dominating diffusion mechanism. For example, at 5 K, the tunneling-assisted coefficient is D = 4.8 · 10−11 cm2 /s, approximately one order of magnitude larger than the classical diffusion coefficient. These results show that the H adatom is mobile down to 5 K on ice Ih and that tunneling plays an important role below 10 K. We attribute this to the relative homogeneity of the ice Ih surface and the low activation energy for hops between available adsorption sites. For ASW, the diffusion coefficient is found to be D = 1.7 · 10−4 cm2 /s at 120 K, for both the classical and the WKB-corrected rate constants. This value is one order of magnitude lower than that of ice Ih, reflecting the difference in average activation energy between the two ice surfaces. Furthermore, as the simulation temperature is lowered down to 25 K, both tunneling and classical diffusion coefficients remain nearly identical with D = 5.8 · 10−11 cm2 /s, which is similar to the value of D for ice Ih at 5 K. It should be noted that the limited size of the model ice surfaces used here may have an effect on the magnitude of the long-timescale diffusion. However, the temperature dependence and the comparison between the crystalline and amorphous surfaces should be less influenced by such finite-size effects. In a previous simulation study of H diffusion on an ASW surface, Al-Halabi et al. 20 estimated the diffusion coefficient to be 1.09 ± 0.04 · 10−5 cm2 /s at 10 K. Therefore, they predict the diffusion coefficient to be several orders of magnitude larger than our esimate at 25 K. Later, it was pointed out by Hama et al. 17 that this value corresponds to an average diffusion barrier of 5 meV and cannot be consistent with their reported average binding energy of 56 meV. This discrepancy may be related to the very short simulation timescale employed in their classical MD simulations, 20 which may prevent the H atom to fully explore the surface and find deep adsorption sites. For simulations at temperature below 25 K, the H atom becomes stuck in the deepest

20


Page 20 of 31

Page 21 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


adsorption sites and the diffusion coefficient drops to nearly zero. It can be concluded that even though tunneling can greatly enhance the low-temperature hop rate between neighboring adsorption sites, as was discussed in connection to Fig. 6 above, the H atom becomes so strongly trapped in deep adsorption sites that any long-distance diffusion is suppressed. This raises the question, what happens if the deepest adsorption sites are blocked, e.g., by an inert species? We simulate this situation by removing the most strongly bound adsorption sites from the KMC simulations. Considering the total number of adsorption sites on the ASW surface, this corresponds to roughly 17% of the sites, being occupied by an inert species. We find that the removal of the deepest adsorption sites greatly increases the H atom mobility, see Fig. 7(b). For instance, the diffusion coefficient is found to be around 3.3 ·10−7 cm2 /s at 25 K, four orders of magnitude larger than when the deepest adsorption sites are not blocked. This trapping mechanism was also observed in a recent, combined theoretical (MD) and experimental, study of bulk diffusion of NH3 , CO2 , CO and H2 CO on low-density ASW. 55 They report that a majority of molecules become trapped in the ice while diffusing and suggest that such trapping phenomena can be reflected by a decrease in the diffusion coefficient. Furthermore, Karssemeijer et al. 56 showed, using adaptive KMC simulations, 57 that CO mobility is strongly dependent on the existence of nanopores on ASW, where the CO becomes trapped in deep adsorption sites at low coverage. If the CO coverage increases, the deep sites become occupied and the diffusion coefficient is increased by several orders of magnitude, in agreement to the observations made here for H.

Conclusions We have performed a simulation study of the adsorption and diffusion of H atoms on both crystalline and amorphous ice surfaces. The latter surface type serves as a model for grain surfaces found in cold and dense interstellar molecular clouds, where H2 formation and various other fundamental hydrogenation reactions take place. The mobility of H atoms tends

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 22 of 31

(b)

Figure 7: Diffusion coefficients obtained from multiple KMC simulations as function of temperature. The coefficients are obtained from simulations employing both classical (HTST) and tunneling-corrected (WKB) rate constants for the various diffusion hops. Results for ice Ih and ASW are shown in panels (a) and (b), respectively. The two additional curves in (b), labeled by ”(w/o 2)” denote results from KMC simulations where the two deepest adsorption sites are blocked. to be the rate limiting factor on such surfaces. Compared to earlier theoretical studies, we present several improvements. First, we located all stable adsorption sites and map out their connectivity and associated minimum energy paths. Using this information, we compute both classical (using harmonic transition state theory) and tunneling-assisted (within the WKB approximation) transition (or hop) rates between the adsorption sites and use these to perform long-timescale KMC simulations for H atom diffusion in a large temperature interval. In this way, both the impact of quantum tunneling and the long-timescale dynamical behavior is described within the same framework. We find that the ASW surface has a wide distribution in adsorption energy (32-77 meV) and activation energy (1-56 meV). Even the crystalline ice surface is found to have a significant, although narrower distributions of adsorption energy (34-43 meV) and activation energy (6-12 meV) because of the proton disorder. We find that the ratio between the activation energy of diffusion and adsorption energy is 0.68 and 0.64 for H atoms on ice Ih and ASW, respectively. These values may find use in coarse-grained astrochemical models. 22


Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For ice Ih the classical and tunneling-assisted diffusion coefficients are observed to be of similar magnitude down to 10 K, with D = 1.6 · 10−7 cm2 /s at 10 K. At lower temperature tunneling takes over as the dominant diffusion mechanism and the H atom is observed to be mobile down to 5 K. In the case of ASW, the diffusion coefficient is found to be 5.8 · 10−11 cm2 /s at 25 K, several orders of magnitude lower than that reported value by Al-Halabi et al. 20 For lower temperature, the H atom becomes trapped in deep adsorption sites on the ASW surface and the diffusion coefficient tends to zero, irrespective of whether classical or tunneling-assisted hopping rates are used in the KMC simulations. Thus, although individual nearest-neighbor hops can be greatly facilitated by tunneling on ASW (see Fig. 5), and the average residence times of H in many energy minima on ASW are more strongly affected by tunneling compared to on the ice Ih surface (see Fig. 6), hops out of deep adsorption sites are rate-limiting and determine the overall diffusion constant. From this we conclude that tunneling is not significant for the long-distance diffusion rate of H atoms on ASW surfaces where deep adsorption sites are present. However, by blocking the deepest adsorption sites on the ASW surface, the diffusion coefficient increases by several orders of magnitude, in agreement to the findings of Karssemeijer et al. 56 for CO adsorption on ASW. This corresponds to a physical picture in which deep adsorption sites are occupied, e.g., by an inert species such as an H2 molecule formed in an earlier formation reaction. Since H2 formation is dominated by the Langmuir-Hinshelwood process in interstellar molecular clouds where the flux of gasphase species is extremely low, 10,53 it seems plausible that most H2 formation reactions take place in deep adsorption sites when a mobile H atom encounters an immobile H atom in a deep site. Some fraction of H2 molecules formed in this way will be able to dissipate the energy released in the exothermic reaction into the substrate and remain adsorbed in the deep site. This picture agrees well with conclusions drawn from recent experimental work by Watanabe and co-workers. 53 At the same time it rationalizes discrepancies observed previously between older experiments where the diffusion of H was reported to be either “fast” or “slow”. 52,54,58,59 At low coverage, H atoms

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

become trapped and slow diffusion is observed. At higher coverage, deep adsorption sites become occupied and the remaining H atoms are able to diffuse rapidly over longer distances. One implication of these observations is that in order to reliably measure H diffusion rates, it is important to carefully control the H atom coverage on the surface. Another more general implication is that understanding all underlying mechanisms that determine the occupation of deep sites, including the probability that H2 molecules formed in deep sites remain adsorbed there after the reaction, is critical in order to estimate reliable diffusion constants of H atoms on interstellar dust grain surfaces. A detailed study of the H2 formation mechanism on both surfaces is underway, where the reaction paths and formation barriers are calculated and the partitioning of the reaction energy is calculated using classical trajectories starting from the transition state, the second step in the two step WKE procedure. 1

Acknowledgements ´ Asgeirsson thanks Christoph Bauer for helpful discussions. This work was supported by the Icelandic Research Fund through Grant No. 141080-052. Computational resources granted through the Swedish National Infrastructure for Computing (SNIC) at the HPC2N center, and through University of Iceland Computer Services are gratefully acknowledged.

References (1) Jónsson, H. Simulation of surface processes. Proc. Natl. Acad. Sci. 2011, 108, 944. (2) Bothma, J. P.; Gilmore, J. B.; McKenzie, R. H. The role of quantum effects in proton transfer reactions in enzymes: quantum tunneling in a noisy environment. New J. Phys. 2010, 12, 055002.

24


Page 24 of 31

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(3) Benoit, M.; Marx, D.; Parrinello, M. Tunnelling and zero-point motion in high-pressure ice. Nature 1998, 392, 258. (4) Li, X. Z.; Walker, B.; Michaelides, A. Quantum nature of the hydrogen bond. Proc. Natl. Acad. Sci. (USA) 2011, 108, 6369. (5) Wikfeldt, K.; Michaelides, A. Communication: Ab initio simulations of hydrogenbonded ferroelectrics: Collective tunneling and the origin of geometrical isotope effect. J. Chem. Phys. 2014, 140, 041103. (6) Lauhon, L.; Ho, W. Direct observation of the quantum tunneling of single hydrogen atoms with a scanning tunneling microscope. Phys. Rev. Lett. 2000, 85, 4566. (7) Jardine, A.; Lee, E.; Ward, D.; Alexandrowicz, G.; Hedgeland, H.; Allison, W.; Ellis, J.; Pollak, E. Determination of the quantum contribution to the activated motion of hydrogen on a metal surface: H/Pt (111). Phys. Rev. Lett. 2010, 105, 136101. (8) Suleimanov, Y. V. Surface diffusion of hydrogen on Ni (100) from ring polymer molecular dynamics. J. Phys. Chem. C 2012, 116, 11141. (9) McIntosh, E. M.; Wikfeldt, K. T.; Ellis, J.; Michaelides, A.; Allison, W. Quantum effects in the diffusion of hydrogen on Ru (0001). J. Phys. Chem. Lett. 2013, 4, 1565. (10) Vidali, G. Cosmic low temperature physics: making molecules on stardust. J. Low Temp. Phys. 2013, 170, 1. (11) Gould, R. J.; Salpeter, E. E. The interstellar abundance of the hydrogen molecule. I. Basic processes. Astrophys. J. 1963, 138, 393. (12) Hollenbach, D.; Salpeter, E. Surface adsorption of light gas atoms. J. Chem. Phys. 1970, 53, 79. (13) Hollenbach, D.; Salpeter, E. Surface recombination of hydrogen molecules. Astrophys. J. 1971, 163, 155. 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) Smoluchowski, R. Rate of H2 formation on amorphous grains. Astrophys. Space Sci. 1981, 75, 353. (15) Smoluchowski, R. Adsorption and mobility on amorphous surfaces. Application to astrophysical problems. J. Phys. Chem. 1983, 87, 4229. (16) Kuwahata, K.; Hama, T.; Kouchi, A.; Watanabe, N. Signatures of quantum-tunneling diffusion of hydrogen atoms on water ice at 10 K. Phys. Rev. Lett. 2015, 115, 133201. (17) Hama, T.; Kuwahata, K.; Watanabe, N.; Kouchi, A.; Kimura, Y.; Chigai, T.; Pirronello, V. The mechanism of surface diffusion of H and D atoms on amorphous solid water: Existence of various potential sites. Astrophys. J. 2012, 757, 185. (18) Watanabe, N.; Kimura, Y.; Kouchi, A.; Chigai, T.; Hama, T.; Pirronello, V. Direct measurements of hydrogen atom diffusion and the spin temperature of nascent H2 molecule on amorphous solid water. Astrophys. J. Lett. 2010, 714, L233. (19) Al-Halabi, A.; Kleyn, A.; Van Dishoeck, E.; Kroes, G. Sticking of hydrogen atoms to crystalline ice surfaces: Dependence on incidence energy and surface temperature. J. Phys. Chem. B 2002, 106, 6515. (20) Al-Halabi, A.; Van Dishoeck, E. F. Hydrogen adsorption and diffusion on amorphous solid water ice. Mon. Not. Royal Astron. Soc. 2007, 382, 1648. (21) Veeraghattam, V. K.; Manrodt, K.; Lewis, S. P.; Stancil, P. The sticking of atomic hydrogen on amorphous water ice. Astrophys. J. 2014, 790, 4. (22) Buch, V.; Czerminski, R. Eigenstates of a quantum-mechanical particle on a topologically disordered surface: H (D) atom physisorbed on an amorphous ice cluster (H2 O)115 . J. Chem. Phys. 1991, 95, 6026.

26


Page 26 of 31

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(23) Senevirathne, B.; Andersson, S.; Dulieu, F.; Nyman, G. Hydrogen atom mobility, kinetic isotope effects and tunneling on interstellar ices (Ih and ASW). Mol. Astrophys., in press. (24) Hayward, J.; Reimers, J. Unit cells for the simulation of hexagonal ice. J. Chem. Phys. 1997, 106, 1518. (25) Pan, D.; Liu, L. M.; Tribello, G. A.; Slater, B.; Michaelides, A.; Wang, E. Surface energy and surface proton order of ice Ih. Phys. Rev. Lett. 2008, 101, 155703. (26) Pedersen, A.; Karssemeijer, L. J.; Cuppen, H.; Jónsson, H. Long-timescale simulations of H2 O admolecule diffusion on ice Ih(0001) surfaces. J. Phys. Chem. C 2015, 119, 16528. (27) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104 (28) Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865 (29) The CP2K Developers Group. Available at: http://www.cp2k.org/ (30) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103. (31) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105 (32) Dyer, K. M.; Perkyns, J. S.; Stell, G.;Montgomery, P. B. Site-renormalised molecular fluid theory: on the utility of a two-site model of water. Mol. Phys. 2009, 107, 423

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(33) Andersson, S.; Al-Halabi, A.; Kroes, G. J.; van Dishoeck, E. F. Molecular-dynamics study of photodissociation of water in crystalline and amorphous ices. J. Chem. Phys. 2006, 124, 064715. (34) Yang, M.; Zhang, D. H.; Collins, M. A.; Lee, S. Y. Ab initio potential-energy surfaces for the reactions OH+ H2 ↔ H2 O+ H. J. Chem. Phys. 2001, 115, 174. (35) G. Mills, H. Jónsson and G. K. Schenter, Reversible work based transition state theory: Application to H2 dissociative adsorption. Surface Science 1995, 324, 305. (36) H. Jónsson, H. G. Mills, G. Jacobsen, K. W. Nudged elastic band method for finding minimum energy paths of transitions. in Classical and Quantum Dynamics in Condensed Phase Simulations, ed. B. J. Berne, G. Ciccotti and D. F. Coker (World Scientific, 1998), page 385. (37) Henkelman, G.; Jónsson, H. Improved tangent estimate in the NEB method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978. (38) Henkelman, G.; Jónsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys. 1999, 111, 7010. (39) Olsen, R.; Kroes, G. J.; Henkelman, G.; Jónsson, H. Comparison of methods for finding saddle points without knowledge of the final states. J. Chem. Phys. 2004, 121, 9776. (40) Benderskii, V.A.; Goldanskii, V. I.; Makarov, D.E. Quantum dynamics in lowtemperature chemistry. Phys. Rep. 1993, 233, 195. (41) Andersson, S.; Nyman, G.; Arnaldsson, A.; Manthe, U.; Jónsson, H. Comparison of quantum dynamics and quantum transition state theory estimates of the H + CH4 reaction rate. J. Phys. Chem. A 2009, 113, 4468.

28


Page 28 of 31

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(42) Gillan, M. Quantum-classical crossover of the transition rate in the damped double well. J. Phys. C: Solid State Phys. 1987, 20, 3621. (43) Chill, S.T.; Welborn, M.; Terrell, R.; Zhang, L.; Berthet, J-C.; Pedersen, A.; Jónsson, H.; Henkelman, G. EON: Software for long time simulations of atomic scale systems. Modelling Simul. Mater. Sci. Eng. 2014, 22, 055002. (44) Bernal, J.D.; Fowler, R.H. A theory of water and ionic soluation, with particular reference to hydrogen and hydroxyl ions. J. Chem. Phys. 1933, 1, 515 (45) Buch, V.; Groenzin, H.; Li, I.; Shultz, M. J.; Tosatti, E. Proton order in the ice crystal surface. Proc. Natl. Acad. Sci. 2008, 105, 5969. (46) Batista, E. R.; Jónsson, H. Diffusion and island formation on the ice Ih basal plane surface. Comp. Mat. Sci. 2001, 20, 325. (47) Batista, E. R.; Jónsson, H. Diffusion and island formation on the ice Ih basal plane surface. Comp. Mat. Sci., Corrigendum 2015, 102, 338. (48) Karssemeijer, L. J.; Pedersen, A.; Jónsson, H.; Cuppen, H. M. Long-timescale simulations of diffusion in molecular solids. Phys. Chem. Chem. Phys. 2012, 14, 10844. (49) Katz, N.; Furman, I.; Biham, O.; Pirronello, V.; Vidali, G. Molecular hydrogen formation on astrophysically relevant surfaces. Astrophys. J. 1999 522(1), 305. (50) Karssemeijer, L. J.; Cuppen, H. M. Diffusion-desorption ratio of adsorbed CO and CO2 on water ice. Astron. Astrophys., 2014, 569, A107. (51) Buch, V.; Zhang, Q. Sticking probability of H and D atoms on amorphous ice-A computational study Astrophys. J. 1991, 379, 647. (52) Matar, E.; Congiu, E.; Dulieu, F.; Momeni, A.; Lemaire, J. Mobility of D atoms on porous amorphous water ice surfaces under interstellar conditions Astron. Astrophys. 2008, 492, L17. 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(53) Hama, T.; Watanabe, N. Surface processes on interstellar amorphous solid water: Adsorption, diffusion, tunneling reactions, and nuclear-spin conversion Chem. Rev. 2013, 113, 8783. (54) Hornekær, L.; Baurichter, A.; Petrunin, V. V.; Field, D.; Luntz, A. C. Importance of surface morphology in interstellar H2 formation Science 2003, 302, 1943. (55) Ghesquière, P.; Mineva, T.; Talbi, D.; Theulé, P.; Noble, J. A.; Chiavassa, T. Diffusion of molecules in the bulk of a low density amorphous ice from molecular dynamics simulations PCCP 2015, 17(17), 11455-11468. (56) Karssemeijer, L. J.; Ioppolo, S.; van Hemert, M. C.; van der Avoird, A.; Allodi, M. A.; Blake, G. A.; Cuppen, H. M. Dynamics of CO in amorphous water-ice environments. Astrophys. J. 2013, 781(1), 16. (57) Henkelman, G.; Jónsson, H. Long time scale kinetic Monte Carlo simulations without lattice approximation and predefined event table. J. Chem. Phys. 2001, 115(21), 96579666. (58) Manico, G.; Raguni, G.; Pirronello, V.; Roser, J. E.; Vidali, G. (2001). Laboratory measurements of molecular hydrogen formation on amorphous water ice. Astrophys. J. Lett., 2001, 548(2), L253. (59) Perets, H. B.; Biham, O.; Manic, G.; Pirronello, V.; Roser, J.; Swords, S.; Vidali, G. Molecular hydrogen formation on ice under interstellar conditions. Astrophys. J., 2005, 627(2), 850.

30


Page 30 of 31

Page 31 of 31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Long-Time Scale Simulations of Tunneling ... - ACS Publications

Recommend Documents