Molecular dynamics simulations of energy dissipation and non-thermal diffusion on amorphous solid water

A. Fredon and H. M. Cuppen *
Radboud University Nijmegen, Institute for Molecules and Materials, Heyendaalseweg 135, 6525 AJ, Nijmegen, The Netherlands. E-mail:

Received 8th September 2017 , Accepted 23rd January 2018

First published on 29th January 2018

Molecules in space are synthesized via a large variety of gas-phase reactions, and reactions on dust-grain surfaces, where the surface acts as a catalyst. Especially, saturated, hydrogen-rich molecules are formed through surface chemistry where the interstellar grains act as a meeting place and absorbing energy. Here we present the results of thousands of molecular dynamics simulations to quantify the outcome of an energy dissipation process. Admolecules on top of an amorphous solid water surface have been given translational energy between 0.5 and 5 eV. Three different surface species are considered, CO2, H2O and CH4, spanning a range in binding energies, number of internal degrees of freedom and molecular weight. The results are compared against a previous study using a crystalline water ice surface. Possible outcomes of a dissipation process are adsorption – possibly after long-range diffusion-, desorption and desorption of a surface molecule. The three admolecules were found to bind at different locations on the surface, particularly in terms of height. Water preferably binds on top of the surface, whereas methane fills the nanopores on the surface. This has direct consequences for desorption, travelled distance, and kick-out probabilities. The admolecules are found to frequently travel several tens of angstroms before stabilizing on a binding site, allowing follow-up reactions en route. We present kick-out probabilities and we have been able to quantify the desorption probability which depends on the binding energy of the species, the translational excitation, and a factor that accounts for difference in binding site height. We provide expressions that can be incorporated in astrochemical models to predict grain surface formation and return into the gas phase of these products.

1 Introduction

Interstellar dust grains in dark molecular clouds are typically covered by ices. These are predominantly made of water in the form of amorphous solid water (ASW) but may also contain other chemical species such as CO or CO2. Many saturated molecules like acetaldehyde and dimethylether have been observed in these environments. Most reactions leading to saturated products need a third body to release their exothermicity and the low number density prevents these from occurring in the gas phase. Surface reactions on the ice mantle of interstellar dust grains play hence a key role in the production of these chemical species. Surface reactions proceed typically through the diffusive Langmuir–Hinshelwood reaction. The limiting step is the diffusion of reactant species adsorbed on the grain. Because of the physical conditions in the ISM, diffusion of species heavier than H and H2 is likely to be prohibitively slow.

In a previous paper,1 we introduced the concept of “non-thermal diffusion” where the energy released during a surface reaction is applied to allow the reactant to travel across the surface. This behaviour was demonstrated by preforming a series of Molecular Dynamics (MD) simulations where an admolecule was given an initial amount of translational energy (between 0.5 and 5 eV). This type of excitation was used as a proxy to study the fate of excited molecules that are formed through surface reaction. The simulations showed that the travelled distance over a crystalline ice surface can be as large as 100 nm when 5 eV is given to a molecule in the form of translational kinetic energy. This concept appeared to be generally applicable, since we found similar radial travelled distances for CO2, H2O, and CH4, the three admolecules studied. The energy can also be applied to desorb the species, i.e., leading to chemical or reactive desorption.2–4 This type of desorption became particularly important after recent astronomical detection of organic molecules, like acetaldehyde and dimethylether, in cold dark molecular clouds.5 These species are thought to have formed through grain surface chemistry, but are observed in the gas phase.6 Without thermal or photon energy present, other desorption routes need to be considered.

The present paper aims to present a similar approach for amorphous ice and to compare the results with the crystalline ice results. Interstellar ices are predominately amorphous and hence such a surface will be a better representative for interstellar surfaces. Because of the corrugated nature of amorphous ice we expect the travelled distance to be shorter due to the presence of stronger binding sites which trap the molecules more easily. We also expect the desorption probability to be higher for CH4 because collisions of the CH4 molecule with dangling protons lead to desorption of the molecule. In contrast, H2O prefers to interact with the dangling protons of the ice hence we expect an opposite trend for the desorption probability. We do not expect the desorption probability to substantially alter for CO2.

