© ESO, 2014

1. Introduction

Diffusion on icy dust grains is a fundamental process in the chemical evolution of molecular clouds (Tielens 2010; Vidali 2012; Herbst 2014). The main reason is that many of the key molecular species are believed to be formed through the diffusive Langmuir-Hinshelwood mechanism on dust grain mantles. At low temperatures, the surface chemistry is dominated by hydrogenation reactions, while at higher temperature, larger species become mobile and lead to the formation of complex organic molecules. In both domains, diffusion is often the rate-limiting step. Furthermore, diffusive processes determine the structure of the ice and lead to trapping of molecular species. Hence, diffusion strongly influences the conditions under which species are released back into the gas phase as the cloud collapses.

Despite its critical importance, diffusion is a poorly understood process in the field of astrochemistry. As the amount of chemical complexity in modern gas-grain simulation codes continues to grow, this lack of knowledge is becoming an increasingly serious bottleneck that limits the accuracy of the models.

The most recent astrochemical models are capable of simultaneously modeling the chemistry in the gas-phase, which is relatively well understood, and on the grain-surfaces, which is treated by a (possibly multilayered) stochastic or rate-equation method (Vasyunin & Herbst 2013; Chang & Herbst 2014; Garrod & Pauly 2011; Taquet et al. 2012). In these models, each species i is assigned a binding energy to treat desorption, E_bind,i, and an energy barrier for grain-surface diffusion, E_diff,i. Usually, the desorption energies are relatively well defined (at least for the stable species) from either experiment or theory, but this is not the case for the diffusion energy barriers. Because reliable data are often not available, most models assume the diffusion energy to be a universal, fixed fraction, f, of the desorption energy: E_diff,i = fE_bind,i.

The use of this fraction is a key limitation for the models, first and foremost because there is no fundamental physical argument for such a universal ratio to even exist. Instead, this ratio depends on the diffusing species, the substrate, the surface coverage, and possibly on temperature. For this reason the fraction f is very poorly constrained and values between 0.3 and 0.8 are used by the modeling community (Hasegawa et al. 1992; Ruffle & Herbst 2002; Cuppen et al. 2009). As shown by Vasyunin & Herbst (2013), however, the value of f seriously influences the outcome of the models. Secondly, from a microscopic point of view, there are also obvious problems with the concept of using a single diffusion and a single desorption barrier because they both vary strongly from site to site, especially in the amorphous ices in the interstellar medium (Karssemeijer et al. 2014b). Inclusion all these local chemical details is not possible for current models, however, and this might not even be necessary because they aim to provide a more macroscopic view. For this purpose, a single diffusion barrier per species might well be sufficient and may even be desirable, in view of simplicity. But then, this should be a well-constrained species- and environment-specific value. Especially when considering the efficiency of reactions with an activation barrier, it is crucially important to know whether the diffusion rate is higher or lower than the reaction rate. Finally, from a practical point of view, an accurate value of f is also important because it affects the conditions under which the accretion limit is reached and thereby whether or not modelers should use stochastic models.

The purpose of this research note is to provide the astrochemical modeling community with more accurate data on the energy barriers for surface diffusion and thermal desorption of two of the key astrochemical molecules, CO and CO₂, on two forms of crystalline water ice. With these data, the accuracy of the models can be improved by making the ratio f species-specific for at least these two molecules.

2. Computational methodology

