Molecular Oxygen Formation in Interstellar Ices Does Not Require Tunneling

The formation of molecular oxygen in and on amorphous ice in the interstellar medium requires oxygen diffusion to take place. Recent experiments suggest that this process involves quantum tunneling of the oxygen atoms at sufficiently low temperatures. Fitting experimental diffusion rates between 6 and 25 K to an expression that accounts for the roughness of the surface yields excellent agreement. The molecular dynamics of adsorbed oxygen is characterized by rapid intrasite dynamics, followed by intersite transitions over distances of ∼10 Å. Explicit simulations using a realistic free-energy surface for oxygen diffusion on amorphous ice down to 10 K show that quantum tunneling is not required for mobility of adsorbed oxygen. This is confirmed by comparing quantum and classical simulations using the same free-energy surface. The ratio of diffusional and desorption energy Edif/Edes = 275/1082 ≈ 0.3 is at the lower end of typically used values but is still consistent with the assumptions made in models for interstellar chemistry. T diffusion of oxygen atoms in and on interstellar ice grains is of fundamental importance for the formation of CO2 and molecular oxygen O2. 1−3 In addition, O diffusion has also been implicated in the formation of ozone in laboratory experiments within a Langmuir−Hinshelwood mechanism (surface diffusion) and has been proposed to be relevant for interstellar grains, eventually leading to the formation of water, although the major water-formation pathway involves mobile hydrogen atoms. In all of these cases the diffusion of atomic oxygen is a primary driver for forming dior triatomic molecules under conditions typical for interstellar environments. These include higher temperatures (T ≈ 50 K) for translucent and diffuse clouds and temperatures T ≈ 10 K for dense molecular clouds. Questions of particular relevance pertain to the diffusional barrier for intersite migration of oxygen atoms on amorphous solid water (ASW) and their desorption energies away from the surface. For this, experiments based on temperature-programmed desorption have been used. For analyzing the diffusion rate (i.e., the probability for one hop between different adsorption sites in s−1, which is referred to as D in ref 4) of atomic oxygen on ASW, either an Arrhenius expression (D(T) ≈ Ae(−βEdif), where Edif is a conformationally averaged diffusional barrier and β = 1/(kBT)) or an expression accounting for tunneling through a square well barrier was used. A diffusion coefficient Ddif (in units of length /time) can be converted into a diffusion rate D (in 1/time) through D = Ddif/a , where a is a typical area sampled by the adsorbate between consecutive hops (see Figure S1). Typically, a is on the order of 10 Å, which corresponds to 10 adsorption sites/ cm2. Using temperature-programmed desorption, the formation of molecular oxygen and ozone from atomic oxygen was studied on amorphous solid water and crystalline water. From monitoring the relative amounts of the two molecular species n(O3)/n(O2) as a function of temperature (between 6 and 25 K), the oxygen diffusion rate as a function of temperature was obtained. These data were analyzed by assuming an Arrhenius law with an activation energy of 450 K and a prefactor of 10 s−1, which was unable to describe the temperature dependence (see blue line in Figure 1). Alternatively, a model that accounts for tunneling through a square barrier was used, which, however, also did not yield a satisfactory fit using a typical activation energy of 500 K and a barrier width of 1 Å. However, it is to be noted that such models assume homogeneous, noncorrugated (i.e., “smooth”) surfaces on which the diffusing particle moves. Once the surface contains variations on atomic length scales (“roughness” of scale ε, see Figure 1), as is typical for ASW, the temperature dependence of the diffusion rate may change. Diffusion in a rough, 1D potential was investigated by Zwanzig, who found that the motion of a particle on a rough potential energy surface (PES) is characterized by an effective diffusion coefficient Ddif * . 9 If the roughness of the surface is random, independent of position x, and Gaussian-distributed in ε (the scale of the roughness), then the effective diffusion coefficient is Ddif * = Ddife −(βε)2, where Ddif is the diffusion coefficient on the surface with zero roughness. On the contrary, for a periodic perturbation ∝ ε cos(qx) of the roughness with periodicity q the effective diffusion coefficient at low temperatures is Ddif * = Ddife −(2βε), that is, an Arrhenius dependence. Received: January 31, 2018 Accepted: March 25, 2018 Published: March 26, 2018 Letter pubs.acs.org/JPCL Cite This: J. Phys. Chem. Lett. 2018, 9, 1822−1826 © 2018 American Chemical Society 1822 DOI: 10.1021/acs.jpclett.8b00328 J. Phys. Chem. Lett. 2018, 9, 1822−1826 Because an Arrhenius expression does not manifestly describe the diffusion at low temperatures correctly, we used an expression D* = D0 + De −(βε) in which D0, D, and ε are adjustable parameters to fit the experimental data from ref 4 Figure 1 demonstrates that such an expression faithfully describes the experiments down to 6 K, whereby the parameter values are D0 = 3.75 s −1, D = 8.83 × 10 s−1, and ε = 93.4 K. For comparison, fits with D0 = 0 and the Arrhenius expression used in ref 4 are also shown. Diffusivities larger than 0 at very low temperatures (≤1 K) have, for example, been observed for atom diffusion in body-centered tungsten. Finite diffusivities at such low temperatures (<1 K) have been primarily associated with zero-point motion. To further corroborate the finding that no major corrections due to tunneling are required, atomistic simulations together with quantum-dynamical calculations for the translational motion of the oxygen atom were carried out. Molecular dynamics (MD) and Monte Carlo (MC) simulations were performed using CHARMM. The initial ASW structure was generated as previously described. Starting from a TIP3P water box (31 × 31 × 50 Å), equilibrated at 300 K, the system was first quenched to 50 K, then equilibrated in the NpT ensemble, followed by further equilibration with NVT using periodic boundary conditions. A cutoff of 13 Å was applied to the nonbonded interactions. All production MD simulations were performed in the microcanonical (NVE) ensemble with a time step of 1 fs. For both simulation temperatures of 10 and 50 K, respectively, heating and equilibration runs were performed for 10 and 50 ps, respectively. At 10 K only a single production simulation was run for 900 ns due to the low mobility of the oxygen atom on the surface. Increasing the temperature to 50 K lead to much wider sampling, and 10 independent simulations, each 900 ns long, were run. Bonds involving hydrogen atoms were constrained using SHAKE. For determining the nonbonded parameters for oxygen, 34 clusters (from an equilibrium MD simulation at 50 K) consisting of 12 water molecules interacting with one oxygen atom were extracted, and total energies at the MP2/aug-ccpVDZ level of theory using Gaussian09 were determined. Next, the oxygen atom van der Waals parameters rmin,O and εO were fitted using a simplex algorithm to best reproduce the target energies, which yielded rmin,O = 1.62 Å and εO = −0.518 kcal/mol. This parametrization gives oxygen desorption energies of 2.15 kcal/mol (1082 K), in good agreement with values of 1380 K derived from temperature-programmed desorption experiments. Using this parametrization, the 2D free-energy surface G(x, y) was determined from Monte Carlo (MC) simulations at 50 K on an equally spaced grid (31 × 31 Å). For each of the 961 grid points, 10 MC steps were run to determine the probability distribution P(x, y) from which G = −RT ln(P) was determined. The interaction of the oxygen atom with the water surface is sufficiently strong to preclude exhaustive sampling of the entire G(x, y) (only 60% of the grid is sampled). Therefore, the same MC simulations were repeated with the original (“orig”) CHARMM36 parametrization for the oxygen atom. For this force field the oxygen desorption energy is reduced to ∼1 kcal/mol, which allows complete and exhaustive sampling of Gorig(x, y) due to the smaller desorption energy, which also lowers the intersite barrier. The two free-energy surfaces can be readily related for the grid points that are sufficiently sampled (60% of all grid points) with a correlation coefficient of 0.925, and hence the less sampled regions using the correct, present parametrization for the oxygen atom can be extrapolated from Gorig(x, y), as shown in Figure S2. In addition, explicit MD simulations were carried out at 10 and 50 K using the refined nonbonded parameters for the oxygen atom. For 10 K the oxygen atom diffuses between neighboring minima for 900 ns; see Figure 2. Increasing the temperature to 50 K, the oxygen atom samples an extended path (length 62.4 Å (minima 1 to 8 and 17 to 26; see Figure 2A)) within 900 ns of explicit MD simulations. Hence, in translucent and diffuse clouds oxygen diffusion is expected to occur. The average value of the free-energy profile is indicated by the dashed red line, whereas the average barrier height (0.55 kcal/mol equivalent to 275 K) is the dashed green line above the average energy. As can be seen, many barriers along the path are below the average barrier height. At 50 K the average dwell time in a local minimum is ∼5 ns, which suggests facile diffusion at such temperatures on comparatively short time scales. To explicitly address the question of whether quantum effects are relevant for oxygen diffusion, the 1D Schrödinger equation for the translational motion of the oxygen atom was solved for part of the free-energy profile from Figure 2 (minima 23, 24, and 25; see Figure 3). For comparison, classical MD simulations solving Newton’s equations of motion were also carried out using the same free-energy profile. Parabolic walls were added at the outer boundaries to confine the motion of the oxygen atom. For following the quantum-mechanical time evolution of the system, the time-dependent Schrödinger equation ħ ∂ ∂ Ψ = ̂ Ψ i t x t H x t ( , ) ( , ) (1) was solved using the Crank−Nicholson (CN) method. Here x is the position, t is time, Ĥ is the Hamiltonian, and Ψ(x, t) is Figure 1. Different models for oxygen diffusion at low temperatures (6−25 K). Experimental points and error bars are those from ref 4. The blue curve is the Arrhenius expression using the parameters from ref 4 (Figure 2, D = Ae(−βε) with A = 10 s−1 and ε = 450 K), which fails to describe the experimental data. Fitting to an expression suggested in ref 9 (D* = De−(ε/T) 2 , D = 2.88 × 10 ± 2.61 × 10 s−1, ε = 78.6 ± 11.0 K, green line) describes the data down to 20 K. Allowing for finite diffusion down to very low temperatures (D* = D0 + De−(ε/T) 2 , D0 = 3.75 ± 0.70 s −1, D = 8.83 × 10 ± 4.82 × 10 s−1, ε = 93.4 ± 6.8 K, red curve) reproduces the experimental data to 6 K. The inset reports the y axis on a logarithmic scale. The Journal of Physical Chemistry Letters Letter DOI: 10.1021/acs.jpclett.8b00328 J. Phys. Chem. Lett. 2018, 9, 1822−1826 1823 the wave function. The CN method is a second-order, symplectic method that conserves total energy and provides high accuracy. Simulations used a time step of δt = 0.2 fs and a grid spacing of δx = 0.01a0. The original width of the groundstate wave function was 2.048 au, which is the width of the harmonic approximation to the potential around the minimum. For the classical molecular dynamics simulations the initial position of all particles was δ(x − x0), where x0 = 2.279 Å is the minimum of the 1D energy profile in region 1 (see Figure 3) and velocities were sampled from a Maxwell−Boltzmann distribution at 50 K. The classical equations for motion were propagated with a 1 fs time step using a velocity Verlet algorithm. Each dynamics was run for 5 ns and averaged over 10000 independent simulations. For solving the 1D timedependent Schrödinger equation using the CN method, 1000 simulations were initialized according to a Maxwell−Boltzmann distribution at 50 K. For each time frame the amplitudes were Boltzmann-averaged to obtain the total probability distribution. Both classical (panel C) and quantum simulations (panel D) reach a steady distribution after 10 (classical) and 20 ps (quantum), respectively, showing a slightly larger probability (10−3 vs 5 × 10−4 for region 3 and 10−2 vs 3 × 10−2 for region 2) if quantum effects are included. At the end (5 ns) of both simulations the probability remains highly localized in region 1, Figure 2. Typical oxygen diffusion path. (A) 2D free-energy surface G(x, y). All of the minima of the free-energy surface are labeled from 1 to 32. They form a continuous closed loop of length 120 Å. The typical region visited during one single 500 ns MD simulation at 50 K is indicated by a dashed, arrowed line and includes minima 1 to 8 and (due to periodic boundaries) 17 to 26. Even for a 900 ns simulation at 10 K (gray dots) several minima (22 to 23) are sampled. (B) 1D projection of the closed loop from panel A with the path starting and ending at minimum number 1. The average free energy is indicated by the red dashed line, and the average barrier height (green dashed line) is 275 K (0.55 kcal/mol). The Journal of Physical Chemistry Letters Letter DOI: 10.1021/acs.jpclett.8b00328 J. Phys. Chem. Lett. 2018, 9, 1822−1826 1824 with just a small portion able to overcome the barriers. This reflects the fact that the system is unable to reach thermal equilibrium at 50 K for both the potential and the free energy because no energy exchange with the environment can take place and only particles with sufficient initial velocity eventually overcome the barrier between the global minimum and the neighboring local minima. Hence the classical and quantum simulations are consistent with one another and suggest that tunneling plays a minor role. Remaining deviations from purely classical behavior, for example, a value of D0 ≠ 0 in the expression for the diffusion rate D* = D0 + De −(ε/T)2, are likely to be related to zero-point motion. Oxygen diffusion on a rough potential energy surface as provided by ASW allows for facile diffusion of the adsorbate down to temperatures relevant to interstellar chemistry. The parametrization for oxygen interacting with ASW used in this work correctly describes the average desorption energy of oxygen from water. Fitting the temperature dependence of the experimentally measured diffusion rates to an expression that accounts for roughness of the underlying PES yields a typical scale of the roughness of ε ≈ 100 K, which compares with typical average barrier heights of Edif = 275 K. Direct comparison of these two quantities should be done with care as they originate from two quite different analyses: ε is an effective roughness from a fit of experimental data sampling a large number of transitions and including the actual dynamics, whereas Edif is an average barrier height experienced by the diffusing oxygen from a few tens of transitions (see Figure 2B). Furthermore, the morphologies of the ASW surfaces from experiments and in the simulations are likely to differ. Nevertheless, both analyses clearly establish the existence of surface roughness, and it is encouraging that the magnitudes of the two effective parameters, ε and Edif, characterizing the surface corrugation are within less than a factor of 3 given the differences in the systems investigated. Explicit quantum and classical simulations for the translational motion of the oxygen atom on the free-energy surface demonstrate that quantum and classical simulations on the 5 ns time scale agree well and tunneling does not play a role. This is particularly true when considering that rough surfaces tend to suppress tunneling effects. The ratio of diffusional versus desorption barrier Edif/Edes = 275/1082 ≈ 0.3 found in the present work (from sampling an exemplary circular path) is somewhat lower than that recently proposed from temperature desorption experiments (0.55) but similar to other typically used ratios of 0.3 to 0.4. This ratio may depend on the surface morphology and change if longer simulations with more extensive sampling are carried out. Furthermore, the diffusion barrier height Edif = 275 K is sufficiently low to yield comparatively rapid mobility at low temperatures. Such a barrier corresponds to ∼10 hops per second (based on an Arrhenius dependence with a frequency factor of 10 s−1) and is consistent with values for D ≈ 10 to 10 s−1 from Figure 1, which suggests that oxygen atoms under such conditions are able to sample appreciable parts of the available surface to eventually find reaction partners. Conversely, the desorption energy Edes = 1000 K is sufficiently high to prevent escape from the surface and leads to long residence times. In conclusion, the analysis of experimental and simulation results shows that oxygen diffusion on ASW down to very low temperatures is possible and does not require quantum tunneling. This insight is of fundamental relevance for the formation of molecules including CO2, water (minor pathway), and other molecules (e.g., OCOH from HCO or OCCHO from HCCO) in laboratory experiments and on grains under conditions characteristic of the interstellar medium. The extension of such considerations to surfaces including graphite and silicon may provide further insight into chemical processing of grains in later stages of their development. ■ ASSOCIATED CONTENT *S Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.8b00328. Figure S1. Explicit diffusion path at 50 K. Figure S2. Additional information about the construction of the 2D free-energy surface for the parametrized model for oxygen interacting with water clusters. Figure S3. Data from Figure 3B on a linear y scale. (PDF) ■ AUTHOR INFORMATION Corresponding Author *E-mail: m.meuwly@unibas.ch.