The chemical evolution of a molecular cloud is often modelled by rate equations including hundreds of species, both in the gas phase and on the grain surface, connected through thousands of reactions. General analytical expressions, often empirically obtained, are used to describe processes such as diffusion, desorption and reaction.7 Here, the aim is to arrive at an expression to describe the fate of the molecules after reaction or photodissociation. It is less straightforward to simulate the actual reaction using classical MD simulations and hence assumptions on the energies involved need to be made. This is however also an advantage since the amount of translational energy can be freely varied to study trends, rather than being limited to a single reaction and a single system. This is beneficial for our aim to arrive at general expressions. Experimentally, this is harder to achieve, since upon changing the chemical species to investigate the role of, for instance, binding energy, the energetics will change as well. In the current study, three different surface species are considered, CO2, H2O and CH4, which are among the most common ice species.8 Moreover, these species span a range in binding energies, number of internal degrees of freedom and molecular weight.

We expect that in general, part of the excess energy will be converted to kinetic energy. The remainder will be rotational, or vibrational energy or electronic excitation. How this energy is spread is highly system dependent;9 it does not only depend on the specific reaction but also on the reaction environment and the configuration in which it occurs. Here we limit ourselves to translational excitation. The chosen range corresponds to typical exothermicities of surface reactions or energy released upon photodissociation. In the latter process, VUV photons electronically excite molecules that then photodissociate and its photofragments fall back into the ground state with excess energies of several eV. At these high energies, we do not expect quantum effects such as zero point energy and tunnelling to dominate. Bäck and Marković10 found for collision simulations resembling a water cluster interacting with graphite, vibrational dissipation to be similar for classical and quantum, wave-packet simulations. Similarly, Glowacki et al.11 used MD simulations to study energy dissipation at surface hotspots under typical CVD conditions.

2 Methodology

The protocol used is similar to the one we used in our previous study.1 We mimic the process in which a reaction product is formed on an amorphous ice surface with an excess energy by giving a specific amount of translational energy to the admolecule in a random direction and following the outcomes. First, the initial ice surface is probed by identifying the different binding sites. Next, the ten strongest binding sites on the surface are selected and from each of them new simulations are issued where the admolecule is given a specific amount of additional translational energy in a random direction.

2.1 Interaction potentials

To allow for the transfer of energy to molecular vibrations, fully flexible interaction potentials are used. For H2O, we use the flexible TIP4P/2005f potential.12 The intramolecular interactions for CO2 and CH4 have been taken from ref. 13 and 14, respectively. The intermolecular H2O–CO2 interaction is taken from Karssemeijer et al.15 For the H2O–CH4 interactions, we fitted a site-site potential from the H2O–CH4 PES2b-CSM potential as obtained by Qu et al.16 All potentials and parameters are described in the section “Interaction potentials” of the previous paper. We observed during our first study that 8.7% of all CH4 trajectories ended up in the “FAIL” category. The “FAIL” simulations are the simulations in which LAMMPS stopped to run for various reasons, mainly due to a too large time step and/or too attractive short range potential. Visual inspection of trajectories at high kinetic energy showed that the carbon atom of CH4 closely approached the atoms of the water substrate. For an amorphous substrate this percentage is much higher. The H2O–CH4 potential uses an undamped Buckingham and this particular fit resulted in a rather broad “hole” at close distance, especially for the OH2O–CCH4 and HH2O–CCH4 interactions. We have therefore exchanged the Buckingham for a Beck potential, which has a damping term:
image file: c7cp06136f-t1.tif(1)
The parameters we used for the Beck potential are listed in Table 1. The newly fitted Beck potential closely follows the original Buckingham potential, except at very close distances.
Table 1 Parameters of the H2O–CH4 potential
A (eV) B (eV Å6) α−1) β−6) a (Å)
OH2O–CCH4 741.538 58.567 2.988 0.0 1.856
HH2O–CCH4 58.373 2.815 3.538 0.0 1.1

For all interactions, the dispersion contribution is cut-off at 11.0 Å and the Coulombic part of the potential is computed using a particle–particle particle–mesh solver17 with a relative RMS error in per-atom force of 10−6.

2.2 Simulation procedure