The diffusion and desorption barriers presented in this note were determined using a computational approach developed by the authors in a recent series of papers. For a detailed description of the methods and force fields, we therefore refer to Karssemeijer et al. (2012, 2014b,a). Summarizing, the calculations are performed using the adaptive kinetic Monte Carlo (AKMC) technique (Henkelman & Jónsson 2001), as implemented in the EON software package (Chill et al. 2014)¹. This is an off-lattice kinetic Monte Carlo (KMC) technique that combines the atomistic detail from molecular dynamics simulations with the ability of probing long timescales of kinetic Monte Carlo. This makes it specifically suitable for studying the details of diffusive processes under astrochemical conditions. To describe the interactions between the molecules, the same set of force fields is used as in the previous studies. The H₂O molecules are described using the TIP4P/2005f (González & Abascal 2011) potential. Interactions between water and CO are described in Karssemeijer et al. (2014b) and interactions with CO₂ in Karssemeijer et al. (2014a). The desorption barrier for each of the surface-binding sites discovered by the AKMC simulations was calculated with respect to the energy of the structurally relaxed isolated substrate and the CO or CO₂ admolecule. Corrections accounting for the quantum mechanical zero point energy contribution to the desorption barriers were applied to all results following the method from Appendix B of Karssemeijer et al. (2014b), to allow for a direct comparison with experimental values.

The surface diffusion of CO and CO₂ was studied on the basal plane of two forms of hexagonal ice (ice Ih). The first is the most common form of crystalline ice, which is characterized by a random hydrogen bond network between the oxygen atoms, which sit on tetrahedrally coordinated lattice sites. This partially disordered structure also leads to a random pattern of dangling OH bonds sticking out of the surface, which are known to strongly influence the local binding energy on the surface (Batista & Jónsson 2001; Sun et al. 2012). The second form of hexagonal ice is the so-called Fletcher phase (Fletcher 1992). This form of ice also has a random hydrogen bond network in the bulk, but has an ordered dangling proton pattern on the surface, with the OH bonds aligned in parallel rows (along the y-axis). This ordered structure is believed to minimize the surface energy (Buch et al. 2008; Pan et al. 2008) and may therefore be the thermodynamically most stable form of hexagonal ice under typical dense cloud conditions. These two samples are referred to as the “disordered” and the “Fletcher” substrate in the remainder of this paper. The surfaces of the substrates are shown in Fig. 1.

Fig. 1 Surfaces of the hexagonal ice substrates used in this work. The proton-ordered Fletcher sample is shown on the left, the right panel shows the disordered substrate. Dangling OH bonds are shown in blue.

Both samples were created following the procedure given in Karssemeijer et al. (2014a) and contain 672 H₂O molecules with a bulk density of 0.94 g cm^-3. They have dimensions of 31 × 31 Ų, with periodic boundary conditions applied along the x,y-plane, parallel to the surface. Because we are only interested in the dynamics of the adsorbed CO or CO₂ molecules, the water molecules were constrained from movement, which prevents the morphology of the ice from evolving during the simulation. Because of this constraint, the substrate cannot fully accommodate the admolecules, which typically leads to slightly lower desorption and surface diffusion barriers (Batista & Jónsson 2001; Karssemeijer et al. 2014b). The magnitude of this effect was evaluated by also performing simulations on the Fletcher substrate where the topmost 448 H₂O molecules were completely free. When the water molecules were allowed to move, we refer to the substrate as “free”, when they were constrained from movement, we use the term “frozen”.

3. Results

Using the AKMC scheme, the surface binding sites of CO and CO₂ were determined on each of the hexagonal ice substrates. The number of binding sites is listed in Table 1 and the corresponding distribution of binding energies is shown in Fig. 2. The distribution of binding energies is the widest on the proton-disordered substrate for both CO and CO₂. This is because the local variation in the arrangement of the dangling OH bonds has a strong influence on the binding energy of small adsorbed molecules. On the Fletcher substrate, the binding sites show more similarity and the distribution of binding energies is sharper. Especially for CO₂, two distinct types of binding sites appear on this substrate. These are discussed in detail in Karssemeijer et al. (2014a). There is still some site-to-site variation, which arises from the disordered hydrogen bond network in the bulk of the ice. As expected, the binding energies on the frozen Fletcher substrate are systematically lower than on the free sample.

Fig. 2 Distribution of binding energies of CO and CO2 on the various crystalline ice substrates. The dashed lines indicate the time-averaged binding energies from the kinetic Monte Carlo simulations for CO and CO2 at T = 25 and T = 70 K.