T he diffusion of oxygen atoms in and on interstellar ice grains is of fundamental importance for the formation of CO 2 and molecular oxygen O 2 . 1−3 In addition, O diffusion has also been implicated in the formation of ozone in laboratory experiments 4,5 within a Langmuir−Hinshelwood mechanism (surface diffusion) and has been proposed to be relevant for interstellar grains, 1 eventually leading to the formation of water, 5 although the major water-formation pathway involves mobile hydrogen atoms.
In all of these cases the diffusion of atomic oxygen is a primary driver for forming di-or triatomic molecules under conditions typical for interstellar environments. These include higher temperatures (T ≈ 50 K) for translucent and diffuse clouds and temperatures T ≈ 10 K for dense molecular clouds. 1,6 Questions of particular relevance pertain to the diffusional barrier for intersite migration of oxygen atoms on amorphous solid water (ASW) and their desorption energies away from the surface. For this, experiments based on temperature-programmed desorption have been used. 4,7 For analyzing the diffusion rate (i.e., the probability for one hop between different adsorption sites in s −1 , which is referred to as D in ref 4) of atomic oxygen on ASW, either an Arrhenius expression (D(T) ≈ Ae (−βE dif ) , where E dif is a conformationally averaged diffusional barrier and β = 1/(k B T)) or an expression accounting for tunneling through a square well barrier was used. 4,7 A diffusion coefficient D dif (in units of length 2 /time) can be converted into a diffusion rate D (in 1/time) through D = D dif /a 2 , where a 2 is a typical area sampled by the adsorbate between consecutive hops (see Figure S1). Typically, a 2 is on the order of 10 Å 2 , which corresponds to 10 15 adsorption sites/ cm 2 . 8 Using temperature-programmed desorption, the formation of molecular oxygen and ozone from atomic oxygen was studied on amorphous solid water and crystalline water. 4 From monitoring the relative amounts of the two molecular species n(O 3 )/n(O 2 ) as a function of temperature (between 6 and 25 K), the oxygen diffusion rate as a function of temperature was obtained. These data were analyzed by assuming an Arrhenius law with an activation energy of 450 K and a prefactor of 10 12 s −1 , which was unable to describe the temperature dependence (see blue line in Figure 1). 4,7 Alternatively, a model that accounts for tunneling through a square barrier was used, which, however, also did not yield a satisfactory fit using a typical activation energy of 500 K and a barrier width of 1 Å. 4 However, it is to be noted that such models assume homogeneous, noncorrugated (i.e., "smooth") surfaces on which the diffusing particle moves. Once the surface contains variations on atomic length scales ("roughness" of scale ϵ, see Figure 1), as is typical for ASW, the temperature dependence of the diffusion rate may change. Diffusion in a rough, 1D potential was investigated by Zwanzig, who found that the motion of a particle on a rough potential energy surface (PES) is characterized by an effective diffusion coefficient D dif * . 9 If the roughness of the surface is random, independent of position x, and Gaussian-distributed in ϵ (the scale of the roughness), then the effective diffusion coefficient is D dif * = D dif e −(βϵ) 2 , where D dif is the diffusion coefficient on the surface with zero roughness. On the contrary, for a periodic perturbation ∝ ϵ cos(qx) of the roughness with periodicity q the effective diffusion coefficient at low temperatures is D dif * = D dif e −(2βϵ) , that is, an Arrhenius dependence.
Because an Arrhenius expression does not manifestly describe the diffusion at low temperatures correctly, we used an expression D* = D 0 + De −(βϵ) 2 in which D 0 , D, and ϵ are adjustable parameters to fit the experimental data from ref 4 Figure 1 demonstrates that such an expression faithfully describes the experiments down to 6 K, whereby the parameter values are D 0 = 3.75 s −1 , D = 8.83 × 10 4 s −1 , and ϵ = 93.4 K. For comparison, fits with D 0 = 0 and the Arrhenius expression used in ref 4 are also shown. Diffusivities larger than 0 at very low temperatures (≤1 K) have, for example, been observed for atom diffusion in body-centered tungsten. 10,11 Finite diffusivities at such low temperatures (<1 K) have been primarily associated with zero-point motion. 11 To further corroborate the finding that no major corrections due to tunneling are required, atomistic simulations together with quantum-dynamical calculations for the translational motion of the oxygen atom were carried out. Molecular dynamics (MD) and Monte Carlo (MC) simulations were performed using CHARMM. 12 The initial ASW structure was generated as previously described. 13 Starting from a TIP3P 14 water box (31 × 31 × 50 Å 3 ), equilibrated at 300 K, the system was first quenched to 50 K, then equilibrated in the NpT ensemble, followed by further equilibration with NVT using periodic boundary conditions. A cutoff of 13 Å was applied to the nonbonded interactions. All production MD simulations were performed in the microcanonical (NVE) ensemble with a time step of 1 fs. For both simulation temperatures of 10 and 50 K, respectively, heating and equilibration runs were performed for 10 and 50 ps, respectively. At 10 K only a single production simulation was run for 900 ns due to the low mobility of the oxygen atom on the surface. Increasing the temperature to 50 K lead to much wider sampling, and 10 independent simulations, each 900 ns long, were run. Bonds involving hydrogen atoms were constrained using SHAKE. 15 For determining the nonbonded parameters for oxygen, 34 clusters (from an equilibrium MD simulation at 50 K) consisting of 12 water molecules interacting with one oxygen atom were extracted, and total energies at the MP2/aug-cc-pVDZ level of theory using Gaussian09 16 were determined. Next, the oxygen atom van der Waals parameters r min,O and ϵ O were fitted using a simplex algorithm to best reproduce the target energies, which yielded r min,O = 1.62 Å and ϵ O = −0.518 kcal/mol. This parametrization gives oxygen desorption energies of 2.15 kcal/mol (1082 K), in good agreement with values of 1380 K derived from temperature-programmed desorption experiments. 7 Using this parametrization, the 2D free-energy surface G(x, y) was determined from Monte Carlo (MC) simulations at 50 K on an equally spaced grid (31 × 31 Å 2 ). For each of the 961 grid points, 10 8 MC steps were run to determine the probability distribution P(x, y) from which G = −RT ln(P) was determined. The interaction of the oxygen atom with the water surface is sufficiently strong to preclude exhaustive sampling of the entire G(x, y) (only 60% of the grid is sampled). Therefore, the same MC simulations were repeated with the original ("orig") CHARMM36 parametrization for the oxygen atom. For this force field the oxygen desorption energy is reduced to ∼1 kcal/mol, which allows complete and exhaustive sampling of G orig (x, y) due to the smaller desorption energy, which also lowers the intersite barrier. The two free-energy surfaces can be readily related for the grid points that are sufficiently sampled (60% of all grid points) with a correlation coefficient of 0.925, and hence the less sampled regions using the correct, present parametrization for the oxygen atom can be extrapolated from G orig (x, y), as shown in Figure S2.
In addition, explicit MD simulations were carried out at 10 and 50 K using the refined nonbonded parameters for the oxygen atom. For 10 K the oxygen atom diffuses between neighboring minima for 900 ns; see Figure 2. Increasing the temperature to 50 K, the oxygen atom samples an extended path (length 62.4 Å (minima 1 to 8 and 17 to 26; see Figure  2A)) within 900 ns of explicit MD simulations. Hence, in translucent and diffuse clouds oxygen diffusion is expected to occur. The average value of the free-energy profile is indicated by the dashed red line, whereas the average barrier height (0.55 kcal/mol equivalent to 275 K) is the dashed green line above the average energy. As can be seen, many barriers along the path are below the average barrier height. At 50 K the average dwell time in a local minimum is ∼5 ns, which suggests facile diffusion at such temperatures on comparatively short time scales.
To explicitly address the question of whether quantum effects are relevant for oxygen diffusion, the 1D Schrodinger equation for the translational motion of the oxygen atom was solved for part of the free-energy profile from Figure 2 (minima 23, 24, and 25; see Figure 3). For comparison, classical MD simulations solving Newton's equations of motion were also carried out using the same free-energy profile. Parabolic walls were added at the outer boundaries to confine the motion of the oxygen atom.
For following the quantum-mechanical time evolution of the system, the time-dependent Schrodinger equation was solved using the Crank−Nicholson (CN) method. 17 Here x is the position, t is time, Ĥis the Hamiltonian, and Ψ(x, t) is the wave function. The CN method is a second-order, symplectic method that conserves total energy and provides high accuracy. Simulations used a time step of δt = 0.2 fs and a grid spacing of δx = 0.01a 0 . The original width of the groundstate wave function was 2.048 au, which is the width of the harmonic approximation to the potential around the minimum.
For the classical molecular dynamics simulations the initial position of all particles was δ(x − x 0 ), where x 0 = 2.279 Å is the minimum of the 1D energy profile in region 1 (see Figure 3) and velocities were sampled from a Maxwell−Boltzmann distribution at 50 K. The classical equations for motion were propagated with a 1 fs time step using a velocity Verlet algorithm. 18 Each dynamics was run for 5 ns and averaged over 10000 independent simulations. For solving the 1D timedependent Schrodinger equation using the CN method, 1000 simulations were initialized according to a Maxwell−Boltzmann distribution at 50 K. For each time frame the amplitudes were Boltzmann-averaged to obtain the total probability distribution.
Both classical (panel C) and quantum simulations (panel D) reach a steady distribution after 10 (classical) and 20 ps (quantum), respectively, showing a slightly larger probability (10 −3 vs 5 × 10 −4 for region 3 and 10 −2 vs 3 × 10 −2 for region 2) if quantum effects are included. At the end (5 ns) of both simulations the probability remains highly localized in region 1, with just a small portion able to overcome the barriers. This reflects the fact that the system is unable to reach thermal equilibrium at 50 K for both the potential and the free energy because no energy exchange with the environment can take place and only particles with sufficient initial velocity eventually overcome the barrier between the global minimum and the neighboring local minima. Hence the classical and quantum simulations are consistent with one another and suggest that tunneling plays a minor role. Remaining deviations from purely classical behavior, for example, a value of D 0 ≠ 0 in the expression for the diffusion rate D* = D 0 + De −(ϵ/T) 2 , are likely to be related to zero-point motion.
Oxygen diffusion on a rough potential energy surface as provided by ASW allows for facile diffusion of the adsorbate down to temperatures relevant to interstellar chemistry. The parametrization for oxygen interacting with ASW used in this work correctly describes the average desorption energy of oxygen from water. Fitting the temperature dependence of the experimentally measured diffusion rates 4 to an expression that accounts for roughness of the underlying PES yields a typical scale of the roughness of ϵ ≈ 100 K, which compares with typical average barrier heights of E dif = 275 K. Direct comparison of these two quantities should be done with care as they originate from two quite different analyses: ϵ is an effective roughness from a fit of experimental data sampling a large number of transitions and including the actual dynamics, whereas E dif is an average barrier height experienced by the diffusing oxygen from a few tens of transitions (see Figure 2B). Furthermore, the morphologies of the ASW surfaces from experiments and in the simulations are likely to differ. Nevertheless, both analyses clearly establish the existence of surface roughness, and it is encouraging that the magnitudes of the two effective parameters, ϵ and E dif , characterizing the surface corrugation are within less than a factor of 3 given the differences in the systems investigated. Explicit quantum and classical simulations for the translational motion of the oxygen atom on the free-energy surface demonstrate that quantum and classical simulations on the 5 ns time scale agree well and tunneling does not play a role. This is particularly true when considering that rough surfaces tend to suppress tunneling effects.
The ratio of diffusional versus desorption barrier E dif /E des = 275/1082 ≈ 0.3 found in the present work (from sampling an exemplary circular path) is somewhat lower than that recently proposed from temperature desorption experiments (0.55) 7 but similar to other typically used ratios of 0.3 to 0.4. 19 This ratio may depend on the surface morphology and change if longer simulations with more extensive sampling are carried out. Furthermore, the diffusion barrier height E dif = 275 K is sufficiently low to yield comparatively rapid mobility at low temperatures. Such a barrier corresponds to ∼10 4 hops per second (based on an Arrhenius dependence with a frequency factor of 10 12 s −1 ) 7 and is consistent with values for D ≈ 10 4 to 10 5 s −1 from Figure 1, which suggests that oxygen atoms under such conditions are able to sample appreciable parts of the available surface to eventually find reaction partners. Conversely, the desorption energy E des = 1000 K is sufficiently high to prevent escape from the surface and leads to long residence times.
In conclusion, the analysis of experimental and simulation results shows that oxygen diffusion on ASW down to very low temperatures is possible and does not require quantum tunneling. This insight is of fundamental relevance for the formation of molecules including CO 2 , water (minor pathway), and other molecules (e.g., OCOH from HCO or OCCHO from HCCO) 3 in laboratory experiments and on grains under conditions characteristic of the interstellar medium. The extension of such considerations to surfaces including graphite and silicon may provide further insight into chemical processing of grains in later stages of their development.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.8b00328. Figure S1. Explicit diffusion path at 50 K. Figure S2. Additional information about the construction of the 2D free-energy surface for the parametrized model for oxygen interacting with water clusters. Figure S3. Data from Figure 3B