All simulations are performed using the LAMMPS package (version 16/02/16).18 An amorphous sample of 360 H2O molecules is used. This is obtained by heating an Ih sample from 10 K to 400 K with a ramp of 78 K ps−1 with fixed volume and subsequently cooling it to 200 K with a ramp of 2 K ps−1. It is then quenched to 10 K and equilibrated during 50 ps at this temperature. The resulting structure is further allowed to relax during a NPT simulation using a Nosé–Hoover thermostat and barostat at a temperature of 10 K with a thermostat parameter of 50 fs and at a pressure of 1 bar with a barostat parameter of 500 fs. The timestep of the simulation is 0.5 fs and after 20 ps the box lengths are equilibrated with the following values: x = 22.239176 Å, y = 22.137243 Å, z = 22.142457 Å. The density of the amorphous ice is 0.988 g cm−3, which is slightly higher than the experimental value for vapour-deposited high-density amorphous ice Iah of 0.94 g cm−3.19 Although the pressure in molecular clouds is on order of 10 cm−3 a pressure of 1 bar was applied to make the integration feasible.

To create a surface along the (x,y)-plane of the simulation box, the box length in the z direction is increased by 100 Å, creating a slab model. This 100 Å prevents the surface and any desorbing molecule to interact through the periodic boundary in the z direction. The surface is further relaxed for 20 ps at 10 K in the NVT ensemble.

To trace the variations in binding sites, 1000 NVT equilibration simulations are performed for each of the three admolecules. For each run, the initial position of the admolecule is randomly chosen in x and y direction and Δz = 7 Å above the surface. The orientation of the admolecule with respect to the surface is also chosen randomly. To promote binding to the surface, the molecule is given a low initial velocity towards the surface. The initial velocities of the substrate H2O molecules were taken from the last step of the surface preparation simulation. Again, the Nosé–Hoover thermostat was used at a temperature of 10 K.

We then select 10 final configurations with the corresponding velocities as starting point for the energy dissipation simulations. The non-bonded interaction energy between the admolecule and the surface is used as a selection criterion where we avoid trajectories that ended in the same binding site (similar (x,y) coordinates). To avoid energy being artificially extracted from the system by the simulation thermostat, we switch the thermostat off and simulate in the NVE ensemble. Thus the total energy of the system remains constant, but the energy can be released from the admolecule to the substrate. The admolecule is given an additional translational kinetic energy between 0.5 and 5 eV in steps of 500 meV. To obtain sufficient statistics for the amorphous, inhomogeneous sample, 2000 simulations per species per initial kinetic energy are performed; as a comparison 1400 simulations were used on the crystalline surface. Two hundreds simulations are performed for each of the ten energy bins and each of the ten binding sites, leading to a total of 10 × 10 × 200 = 20[thin space (1/6-em)]000 energy dissipation simulations per admolecule. The direction of the additional velocity is chosen randomly and is uniformly distributed. The simulations run with a maximum of 20 ps, but can be terminated earlier if they meet one of the following criteria:

Desorption. The admolecule is outside the cutoff range of the potential, and hence no longer feels the surface.
Penetration into the surface. The admolecule is below a certain threshold for at least 4 ps.
Adsorption. The admolecule has a total energy smaller than −0.1 eV (CO2), −0.4 eV (H2O) or 0.0 eV (CH4) for at least 4 ps. These thresholds are the upper limits of the total energy distribution.

In summary, we perform (1) an NPT simulation for the ice bulk, (2) followed by an NVT simulation to obtain an ice surface and (3) an NVT simulation for the ice + admolecule followed by an NVE simulation where the admolecule has been given additional translational energy.

3 Results

3.1 Binding energy of the admolecule to the amorphous solid water surface