The mobility of CO and CO₂ was studied by generating KMC trajectories at temperatures between 15 and 50 K for CO, and between 50 and 100 K for CO₂. From these trajectories, the diffusion constants, D, were extracted from the mean squared displacement of the respective admolecule as a function of time (Frenkel & Smit 2002): $\mathit{D}\mathrm{=}\underset{\mathit{t}\mathrm{\to }\mathrm{\infty }}{\lim }\frac{\mathrm{1}}{\mathrm{2}\mathit{dt}}\mathrm{⟨}\mathrm{\left|}{r}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{r}\mathrm{(}\mathrm{0}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}\mathit{.}$(1)Here, d refers to the dimensionality, which is equal to 2 in the case of surface diffusion. As shown in Fig. 3, the diffusion constants show an Arrhenius behavior as a function of temperature: $\mathit{D}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{0}}\exp \left(\mathrm{-}\frac{{\mathit{E}}_{\mathrm{diff}}}{{\mathit{k}}_{\mathit{B}}\mathit{T}}\right)\mathrm{·}$(2)By fitting this expression to the data, the effective activation barrier for diffusion, E_diff, and the pre-exponential factor D₀ were determined (see Table 1).

Fig. 3 Surface diffusion constants for CO and CO2 on the various crystalline ice substrates.

The KMC trajectories were also used to determine the binding energies. These were time-averaged over the whole simulation at temperatures of 25 K for CO and 70 K for CO₂. Because most time is spent in the strongest binding sites, the binding energies are at the edges of the distributions in Fig. 2. On the frozen substrates, the strongest binding energies are found on the proton-disordered ice for both CO and CO₂. Using these binding energies, the ratio E_diff/E_bind was determined for each adsorbate/surface combination. For CO, this ratio is 0.31 on average, whereas for CO₂ a value of 0.39 is found.

The simulations themselves provide very well constrained, reproducible energy barriers for diffusion and desorption, which results in only small uncertainties on the ratio f for a given substrate. When comparing the different crystalline substrates, however, variations in f on the order of 10% are observed because of the different dangling-proton patterns. Compared with the values between 0.3 and 0.8 used in chemical models, however, these variations remain small.

The ordered surface of the Fletcher phase leads to a highly anisotropic surface diffusion. CO and CO₂ are both found to diffuse more rapidly along the y-direction, parallel to the dangling OH bonds, than along the perpendicular, x, direction. To study this effect, separate analyses were made on the one-dimensional diffusion along these two directions by calculating the one-dimensional mean squared displacements and using Eq. (2). The diffusion is found to follow an Arrhenius behavior in one dimension as well, and the activation barriers for the x- and y-directions are also listed in Table 1. For both adsorbates, the activation energies are higher along the x- than along the y-direction. Similar to the effect on the binding energies, however, the effect is significantly stronger for CO₂ than for CO. For completeness, this analysis was also performed on the disordered structure, where the difference between the energy barriers between the two direction is less than half of the difference on the Fletcher-phase surface. The anisotropy on the disordered substrate arises mainly because of the finite size of the sample, which still introduces some inequality between the x- and y-directions.

4. Discussion

Surface diffusion and desorption are complex processes on the microscopic scale that strongly depend on the local molecular environment. To describe them efficiently in grain-surface chemical codes, however, a simplified treatment is needed. In present-day codes, this simplification is achieved through a global ratio f = E_diff/E_bind that is used to define the diffusion energy barrier based on the binding energy of the specific molecular species. From our calculations, we have determined f for CO and CO₂ based on the absolute values of E_diff and E_bind on crystalline water ice. Although using these species-specific values arguably is an a priori improvement over using a single value for all species, we discuss the validity and assumptions of the calculations below.