Fig. 1 shows a top view of the ASW substrate with the positions of the oxygen and hydrogen atoms indicated by the dots. The colour coding represents the relative height of the atoms which spans over 7 Å. The stars in this figure represent the centres of mass (COMs) of the ten selected binding sites for the three admolecules. These, hence, represent the strongest binding sites. What is immediately apparent is that the (x,y) location of the strongest binding sites on substrate is highly species dependent. We will return to this point later. But also the z coordinate is different for the different species. Water admolecules seem to favour the binding sites with the high z values, i.e., on top of the substrate. We believe that this is because the deposition of H2O molecules is mostly ballistic. Approaching H2O molecules tend to form hydrogen bonds as soon as there is another water molecule in their vicinity. Because of their low incoming velocity, they will stick here, which is typically at high z values. Carbon dioxide has a large range of different heights; some are very close in terms of (x,y) coordinates. The heights depend on whether the molecule is parallel to the surface or perpendicular. In the latter case, the C atom will be higher above the surface. Finally, the methane sites are rather deep into the substrate. We think that this is because of its rather low binding energy, which allows the admolecule to scan the surface and reach a pore. In the ESI, there is movie available which shows the adsorption of a CH4 molecule followed by a kick. The first part of the movie clearly shows that the methane molecule initially finds a higher lying binding site and then slowly moves further to find a deep binding site.
image file: c7cp06136f-f1.tif
Fig. 1 Top view of H2O substrate and the binding sites for CO2, H2O, and CH4, from top to bottom, after the equilibration simulations. The atoms of substrate molecules represented by circles, connected by lines. The stars represent the COMs of the admolecules in the 10 strongest binding sites existing on the surface. The colour coding indicates the height of the atoms/COMs.

Fig. 2 shows the interaction energies for the different binding sites. The black dots represent the oxygen atom of the substrate; the coloured dots represent now the ending positions of the COMs of the different admolecules after energy dissipation in the NVE ensemble. The colour coding indicates the interaction energy. The interaction energy is here defined as the sum of the Coulombic and van der Waals interactions between the admolecule and the surface. For all three admolecules this spans a large range: between −0.6 and −0.2 eV, −1.2 and −0.4 eV, and −0.4 and 0.0 eV, for CO2, H2O, and CH4, respectively. Again the different (x,y) locations of the strongest binding sites for the three admolecules is clearly visible. The stars indicate the ten selected binding sites.

image file: c7cp06136f-f2.tif
Fig. 2 Top view of H2O substrate and the binding sites for CO2, H2O, and CH4, from top to bottom, after the energy dissipation simulations. The COMs of substrate molecules are represented by the black dots. The stars are the ten selected sites for the follow-up simulations (cf.Fig. 1) and the coloured dots represent the COMs of the admolecules at the end of the energy dissipation simulations. The colour coding indicates the interaction energy of the admolecules with the substrate.

Experimentally, the binding energies of CO2, H2O, and CH4 on amorphous solid water are −0.195 eV,20 −0.414 eV,21 and −0.118 eV,22 respectively, as determined by temperature programmed desorption (TPD) experiments. This is, in all cases, on the lower end of our energy range. One should realise that our simulation results are for very low coverage, where our admolecule can always take up the strongest binding site. In the experiments the coverage is significantly higher and molecules are forced to occupy less favourable binding sites. Moreover, TPD experiments are more sensitive to the more weakly bound species.

3.2 Outcome classification

Fig. 3 and 4 show the distribution of the trajectories over the different categories as a function of the initial translational energy given to the three admolecules. Contrary to the simulations on a crystalline surface, there are several trajectories that lead to a “kick-out” of a water substrate molecule. We have hence added three more categories with respect to our previous work: desorption/adsorption/penetration with desorption of a H2O substrate molecule. The distributions over the different outcome categories are presented in two histograms for legibility: Fig. 3 shows the dominant outcomes (desorption, adsorption, and penetration) with an additional “others” category which hold the less frequently occurring outcomes. These are presented in Fig. 4 (the three kick-out categories, time limit, and FAIL).
image file: c7cp06136f-f3.tif
Fig. 3 Probabilities of occurrence for the main categories as a function of the initial translational kinetic energy given to the admolecules.

image file: c7cp06136f-f4.tif
Fig. 4 Probabilities of occurrence for the minority categories as a function of the initial translational kinetic energy given to the admolecules.

Fig. 3 clearly shows, like for the crystalline surface, that the desorption probability increases with increasing kick energies. As expected the stronger the kicking energy, the higher the probability to desorb from the surface. However it appears that the desorption probability reaches an asymptotic value after a certain kicking energy. This is especially true for CH4. Methane desorbs from crystalline ice with nearly 100% probability with a kicking energy higher or equal to 2 eV; in the case of ASW, the CH4 admolecule hardly reaches 53% of desorption with a kicking energy of 5 eV.

Fig. 4 shows that the “FAIL” outcome does not occur any longer in the case of CH4 meaning that the modification of the pair potential is successful. In the case of H2O the “FAIL” outcome accounts for a substantial fraction of the trajectories at high energy (7.9% for a kicking energy of 5 eV). Inspection of some of these trajectories shows that this occurs when the translational energy is given along an hydrogen bond and the oxygen of one molecule and the hydrogen of another one come too close. Unfortunately, this is intrinsic to the TIP4P potentials, since a smaller timestep cannot directly solve this problem. We restarted the 361 “FAIL” simulations with a timestep of 0.1 fs (first set) and 0.01 fs (second set). For the first set, only one simulation finished normally and for the 2nd set, ten simulations.

3.3 Desorption

The present section describes trajectories leading to desorption of the admolecule. This includes the “desorption” and the “desorption + H2O desorption” outcomes. Fig. 5 gives the desorption probability as a function of the initial translational energy. It is compared against crystalline surface results. For all species the desorption probability is reduced on an ASW surface as compared to a crystalline surface. The effect is the largest for CH4, with nearly a factor of 2 reduction, and fairly small for H2O. This is different from what we initially expected. For CH4 on a crystalline surface, the desorption was mainly triggered by CH4 colliding with a dangling proton on the surface. We had, therefore, expected that on a substrate with more protrusions, like ASW, the desorption would be even more frequently occurring. The desorption of H2O on the other hand was found to be hampered by dangling protons since they would allow the admolecule to form strong hydrogen bonds. We had expected this to become even more important on ASW, which would strongly reduce the desorption. This is not what the results show.
image file: c7cp06136f-f5.tif
Fig. 5 Desorption probability as a function of the initial translational energy given to the admolecule. The dots represent the data from our simulations with an error bar of 3σ. The dashed lines represent our models for desorption from crystalline ice (eqn (3)) in black and from ASW (eqn (4)) in blue.

We now believe that this is caused by the location of the binding sites for the three different admolecules on the surface. As shown in Section 3.1, CH4 preferably resides in lower lying pockets in the substrate, from which it has to escape. Only specific directions, those pointing upwards, result in desorption and hence the desorption is only possible for a fraction of the trajectories. H2O, on the other hand, preferably resides on top of substrate, after ballistic deposition. Its desorption behaviour hence resembles much more the crystalline case.

The desorption probability presented in Fig. 5 could be a measure of chemical desorption where some of the excess energy is released in the form of kinetic energy. As mentioned before this type of non-thermal desorption has gained in attention recently and is included in several astrochemical models. These models would benefit from a general law that can be applied to predict the efficiency of chemical desorption with the excess energy and the binding energy of the reaction product as input parameters. Here, we have access to the latter. We do, however, not know how the reaction enthalpy and the translational excitation energy are related. This is heavily dependent on the specific reaction and reaction environment.9 For simplicity, it is often assumed that the excitation energy spreads equally spread over all 3N degrees of freedom. We will also use this approximation

image file: c7cp06136f-t2.tif(2)
Here we ignore the Δ(PV) difference. For the crystalline substrate, we found the following empirical relation between this quantity and the desorption probability
image file: c7cp06136f-t3.tif(3)
with ω = 0.3 for all three admolecules. This value of ω can be reasoned by assuming that only the projection of the Etransreact on the z axis, both positive and negative values, contributes to the desorption. This projection could effectively be a third of the energy integrate over all directions. The initial kinetic energy is used as a proxy for Etransreact. The black dashed line in Fig. 5 represents this general expression. To compute the binding energies of the three species, the change in kinetic energy upon desorption was plotted against the initial kinetic energy for all “desorption” simulations with a kick angle θ between 0 and 20° with respect to the surface normal. In these simulations, the admolecule escapes the surface without a collision. The binding energy was then obtained by extrapolating with a linear fit to the energy required to leave the surface at zero kinetic energy. The binding energies thus obtained are 0.495, 0.839, and 0.160 eV for CO2, H2O, and CH4, respectively. These agree favourably with the interaction energies at the initial binding sites.