The first point to address is the crystalline nature of the water ice substrates we used. Because interstellar ices are believed to be mostly amorphous, the question naturally arises what effect this has on the diffusion to desorption energy ratio. From our previous work on CO dynamics on amorphous solid water (ASW) ice (Karssemeijer et al. 2014b), we know that the mobility of CO is strongly related to the presence of strong binding nanopores on the amorphous ice surface. These sites can have binding energies in excess of 200 meV and can effectively immobilize adsorbed CO molecules. However, these pores do not only increase the effective diffusion energy barrier, but also the average binding energy. Thus, the effect on f of the amorphicity of the substrate will be weaker than the effect on the binding and diffusion energies themselves, compared with crystalline ice. Although a direct calculation of f on ASW at T = 25 K is computationally not feasible, we calculated the ratio to be 0.42 and 0.36 for CO on amorphous ice where the surface pores are partially occupied with either three or six additional CO molecules (see Table 1).

These values for ASW with an increased adsorbate coverage are not only similar to those on crystalline ice, they may also be more relevant for the interstellar medium than values on the bare ASW substrates. The reason is that the ASW nanopores in molecular clouds are likely to be occupied by species like molecular hydrogen, which is far more abundant and diffuses more rapidly than CO or CO₂. In addition, there might not be so many nanopores in interstellar ices because they may simply not survive the long timescales because of pore collapse (Bossa et al. 2012) or because of energy release by exothermic reactions in the pore sites.

The second point is that in most of the calculations we have made use of frozen substrates, where the water ice itself cannot evolve in time. This simplification effectively lowers the energy barriers for both diffusion and desorption. The diffusion-to-desorption ratio, however, remains largely unaffected, as is clear from comparing the results on the Fletcher substrate, which we simulated in both the free and frozen form.

This discussion shows that chemical models with a relatively fast diffusion (f between 0.3 and 0.5) are most realistic for describing the reaction rates with CO and CO₂. Diffusion of these species is underestimated in models using higher ratios (0.5−0.8). This is hardly surprising because these high ratios are based on the experiments on H₂ formation through H atom recombination (Katz et al. 1999; Perets et al. 2005). For these light atoms, the mechanisms underlying diffusion and desorption may be very different from the small, van der Waals-bonded molecules considered here; especially given the discussion on the possibility of H atom tunneling (Cazaux & Tielens 2004; Hama & Watanabe 2013). Furthermore, the value of E_bind from Katz et al. (1999) was given as a lower bound, which means that the ratio f based on these experiments (0.77) is actually an upper bound.

As long as no additional data are available, we recommend using the lower ratios presented here for other small, neutral molecules on water-dominated substrates as well. For radical,

light atomic, and charged species, however, the diffusion and desorption mechanisms are most likely very different, and the low ratios we found here may not be appropriate at all. The same holds for diffusion on substrates other than water ice. Here, more research is clearly needed.

5. Conclusions

We have presented simulation results on the surface diffusion and desorption energy barriers for CO and CO₂ on crystalline water ice. The diffusion-to-desorption barrier ratio is one of the crucial ingredients in current astrochemical gas-grain models, and the accuracy of the models can be significantly improved by using species-specific ratios of 0.31 for CO and 0.39 for CO₂.

These new data close but a small gap of the missing information, and it is up to the astrochemical community to constrain these ratios as well as possible for the other key species and substrates.

Acknowledgments

This work has been funded by the European Research Council (ERC-2010-StG, Grant Agreement no. 259510-KISMOL), and the VIDI research program 700.10.427, financed by The Netherlands Organization for Scientific Research (NWO). We are thankful for the stimulating discussions on the so-called X-factor (f) among the participants of the Lorentz workshop on “Grain-Surface Networks and Data for Astrochemistry” (July 2014, Leiden, The Netherlands). These inspired us to write this note.

References

All Tables

All Figures

	Fig. 1 Surfaces of the hexagonal ice substrates used in this work. The proton-ordered Fletcher sample is shown on the left, the right panel shows the disordered substrate. Dangling OH bonds are shown in blue.
In the text

	Fig. 2 Distribution of binding energies of CO and CO2 on the various crystalline ice substrates. The dashed lines indicate the time-averaged binding energies from the kinetic Monte Carlo simulations for CO and CO2 at T = 25 and T = 70 K.
In the text

	Fig. 3 Surface diffusion constants for CO and CO2 on the various crystalline ice substrates.
In the text