For the ASW results, we can observe a similar trend, but as mentioned earlier, the curves do not approach unity but some other, smaller, asymptotic value. The blue dashed curve plots a new relation

image file: c7cp06136f-t4.tif(4)
with ω = 0.3 for all three admolecules and f an admolecule-dependent parameter. We use f = 0.6, 0.8, and 0.5 for CO2, H2O, and CH4, respectively, to fit the results. For the binding energies we use −0.468 eV, −0.746 eV, and −0.280 eV for CO2, H2O, and CH4, respectively. The binding energies are obtained in the same way as described above. Both for f and |Ebind| we see the largest effect for CH4 with respect to the crystalline results. Again, we think this is due to environment of the binding sites.

As mentioned earlier, it is more difficult to arrive at these general expressions experimentally. Minissale et al.23 have measured the chemical desorption using Temperature Programmed Desorption (TPD) for several reactions using several substrates, mostly HOPG and silicates. TPD is not particularly sensitive for this and their values for non-porous ASW have either large uncertainties or are upper limits. For the reaction

H + OH → H2O.(5)
they found a desorption efficiency of 0.30 ± 0.15 on non-porous ASW. We will use this number to test the validity of our expression. The enthalpy of reaction is 4.89 eV. If we assume an equipartition of the energy over all 9 degrees of freedom (Etransreact = 3 × 1/9 × 4.89 = 1.63 eV), the desorption probability from ASW according to eqn (4) is
image file: c7cp06136f-t5.tif(6)
This is in good agreement with the experimental value, considering the large uncertainties.

3.4 Kick-out of a substrate molecule

A substantial fraction of trajectories results in the kick-out of a substrate molecule. This desorption mechanism was first reported in molecular dynamics simulations of photodissociation of water ice where the photofragments of a water molecule in the lower layers were given excess energy.24,25 It was later also experimentally shown to occur.26 We did not observe this in our crystalline ice study. The setup where molecules on top of the surface were given excess energy, most likely prevented this kick-out mechanism from occurring. In the current setup, this kick-out mechanism is observed. Fig. 6 shows the total fraction of trajectories resulting in desorption of a substrate H2O. These include desorption, adsorption and penetration of the admolecule.
image file: c7cp06136f-f6.tif
Fig. 6 Probability for a substrate molecule to desorb through the “kick-out” mechanism as a function of the initial translational energy given to the admolecule.

Methane can already kick out substrate molecules for a moderate kicking energy of ≥1.5 eV and it remains the admolecule which can most effectively kick out substrate molecules over the whole energy range. This is again due to the fact that the binding sites of CH4 are more encapsulated. The results by Andersson and Arasa et al.27,28 also showed that this mechanism is more effective when the excited species is subsurface. The movie in the ESI, shows a trajectory where a CH4 is given a kick downwards to right hand side of the frame. It kicks out a water molecule on the other side of the periodic box. The methane molecule leaves the surface highly rotationally excited. The H2O molecule is not rotationally excited.

Water as an admolecule is surprisingly efficient in kicking out substrate molecules, especially in comparison to CO2. Its binding sites are located more on top of the surface and it is much lighter. Here most likely collisions with protruding water molecules on the substrate result in kick out. These protruding molecules are less strongly bound.

3.5 Adsorption following non-thermal diffusion

The present section describes trajectories leading to adsorption of the admolecule. Fig. 7 indicates the travelled distance and direction of the admolecules for each “adsorption” trajectory. The figures shows multiple copies of the same simulation cell stacked in both x and y direction; the individual cells are indicated by the grid. All admolecules have started from the central simulation box and they moved to neighboring replicas (in some cases multiple times) to their final position indicated by the dots. The dots are colour coded according to their initial kicking energy. H2O can travel the largest distance with many trajectories resulting in a travelled distance around 75 Å. This is closely followed by CO2. The travelled distance by CH4 is substantially less. For all species, the travelled distance is significantly reduced (more than a factor of 2) with respect to the crystalline case. Again, the largest difference can be observed for CH4.
image file: c7cp06136f-f7.tif
Fig. 7 Top view of 8 × 9 replicas of the simulation cell indicated by the grid. Each coloured dot represents the ending position of the admolecule (CO2 (top), H2O (middle), and CH4 (bottom)) in one simulation taking into account the periodic boundary conditions. The admolecule started in the central cell. The figure hence displays the travelled distance and direction of the admolecules for each “adsorption” trajectory. The colour coding represents the initial translational energy of the admolecules.

Intuitively, one would expect the dots to show a colour gradient in radial distance, with the highest initial translational energy resulting in the largest travelled distance. This can not be observed in Fig. 7. For all three admolecules, the distribution in travelled distance is rather independent of the initial energy.

Fig. 8 shows the overall radial distribution of the final location. Each dot represents the relative count of the number of simulations in which the travelled distance of the admolecule falls within a certain distance bin. Each bin is 1 Å wide. The dash-dotted curve represents this cumulative probability (right axis scale) and follows some sort of exponential dependence. We showed in our previous paper how this cumulative probability P(ri) can be applied to estimate the probability for follow-up reactions, where a newly formed product can move across the surface by non-thermal diffusion and meet another species for a subsequent reaction. We derived that the probability for the product S1 to reach new species S2 is

image file: c7cp06136f-t6.tif(7)
where a is the distance between two sites and 〈d〉 the average distance between S1 and S2. This can be estimated by image file: c7cp06136f-t7.tif with N(S2) the number of S2 species on a grain and s the surface area of a grain. From Fig. 8 we can obtain an empirical estimate for the cumulative probability P(ri)
P(ri) = exp(−(Ari)2/3)(8)
where A can be applied as a fitting parameter. The best fits are represented by the dashed curves and correspond to A = 0.349, 0.310, and 0.583 Å−1 for CO2, H2O, and CH4 respectively.

image file: c7cp06136f-f8.tif
Fig. 8 Distance travelled distribution of the admolecule over the surface (solid curves and left axis scale). The dash-dotted curves (right axis scale) represent the cumulative distribution which are fitted by eqn (8) (dashed lines).

4 Summary and conclusions

In conclusion, our study showed that the energy dissipation of translationally excited molecules on water ice surface is not immediate. Translational energy can be applied for diffusion and desorption of the reaction products. Our conclusions can be summarized as follows:

(1) H2O was found to have binding sites on top of the ASW substrate, CH4 fills the nanopores on the substrate, and CO2 has intermediate binding sites. The effect of the surface morphology on desorption and diffusion is strongest for CH4, whereas H2O is least affected.

(2) The exact location in terms of height of the admolecules is important for its desorption. This is in contrast with a crystalline ice surface, where the binding energy of the admolecule to the surface is the main determining parameter in the desorption probability. The desorption behaviour of all three admolecules were found to follow

image file: c7cp06136f-t8.tif(9)
with ω = 0.3 and Etransreact the translational component of the reaction energy. The effect of the binding site location is captured by the f factor, with f = 0.6, 0.8, and 0.5 for CO2, H2O, and CH4, respectively.

(3) As expected the travelled distances are shorter on ASW than on crystalline ice because of the depth of the initial binding sites. The probability P(ri) for a molecule to travel more than distance ri, i.e., the probability to find an admolecule in the range [ri;∞) follows the empirical relation

P(ri) = exp(−(Ari)2/3)(10)
with A = 0.349, 0.310, and 0.583 Å−1 for CO2, H2O, and CH4, respectively, on top of an amorphous water surface.

Finally, the current work only includes translational excitation. We used this type of excitation as a proxy to study the fate of excited molecules that are formed through surface reaction. We expect that, in general, part of this excitation energy will be converted to kinetic energy. Our aim for future study is to also treat vibrational and rotational excitation. We expect these types of excitation to be less effective for the desorption of the admolecule as well as the kick-out mechanism and the distance travelled by the admolecules.

Conflicts of interest

There are no conflicts of interest to declare.


AF and HMC acknowledge NWO for financial support (VIDI research program 700.10.427 and ECHO 712.014.004).


  1. A. Fredon, T. Lamberts and H. M. Cuppen, Astrophys. J., 2017, 849, 125 CrossRef.
  2. R. Garrod, I. H. Park, P. Caselli and E. Herbst, Faraday Discuss., 2006, 133, 51 RSC.
  3. R. T. Garrod, V. Wakelam and E. Herbst, Astron. Astrophys., 2007, 467, 1103–1115 CrossRef CAS.
  4. M. Minissale and F. Dulieu, J. Chem. Phys., 2014, 141, 014304 CrossRef CAS PubMed.
  5. A. Bacmann, V. Taquet, A. Faure, C. Kahane and C. Ceccarelli, Astron. Astrophys., 2012, 541, L12 CrossRef.
  6. G. Fedoseev, H. M. Cuppen, S. Ioppolo, T. Lamberts and H. Linnartz, Mon. Not. R. Astron. Soc., 2015, 448, 1288–1297 CrossRef CAS.
  7. H. M. Cuppen, C. Walsh, T. Lamberts, D. Semenov, R. T. Garrod, E. M. Penteado and S. Ioppolo, Space Sci. Rev., 2017, 212, 1 CrossRef.
  8. A. A. Boogert, P. A. Gerakines and D. C. Whittet, Annu. Rev. Astron. Astrophys., 2015, 53, 541 CrossRef CAS.
  9. J. C. Polanyi, Noble Lecture, 1986, Chemistry.
  10. A. Bäck and N. Marković, J. Chem. Phys., 2005, 122, 144711 CrossRef PubMed.
  11. D. R. Glowacki, W. J. Rodgers, R. Shannon, S. H. Robertson and J. N. Harvey, Philos. Trans. R. Soc., A, 2017, 375, 20160206 CrossRef PubMed.
  12. M. A. González and J. L. F. Abascal, J. Chem. Phys., 2011, 135, 224516 CrossRef PubMed.
  13. S. Zhu and G. Robinson, Comput. Phys. Commun., 1989, 52, 317 CrossRef CAS.
  14. T. J. Lee, J. M. L. Martin and P. R. Taylor, J. Chem. Phys., 1995, 102, 254 CrossRef CAS.
  15. L. J. Karssemeijer, G. A. de Wijs and H. M. Cuppen, Phys. Chem. Chem. Phys., 2014, 16, 15630–15639 RSC.
  16. C. Qu, R. Conte, P. L. Houston and J. M. Bowman, Phys. Chem. Chem. Phys., 2015, 17, 8172 RSC.
  17. R. W. Hockney and J. W. Eastwood, Computer simulation using particles, Hilger, Bristol, 1988 Search PubMed.
  18. S. Plimpton, J. Comput. Phys., 1995, 117, 1 CrossRef CAS.
  19. P. Jenniskens, D. F. Blake, M. A. Wilson and A. Pohorille, Astrophys. J., 1995, 455, 389–401 CrossRef CAS.
  20. J. A. Noble, E. Congiu, F. Dulieu and H. J. Fraser, Mon. Not. R. Astron. Soc., 2012, 421, 768–779 CAS.
  21. F. Dulieu, E. Congiu, J. Noble, S. Baouche, H. Chaabouni, A. Moudens, M. Minissale and S. Cazaux, Sci. Rep., 2013, 3, 1338 CrossRef PubMed.
  22. R. S. Smith, R. A. May and B. D. Kay, J. Phys. Chem. B, 2016, 120, 1979 CrossRef CAS PubMed.
  23. M. Minissale, F. Dulieu, S. Cazaux and S. Hocuk, Astron. Astrophys., 2016, 585, A24 CrossRef.
  24. J. Crouse, H. Loock and N. M. Cann, J. Chem. Phys., 2015, 143, 034502 CrossRef CAS PubMed.
  25. S. Andersson, A. Al-Halabi, G.-J. Kroes and E. F. van Dishoeck, J. Chem. Phys., 2006, 124, 064715 CrossRef PubMed.
  26. A. Yabushita, T. Hama, M. Yokoyama, M. Kawasaki, S. Andersson, R. N. Dixon, M. N. R. Ashfold and N. Watanabe, Astrophys. J., Lett., 2009, 699, L80–L83 CrossRef CAS.
  27. S. Andersson and E. F. van Dishoeck, Astron. Astrophys., 2008, 491, 907–916 CrossRef CAS.
  28. C. Arasa, S. Andersson, H. M. Cuppen, E. F. van Dishoeck and G. Kroes, J. Chem. Phys., 2010, 132, 184510 CrossRef.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c7cp06136f

This journal is © the Owner Societies 2018