The Hidden!: The Kinetic Monte Carlo Method as a Way To Solve the Master Equation for Interstellar Grain Chemistry

Herma Cuppen studied Chemistry at the Radboud University Nijmegen (then named University of Nijmegen), The Netherlands, where she earned her M.Sc. degree in 2000 and her Ph.D. in 2005 (with highest honors) on theory and simulations of crystal growth. During this time, she was introduced to the kinetic Monte Carlo technique. She then moved to The Ohio State University as a postdoctoral fellow in the group of Eric Herbst where she applied KMC simulations to astrochemically relevant surface processes. In 2006, she moved back to The Netherlands to join the Sackler Laboratory of Astrophysics in Leiden where she led the theoretical interpretation of experimental data and bridged the gap between laboratory and space conditions. She was offered an assistant professorship in the theoretical chemistry at the Radboud University Nijmegen in 2011, where she studies both interstellar ices and molecular crystals using a range of different computational chemistry techniques, including KMC and molecular dynamics. Herma Cuppen is recipient of a NWO VENI grant (2006–2009), an ERC Starting Grant (2011–2016), a NWO VIDI grant (2011–2016), and the 2005 Crystal Growth Award for the best doctoral thesis from the Dutch Crystal Growth Society.BIOGRAPHY

Leendertjan Karssemeijer was born in The Netherlands and obtained his B.Sc. and M.Sc. degrees in physics at the Radboud University Nijmegen. He is currently a Ph.D. candidate under the supervision of Dr. Herma Cuppen at the Institute of Molecules and Materials of the Radboud University Nijmegen. His project involves atomistic computational modeling, applying mostly AKMC and molecular dynamics, of long time scale kinetic processes, such as diffusion and segregation, in interstellar ice mantles.BIOGRAPHY

Thanja Lamberts studied chemistry at the Radboud University Nijmegen, where she obtained her B.Sc. and M.Sc. degrees (with honors). She is currently a Ph.D. candidate with a shared position between the Institute of Molecules and Materials at the Radboud University Nijmegen under the supervision of Dr. Herma Cuppen and Leiden Observatory at Leiden University under the supervision of Prof. Dr. Harold Linnartz. Her project involves the simulation of surface chemistry of water formation in interstellar medium, both experimentally on grain analogues (Leiden) and by grain surface models (Nijmegen). The latter aims at understanding the relevant experiments as well as modeling chemistry in astrochemical environments using lattice-gas KMC.

1 Introduction

ARTICLE SECTIONS

Jump To

The interstellar medium (ISM) is far from empty; rather it contains large molecular clouds consisting of dust and gas. It is in these clouds that new stars form—possibly with habitable planets—with molecules playing a crucial role.(1, 2) At the moment, over 170 different molecules have been detected: stable molecules, radicals, cations, and anions. Many molecules only consist of a few atoms and are familiar to us from a terrestrial point of view, such as water, molecular hydrogen, methanol, formaldehyde, formic acid, and dimethyl ether. Other molecules are more exotic and consist, for instance, of very unsaturated carbon chains. Gas phase chemical models show that these latter molecules can be quite easily formed through gas phase reactions.(3) The saturated molecules are however only formed in small quantities through gas phase chemistry, since they require either high temperatures to be formed or three-body reactions, which are both not available in cold molecular clouds. Dust grains can act as a third body, facilitating addition reactions that lead to the formation of saturated molecules. Already in 1949 de Hulst(4) realized the importance of grain surface chemistry and he suggested surface formation routes for H2 from two hydrogen atoms and the formation of H2O from O2 with hydrogen peroxide as intermediate.

From various sources of information, we know that approximately 1% of the mass of the material in the ISM, and 10–10% in terms of numbers, consists of small, nanosized silicate and carbonaceous particles called dust grains. Although this seems a small portion of the total amount of matter, dust grains play an important role in the ISM, not only by acting as catalytic sites for molecule formation but also by shielding molecules from dissociating UV radiation. Depending on the physical conditions, the dust grains can be coated with a layer of “dirty” ice consisting of a mixture of different species. Their main components are water, carbon monoxide, and carbon dioxide, but also more complex molecules such as methane, methanol, formaldehyde, and ammonia may be present. Some of the ice species accrete from the gas phase, but as mentioned before, many of the simple important molecules such as H2 and H2O and also several complex organic molecules are not formed in the gas phase, but rather on the grain surfaces themselves. As a result, many modelers have begun to include grain surface processes in their models to describe chemical evolution during, for instance, star formation.(1, 5, 6) Originally, the development of grain surface models was mostly driven from the astrochemical modeling community(7-9) without much input from specialists in ice chemistry. Gas phase models were extended with a surface phase and the grain surface chemistry was treated in a way similar to that of gas phase chemistry. This was found to lead to problems in the regime where the surface abundance is low—the so-called “accretion limit”—where formation rates could be orders of magnitude off. Since then, different modeling methods have been applied to tackle this and other modeling problems.

From the physicochemical point of view, refinements have been made as well. Many of the initial expressions and input parameters originated from surface science or were extensions from gas phase experiments. In surface science, mostly metal or semiconductor surfaces are studied on which absorbates are deposited within the monolayer regime, around room temperature. This is quite different from the rough grain surface on which inhomogeneous ice layers of tens of monolayers are formed at very low temperatures. Fortunately, technological developments have now enabled the experimental verification of grain surface chemistry. Many of the assumptions within the models are now put to the test by applying to experimental surface science techniques to real interstellar analogues. With this new experimental information comes the need for detailed surface models to understand the physicochemical processes behind the experimental results and to translate the experimental findings to astrochemical time scales.

Different modeling techniques have been applied to grain surface astrochemistry, covering a large range of time and length scales. More molecular detail comes at the expense of more CPU time to cover the same evolution in real time. How different techniques relate to one other in this respect is indicated in Figure 1. Molecular dynamics simulations trace the exact location and orientation of the molecules including lattice vibrations and in some cases even the intramolecular movement, but they typically stay within the picosecond to nanosecond time frame. Rate equations on the other hand can easily handle 108 years—much longer than the lifetime of a molecular cloud—but adsorbed molecules are only treated in terms of numbers and their exact locations are not known. The present review aims to discuss the application of the kinetic Monte Carlo technique (KMC) to grain surface chemistry. This method was initially introduced into astrochemistry as a solution to the “accretion” limit, but it is now more and more applied to gain insight into the physicochemical surface processes. The method has, in principle, not a real restriction in how molecules are represented: in terms of number densities or with their exact location and orientation. The different implementations of kinetic Monte Carlo can therefore cover a huge time and length scale range.

Figure 1

Figure 1. Overview of the different simulation methods mentioned in the present review.

The kinetic Monte Carlo technique is a way to solve the master equation, and this review will start by deriving this expression in section 2. Several modeling methods have been used to simulate grain surface chemistry, and section 3 gives a historical overview of the introduction of these techniques into astrochemistry and their main advantages and disadvantages. The reason to include a short discussion on these other techniques is that this allows us to put the results on kinetic Monte Carlo into a better perspective. The review will then continue by discussing in section 4 some of the technical aspects involved in constructing a KMC model and the input parameters that are needed for such a model. Sections 5 and 6 summarize the application of the kinetic Monte Carlo technique in astrochemistry. Here the distinction is made between simulations that try to reproduce astrochemically relevant experimental results and simulations that attempt to predict astrochemical environments. Here, only models that treat the grain with some microscopic detail will be discussed, starting from the first model in 2005 until the beginning of 2013. Finally, some recent advances in the application of kinetic Monte Carlo are discussed in section 7. We will go into detail on the challenges and opportunities involved.

2 Theoretical Background

ARTICLE SECTIONS

Jump To

Kinetic Monte Carlo is a method to solve the master equation. We will give a brief derivation of this equation in order to discuss the underlying assumptions to the derivation of this expression, the different ways in which this equation is solved, the approximations involved, and their advantages and disadvantages. A more thorough derivation can be found in textbooks such as refs 10−12.

Many processes occurring on the surfaces of grains are in principle stochastic processes. We can describe these systems by a state vector x and a time coordinate t. In this section, x is an abstract quantity which may be some local minimum on the potential energy surface (PES), which is schematically drawn in the left panel of Figure 2. Here the black circles represent different states. In many grain surface models a PES is however not explicitly considered and in those cases a state x can be represented by, e.g., the species on the grain, their position, the temperature of the grain, etc., or x can simply contain the coordinates of all particles on the grain.

Figure 2

Figure 2. Schematic representation of a potential energy surface (PES). The states, local minima on the PES, are indicated by black solid spheres. (left) The system moves from state to state, or (right) the vibrational movement is included as well and the trajectories spend most of their time around the local minimum.

Let us call p(x, t) the probability density for the system to be in state x at time t. The probability density can be written as

(1)where P(x, t) is the probability that the system is in an infinitesimal volume element dx around x at time t. Since the probability to find the system in any x at a given time t is unity, we have

(2)The definition of the probability density, eq 1, can be extended to an nth order probability density function p(xn, tn; xn–1, tn–1; ...; x1, t1) which gives the probability density for the system to be in x1 at t1and in x2 at t2, etc. To describe the time evolution of the system, the conditional probability density is introduced

(3)which gives the probability density to be in states xn+1 through xm at times tn+1 through tmgiven that it is in states x1 through xn at times t1 through tn. The higher order conditional probability densities are complicated functions, and therefore one usually assumes a Markov chain. To define a Markovian system, we imply time ordering (t3 > t2 > t1) and impose the following restriction on the conditional probability density functions:

(4)This means that the probability of arriving in xn, tn is only dependent on xn–1, tn–1 and not on the states before tn–1. A Markovian system is completely defined by an initial state, say x1, t1, and the first order conditional probability densities since any other probability density function can be constructed from these. For example

(5)where we used eqs 3 and 4.

The Markov chain assumption is valid for memoryless processes. When simulating grain surface chemistry, we are typically interested in rare events, such as diffusion and desorption. This type of events occurs on much longer time scales than the (lattice) vibrations, and the event time scale and the vibrational time scale become decoupled. This makes rare events effectively memoryless, since all history about which states were previously visited is lost due to the vibrations between two transitions. This is schematically depicted in Figure 2. The residence time of the system in a local minimum of the potential energy surface is several orders of magnitude larger than the vibrational time scale. The system will therefore typically proceed according to a trajectory like the one drawn in the right panel of Figure 2. At the time that the system leaves the potential well, all information about the direction from which it initially entered this potential well is lost. In kinetic Monte Carlo this is approximated by a trajectory as drawn in the left panel of Figure 2. In some cases, when the residence time in a state is short as compared to the vibrational time scale, this approximation breaks down, as we will discuss in sections 5.1 and 5.2. Examples of these are H atoms diffusing on a graphite surface without feeling a potential (section 5.1) and the “hot” reaction products which move away from each other with a certain momentum upon reaction (section 5.2). For these cases, molecular dynamics simulations, which use Newton’s equations of motion to determine the molecular trajectories, are better suited. These types of simulations would result in trajectories plotted in the panel on the right in Figure 2. A drawback of this technique is however that their simulation time scales do not match astrochemical time scales—they are shorter by roughly 20 orders of magnitude—which makes them less suitable for astrochemical simulations. Individual processes with astrochemical relevance can however be treated with molecular dynamics.

The transition probabilities of Markov processes obey the Chapman–Kolmogorov equation:

(6)Equation 6 describes the transition from x1, t1 to x3, t3 through all possible values of x2 at t2. Let us now assume that the conditional probability densities are time independent, so we can write

(7)because there is now only a dependence on the time difference Δt. The expression pΔt(x2|x1) denotes the transition probability density from state x1 to state x2 within the time interval Δt. This can be expressed as a rate, k, by

(8)Using the Chapman–Kolmogorov equation, eq 6, for pΔt, we get

(9)After some manipulation of this expression, we arrive at the master equation for the transition probability densities

(10)If we want to obtain a master equation for the state probability density itself, we can multiply this equation by p(x1, t) and integrate over x1. This will ultimately lead to the master equation

(11)On a discrete state space, which is usually used for grain simulations, the equation is written in terms of probabilities in the following way:

(12)The first part of the sum in eq 12 represents the increase in P(xi, t) because of events that enter state xi from any other state xj; the second part is the decrease in P(xi, t) due to all events leaving xi. Most chemists will construct this expression intuitively. Here we made the formal derivation of the master equation to show the assumptions involved at the root of this expression which could be a severe restriction, for instance the Markov chain assumption.

3 Historical and Theoretical Overview of Grain Modeling Techniques

ARTICLE SECTIONS

Jump To

3.1 Rate Equations

One of the most applied methods in astrochemistry is solving rate equations numerically. These equations model the abundances of species and can be seen as the deterministic counterpart of the stochastic master equation. They do not describe the probabilities, P, to be in individual states, but rather the surface abundances, ns, of the different species on the grain. Since a similar quantity is used to describe the gas phase chemistry (the gas phase abundance ng), these types of equations can be easily coupled to gas phase chemistry networks. The first gas–grain simulations were therefore rate equation based.(-8, 13-15) Rate equations are often applied in chemistry to describe macroscopic (experimental) effects and account for many-body effects with a mean-field approach. A set of rate equations to describe surface reactions in astrochemistry was first constructed by Pickles and Williams(16) to describe the formation of molecular hydrogen from H atoms and water from O2 and H through HO2 and H2O2.

Grain species can generally undergo four types of processes: accretion onto the surface, desorption, diffusion, and reaction. This leads to the following expression for the change in surface abundance of species A at time t

(13)where Nsites is the number of sites on the grain. Here the subscripts “s” and “g” of n represent the surface and gas phase abundances of the species, respectively. The first two terms in this expression denote the gain and loss of species A due to reactions which form and destroy A. These are described in terms of reaction rates kreact and hopping rates khop. The third term accounts for the accretion of new species from the gas phase, in terms of the accretion rate kacc. The final term represents the desorption of the species into the gas phase. This can either be due to thermal desorption (kevap) when the temperature is high enough to overcome the binding energy or be due to nonthermal effects (knon-th) such as photodesorption,(17-20) sputtering by cosmic rays, or grain heating by cosmic rays.(21, 22) Here first order kinetics are assumed for the desorption, although experimentally zero order kinetics are often observed for the desorption of ices in the multilayer regime. We will come back to this in section 4.2.2.

The reactions included in eq 13 are assumed to occur through the Langmuir–Hinshelwood mechanism. Generally, the distinction is made between three types of surface reactions that are schematically depicted in Figure 3: the Langmuir–Hinshelwood (LH) mechanism, which is a diffusive mechanism in which the species move over the surface and try to react upon meeting; the Eley–Rideal (ER) mechanism, where one (stationary) reactant is hit by another species from the gas phase; and the hot-atom mechanism, which is a combination of both mechanisms where nonthermalized species travel some distance over the surface finding reactants on their way. In all cases, the excess heat which becomes available through the reaction can be used for desorption of the products. Since in dense clouds the temperature of the gas phase is similar to the grain temperature, the hot-atom mechanism is generally not considered. It can be important in experimental conditions or for instance in shock regions where the gas phase is often much warmer than the grain. The diffusive LH mechanism leads to second order rate kinetics as can be seen from the first two terms of eq 13.

Figure 3

Figure 3. Schematic representation of the three mechanisms for surface reaction: the Langmuir–Hinshelwood mechanism, which is a diffusive mechanism in which the species move over the surface and try to react upon meeting; the Eley–Rideal mechanism, where one (stationary) reactant is hit by another species from the gas phase; and the hot-atom mechanism, which is a combination of both mechanisms where nonthermalized species travel some distance over the surface finding reactants on their way.

Eley–Rideal is generally considered to be only important for high surface coverages or low surface mobility(23) and is therefore not explicitly considered in most models. However, catastrophic freeze-out of CO in prestellar cores results in a layer of CO ice,(24) i.e., a high surface coverage of reactive species. Under these circumstances, Eley–Rideal could be an important mechanism in the formation of methanol. Rate equations that include surface reactions following the Eley–Rideal mechanism contain an additional term that depends on the gas phase abundance of one reactant, the surface abundance of the other reactant, and the rate of reaction between them.(25) The Langmuir–Hinshelwood expression for the same reaction contains the surface abundance of both reactions, the reaction rate, and a meeting rate, as can be seen in eq 13. The meeting rate due to the movement of species A is given by khop,Ans(A,t)/Nsites. The hopping rate khop,A used here is the rate of an individual hop of species A. Often, however, the hopping rate quoted in the literature dealing with rate equations is the hopping rate of an individual hop divided by the number of sites on the grain. It is worth mentioning that the reaction rate is different for Eley–Rideal and Langmuir–Hinshelwood mechanisms: for ER the two reactants have one attempt to cross the reaction barrier, whereas for LH the two reactants will remain adsorbed in close vicinity until they react or one of them diffuses or desorbs. Reaction is therefore in competition with diffusion and desorption.(26) Rate equation models, however, often ignore this.

To relate the deterministic rate equations to the stochastic master equation, one can think of the abundance of the species A in terms of the expectation values of the number of particles of type A:

(14)where the brackets denote the expectation value and NA(t) is the number of particles of type A on the grain at time t. On the right of expression 14 we have the probability of the state vector xA[m] which describes a state with m particles of species A. An obvious disadvantage of the rate equation method is that it does not treat statistical fluctuations. One of the assumptions underpinning the rate equations leads to following assumption:

(15)When the number of species on the surface is small, this assumption breaks down and the system reaches the accretion limit as it was called by Charnley et al.(9) In this case the rate limiting step for reaction is not diffusion but the arrival of species on the surface from the gas phase. Figure 4 shows the evolution of several surface species as calculated by Stantcheva and Herbst(27) using the rate equation method and the “master equation method”, which will be discussed in section 3.2. The main difference between the two methods is that the rate equations use expectation values whereas the master equation method does not. We show these models here to illustrate the effect of the accretion limit. The results from Stantcheva and Herbst show that the abundances of the probabilistic species, such as H atoms and other reactive species, are highly overestimated by rate equations. The most extreme case is for atomic H, shown in the top panel of Figure 4, where the difference is many orders of magnitude. Since these reactive species are minority species which cannot be observed directly, we are usually more interested in stable reaction products as plotted in the bottom panel of Figure 4 for the CO hydrogenation route leading to methanol. Also here the effect is striking, with CH3OH being underproduced by orders of magnitude depending on the age of the cloud. In the master equation model, all CO molecules appear to be converted almost immediately to methanol, while in the rate equation model the conversion occurs much later. This discrepancy arises from the huge difference in surface abundance of atomic H and the inefficiency of the H + CO reaction, which proceeds with a barrier.

Figure 4

Figure 4. (top) Absolute surface abundance of atomic H at 10, 15, and 20 K and (bottom) fractional abundances (relative to H2) for CO, H2CO, CH3OH, and CO2 at 10 K as calculated by master equation (meq) and rate equation (req) methods. Reprinted from ref 27. Credit: Stantcheva, Herbst, Astron. Astrophys. 2004, 423, 241, reproduced with permission © ESO.

To account for the accretion limit problem in the rate equations, alternative methods have been proposed which either adjust the rate equations in the accretion limit regime to result in the correct rates or do not use the expectation value assumption and are so-called stochastic methods. The master equation method, applied in Figure 4, is an example of the latter. Caselli et al.(28) were the first to propose a way to adjust the rate equations. They applied a semiempirical approach to scale down the reaction rates for those cases where the surface migration of atomic hydrogen is significantly faster than its accretion rate onto grains. This method appeared to give a good agreement with stochastic methods for a limited number of cases. Unfortunately, due to the ad hoc way in which the rates were scaled down, the validity of the results outside the tested regime cannot be guaranteed.

More recently, a new modified rate approach was suggested by Garrod.(29) This method is more rigorously derived than the original modified rate approach. A modified rate is determined by the accretion rate of one of the reactants multiplied by the surface abundance of the other reactant times an efficiency term that takes into account the competition between diffusion across the surface for the reactants and the desorption of the reactants. The expressions for the reactions are adjusted if the modified reaction rate is smaller than the classical diffusive rate. The main advantages of these modified rate equation models are that they are computationally inexpensive as compared to the stochastic methods and that they automatically switch to the normal rate equations in the regime where those are still valid.

3.2 Master Equation Method

The accretion limit problem led to the realization that grain surface processes can be better described by a master equation than by rate equations. This resulted in the introduction of several techniques that solve this master equation. One of these constructs differential equations for the probabilities to have a certain number of species on the surface and uses a numerical integrator to obtain the average number of species on the grain as a function of time. This method was introduced into astrochemistry to circumvent the accretion limit problem(30-32) and is usually referred to as the “master equation method”. This can be confusing since the kinetic Monte Carlo approach, which is discussed in section 3.3, also solves the master equation. A better description, although rather lengthy, would perhaps be “numerical integration of master equations”, but in line with previous authors we will stick to using the term “master equation method” in the remainder of this review.

Since most gas phase models also use numerical integration, the master equation approach can again be easily integrated within a gas phase astrochemical code. If only two types of species, A and B, are considered, the expressions that have to be integrated become of the type

Here P(xA,B[k, l], t) denotes the probability that k species A and l species B reside on the grain at time t. One can again recognize the reaction, accretion, and desorption terms. The factors (k + 1) and (l + 1) in the first reaction term reflect the number of possible combinations to form new molecule AB of k + 1 species A and l + 1 species B and leave k species A and l species B behind. In principle, one should construct these expressions for all possible combinations of species A and B (P(xA,B[0, 0], t), P(xA,B[1, 0], t), P(xA,B[0, 1], t), P(xA,B[1, 1],t), etc.). Obviously, these expressions become much more involved if more types of species are included in the model. Theoretically the number of expressions to solve is infinite. However, cutoffs can be used to exclude the higher order terms, if their probability of occurrence becomes very small. In other words, one can set upper limits on the maximum number of particles per species that can reside on the grain, Nmax(A). Stantcheva et al.(32) used this approach for a small network containing 11 species connected through 10 surface reactions. Only the five minority species O, H, OH, HCO, and H3CO were treated probabilistically; the other species were described using rate equations. For the species O and H the maximum number on the grain considered never exceeded five. The resulting expressions, both master equations and rate equations, were integrated numerically using a numerical integrator. In this way, they were able to simulate a small network of species under conditions where the abundance of these minority species remains low.

If the number of particles expands, either by a change in physical conditions (Nmax grows) or by increasing the number of considered species in the model, the number of equations blows up rapidly. Algorithms have been suggested to make the effect less severe and to extend the range in which the master equation method is applicable.(33, 34)

The lowest possible cutoff would be Nmax = 2 for species that form homonuclear diatomic molecules (H2, O2, etc.) and Nmax = 1 for all other species, but obviously this would lead to a large error. To make this less dramatic, one can convert the original master equations into moment equations.(34-36) Let us consider moments of the distribution P. The kth moment is defined as

(17)and the expectation value used in the rate equation (eq 14) is therefore simply the first moment. Differently from the master equation approach, which is expressed in differential equations of the probabilities P(xA,B[k, l], t), the moment equations are differential equations of the different moments ⟨NAkNBl⟩. These equations can be directly obtained by expanding the time derivatives of the moments using eq 17 and inserting the expression for P(xA[m]) from the master equation. Obviously, still an infinite number of equations can be composed, with the order of the moment running to infinity instead of the number of species. But the expressions for the first moments, the population sizes we are interested in, are completely determined by first and second moments, and we could use the second order as a upper cutoff. The second moment equations, however, are determined by first, second, and then unknown third moments. In order to close the set of moment equations, Barzel and Biham(34) imposed the following constraint on the moment equation: at any given time, at most, two atoms or molecules can be adsorbed simultaneously on the surface. Furthermore, these two atoms or molecules must be from species that react with each other. In a later paper, Barzel and Biham(37) introduced a diagrammatic way to visualize a surface reaction network and to aid in the construction of the moment equations to make the method more applicable for grain surface chemistry in general. This formulation was applied in ref 38 to incorporate stochastic chemistry in the form of moment equations into an astrochemical code to simulate photon dominated regions. The specific focus was on the formation of H2, D2, and HD.

While the constraint that at any given time, at most, two atoms or molecules can be adsorbed simultaneously on the surface can work for a small network under specific physical conditions where the surface density is low, this constraint cannot be obeyed in the general case and application of moment equations to large grain surface networks could result in undesirable effects. For this reason, Du and Parise(39) introduced a hybrid method that combines moment equations and rate equations. In general, both gas phase and surface reactions are treated by rate equations, but if the average population of a certain surface species becomes smaller than one per grain, the algorithm switches to moment equations.

3.3 Macroscopic Kinetic Monte Carlo

The first astrochemical application of simulations based on the Monte Carlo method was performed by Tielens and Hagen.(7) They used the method to simulate the evolution of the grain mantle by one-by-one accretion of gas phase species and calculating their fate (reaction, desorption, trapping) before a new species would land on the surface. In this way, the time evolution is determined by the accretion rate and the Monte Carlo aspect determines the fate of accreted species. Steady-state gas phase abundances were used to determine the relative accretion rates for each of the gas phase species. These were obtained from a pure gas phase model and gas–grain interactions were not directly included in this calculation other than using a depletion parameter which decreases the relative elemental abundances of the heavier species such as O, N, and C toward clouds with higher densities. Using the resulting gas phase abundances, the mantle composition was calculated as a function of visual extinction AV or radiation field, cloud density, and depletion factor. Since the results of the simulations were presented as the steady-state fractional composition of the grain mantles, the time evolution of the grain was not discussed.

Charnley(40, 41) was the first to introduce the kinetic Monte Carlo method to determine the chemical evolution of a cloud, including both gas phase and grain surface. This was done in a more rigorous way than by Tielens and Hagen(7) and built strongly on the work by Gillespie,(42) who used the KMC method to numerically simulate the stochastic time evolution of coupled chemical reactions. Both Charnley and Gillespie started from the master equation (eq 12). Using this equation, one can derive the time interval between two events, which is one of the fundamental expressions in kinetic Monte Carlo, since it describes how time progresses during a KMC simulation. If one starts from state i at time t (so P(xj, t) = 0 ∀ j ≠ i), then for small Δt we have

(18)from which we can calculate the probability that the system has not left state xi at time t + Δt

(19)where ktot(xi, t) is the total rate leaving state xi. If we assume stationary rates, then eq 19 becomes

(20)Stationary rates are a reasonable assumption for grains of constant temperature. If the surface temperature fluctuates, for instance due to stochastic heating of the grains by incoming UV photons or cosmic rays(21, 22, 43, 44) or in a temperature programmed desorption experiment, the rates change over time and the time interval Δt should be determined differently.(45) This will be discussed in section 4.

In practice, this means a KMC cycle starting from state xi consists of three steps if all possible events are known: first a final state xf is picked with probability

(21)using a random number. Next, the time is advanced where the values for Δt are chosen such that they follow the distribution dictated by eq 20. This can be numerically achieved by generating a uniform random number Z in the range (0, 1] and equating this to the probability that the reaction has not yet occurred:

(22)which leads to

(23)for stationary rates. Finally, the picked transition is evaluated and the transition rates leaving this new state, xf, are determined. The consequence of eqs 21 and 23 is that at each given time all possible transition events leaving the current state and their corresponding rates should be known. This combination of events and rates is often referred to as “table of events”. For the systems discussed so far, the table of events can be easily constructed, since as in the rate equation approach and numerical integration of the master equation approach all possible events are known at the start of the simulation. If we want to look at grain surface chemistry in more microscopic detail, this can become problematic as we will discuss in sections 4 and 7.

The KMC realization of the master equation is not as easily implemented in a gas–grain astrochemical code as the rate equation and master equation methods. The reason for this is that for the gas phase usually rate equations are used and the Monte Carlo algorithm and numerical integrators cannot be easily coupled. Chang et al.(46) presented an iterative solution where the grain surface chemistry is treated by a microscopic KMC (see section 3.4) method and the gas phase is treated by means of rate equations. Attempts to couple gas phase models to grain chemistry models are discussed in more detail in section 7.

3.4 Microscopic Kinetic Monte Carlo

The methods discussed so far only considered the abundances or numbers of the different surface species and not their exact location and exact movement on the grain. The diffusion of the different species is included in some average way. In microscopic KMC simulations, the location of the particles is known. They undergo a random walk and they can revisit sites multiple times. This revisiting is often referred to as back diffusion. Although it is possible to express this effect in rate equations,(47) it is included in a more straightforward way in microscopic KMC simulations. In these kinds of simulations, accretion, desorption, diffusion, and reaction are again considered, with rates similar to those in the macroscopic KMC method, but with the difference that for each individual atom or molecule on the surface its specific position on the grain is monitored. Diffusion now occurs through individual hops of the species from one site to the next. This includes hopping back to where it came from. The information of the position of the atoms/molecules on the surface is an advantage over macroscopic KMC.

Another advantage is that the surface structure can be included. The rate of diffusion and desorption can be made site-specific, mimicking steps on the surface, amorphicity of the grain, and/or inhomogeneity of the grain material. The species might undergo a hindered random walk that does not result in an exact second order for reaction. Since, in microscopic KMC, only species either in neighboring sites or occupying the same site are allowed to react and the reaction order is not included specifically beforehand, this deviation from second order will come out naturally. In the same way, crossing a reaction barrier is treated in competition with diffusion and desorption automatically in microscopic KMC.(48)

One can imagine that, while the grain mantle grows, species can become trapped in the lower layers and are not able to participate in the chemistry that occurs on the surface of the formed ice, or that these species will not be able to desorb. Again this will come out naturally from the microscopic Monte Carlo simulations. We will come back to this point at a later stage (sections 4, 5, and 7).

The first implementation of microscopic KMC in astrochemistry was performed by Chang et al. and focused on H2 formation and was called continuous-time, random-walk (CTRW) Monte Carlo(49) following the algorithm by Montroll and Weiss.(50) We will refer to this application as microscopic KMC throughout this review to set it apart from macroscopic KMC. When referring to the algorithm, we will still use the term continuous-time, random-walk KMC. In their study, Chang et al. used the microscopic aspect of the method to study the effect of back diffusion and of site-specific diffusion and desorption rates. The site-specific diffusion and desorption rates were applied to mimic either the amorphous character of the grain or the effect of composite grains of carbonaceous and silicate material. The simulations start with an empty grain, but depending on the temperature and surface structure of the grain, the H-atom surface abundance slowly builds up until a steady-state value has been reached. Chang et al. determined the steady-state recombination efficiency

(24)where NH is the number of arriving H atoms and NH2 is the number of formed H2 molecules as a measure for the formation rate of H2. Back diffusion was found to be of importance, through a comparison of the η obtained with a rate equation model, a master equation model (both from Biham et al.(30)), and two microscopic KMC models differing in the amount of nearest neighbors (Figure 5). Since the likelihood of revisiting a site decreases with the number of nearest neighbors, more H2 is formed for surfaces with a higher coordination number, although still less than predicted by the rate equation and master equation models. Since the latter models mimic the limit in which no back diffusion is present, they easily overestimate the efficiency. Figure 5 shows an overestimation of a factor of 2–3 at 10 K. Additional simulations were performed with site-specific rates. The inhomogeneity of the diffusion and desorption rates was found to have a strong effect on the temperature window, where H2 can be efficiently formed, much stronger than back diffusion. Molecular hydrogen could be produced up to higher temperatures than in the case of homogeneous rates.

Figure 5

Figure 5. Hydrogen recombination efficiency as calculated by the microscopic KMC approach for homogeneous olivine and amorphous carbon as a function of surface temperature. For olivine, the circles represent a lattice with four nearest neighbors while the squares represent a lattice with six. For carbon, the triangles pointing upward and downward represent, respectively, lattices with four and six nearest neighbors. The master equation/rate equation results of Biham et al.(30) for olivine are shown as diamonds. Reproduced from ref 49. Credit: Chang, Cuppen, Herbst, Astron. Astrophys. 2005, 434, 599, reproduced with permission © ESO.

4 Technical Aspects of KMC

ARTICLE SECTIONS

Jump To

This section will discuss the technical aspects of the kinetic Monte Carlo technique applied to surface astrochemistry. Sections 2 and 3 discussed the theory behind KMC. In this section, we will go into more detail on the implementation of the theory and its technical limitations. As discussed in section 2, the theory behind KMC demands the formation of a table of events. We will discuss the type of information that is needed to construct such a table, how this is obtained, and what the common assumptions are. Different applications of the KMC technique require different algorithms, and therefore a wide selection of different algorithms has been developed, usually with only slight variations. Here, three algorithms that are most commonly used within astrochemistry will be discussed, together with their advantages and disadvantages. Additionally, an algorithm which is able to handle nonstationary rates will be discussed as well. Let us start by discussing the representation of the grain in the simulations.

4.1 Representation of the Grain

As discussed in section 2, KMC follows the evolution of a Markov chain from one state, or local minimum on the potential energy surface, to the next. These states are usually put on a regular grid represented by a two- or three-dimensional array where for each grid point the occupancy is monitored. A reaction or a diffusion from one site to another will be modeled as a sudden change in the occupancy of the sites. A model that puts sites on a grid and only registers their occupancies is often referred to as a lattice-gas model.

One of the largest disadvantages of using a lattice model is that a large part of the molecular detail is lost in these models. Because of the amorphous character of the grain and ice, the local environment of the sites probably varies across the surface and the states are most likely not at exact equidistant locations. Nevertheless, the bare grain or its ice mantle is usually considered to have well-defined positions on its surfaces where other species can adsorb and these sites are probably distributed with some limitations in terms of number of neighbors. In this sense, there is some regularity, even for amorphous grain mantles, which would justify the use of lattice-gas models. Deviations from the regular character of the grid can be accounted for by using site-specific rates as was done in ref 49.

On the other hand, since the molecules are confined to lattice positions, using realistic molecular potentials can become problematic. Furthermore, making guesses on the possible events is far from straightforward and we can easily miss important mechanisms. Therefore, ideally, one wants to move away from lattice models. In computational chemistry and statistical physics so-called off-lattice models are more commonly used for Monte Carlo simulations. These Monte Carlo simulations are however not kinetic Monte Carlo simulations, but what we will refer to as Metropolis Monte Carlo simulations based on the early work by Metropolis and Ulam(51) and Metropolis et al.(52) These types of simulations do not follow the time evolution of a system but rather obtain an ensemble average of certain properties of the system, and their algorithms only involve one trial configuration at each point in the simulation, which does not necessarily need to be a minimum on the PES. Kinetic Monte Carlo on the other hand allows simulation of the evolution of a system in time, but the drawback is that the simulation runs from state to state and all accessible processes at a given cycle in the simulation need to be known. This difference in approach has a direct effect on the effort it takes to use an off-lattice method. In Metropolis Monte Carlo this can be achieved fairly easily, whereas for kinetic Monte Carlo this is not so obvious. Examples of off-lattice kinetic Monte Carlo techniques will be discussed in detail in section 7.3. For now, we will limit ourselves to lattice-gas simulations, since these are most commonly used in astrochemistry.

A clear advantage of a lattice-gas model is that one can work with a predefined event table. Since the molecules are confined to a lattice, only a limited amount of events is possible of which the rates can be determined at the start of the simulation. The limited amount of possible events in lattice-gas models allows one to cover very large time scales. The simulations can really mimic laboratory experiments of ices, covering hours of simulation time and tracing similar properties.

4.2 Input Parameters in Grain Models: Filling the Table of Events

A requirement for KMC is to have a table of events with corresponding rates. For macroscopic KMC this table contains processes and input parameters similar to those in rate equation and master equation models. These include sticking fractions, diffusion barriers, binding energies, and activation energies for reaction. The sticking fraction is the fraction of impingements that leads to adsorption. Fractions lower than 1 effectively decrease the flux of incoming particles. The binding energies determine the temperature regime in which the species are on the surface of the grain and are available for reaction. Diffusion barriers are crucial in determining the rate of reaction, since most reactions occur through the Langmuir–Hinshelwood mechanism and should also be given for each species. Since these are intrinsically hard to determine experimentally, often a fixed fraction of the binding energy is taken. Originally H atoms were assumed to diffuse through quantum tunneling. Nowadays models favor thermal hopping after experimental studies that showed that the diffusion of H atoms is rather slow.(53-56) Finally, the activation barrier for reaction should be given for all reactions in the network. These are often guessed from analogous reactions or taken from gas phase data, ignoring the influence of the grain. This ignores the catalytic effect of the grain, however, which can lower the barrier and change the branching ratios of reaction products. Often quantum tunneling is assumed using an approximation with tunneling through a square barrier.(8)

In summary, macroscopic KMC requires at least two input parameters per species (diffusion barrier and binding energy) and one per reaction (reaction barrier). If other types of desorption are considered, such as photodesorption,(17-20) sputtering by cosmic rays, grain heating by cosmic rays,(21, 22) or desorption upon reaction using the excess energy,(2, 57) more input parameters are needed. For microscopic KMC the list is much longer, since the location of each species is known and the exact environment can influence binding energies and barriers.

4.2.1 Sticking Fractions

The sticking fraction SA of a species A to the surface depends on both the grain and gas temperature and is one of the quantities determining the accretion rate of said species on the grain:

(25)where rgrain is the average radius of an interstellar grain (∼0.1 μm) with number density ngrain, and vA is the average gas phase speed

(26)which is determined by the gas phase temperature Tgas, the mass of the species mA, and the Boltzmann constant k. Sticking fractions are often assumed to be equal to unity, which is a good approximation for species other than hydrogen at low gas and grain temperatures.(58-60) They can deviate from 1 if the incoming species cannot convert their momentum into phonons or if a barrier for sticking, typically restricted to chemisorption, exists. Sticking fractions can be determined by molecular dynamics,(58-62) perturbation and effective Hamiltonian theories, close coupling wave packet, and reduced density matrix approaches,(63) or by the much more approximate soft-cube method.(64, 65) Experimental studies on the sticking of H, H2, D, and D2 have been performed.(66, 67) It was found that, at low coverage, the sticking coefficient of H2 increases linearly with the number of deuterium molecules already adsorbed on the surface. Similar trends were found for the hydrogenation of CO, where the KMC simulations could only explain the experimental trends at intermediate temperature if a H2-coverage-dependent sticking fraction of H atoms was assumed.(68)

4.2.2 Thermal Desorption Rates

The residence time of a species A on the surface is predominantly determined by their desorption rates. The thermal desorption rate, in turn, depends on the binding energy Ebind,A

(27)where ν is an attempt frequency. In an off-lattice simulation, one should be able to obtain the binding energy from the potential energy surface; in lattice-gas models this is not possible and the binding energies should be taken from some external source. One of these sources is temperature programmed desorption (TPD) experiments. Since this is such an important technique for the determination of input parameters and since also KMC simulations of TPD experiments have been reported, we will here briefly explain the technique. See ref 69 for a review. These experiments are usually performed in an ultra-high-vacuum setup equipped with a quadrupole mass spectrometer. The temperature of the substrate of interest can be carefully controlled. The experiment consists of two phases: first, the substrate is kept at a constant low temperature and a known quantity of one or multiple species is deposited; during the second phase the temperature is linearly increased and the desorption of surface species is monitored using the mass spectrometer. The measured desorption rate of the species is then plotted as a function of temperature.

One can perform a series of these experiments, varying the initial deposited amount, the deposition temperature, and/or the heating ramp. Often, the resulting TPD spectra are then fitted by the Polanyi–Wigner equation:

(28)where o is the order of the desorption process and β is the heating rate. Generally the order is assumed and the prefactor and binding energy are used as fitting parameters. Unfortunately, these two parameters are correlated and it is not straightforward to get accurate values for both parameters at same time. For this reason, sometimes the prefactor is assumed as well and only the binding energy is obtained from the fit.

Figure 6 shows synthetic TPD spectra using this expression with different desorption orders. The left panel of Figure 6 shows three TPD spectra of zero order desorption with different initial deposited amounts. The leading edges of these three curves coincide. Zero order desorption generally occurs when multiple layers of the same species are deposited. During the desorption, the number of surface species available for desorption (the top layer) remains the same. The middle panel of Figure 6 shows first order desorption. Here, the temperature for which the desorption rate peaks is independent of the deposited amount. First order desorption generally occurs in the monolayer regime, where all surface species are also available for desorption. Please notice that most rate equation, master equation, and macroscopic KMC models assume first order desorption in all cases, even in the multilayer regime. Section 7 discusses some more recent models that move away from this assumption.

Figure

Figure 6. Synthetic TPD spectra using eq 28 for zero (left), first (middle), and second (right) order desorption. Three different initial deposition amounts have been used.

Finally, the right panel of Figure 6 shows second order desorption. Here the trailing edges coincide. Second order desorption is mainly due to two reasons: either desorption of reaction products that were formed in a second order surface reaction, or the surface exhibits a distribution of binding sites. In the latter case, the stronger binding sites are filled first and the average binding energy decreases with coverage, resulting in TPD curves that resemble second order behavior. Second order behavior is also found for the desorption of HD that is formed by reaction of H with D. Katz et al.(56) obtained binding energies and diffusion barriers for H2 and H by fitting a rate equation model to the experimental data where the binding energies and diffusion barriers are used as fitting parameters. A diffusive model was applied to fit the results of Pirronello et al.,(53-55) which could only give a lower limit for the desorption energy of H atoms. Higher values for this energy were found to give degenerate results. In general, one could state that the desorption energies of stable species such as H2 and D2 are relatively straightforward to determine whereas desorption energies of unstable species such as H and D are much harder to determine. For most cold dark cloud simulations, this is however not a problem for most species other than H and D, since desorption of these species hardly occurs because of the low temperature. Unfortunately, binding energies of H and D are especially important in the determination of the temperature window in which efficient hydrogenation can occur. For the diffusion barrier, the situation is the reverse: they are easier to obtain for reactive species than for nonreactive species. Unfortunately, the ratio between the found diffusion barrier and the lower limit of the binding energy of H is often generalized for other species as well. This is however an upper limit.

The binding energies of a wide collection of stable species have been determined using the TPD technique. Examples are N2,(70) CO,(70) O2,(71) H2O,(72, 73) and CH3OH.(74) The binding energies have been mostly determined for the desorption of pure ices from different substrates. The differences between the different substrates are rather small and become negligible in the multilayer regime.(74) The experiments show that the molecules indeed desorb with a (near) zero order rate in the multilayer regime whereas they desorb with a (near) first order rate in the monolayer regime, as explained above.(73)

Since interstellar ices are not homogeneous, the desorption of mixed layers is more relevant for astrochemical modeling. However, the introduction of more species in the ice makes the desorption process immediately much more complex. First, the binding energy of a particle can change depending on its surrounding material. Second, the dominant mantle species can prevent other species from desorbing. Collings et al.(75) showed, for instance, that molecules such as CO and CO2 can become trapped in an ice mantle which consists predominantly of water ice since the desorption of water occurs at much higher temperatures than those for CO and CO2. However, at the long time scales available in the ISM, some of these trapped species might be able to escape because of a segregation process where the two main fractions of the mantle slowly separate.(76) This process depends on a large number of parameters including surface temperature, ice composition, and mixing ratio. The effect of trapping due to layering can be easily simulated using microscopic KMC simulations; segregation is harder to model, since the exact mechanism is not yet understood.(76)

4.2.3 Diffusion Rates

TPD experiments can also be applied to obtain information on diffusion barriers. This information is however indirect, and the result is rather sensitive to the formation mechanism/reaction order that is implicitly assumed in the model used for fitting. For example, Katz et al.,(56) Cazaux and Tielens,(77) and Cuppen and Herbst(78) were able to fit the same experimental results from refs 53−55 by applying different models. Katz et al. had a homogeneous model with one binding site; Cazaux and Tielens applied a binary model with both physisorption and chemisorption binding sites, as well as quantum mechanical tunneling and thermal hopping between sites. Cuppen and Herbst introduced surface roughness in their model.

Recently, TPD experiments were applied to obtain the diffusion barrier of D atoms on porous amorphous water.(79) The barrier of 22 ± 2 meV is in agreement with earlier experimental(80) and theoretical results.(58) The determination of the diffusion barrier remains however experimentally challenging and is one of the largest reasons for uncertainty in the modeling of surface chemistry in general. Experimental information about site-specific barriers is practically nonexistent. Computational efforts to obtain more insight have been reported.(-81-83) A few of these studies applied KMC to obtain diffusion rates, and these will be discussed in section 7.

A requirement for every rate is that it should fulfill macroscopic reversibility. In equilibrium, the master equation should result in a steady state

(29)where k̃(xf|xi) denotes the equilibrium rates, which are equal to the transition probabilities times some unit of time. The easiest way to fulfill this equation is to require detailed balance or microscopic reversibility, i.e., the net flux between every pair of states is zero:

(30)For a canonical ensemble

(31)For KMC models with homogeneous diffusion and desorption rates, this requirement is met naturally. When different types of binding sites are introduced, one should be careful to check whether especially the diffusion rates obey this detailed balance criterion. This requirement is not met for all models with inhomogeneous binding. For instance, the model with the distribution of diffusion barriers in ref 49 or the work by Dumont et al.(84) does not obey detailed balance. For some models it not clear whether they follow microscopic reversibility since the barriers are not mentioned in the paper explicitly.(85, 86) Fortunately, most models do, and in general two types of expressions for the diffusion barriers have been applied. The first expression for the hopping barrier from site i to site f is

(32)where Ehop is an input parameter determining the diffusion barrier on a flat surface, i.e., between sites of equal binding energy, Ebindi is the binding energy at site i, and Ebind0 is some reference binding energy, usually the binding energy of the absorbate without lateral interactions. The second is

(33)Both expressions 32 and 33 hold microscopic reversibility, but the exact choice influences the final results. Figure 7 plots the recombination efficiency for H2 formation. This is the fraction of H atoms that hits a grain and leaves again as H2. The top panel of Figure 7 plots the results as presented in ref 78, where the four surfaces mentioned in the legend have different degrees roughness, starting from completely flat (surface a) to a very rough surface (surface d). Here eq 32 is used to describe the diffusion barrier. The bottom panel of Figure 7 shows the results for the other choice of diffusion barrier. The choice of barrier clearly has an influence on the exact value of the efficiency, but more importantly it influences the exact temperature range in which H2 can be efficiently formed, with eq 33 leading to a wider range.

Figure 7

Figure 7. Molecular hydrogen recombination efficiency as a function of temperature on olivine substrates of different surface roughness. The top panel is reproduced from Figure 6a in ref 78 with permission of Oxford University Press on behalf of the Royal Astronomical Society. It uses eq 32 to calculate the diffusion barriers. Bottom panel uses eq 33.

Expression 32 was introduced by Cuppen and Herbst(78) to introduce hampered movement and desorption because of surface roughness. Equation 33 was later introduced based on a study of hydrogen atom diffusion on graphite(87) (see Figure 8). Here the barrier for diffusion was found to be linearly dependent on the energy difference between the initial state and the final state. For H on graphite, the binding energy of an atom is affected by the specific arrangement of other hydrogen atoms chemisorbed in its vicinity. The same dependence of Ediff and ΔE is found for dimer structures, where two H atoms are involved, and for tetramer structures with four atoms. This very nice correlation is probably an effect of the isotropic nature of the system and is possibly due to the strong chemisorption interactions that are at play. For diffusion of physisorbed CO on a hexagonal water ice surface, no correlation was found as can be seen from the top panel of Figure 9. Again the energy difference between the initial and the final state is plotted versus the diffusion barrier. These barriers are obtained in off-lattice KMC simulations(83) and will be discussed in more detail in section 7. Even less correlation is found for diffusion on an amorphous water substrate (bottom panel of Figure 9 based on unpublished results). Checks for the validity of eq 32 for diffusion on interstellar ice show a similar lack of structure.

Figure 8

Figure 8. Diffusion barrier for chemisorbed H atoms on graphite as a function of the energy difference between the initial and final configurations. The diffusion barriers between dimer and tetramer configurations follow the same linear dependences (eq 33). Figure reproduced with permission from ref 87. Copyright 2008 AIP Publishing LLC.

Figure 9

Figure 9. Similar to Figure 8 for CO diffusion on crystalline (top) and amorphous (bottom) water ice. The black solid line indicates the minimum barrier possible to ensure microscopic reversibility (ΔE). Top panel constructed on the basis of data produced for ref 83.

For both eqs 32 and 33, a choice for Ehop needs to be made. This is often chosen as a fraction of the binding energy, because of lack of information on the exact barrier. This fraction is preferably chosen to be rather high, since a low diffusion barrier slows down the simulations substantially. In that case, the diffusion time scale becomes very small as compared to, for instance, accretion or desorption rates and the simulation spends most of its CPU time on diffusion steps that do not lead to a fast evolution of simulated time. On first instance, one would expect the diffusion rate to play a crucial role in the overall evolution of the systems. However, this depends strongly on the circumstances. For H-atom bombardment of a CO ice under laboratory conditions, the effect was found to be rather small(68) since a barrier exists for the reaction between CO and H. A faster diffusion rate leads to more frequent meetings between H and CO, but per encounter the time to react is shorter with respect to lower diffusion rates. Under interstellar conditions the effect of the diffusion/desorption ratio on the CH3OH formation by CO hydrogenation was found to be somewhat stronger. Chang and Herbst(88) developed a macroscopic gas–microscopic grain Monte Carlo model which treats a moderately complex grain surface chemical network that leads to the formation of the stable molecules CH4, CO2, CH3OH, and H2O. They found that the diffusion ratio has a moderate influence on the abundance of most surface species; CO and CH3OH were most affected. The main difference between their simulations and the laboratory simulations,(68) apart from the very different fluxes, is that CO is a minority species in the ice. The laboratory experiments, on the other hand, start with a pristine CO ice surface and only after some time when H2CO and CH3OH are formed, become CO and H encounters more rare. This difference in CO distribution could explain the difference in sensitivity toward the height of the diffusion barrier between the laboratory and ISM simulations.

But even in simulations of the hydrogenation of pure ices, diffusion rates can have a significant effect on the final yield of many of the species, as is the case for O2 hydrogenation under laboratory conditions.(89) For this particular system, bulk diffusion is crucial in determining the final yield.(90) Simulating bulk diffusion using lattice-gas KMC simulations is however hard, since the exact mechanism is unclear. The discrepancy between simulations and experiments in the cases of both CO(68) and O2(89) hydrogenation is probably mainly due to a missing bulk diffusion mechanism. KMC simulations on segregation of ice mixtures, a direct measure of bulk diffusion, showed that experimental results can only be qualitatively reproduced if both surface diffusion and bulk diffusion by swapping of molecules between sites is included.(76) The exact molecular mechanism behind this process is unclear, though.

4.2.4 Theoretical Determination of Rates

Results from theoretical calculations are also often used as a source of rate constants, especially when experimental data is scarcely available. Binding energies can be determined in a relatively straightforward way using different levels of theory. The density functional theory calculations mentioned in section 5.1 are an example. Wave packet propagation and molecular dynamics techniques can be applied to determine rates or branching ratios. A review of these applications is however beyond the scope of this review. Here, we will focus on reaction path finding methods since these are—in combination with transition state theory (TST)—frequently applied to determine rates for KMC simulations. With the use of TST it is possible to calculate process rates if a suitable dividing surface (the transition state) between the initial state and the final state of the process is known. The identification of this transition state is the crucial step in the procedure and can be handled using a variety of techniques.(91) In the often-used classical harmonic approximation to TST, the transition state is simply a first order saddle point on the potential energy surface. In this case, the nudged elastic band (NEB) method(92) is one of the most applied methods for locating the transition state, if the initial and final states of the reaction are known. The NEB algorithm can determine the minimum energy path between the two minima on the PES. A NEB calculation starts from a sequence of images of the system, initially placed along some pathway between the two minima, usually just along a straight line connecting the two. Every image is then connected to its two neighboring systems by fictitious springs, and the energy of this constrained system is then minimized using any minimization algorithm. Analysis of the coordinates and energies of the highest energy image, after minimization, then gives an estimate of the process energy barrier and the transition state itself. With the improved, climbing image NEB(93) (CI-NEB), algorithm, it is also possible to find the location of the saddle point with arbitrary precision.

Within an astrochemical context, the NEB method has been used, for example, to determine the diffusion rate of atomic hydrogen in a CO ice environment.(68) These rates were then used as input in microscopic KMC simulations. One of the major advantages of using theoretical calculations to obtain rate constants is that the environment can be precisely controlled, which can greatly add to the physical accuracy of the KMC simulations. Fuchs et al.(68) for example found a strong dependence of the barrier for a CO and H to swap on the depth from the surface where the process occurs.

The determination of diffusion barriers is a typical case where NEB calculations can be of use because, as said before, experimental data for these processes is often hard to obtain. An example where this has been employed is the work by Batista and Jónsson.(81) Here diffusion rates for H2O on hexagonal ice were determined using NEB and subsequently included in a lattice-gas KMC model to determine the diffusion constant which is in agreement with the experimental diffusivity at the same temperature. The NEB method can in principle be used for any kind of process however, including chemical reactions. The only prerequisites are that the initial and final states of the process are known and that an adequate theory is available to describe interactions. In this spirit, the studies by Goumans et al.(94, 95) are certainly worth mentioning. They used an embedded cluster method to study the reaction pathways of H2O and H2 formation on olivine surfaces using the CI-NEB method. Also here, the atomistically detailed information about the environment of the reaction revealed important aspects of the reaction pathway.

Summarizing, the NEB method is a powerful theoretical tool which can be of great value when creating a rate table for KMC simulations. Also, when accurate experimental data is available, NEB calculations can often still reveal important aspects of the processes under consideration which can then also be included in the KMC simulation.

4.3 Kinetic Monte Carlo Algorithms

The theory behind KMC can be implemented into an algorithm in several ways. Here, we will discuss a few algorithms that are most commonly used for microscopic simulations in astrochemistry. Figure 10 shows a flow diagram of the KMC algorithm as proposed by Nordmeyer and Zeara(96, 97) which is based on the original paper by Gillespie.(42) This algorithm is mainly developed to study the kinetics of gas adsorption on solid surfaces and it is used to simulate homogeneous surfaces; in other words, the desorption and diffusion rates are site independent and also only one type of species is considered. The algorithm starts by inputing all rates, i.e., creating the table of events. Then, the surface is initialized by creating an empty matrix, M, where each entry represents a site and the occupancy of this site can be changed by changing the number. In the next step the desorption, diffusion, and deposition rate are calculated. Since a homogeneous surface is assumed this corresponds to

(34)

(35)

(36)

(37)where N is the number of species on the grain, Nsites is the number of sites (the size of matrix M), and kx is the rate for an individual atom or molecule to undergo process x. Then one of the three possible events is chosen by picking a random number between 0 and ktot, a second random number is picked to choose the site where the event will occur, and a third one is chosen to increment the simulated time according to eq 23, and if the termination criterion is met, the simulation will finish; otherwise the loop is entered again.

Figure 10

The termination criterion can be in terms of the number of iterations or of a certain final time. In astrochemistry the time steps are often very inhomogeneous: there can be a long time span between two depositions but once, for instance, a hydrogen atom lands on the surface it can make many fast hops before it reacts or desorbs. Setting a termination time is therefore a better criterion than simply the number of iterations, although this is traditionally more commonly used in, for instance, catalysis. The same applies to the sampling criterion, which is checked of the end of every loop in Figure 10. One would like to take some averages, take some snapshots, or write the current surface abundance at regular intervals. These intervals can again be defined in terms of simulation time or number of iterations. When taking averages of, for instance, the average production rate under certain conditions, one should be careful that steady state has been reached and that the time between two samples is long enough for them to be independent.

Upon reentering the simulation loop, the following iterations simply start by recalculating the desorption, diffusion, and deposition rates. For homogeneous surfaces this can be done in a straightforward way and this step does not require much effort. Extending the model to systems with more species would already be more involved and if one would like to treat inhomogeneous surfaces, smart updating routines should be applied to keep the computational cost of this step limited. In that case, however, separating the determination of the overall event and the specific location becomes unnecessarily cumbersome and a different algorithm would be preferred, for instance, the commonly used n-fold way algorithm.(98)

Figure 11 shows a flow diagram based on the n-fold way KMC algorithm.(98) The main difference between this algorithm and the algorithm of Figure 10 is that the event and the location are chosen in one step. This makes this algorithm more appropriate for systems where the surface is inhomogeneous or where the adsorbed species interact with each other, since not all sites are equally likely to have an event. Only two random numbers are needed per iteration in this routine, but the determination of which event will occur is more computationally intense as well as determining the total rate, since the table of events is longer. Updating for the next iteration can be optimized by only recalculating those entries in the table of events that refer to the changed species or, in the case of interacting species, to its neighbors. The total rate can be determined by subtracting the old rates at these entries and adding the new ones. This works well for systems where the rates are rather homogeneous. For systems where certain events can result in large changes of the rates, numerical errors can occur. An example would be a case where the reactivity changes dramatically between sites. If an atom or molecule would hop up and down between this reactive site and a less reactive site, large rates would be constantly added and subtracted, resulting in a numerical error of the same order as the rate of the slow process. In the selection process, the wrong event could then be selected.

Figure 11

Figure 11. Flow diagram of the n-fold way KMC algorithm derived from ref 98 and adjusted to simulate grain surface chemistry.

These numerical issues do not occur in the continuous-time, random-walk algorithm.(50) The first microscopic KMC simulation in astrochemistry used this particular algorithm, and its flow diagram is depicted in Figure 12. This algorithm is designed to describe noninteracting random walkers on a surface. The main difference between the previous two algorithms is that a separate time line for each individual random walker is determined. The time increment is therefore not determined by the total rate for the whole system, but is determined by the total rate for an individual random walker. The event times for each walker are then compared, together with the next deposition time, and the first event is then executed. Only one random number per iteration is needed, two in the case of deposition. In this way when updating the system for the next iteration only the random walker which moved in the previous iteration has to be considered. This random walker will get a new event time and will be placed in the species list which is sorted according to its event time. This algorithm works very well for noninteracting species. Obviously when species interact, a move of one random walker can affect the total rate of a neighboring random walker and, in turn, its event time. A good example would be when a stationary species, which has an event time well above the termination time of the simulation, gets a mobile neighbor with which it can undergo reaction. The event time of the stationary species changes drastically. A way to adjust for this is to store, for each species, the random number Zi that was drawn to determine its initial event time, its total rate, and the time at which its last event occurred, tevent,oldi. The adjusted event time, tevent,adjusti, can then be calculated using

(38)Again, numerical errors can occur in adjusting this event time, if large time differences are subtracted. In this case, it is best to determine the time according to the usual scheme since the first term in eq 38 will probably be small anyway.

Figure 12

Figure 12. Flow diagram of the continuous-time, random-walk KMC algorithm derived from ref 50 and adjusted to simulate grain surface chemistry.

Algorithms that simulate the master equation in an approximate manner have been suggested as well with the aim of speeding up the simulation. One such suggestion is by splitting by physical processes where each of the basic processes—deposition, diffusion, reaction, and desorption—is considered as an independent Markov chain.(99) This algorithm was tested for H2 formation on interstellar grains, and similar results were obtained as in refs 30 and 49. Exact timing differences with other algorithms are not given; it is therefore unclear how large the final speedup will be.

4.4 Simulations with Varying Rates

Temperature programmed desorption is one of the main techniques within solid state laboratory astrophysics. As such it is interesting to simulate these experiments using KMC. There are also advantages of using KMC over the Polanyi–Wigner equation (eq 28) for simulating TPD spectra. First, the order of desorption might change during the experiment, for instance if one goes from multilayer to monolayer desorption. KMC will automatically take this into account. Second, there might be lateral interactions between the adsorbates, either repulsive or attractive, and this might result in nonlinear effects that the Polanyi–Wigner equation cannot account for, but which can be included in a KMC model. Finally, the Polanyi–Wigner equation provides an empirical description of the data and does not give real microscopic insight into the mechanism behind the desorption.

As described in section 4, the two main equations on which KMC algorithms rely are eqs 19 and 21. While these are relatively straightforward for stationary rates, they become more complex in the case of TPD experiments where the changing surface temperature causes the rates to change in time. The good news is however that we know in advance how the temperature will change. For TPD experiments this is usually T = T0 + βt, where T0 and β are constants. Jansen(45) has described a way to account for this time dependence and to obtain values for Δt that follow the correct distribution. Both eqs 19 and 21 have the term

(39)For thermal rates and assuming a linear increase in temperature, this becomes

(40)with

(41)Here E2 is an exponential integral.(100) Since tk–1 is known, it is the current time, and β and T0 are experimental parameters, i.e., the heating ramp and the initial deposition temperature, respectively. A value for tk can be obtained by solving

(42)This effectively means finding the root of

(43)This can be most conveniently done with the Newton–Raphson method.(101) Since Ω and dΩ/dt are both monotonically increasing, the method always succeeds.

The derivation here used a linear increase of the surface temperature. Different time dependences are also possible. The numerical effort depends on the integration of ∫0t ν exp(−E/kT(t)) dt′. One such application could be the thermalization of species upon reaction or deposition. Here the surface temperature does not really change, but one could consider that each species has its own internal “temperature” which starts at some high value and then decreases until it reaches the surface temperature. Cuppen and Hornekaer(87) applied this approach to account for the thermalization of H atoms landing on graphite. A discussion of the results of these simulations is provided in section 5. In the experiments that they aimed to reproduce, the H atoms arrive at the graphite surface at a temperature around 2000 K, which is much higher than the surface temperature. Furthermore, the atoms will gain energy due to the high binding energy if they chemisorb. The atoms will not be thermalized instantaneously, but will most likely gradually lose their energy to the substrate. The dissipation of the excess energy into the substrate is expected to be exponential.(102) When using an exponential function, the expression for Ω in eq 43 becomes rather complex and computationally expensive; the time dependence of the energy loss was therefore emulated by a simpler expression. The following function was used to describe the temperature:

(44)where Ts and Tstart are the surface and starting temperature, respectively, ta is the time at which the atom has adsorbed on the surface, and B is some thermalization parameter. A value for tk can now be obtained using

(45)where “erf” is the error function. Several functions of different order in t or with different functional forms were tested, but this appeared to have a minor effect on the results. For higher orders in t, solving Ω becomes more computationally expensive.

5 Simulations of Experiments

ARTICLE SECTIONS

Jump To

An advantage of KMC over, for instance, molecular dynamics, is the relatively long time scales that can be simulated. Laboratory surface science experiments of several hours can be easily reached with the possibility to include atomistic detail in the simulations. For this reason, KMC is often applied in surface science and catalysis to simulate experimental conditions and determine experimentally measurable quantities. These can then be directly compared to experiments. Section 4 made it very clear that lattice-gas KMC is far from being an ab initio method, in the sense that all input parameters need to be given. This can also be an advantage, however, since the effect of these parameters on the output can be tested and different mechanisms can be switched on and off by changing these parameters. By testing different mechanistic scenarios in this way and comparing the results with experiments, insight into the mechanisms can be obtained. This section will discuss the application of the KMC technique to gain more insight into experimental observations. Since this is a huge area of research with many different applications, we will limit ourselves to the simulation of experiments that are astrochemically relevant. All simulations discussed in this section will use the KMC technique either to verify the experimental mechanism or to parametrize the KMC model using experimental constraints. In section 6, we will discuss KMC simulations under astrochemical conditions: low fluxes and long time scales. The simulations in both sections will generally involve the same systems—H2 formation on different substrates or ice chemistry—but in section 6 the model and mechanism are considered a given and are not tested against experiments; they are rather applied to predict abundances in the ISM.

5.1 H on Graphite

The formation of the most abundant molecule in the universe, H2, occurs mainly through surface reactions, and as such it has received considerable attention. These reactions have been the topic of several astrochemical modeling studies as mentioned in previous sections, but also of experiments and quantum chemical calculations. A model system for H2 formation on carbonaceous grains is the deposition of H atoms on a graphite surface or a graphene layer. Figure 13 shows a schematic representation of the honeycomb lattice of a graphene layer. Hydrogen atoms can be chemically attached to the carbon atoms in this hexagonal lattice. The gray circle represents such an atom. Molecular hydrogen can be formed by a process called “abstraction”, where the atom chemisorbed to the surface reacts with another H atom to form H2. Hydrogen abstraction from graphite has been studied extensively both experimentally(103, 104) and by density functional theory (DFT) and quantum wave packet calculations.(105-110)

Figure 13

Figure 13. Schematic representation of the honeycomb lattice of a graphene layer with a hydrogen atom attached (gray circle). Stable dimer configurations can be formed by attachment of a second hydrogen atom at the o(rtho) or p(ara) position.

Experimental(111, 112) and theoretical(111, 113-116) studies show that hydrogen atoms preferentially adsorb in clusters on graphite. Particularly preferred clusters are so-called ortho and para dimers, where two hydrogen atoms adsorb to carbon atoms that are within the same hexagonal ring in the graphite structure at the ortho and para positions, indicated in Figure 13 by an “o” and a “p”, respectively. DFT calculations show that not only the binding energy of these dimers is higher than for two individual monomers, but also that the sticking into dimer configurations is more favorable. More complex structures consisting of three or more hydrogen atoms are found to be favorable as well.(115, 117) A recent study has applied ab initio molecular dynamics to study the competition between the Eley–Rideal abstraction reaction and dimer formation.(118)

Performing these types of competition calculations at a high level of theory is however computationally demanding and could, therefore, only recently be investigated. Typically, either no C atoms or only a single C atom on the simulated graphite surface is allowed to relax during the reaction, which can artificially close other low barrier reaction channels. Another approach is to apply kinetic Monte Carlo models to simulate the experimental measurements and to connect the atomistic results from the quantum calculations to the more macroscopic experimental observations. Cuppen and Hornekaer have used such a hybrid approach by using binding energies and barriers derived from DFT and wave packet calculations as input parameters for Monte Carlo simulation and allowing different mechanisms for H2 formation. Independent of the exact values of the input parameters, they found that certain mechanisms could better explain the experiments than others. Fast diffusing physisorbed atoms(119) could explain the high percentage of H atoms present on the surface as part of a cluster (dimer, trimer, etc.).(115) The fast-moving physisorbed atoms can scan a large part of the surface, even considering their low binding energy, and find chemisorbed atoms to form a cluster. Some excess energy will probably be released in the binding process of this atom. The experimental data of Zecho et al.(103) could be explained by an abstraction mechanism in which part of this excess energy is used to overcome the barrier for abstraction.

The combination of these mechanisms together with the parameter settings that reproduced the experiments best were subsequently used by Gavardi et al.(120) to study the desorption behavior of H2 from graphite, testing the KMC simulations against experimental temperature programmed desorption curves.(104) Again the various clustered hydrogen surface structures were taken into account, using binding energies based on DFT calculations and accompanying estimations. Experimentally two desorption peaks are observed as well as a shoulder on the trailing edge of the first peak. The top panel of Figure 14 plots the simulated TPD spectrum at a coverage of 0.1 monolayer (ML) and shows different contributions to the TPD spectrum. This graph is a different representation of the data from ref 120 to allow for better comparison with ref 84. The main graph shows the coverage loss of the different configurations in time, whereas the inset shows the desorption of H2 and the configurations that are directly responsible for this. From the simulations it becomes apparent that the first peak arises from H2 desorption from the para dimer or larger clusters containing a para configuration as a substructure. The second peak was attributed to the two-step diffusion of hydrogen atoms in an ortho dimer into a para dimer. This para dimer can then desorb again in the form of H2, resulting in the second peak. More complex structures were found to be responsible for the third peak, experimentally observed at high coverage.(120)

Figure 14

Figure 14. KMC TPD simulations of coverage loss rate from graphite as determined by Gavardi et al.(120) (top) and Dumont et al.(84) (bottom). The inset in the top panel indicates the H2 desorption rate and the different contributions. The bottom panel is reproduced with permission from ref 84. Copyright 2008 American Physical Society.

Gavardi et al.(120) also recognized the possible existence of a directional bias for the diffusion of physisorbed H atoms. Due to their fast mobility, the Markov chain assumption of memory loss is no longer valid. It is in fact more likely that an atom proceeds in the same direction for some time. Therefore, the six possible directions of diffusion on the honeycomb structure were assigned different weights, with the highest probability for the atom to continue in its original direction. This directional bias was found to have very little effect on the final results. Only when no deviations from a straight line were allowed was a lower H2 formation rate obtained.

Dumont et al.(84) performed a similar study in which they simulated the TPD experiment by Zecho et al.(104) using a kinetic Monte Carlo algorithm. It is not clear how they accounted for the temperature dependence of the rates. Their explanation for the two TPD peaks agrees with that of Gavardi et al., despite some differences in implementation of the H–graphite system in terms of algorithm. Dumont et al. allowed absorption, desorption, and dimerization, but they did not include the formation of clusters larger than two atoms and diffusion of physisorbed H atoms on the surface. The latter was found to be a dominant mechanism for the formation of clusters in the study by Cuppen and Hornekaer. Consequently, the surface snapshots presented by the Dumont group show only dimer configurations. Their results are reproduced in the bottom panel of Figure 14. The first peak is again attributed to para-dimer desorption, whereas the second peak is due to the disappearance of the ortho dimer. The latter can proceed through different routes: through the evaporation of individual H atoms in the dimer and through hydrogen diffusion into different dimer structures.

5.2 Ice Chemistry

In dark regions of the interstellar medium ices have been abundantly observed. The name “ices” in astrochemistry is not reserved only for water ice, but is used for all solids that consist of volatiles (molecules that are liquid or gaseous at room temperature). Since frozen molecules are not isolated when they react with each other, gas phase reaction dynamics cannot be directly extended to the solid phase, although the initially proposed surface reaction networks were largely constructed based on analogies with gas phase data.(7) In the past few decades, coinciding with improvements in (ultra) high vacuum techniques, more and more ice processes have been studied experimentally, to test and improve the postulated grain surface reaction networks. These include both energetic processing initiated by Hagen et al.(121) and lower-energetic processes such as surface formation of interstellar relevant molecules by atomic addition pioneered by Hiraoka et al.(122) To translate the experimental results into reaction rates or barriers, however, is far from trivial. The formation of new species is the result of an interplay between diffusion, desorption, and reaction. This is where KMC comes into play as it can explicitly take the microscopic structure into account and reactions and diffusion can be evaluated in time.

In the field of experimental ice chemistry, a large variety of processes can be explored. Here we discuss only thermal processing and atom addition reactions since they have resulted in follow-up studies involving KMC. Thermal processing consists of probing a reaction and/or diffusive behavior that is temperature driven without other energetic input. Experimentally, there are two types of atomic addition experiments: sequential and codeposition. For the former type, the ice under investigation is grown to obtain the desired number of monolayers (ML). Subsequently, the ice can be bombarded with atoms. This is beneficial for studying stable final products. For a codeposition experiment, however, all species are allowed to freeze-out simultaneously and can therefore be isolated in the ice. Hence, it is possible to inspect also radicals and intermediate species.

Let us first consider the thermally induced process of segregation in mixed H2O and CO2 ices, where the key parameters are bulk and surface diffusion and the difference in binding energy between a water-rich phase and a carbon dioxide rich phase. The latter drives the segregation process; the first determines the rate by which this occurs. Astronomical observations show the presence of pure CO2 ice in protostellar envelopes, whereas CO2 is predominantly mixed with H2O ice in the precursor clouds.(123) This has been attributed to the higher temperatures in protostellar envelopes that would allow segregation of H2O and CO2 ice mixtures to occur. In order to be able to use the fraction of pure CO2 ice as a probe for the thermal history toward a protostar, mechanistic information about how this process occurs as well as quantitative information on the barriers involved are necessary. This motivated a combined experimental and computational study, where three sets of experiments concerning the segregation in H2O:CO2 and H2O:CO mixtures under different vacuum conditions were complemented by KMC simulations.(76) The experimental data could be fitted with two exponential functions, suggesting the existence of two first order mechanisms: a fast one and a slow one. Subsequently, KMC was applied to test different segregation mechanisms. The lattice-gas model included desorption, hopping, and swapping. The interactions between molecular species of similar kinds were set to be higher than the interactions between H2O and CO2, partially corresponding to values found in TPD experiments. Hopping could occur from site to site, and the rate followed eq 33. Finally, swapping is a process in which molecules exchange position, accounting for the bulk diffusion. Also, for this, eq 33 was applied but with a barrier Eswap instead of Ehop, much higher (factor ∼3–7) than for surface diffusion. It was found that hopping is responsible for segregation at the surface, which is described by an exponential function and reaches steady state. Swapping results in some surface segregation which is slow and does not reach steady state (at the simulated time scales). When both mechanisms are included, the experimental temperature dependence of the segregation rate is correctly reproduced as well as the fact that, at higher temperatures, the bulk segregation seemed to play a larger role. In a following paper, the role of segregation was further investigated through TPD studies of the desorption of mixed ices, where segregation determines the amount of more volatile species to be trapped. A three-phase (gas–ice mantle–bulk ice) rate equation model was constructed with which good experimental agreement was obtained.(124)

Vasyunin and Herbst(125) used these TPD experiments to calibrate their multilayer macroscopic KMC model MONACO. In this model only the outermost layers are chemical active; i.e., chemical reactions and desorption can only occur from these layers. The exact number of layers where these processes can occur is determined by comparing simulations and the TPD experiments in ref 124. They found that MONACO treats zero and first order desorption correctly and further concluded that the experimental curves are well reproduced assuming four monolayers to be chemically active. This number is probably larger than in reality since the authors did not allow bulk diffusion to occur.

Examples of surface reaction experiments that are modeled using KMC are the hydrogenation of solid CO

(46)and the hydrogenation of molecular oxygen

(47)The CO hydrogenation scheme was the first ice reaction route to be studied extensively by a number of groups(68, 126, 127) and leads to the formation of methanol, which is an important precursor molecule for more complex organics.(128, 129) The latter route has also been studied extensively experimentally,(79, 90, 130-133) but the exact formation route was found to be more involved than the simple three-step process mentioned here.(133)

For the KMC modeling of the experimental CO hydrogenation, temperature-dependent reaction rates for the first and third reaction steps were used as fitting parameters to fit the experiments.(68) Quantum chemical calculations have shown that these two reaction steps possess a barrier for reaction(134) and that these barriers can probably be crossed through quantum tunneling.(135) The temperature dependence of the rates found in the KMC study was indeed in agreement with this. Experimentally, surface abundances of CO, H2CO, and CH3OH were followed in time by means of RAIRS in a sequential experiment: CO was first deposited and subsequently hydrogenated for 3 h. Good agreement was obtained between KMC simulations and experiments for low surface temperatures where only the top few monolayers of the deposited CO ice were chemically altered as can be seen in the top left panel of Figure 15. At higher temperatures, the hydrogen atoms were found to be able to penetrate the CO ice deeper, but this penetration effect was not fully reproduced by the KMC simulations, which resulted in a lower final yield in the simulations. The initial reaction rates of the pristine CO ice, when penetration does not play a role yet, was however in good agreement with experiment.

Figure 15

Figure 15. Time evolution of the surface abundance of CO, H2CO, and CH3OH during H-atom bombardment of a predeposited CO ice with at four different surface temperatures. Experimental data (symbols) and Monte Carlo simulation results (solid lines) are shown. Reproduced from ref 68. Credit: Fuchs, Cuppen, Ioppolo, Romanzin, Bisschop, Andersson, van Dishoeck, Linnartz. Astron. Astrophys. 2009, 505, 629 reproduced with permission © ESO.

Considering simulations of water formation under experimental conditions, Lamberts et al.(89) used the reaction network proposed in ref 133 to reproduce the sequential and codeposition experiments reported in refs 90 and 133. Since this reaction network contains many more reactions than in the case of CO hydrogenation, constraining the rates of these reactions proved much harder. They therefore performed a sensitivity analysis to understand which reactions had the largest effect on the experimental results and how well these reactions could be constrained. A reasonable agreement between KMC simulations and the different types of experiments was obtained. Experimentally it was found that hydrogen atoms can penetrate into tens of monolayers of O2, unlike the case of CO where only a maximum of four monolayers is affected. Again this penetration mechanism proved hard to reproduce in the KMC simulations. Only a fast H diffusion combined with access to interstitial sites was capable of reproducing some kind of penetration. The agreement with the codeposition experiments where O2 is constantly deposited was much better. For astrochemical purposes, these codeposition experiments are in closer agreement with the conditions in the interstellar medium where all species land on the grains simultaneously and penetration effects are of less importance. However, the diffusion processes covering the penetration of reactive species into the ice can be very similar to the ones that drive segregation of ices, discussed earlier, and as such they are interesting to study.

One aspect of diffusion mechanism that is missing in both atomic addition studies, mentioned above, is the directionality of the movement of hot species upon formation. If hot reaction products are formed in the ice with some momentum, the Markov chain assumption breaks down, which results in an effective directionality of diffusion. Gavardi et al.(120) discussed something similar for the diffusion of physisorbed H atoms on a graphite surface as mentioned in section 5.1. An example for ice chemistry is the formation of two OH radicals from H + HO2. The two fragments are expected to move away from each other, and therefore, recombination is less likely. We refer to ref 89 for a more in-depth discussion. Currently, this directionality is not included in the models. All ice KMC simulations, reactive and nonreactive, have shown that bulk diffusion is important in reproducing the experimental results, but that this is still an ill-understood process. Whether directionality and swapping are good handles on this complicated process is yet to be determined by future theoretical and experimental studies.

6 Applications of the KMC Technique to Astrochemistry

ARTICLE SECTIONS

Jump To

Section 5 discussed KMC simulations which aim at modeling (experimentally observed) behavior in order to obtain relevant parameters, e.g., activation barriers. In this section, kinetic Monte Carlo simulations under astrochemically relevant conditions will be passed in review. Both H2 formation in the monolayer regime and ice mantle growth and accompanying chemistry will be discussed. Typical astrochemical environments, e.g., dark or diffuse clouds, are simulated to predict or reproduce molecular abundances.

6.1 H2 Formation

Hydrogen is the most abundant element and molecular hydrogen the most abundant molecule in the interstellar medium (ISM). H2 is largely responsible for shielding of the denser regions from the interstellar radiation field, thus allowing more complex molecules to survive. It furthermore plays an important role in the cooling of the interstellar medium. Molecular hydrogen is clearly a critical molecule in astrochemistry, and its formation routes have therefore been intensely studied. Its most obvious gas phase formation route, H + H → H2, is not efficient under interstellar conditions, since there is no third body to act as a heat sink and radiative association should occur through a forbidden transition.(136) Because of the high observed abundances of H2, formation must be driven by another process, such as formation on grain surfaces.(4) H2 is assumed to be formed through the Langmuir–Hinshelwood mechanism, where two physisorbed hydrogen atoms scan the surface through thermal hopping or tunneling. Molecular hydrogen can be formed barrierlessly when the two atoms reach the same site. Commonly it is assumed that the new hydrogen molecule immediately desorbs back into the gas phase. Although this mechanism was already postulated in 1949, laboratory experiments were only performed at the end of the 20th century.(53-55) Although the mechanism works at 10 K and lower temperatures, a model based on these experimental results showed that the formation efficiency was too low to account for the observed abundances at typical grain temperatures of roughly 20 K in diffuse clouds.(56) Most of the H2 formation studies are devoted to understanding this discrepancy. Inhomogeneity in terms of binding sites was introduced to increase the temperature range in which H2 can be formed toward higher temperatures. This was initially done by introducing chemisorption sites using a rate equation model,(77) but since microscopic KMC is a more obvious way of introducing site-specific information, later mostly this technique was used.

6.1.1 Effect of Inhomogeneity

The absence of substructure in the observed 9.7 μm band leads to the belief that interstellar silicates are amorphous rather than crystalline (see ref 44 and references therein). As such, it is unlikely that the adsorbate behavior on these grains can be modeled by a single diffusion barrier and binding energy. Inhomogeneous surfaces can be modeled in different ways. Here we will discuss four implementations of surface inhomogeneity. The first is due to amorphicity. Amorphous surfaces have binding energies and barriers that can be thought of as following a (yet unknown) distribution, with slightly changing barriers depending on the small differences in local environment due to the lack of long-range order. Chang et al.(49) studied the Langmuir–Hinshelwood mechanism for physisorbed hydrogen atoms on (i) flat, homogeneous carbonaceous and olivine surfaces with a single binding energy and diffusion barrier, (ii) “inhomogeneous surfaces” of amorphous carbon and olivine in which binding energies and diffusion barriers were taken as continuous distributions (Gaussian and exponential), and (iii) “mixed” grains, consisting of both amorphous carbon and olivine. In agreement with earlier studies,(30, 56) perfectly flat homogeneous surfaces (both olivine and carbon) result in narrow temperature ranges within which H2 formation is efficient. These temperatures were too low to explain the observed H2 formation rate by Jura.(137) Molecular hydrogen was found to form at 10 K and below, whereas the dust in diffuse clouds, where H2 formation is observed, is closer to 20 K. For the inhomogeneous systems the sites with higher binding energy act as sinks in which H atoms can reside up to higher temperatures, thus increasing the surface abundance of atomic hydrogen. The net result is an increase in reaction efficiency. Especially amorphous carbon appears a good candidate for surface H2 formation where recombination efficiencies of at least 0.1 were obtained for temperatures as high as 21.8 K, high enough to cover the grain temperatures of interstellar dust in diffuse clouds.

Since amorphous carbon is not the main component of interstellar dust, it is useful to consider mixed grains consisting of both olivine and amorphous carbon (<40%). In all mixed cases the temperature window for adequate recombination is expanded. Forty percent of clustered, inhomogeneous carbon was found to give an upper bound temperature of 23.8 K.

Inhomogeneities can also occur due to surface roughness. Even for crystalline systems it is well-known that surfaces are usually not flat but exhibit islands with steps and kinks. Cuppen and Herbst studied this effect by using surfaces of different roughnesses.(78) In this case, one can distinguish between vertical interactions, Ev, of a deposited hydrogen atom to the substrate and lateral bonds, El, in case the atom is located near a protrusion. The total binding energy, EB, can then be expressed as

(48)where i is the number of lateral bonds. The top panel of Figure 7 shows the results for four different surfaces with increasing roughness. The surfaces with high surface roughness have more sites with higher values of i and thus more H atoms stay adsorbed which can react with newly deposited species, leading to an extended temperature range in which molecular hydrogen can be efficiently formed. For these simulations, the Eley–Rideal mechanism is included as well; this is especially important at low temperature where the surface coverage is high and the atom mobility is low. The lateral interactions in these simulations are chosen to be as equally strong as the vertical interaction. In reality, these interactions are probably weaker. Since the temperature window for efficient hydrogen recombination scales with El, the assumption for the exact value of El is of key importance. For olivine surface d, still an efficiency of 50% could be obtained for El = 0.4ED at grain temperatures representative for diffuse clouds.

A fourth way to model inhomogeneity is to introduce a second type of binding site with higher binding energy than standard, instead of using a continuous distribution or explicit surface characteristics. Wolff et al. constructed such a system with two binding sites.(138) If the energy difference between the two types of binding sites is chosen to be large, two temperature regimes exist in which H2 is efficiently formed. The authors were especially interested in circumstances where the presence of the deep sites extends the temperature regime belonging to the shallow sides until the two regimes are bridged. They obtained a rate equation model which describes their KMC results with a high accuracy. A strong effect of the formation efficiency on the arrangement of the deep sites was found. The border length which separates the deep and shallow binding sites appeared to be crucial in determining the final H2 formation rate. Surfaces with randomly placed deep sites will have a longer border length than surfaces with one island of deep sites, which results in more efficient H2 formation.

Cuppen and Garrod(139) constructed a similar KMC model with two binding sites. The aim was to develop a method that can grasp some of the surface inhomogeneity that is most likely present on interstellar grains but that is, at the same time, simple enough to describe by modified rate equations that can then be applied to simulate large grain chemistry networks over large time scales and a wide range of different circumstances. A different choice for the diffusion barrier was used—eq 33 instead of eq 32 used in ref 138 —and the focus lay on a different regime in parameter space—low surface coverage, where standard rate equations do not apply. Wolff et al.(138) considered the regime with high efficiencies and surface coverage. As mentioned earlier, traditional rate equation models incorporating the surface formation of H2 failed to take into account the accretion limit problem, first ventilated by Tielens and Hagen.(7) The size limit for this to occur was found to be between 100 and 10 000 sites per grain at temperatures from 10 to 22 K. The exact limit depends on the composition of the grain and the amount of surface inhomogeneity, likely due to the diversity of surface diffusion and binding energies.(49) Note that generally the lower bound of the grain-size distribution is considered to be 0.005 μm,(140, 141) which corresponds to roughly 3000 sites.

Cuppen & Garrod found that the H2 formation efficiency increases with increasing binding energy, fraction of the second binding site, and decreasing temperature. They furthermore managed to construct a modified rate equation method that was in reasonable agreement with the KMC results, if the second type of binding sites were randomly distributed. For distributions other than random (clusters with different aspect ratios, clusters of different sizes, or dispersed in lines), the agreement was not obtained and the results varied by as much as 4 orders of magnitude. As stated before, surface inhomogeneity is most likely due to the amorphous character of the grain and the irregular surface topology of the grain. Since in both cases the strong binding sites will be randomly distributed across the grain, the authors recommend that the modified rate equation model with the introduction of the second type of binding site can be used.

6.1.2 Models Including Chemisorption Sites

The existence of chemisorption wells at the surface can also be viewed as inhomogeneity. Depending on the surface, there can be either a barrier into this well, or not. Cazaux and Tielens(77) showed that this type of surface inhomogeneity can play a role in the formation of H2 and that the presence of chemisorption sites can extend the temperature range in which reaction occurs to quite high temperatures, depending on the height and width of the barrier into the well. The study was initially performed using only rate equations, but in follow-up studies, KMC simulations were applied to check the validity of the rate equation model. In ref 143 the role of carbon grains in the deuteration of H2 was studied using both rate equations and KMC simulations. A model with a high barrier for chemisorption was applied to simulate polycyclic aromatic hydrocarbons (PAHs), which represent the small grains in the grain size distribution. An amorphous carbon model with no barrier to chemisorption was used to describe large grains, following experimental results on the hydrogenation of nanosized carbon grains.(144) Small grains were found to play an important role for dust temperatures below 25 K and large grains the same for the high temperature range. An empirical formula was obtained to describe the efficiency as a function of grain size and D/H ratio. Since the interstellar D/H ratio is roughly 10–5, it is very hard to obtain enough statistics on the formation of HD and D2, since most of the simulation time is spend on H deposition and H2 formation. For this reason the D/H ratio is increased by a factor of 100 or 1000 and the results are then extrapolated.

In a follow-up study,(145) the model was used to study H2 and HD formation at high red shift where the dust abundance is low and gas phase routes compete with the surface formation of these molecules. Three different grain types were considered: amorphous carbon, silicates, and graphite to represent PAHs. The latter model included preferred sticking at para sites as discussed in section 5.1, although a square lattice was used instead of an hexagonal grid. Rate equation results were compared with microscopic KMC results. A grain size dependence of the formation efficiencies was only found for large D/H, which is not so astrochemically relevant. For this reason, rate equations were used in the remainder of the paper to construct fitted expressions that can be used in gas phase models and give the H2 and HD formation rate as a function of red shift, gas temperature, and dust temperature. Even for high red shift where little dust is present, an enhancement of H2 was obtained by surface chemistry with respect to models without surface chemistry. HD was predominantly formed through gas phase route D+ + H2, and therefore the enhancement of H2 led indirectly to an enhancement of HD as well.

The model by Cazaux and Tielens was further extended to olivine by Iqbal et al., who used chemisorption sites for this type of grain as well.(146) Again, both rate equations and KMC were used to describe this system. They performed simulations with different heights and widths for the physisorption to chemisorption barrier for both olivine and carbonaceous material and obtained efficient H2 formation over a wide temperature range. They further found that the rate equations overestimate the KMC results over the entire investigated regime, since for all cases the number of atoms on the surface is less than one. They concluded that coverage can be misleading as an indicator for efficient H2 formation, because low coverage can also mean that all H atoms on the surface immediately convert to H2, resulting in a very high efficiency.

The hydrogen abstraction model described in section 5.1 using a honeycomb graphite structure(87, 120) was used to study the onset of H2 re-formation in postshock regions.(142) During J-type shocks the velocities and temperatures are high enough to dissociate molecular hydrogen and cooling takes places through, subsequently, H, O, and H2 as soon as it is re-formed. This re-formation takes place on the surfaces of dust grains, but the exact conditions under which this occurs are important components in the interpretation of observational data. Observations of O and H2 emission are complementary, and they can be translated into physical conditions using accurate predictions and models. One of the ingredients in such models is the H2 re-formation rate. Using KMC, the product of the sticking efficiency and the H2 formation efficiency (Sη) was determined on grains at a constant temperature of 15 K. The sticking efficiency was first calculated for bare grains, incorporating the effect of chemisorption, leading to the enhancement of Sη at higher temperatures (>500 K) with respect to the original model. Running a shock model including these new results leads to H2 re-formation at earlier times since chemisorption is efficient at high temperatures, i.e., on short time scales after the shock. The higher rate has a substantial effect on the predicted line fluxes of O, and a better agreement was found between the model and observations.

6.1.3 Temperature Fluctuations

Most studies consider a grain at constant temperature, which is a good approximation for large grains and in dark regions. However, since the grain size distribution follows a power law(140)

(49)small grains are responsible for a considerable fraction of the total grain surface area. These grains can be easily heated by UV photons or by cosmic rays. Like grain surface chemistry, this grain heating is a stochastic process and can therefore be perfectly studied by means of KMC simulations. Cuppen et al.(43) heating studied the effect of photon heating on H2 formation in this way, whereas Herbst and Cuppen(22) investigated the effect of cosmic ray heating on the desorption of interstellar grain mantles of different compositions. Since the latter considers ices, it will be discussed in section 6.2. UV photons hit the grains in diffuse clouds with a frequency similar to that of H atoms. If these hits lead to heating, they have a strong effect on the surface chemistry.

When the temperature of the grain fluctuates, the rates become time dependent as in the case of the TPD experiments. Section 4 discussed an elegant algorithm for such a case, using the fact that the time dependence of the rate is known. Unfortunately, this information is not known in the case of stochastic fluctuations. Now, species can land on the grain when it is in its cold phase, at which time the event time is determined using eq 23. This time can easily be longer than the inverse frequency with which the grain is heated. When the grain is heated by cosmic rays or photons, these times change from tevent,old to tevent,new according to

(50)where t is the current time and Told and Tnew are the temperatures before and after heating, respectively. During the cooling the times are changed, using the same equation. Cooling occurs through radiative cooling applying a Debye model for the heat capacity and empirical emission coefficients.

Stochastic heating was found to have the strongest effect for smaller grains: they have fewer incoming photons per time interval and the temperature increase per hit is larger since it is distributed over less volume. The temperature fluctuations can be best described in terms of modal temperature, the temperature at which the grain is most of the time, and the temperature range. For photon heating, the modal temperature is low for small grains and larger values of AV and peaks around r = 0.02 μm and AV = 0.(43) The variation of temperature is largest for the smallest grains. Adding H2 formation to the equation shows that the temperature fluctuations decrease the residence times of H atoms and thus decrease the efficiency. The efficiency is substantially less than the 100% formation efficiency that is often assumed to be required to account for the observational H2 formation rate found by Jura.(137) This assumption however is based on a very narrow grain size distribution of 0.1 μm. Integration of the recombination efficiency over the full size distribution shows a less dramatic effect, since the small grains have a considerably larger surface area. It was found that Jura’s rate could be reproduced as long as grain surfaces of at least moderate roughness were used.(43)

6.2 Ice Chemistry

The astrochemical simulations presented in section 6.1 are relevant for environments where there is either too much radiation for ices to build up or it is too warm for ices to sustain. Dense molecular clouds, on the other hand, provide a surrounding with low temperature, high density, high H2 abundance, and typically high values of AV. The latter is especially ideal for the formation of icy mantles reaching thicknesses up to 100 ML. The high extinction value shields the inner regions of the cloud from UV light, and therefore, atom addition reactions are the prevailing process, together with cosmic ray induced mechanisms. Typical values for diffuse and dense clouds are summarized in Table 1.

Table 1. Typical Physical Conditions for Interstellar Clouds

AV (mag)Tgrain (K)Tgas (K)nH (cm–3)H/H2diffuse ≤1 15–20 60–100 ∼102 H

dense ≥5 10–15 10–30 ∼104 H2

The main problem when simulating grain surface chemistry in these dense conditions is that one leaves the submonolayer regime and ice layers build up. The grain is therefore constantly changing, and the binding energies and hopping barriers change consequently. In diffuse clouds, the desorption and diffusion behavior of the hydrogen atoms is mainly determined by the underlying grain and the interaction between absorbates occurs through reaction. Now, in dense clouds newly formed species determine the binding energy, barriers, and reaction possibilities of their neighbors. Many more interaction energies are thus needed to describe the system, because there is interaction not only between the different species with the bare substrate, but also with all other adsorbed species. The problem is that many of these interaction energies are unknown or at best poorly constrained. We are aware of 12 KMC studies of interstellar ices, treating the ice at least in terms of layers, but often with more microscopic detail. A short summary of these studies is given in Table 2. We will discuss them separately in the remainder of this section.

Table 2. Overview of KMC Ice Studies with Some Microscopic Detail

refno. of reactionsareaction routesgas phase?bD/H?photoprocesses?commentsHerbst and Cuppen(22) 0 – N N N CR desorption

Chang et al.(46) 12 CH3OH, H2O, CO2 Y N N

Cuppen and Herbst(147) 10 H2O2, H2O, O3 N N Y

Muralidharan et al.(148) 0 – N N N water chemisorption model

Fuchs et al.(68) 5 CH3OH N N N sequential accretion of CO and H

Cuppen et al.(48) 5 (12) CH3OH, (CO2, H2O) N N N one simulation with larger network

Cazaux et al.(149) 101 H2O2, H2O, O3 N Y Y

Das et al.(85) 10 CH3OH, H2O, CO2 N N N

Das and Chakrabarti(86) 11 CH3OH, H2O, CO2 N N Y

Marseille and Cazaux(150) 0 – I N N water accretion time scale

Chang and Herbst(88) 29 CH3OH, H2O, CO2, CH4 Y N Y only photodesorption is included

Vasyunin and Herbst(125) full gas–grain network Y N Y only layers

Meijerink et al.(151) 101 H2O2, H2O, O3 I Y Y included in XDR model

Photodissociation reactions are not included in this count.

“I” stands for “indirectly” if the results are later fed to another simulation program treating the gas phase.

The first paper to simulate interstellar ices using microscopic KMC studied the desorption of these ices by cosmic ray heating of grains.(22) In dense clouds, the UV photon flux is substantially reduced with respect to diffuse clouds and cosmic ray heating becomes an important desorption mechanism. Since especially the Fe cosmic rays can lead to a temperature increase of several hundreds of kelvins, cosmic ray hits will be able to lead to substantial desorption of grain mantles. The applied heating model is similar to the one reported in ref 43, but the deposited energy per hit is different and the grains not only cool radiatively but also cool by desorption, which extracts energy from the grain. Every desorption the grain loses Ebind in energy. N2, CO, and H2O ices of different thicknesses and underlying substrates were studied. A single monolayer of water was found to efficiently desorb for all but the largest grains studied as also did the monolayers of CO and N2. Desorption of five monolayers of water ice was found to proceed because of the strong bonds between water molecules.

One of the first papers to use microscopic KMC to study ice buildup was by Cuppen and Herbst.(147) The paper concentrated on the formation of water ice, and it aimed to understand the threshold value of infrared water ice detection at AV = 3.3 ± 0.1 mag.(152) A range of different conditions was simulated. Density, grain and gas temperatures, and visual extinction were varied, representative for diffuse, translucent, and dense interstellar clouds. Water ice was formed under all circumstances, but in diffuse clouds the amount stayed within the submonolayer regime, well below the detection limit. The threshold value for detection was reproduced, taking into account this detection limit. In translucent clouds, the ice was heavily processed by UV radiation and the abundances of oxygen-rich species, such as O2 and O3, were found to increase at later times due to the loss of hydrogen by photodesorption. The route to water formation also changed depending on the conditions—atomic or molecular gas. In diffuse clouds OH + H → H2O was dominant, whereas in dense clouds OH + H2 → H2O + H and also H2O2 + H → H2O + OH became more important. Hydrogen peroxide in the latter route is formed through O2 hydrogenation, and the finding that this route can account for 20% of the water production initiated the experimental study of this formation route.(79, 130, 131)

Another early model focusing on ice formation was by Chang et al.(46) They used a coupled gas–grain code where the gas was described by rate equations and the grain chemistry by microscopic Monte Carlo simulations to simulate the evolution of a grain. Since the main aim of the paper was to show the coupling between the gas and grain chemistries, only a small network of 12 reactions was employed, leading to the formation of CH3OH, H2O, and CO2. The results were compared to a master equation model. Rough and flat surfaces were used at T = 10 and 15 K, and lateral interactions were switched either on or off. A reasonable agreement with ice observations was obtained, better than for the master equation model, and this agreement improved with the lateral interactions switched on. The formation of methanol was found to be less efficient at 15 K than at 10 K, in agreement with experiments,(68) because of a competition between thermal hopping and tunneling for reaction. CO2 was found to be formed through the reaction CO + OH with the OH formed at the site next to the CO. Starting from a rough surface, the grain was found to be slowly smoothed during the buildup of the ice layer. This is caused by the lateral interactions, which lead to filling the strong binding sites, i.e., the valleys and holes.

Section 5.2 discussed a CO hydrogenation KMC model which was compared against experiments.(68) In this paper, also a few simulations on CO hydrogenation under interstellar conditions were performed. In these simulations, a formed CO ice was exposed to hydrogen atoms. The simulations showed that the abundances do not scale with flux, going from the experimental flux to interstellar fluxes over all these orders of magnitude. This can be seen in Figure 16, where the simulated experimental and interstellar ice evolution is plotted as a function of H atom fluence. One reason for this discrepancy is the high surface H abundance in the experimental conditions, leading to more H + H reactions as compared to the CO hydrogenation channel. A second effect is the difference in sticking probably: for hydrogen atoms this is rather temperature dependent and it changes from roughly 0.1 for room temperature H atoms (laboratory) to 1 for molecular cloud conditions.(61) This example shows the need for applying models to bridge the gap between experimental and interstellar conditions.

Figure 16

Figure 16. KMC simulations of CO ice hydrogenation at 12.0 (top) and 16.5 K (bottom) as a function of fluence. The thick lines are at a constant atomic hydrogen gas phase density of 10 cm3 and a gas temperature of 20 K. The thin lines represent experimental conditions with a particle flux of 5 × 1013 atoms cm–2 at room temperature. The results are shown in terms of column density. Figure adapted from ref 68. Credit: Fuchs, Cuppen, Ioppolo, Romanzin, Bisschop, Andersson, van Dishoeck, Linnartz. Astron. Astrophys. 2009, 505, 629 adapted with permission © ESO.

The final H2CO/CH3OH ratio that was obtained in the interstellar simulations showed only a reasonable agreement with observations. The authors argued that a proper model of interstellar CH3OH formation should have simultaneous deposition of CO and H. These types of simulations were performed in a follow-up paper where the model was used to simulate the surface formation of H2CO and CH3OH under different conditions, including density and temperature.(48) The applied H flux was constant in time, whereas the CO flux gradually diminished because of CO freeze-out. The obtained CO/CH3OH and H2CO/CH3OH abundance ratios were found to be predominantly determined by the H/CO flux ratio and the surface temperature. This resulted in a layered structure of simulated grain mantle with CH3OH on top. The species CO and H2CO were found to exist mainly in the lower layers of ice mantles, which are formed in early times when the CO flux is relatively high, and these layers are not available for hydrogenation at late times. This finding discredits many gas–grain models, which do not take into account the layering of the ice. The extent of effect of layering can be seen in Figure 17. The right panel of Figure 17 shows the evolution of the CO, H2CO, and CH3OH ice abundances for a dense cloud of density of nH = 105 cm–3 and grain surface temperature of 12.5 K simulated with a microscopic KMC simulation. The left panel of Figure 17 shows the abundance evolution for the same species obtained with a rate equation model with similar input parameters and conditions. One can clearly see that whereas the initial ice composition is very similar, at later times the ice composition starts to deviate, since in the rate equation model also the lower layers are still available for hydrogenation, whereas in reality they are blocked by the top layers, which is included in the microscopic KMC model.

Figure 17

Figure 17. Evolution of CO, H2CO, and CH3OH ice abundance for a dense cloud of density of nH = 105 cm–3 and grain surface temperature of 12.5 K obtained by (left) microscopic KMC simulations and (right) a rate equation model. Adapted from ref 48. Credit: Cuppen, van Dishoeck, Herbst, Tielens. Astron. Astrophys. 2009, 508, 275 adapted with permission © ESO.

The obtained model results in ref 48 agree reasonably well with observations of solid H2CO/CH3OH and CO/CH3OH abundance ratios in the outer envelopes of an assortment of young stellar objects. The results suggest that the large range of CH3OH/H2O observed abundance ratios is due to different evolutionary stages of the objects. Extension of the limited network to include competing reactions that also lead to the formation of H2O did not alter the overall conclusions concerning H2CO and CH3OH formation.

Cazaux et al.(149) developed a model to form water and its deuterated forms including many different reaction routes. This model simulated the initial stages of water ice formation on carbonaceous grains in the submonolayer regime and builds on earlier models for H2 formation.(143) The model does not allow for a second layer to build up, and reaction products can desorb using the excess energy from the reaction. The aim of this model was to study the effect of surface chemistry on the gas phase composition. Different physical conditions were simulated using a microscopic KMC model to study this effect. Return to the gas phase occurs through photodesorption and chemical desorption. To model these energetic processes, molecular dynamics would probably be better suited, although the required time scales cannot be reached. At a low dust temperature of T = 10 K, hydrogenation reactions leading to species such as H2 and H2O are found to dominate; for higher temperatures (T = 30 K) oxygenation starts to take over, leading to the formation of O2 and O3. Different water formation routes were found to be important in different conditions: in diffuse clouds, formation is mainly through H + O → H2O in agreement with ref 147. Since this is a highly exothermic reaction and the chemical desorption probability is assumed to scale with the exothermicity, a large fraction of the formed water molecules is released into the gas phase. In dense clouds H2 + O → OH + H dominates (this reaction was not included in ref 147, and there is no experimental evidence for it to occur(153)), and since this reaction is less exothermic (it is in fact endothermic(154)), a smaller fraction desorbs into the gas phase. In photon dominated regions, more oxygenated species are released and also the water formation routes through O2 + H and O3 + H increase in importance.

Das et al.(85) developed a discrete-time, random-walk KMC model to look at the formation of CH3OH, H2O, and CO2 under different physical conditions. In discrete-time simulations, time is advanced using

(51)instead of eq 23 and the time steps do not follow a Poisson distribution. Very thick ices were obtained in this way (60–500 ML) and a range of physical conditions was determined where observations were best reproduced: the so-called favorable zone. Ice abundances of CO2 were underproduced in all conditions, probably since the OH + CO reaction was not included in their model. Only a small network containing 10 reactions was used; reactions involving H2 are the most striking omission, especially since both Cuppen and Herbst(147) and Cazaux et al.(149) concluded that reactions with H2 are the most important formation route for the production of H2O in dense cloud conditions. Also the H + O2 route is omitted, whereas large amounts of O2 are formed on the grain. In a follow-up paper, the reaction network was extended with six more addition reactions, now also including the H + O2, H + O3, and H2 + OH reaction channels,(86) as well as with 11 photodissociation reactions. Molecular hydrogen was still not allowed to land on the surface—probably to keep the simulation time within reasonable limits—and the H2 + OH reaction could only proceed through Eley–Rideal reactions, which seems a strict limitation. The aim of the paper was to study the effect of photodissociation on the evolution of the grain mantle. Photodesorption is only allowed to occur in the topmost layer; photodissociation is allowed deeper into the mantle. Physical conditions within the favorable zone were chosen. The percentage of photodissociated molecules as a function of AV and number of UV-exposed layers was studied. A high percentage was found for low AV, as could be expected, and this percentage increases with number of UV-exposed layers, until it saturates beyond 50 ML. The mantle was found to be more oxygenated for low visual extinction, in agreement with the work of Cuppen and Herbst.(147)

Several microscopic KMC studies investigate the so-called “snow line”. This is the line around a protostar that separates the inner area where bare grains exist from the outer area where grains are ice covered. Earth has formed within the snow line and should therefore be rather dry. Several hypotheses exist on how water arrived in the inner solar system. One of these hypotheses is put forward by Drake(155) and Stimpfl et al.,(156) who showed that chemisorption sites exist on forsterite grains where water can absorb. Grains can retain a submonolayer of water well beyond the snow line in this way. Muralidharan et al.(148) built on this work by applying a lattice-gas KMC simulation. The adsorption/desorption energies were mapped on three different fosterite surfaces ({100}, {010}, and {110}) by constructing a fine grid, placing one water molecule at a time on a grid point, and performing a geometry optimization to obtain the minimum energy. The resulting map for the {100} surface is shown in Figure 18. This map was then fed to a KMC model that determined the surface coverage of a grain as a function of grain temperature. Hopping rates were obtained from the diffusion constant that was used as an input parameter. Due to the many assumptions used in the model, the authors argue that the results only give a conservative lower limit to the amount of water that can be absorbed onto grains prior to accretion into planets. Their results show that this amount can indeed be substantial.

Figure 18

Marseille and Cazaux(150) constructed a water-accretion model where water molecules are physisorbed on carbonaceous grains and can cluster together, leading to a stronger binding, with a maximum of six water neighbors. Time scales for condensation of water ice below 100 K were determined in this way, and these were found to be rather long compared to time scales of grain mixing, caused by turbulence and infall. This results in a snow border, rather than a snow line, a region where bare and icy grains can coexist. This region can be significant compared to the size of the object. For a massive dense core this region can take up a fourth of the total system size.

7 New Directions

ARTICLE SECTIONS

Jump To

This section discusses the new directions toward which the field of grain surface chemistry modeling is moving. Because of the increased computer power, one can see an increase in complexity of the Monte Carlo models that are being reported. On the one hand, more molecular detail is introduced. This was already visible in the H on graphite simulations where more and more structural and system specific information was introduced. This evolves even further in off-lattice KMC models. On the other hand, models are developed in which the complexity of the chemical network is increased or in which gas and grain surface chemistries are combined. In the latter case, one can both notice existing gas–grain routines to get more detail in their grain surface treatment and surface models being extended with a gas phase treatment.

7.1 Unified Gas–Grain KMC Models

Several unified gas–grain KMC models have been recently developed.(88, 125, 157) The grain surface is treated at different levels of detail in these models. A first attempt to unify a microscopic KMC model of the grain surface with a method that treats the gas phase was by Chang et al.(46) They used rate equations to calculate the time-dependent chemical evolution of the gas phase including accretion onto grains, but no surface chemistry or evaporation was included. From this information, temporal fluxes were determined to use as input in a microscopic KMC model to treat the surface chemistry. Desorption of species was then included again in the rate equation treatment of the gas phase chemistry, and in this way the chemical evolution of both gas and grain mantle was obtained iteratively. For dense clouds, convergence was reached in two iterations. The success of this method comes from the fact that, for the conditions chosen in their paper, dense cloud conditions, there is very little return of grain surface species into the gas phase and the influence of the grain surface chemistry on the gas phase is therefore limited. Under these conditions, only atomic hydrogen evaporates to a significant extent. Although it has little effect on other gas phase species, the evaporation of atomic hydrogen changes its gas phase abundance, which in turn changes the flux of atomic hydrogen onto grains. For astronomical objects where more grain mantle material is evacuated into the gas phase, such as hot cores or protoplanetary disks, this influence becomes clearly larger and many more iterations will be needed to reach convergence.

Another approach to coupling gas phase chemistry and grain surface chemistry is to use the KMC approach for both phases. Early attempts did not result in the required resolution of the gas phase abundances of rare species,(158) but Vasyunin et al.(6) managed a way to circumvent this by cumulative averaging of the gas phase abundance in time. They demonstrated this by simulating molecular clouds at different temperatures and densities and comparing this to rate equation and master equation results for a large gas–grain chemistry network of roughly 6000 gas phase and 200 surface reactions.

Recently, a new unified KMC model was introduced in which macroscopic KMC is used to describe the gas phase chemistry, whereas microscopic KMC treats the surface chemistry.(88) A medium-size surface network of 29 reactions was used, leading to the formation of stable molecules such as CO2, CH4, CH3OH, and H2O. Although CH4 is overproduced in the simulations, the agreement with ice observations(159) for other carbon-bearing ices such as CO, CO2, and CH3OH is much more reasonable. In the analysis of the simulated ice mantles, an anticorrelation between CO and CO2 was found, in agreement with the observations. The observed CO-rich nonpolar ice segment is however not reproduced. In a few of the models, CO is formed early in the evolution of the grain mantle by reaction of C with O. Increase in the initial H/H2 ratio destroys this route, however, since then most of atomic oxygen reacts with H to form OH. Part of the simulated CO2 ice abundance is obtained through a “chain reaction” in which OH is formed by hydrogenation of O in the vicinity of a CO molecule. CO and OH can then react without much diffusion involved. In this way, the mobility leading to the formation of CO2 is done by the H atom that performed the initial reaction. Detailed mechanisms like this show the value of microscopic models when simulating grain chemistry.

Macroscopic models include more and more “microscopic” information in their models. This is achieved by including layering into the macroscopic model, where the composition of each layer is followed in time. Other models only distinguish between the surface and the bulk of the ice mantle: the so-called three-phase model.(160) Examples of gas–grain models that treat the ice mantle in at least two phases can be found in refs 124, 125, and 160−162. Since ref 125 is the only model that makes use of KMC, we will restrict ourselves to the discussion of this model.

This macroscopic model tracks the species per layer. It is assumed that the formation of a new monolayer starts when the previous layer is completely filled. Chemical reactions and desorption are only allowed to occur in the top monolayers. The exact number of chemically active layers was calibrated against desorption experiments as discussed in section 5.2. Here the assumption is made that chemical reactions only occur in the same layer as where desorption occurs. This requirement is probably too strict. High vacuum experimental studies of energetic processing clearly show that reactions in the bulk can occur. The model was employed to study the chemical evolution starting from a molecular cloud through the collapse and warm-up phases leading to a hot core/corino. The main motivation to use a layer-by-layer model was that it treats the desorption correctly, which is important in the warm-up phase. In general, a good agreement with ice and gas phase observations is obtained for all stages, better than when using a two-phase model: in particular, the formation of polar and apolar ice fractions in the cold stage and the observed jump abundances in H2CO in the hot corino phase. On the other hand, species such as HCOOCH3 which have photoproducts as precursors are underproduced. This is most likely because photochemistry and bulk diffusion should be active in the inner layers as well.

7.2 Combining Different Models

Introducing more microscopic detail into macroscopic models can also be done in a more pragmatic way, where dedicated, detailed microscopic simulation programs are used to obtain information about reaction efficiencies/rates as a function of density, temperature, etc. This is then given to the macroscopic model in terms of lookup tables or fitted analytical expressions. Chang et al.(163) give such lookup tables for H2 formation rates for different H fluxes, grain structure—rough vs flat surface and olivine vs carbonaceous—and surface temperatures. Separate tables are provided for the use in pure gas phase models and gas–grain models. In the latter case, effective evaporation rates and diffusion rates are provided. Reference 147 provided an analytical expression for the H2 formation rate on graphitic grains, resembling PAHs, with the density and gas temperature as input parameters. The dust temperature was kept constant at 15 K.

Something similar was done for the water formation in X-ray exposed environments (XDRs).(151) Here a rate equation model was constructed to determine the formation rates of OH and H2O and the release into the gas phase and analytical expressions for these quantities were then obtained. These were then included in a XDR model. KMC simulations using the model presented in ref 149 were used as a consistency check for the rate equation model. At dust temperatures below 35–40 K, the KMC and rate equation method yield the same results within the uncertainties. Since the simulations stay within the submonolayer regime, rate equations are relatively simple to construct and layering does not need to be considered. At higher dust temperatures, the system enters the stochastic regime and the rate equation method systematically obtains too high efficiencies. OH formation on grains and release to the gas phase is dominated by O + H → OH for low temperatures and by O3 + H → OH + O2 at higher temperatures, consistent with the earlier findings in ref 149. H2O formation on grains and release to gas phase is dominated by OH + H → H2O. The fraction of H2O desorbing into the gas phase is further increased under the influence of strong UV radiation fields by formation–dissociation–formation loops. These additional formation routes for warm gas phase water could help explain unusually strong water lines observed for Mrk 231, which has an accreting supermassive black hole with X-ray luminosity LX ∼ 1044 erg s–1. Typical star-forming environments, without an X-ray source, do not show these strong water lines.

7.3 Off-Lattice Models

One of the main restrictions of traditional KMC methods is the confinement of the system to a predefined lattice. The restriction of the atomic coordinates limits the amount of physical detail contained by the simulation and makes the use of realistic interaction potentials far from straightforward. In the solid state, the lattice approximation is often well justified for simple crystalline systems, but the situation already becomes more complicated for molecular crystals which often show some degree of disorder, such as many hydrates and the hydrogen bond network in H2O ices. Small site-to-site differences in these systems affect the rate constants and make the generation of the table of events a tedious task. For amorphous systems, the situation obviously becomes even worse and lattice-based KMC simulations of these systems should be more seen as toy models. Needless to say, the same is true for simulations of systems in the liquid phase.

A second problem with traditional KMC simulations is the need to define the table of events before the start of the simulation. As with the limitation to a lattice, this can be done for simple crystalline systems but quickly becomes more challenging as the complexity of the system increases and the transition mechanisms become less intuitive. A prime example of such a mechanism is the diffusion of H atoms into a CO ice, which proceeds through an exchange mechanism rather than the H atom taking up interstitial positions.(68) Ideally, one would like to evaluate every process individually, as accurately as possible. This can be done in off-lattice KMC by applying force fields(-83, 164-166) or even electronic structure methods(167) to calculate individual process rates. As every process is considered to be unique, this evaluation needs to be done on-the-fly, which considerably increases the computational costs. The advantage of not having any a priori assumptions about processes in the system is, however, often worth this additional effort.

In sections 7.3.1 and 7.3.2 we will describe two different approaches which take KMC simulations beyond the lattice approximation. The continuum kinetic Monte Carlo method(168) is a particularly interesting off-lattice KMC method which incorporates diffusion of particles analytically into a continuous time KMC scheme. Reaction rates and diffusion constants still need to be specified before the simulation starts, though. The adaptive kinetic Monte Carlo(169) method is an on-the-fly, off-lattice, KMC method which allows for very high atomistic detail, requiring no other input than an interatomic interaction potential.

7.3.1 Continuum Kinetic Monte Carlo

The continuum kinetic Monte Carlo method(168) was designed for simulating reactions in solution, a typical case where one would like to go beyond the lattice approximation. The challenge lies in the fact that conformational changes, reorientations, and individual diffusion steps often occur on much shorter time scales than the chemistry itself. In continuum KMC, this problem is solved by coarse graining the master equation, which allows for diffusion to be treated in an analytic way, whereas the chemical reactions are simulated in a variation on the CTRW algorithm.

A continuum KMC cycle starts, just like CTRW KMC, by constructing a list of possible reactions and their reaction rates. In CTRW the event time is then determined for each particle following eq 23. In continuum KMC the event time for each possible reaction is calculated. This calculation of reaction times, though, is somewhat more involved than in the case of a lattice gas. The rate constants depend on, besides the reaction barrier, also the diffusion constants of the reactants, their initial positions, and time itself. Once the list is complete, the first reaction is executed and the positions of the reacting particles are updated. These new positions are determined by assuming a Gaussian distribution for the diffusing particles by generating a random number from the appropriate probability distribution. After this, the reaction list is updated and a new cycle starts.

The advantage of the continuum KMC method lies in the fact that it treats diffusion analytically, enabling the simulation of rather large systems for long time scales. Predefined reaction rates and diffusion constants keep the computational effort down to a very moderate level but are at the same time an assumption in the model. For dilute systems, where the reactant diffusion is limited by the inert solvent, this is likely well justified, but at higher concentrations this could become a problem. Until now, the applications of the method remain limited(168, 170) but the approach certainly offers appealing features which may possibly also be useful for simulating astrochemically relevant processes such as reactions between minority species in ice mantles in molecular clouds.

7.3.2 Adaptive Kinetic Monte Carlo (AKMC)

Predefined event tables (process rates) in KMC simations are not only severe assumptions in the model; they can also be very tedious to construct. During its construction, one should always ensure that all important, high rate, processes have been included, because a missing process in the event table will not only close certain evolution routes; it will also cause the simulation time to progress too slowly. For complex systems it is therefore desirable to have an unbiased method for finding transition mechanisms, eradicating possible human errors. This is the essence of the adaptive kinetic Monte Carlo method.(169) With this method, it is possible to run off-lattice kinetic Monte Carlo simulations without having to define the event table beforehand. The two most severe limitations of traditional, lattice-gas KMC are therefore not present when using AKMC.

Since the AKMC method is an off-lattice method, a state is defined as a local minimum in the potential energy surface. At the start of a simulation, only one of these minima is known: the initial state. As a first step, the AKMC program needs to determine the table of events for all processes leaving this initial state. This requires the exploration of the PES in the vicinity of the current state to search for other local minima that can be reached via a transition state and to determine the corresponding rate for this process. To locate the relevant transition states out of this minimum, single-ended search methods are used. Double-ended methods such as NEB cannot be applied to find transition states, since they require knowledge of both the initial and final states, whereas the AKMC simulation has information only about the initial state. Several single-ended methods exist,(171, 172) but the minimum-mode-following method(173) is the most frequently used. This method can be used when the transition state is assumed to be a single point on the potential energy surface (a first order saddle point). Olsen et al.(171) have performed a comparative study of different methods for finding saddle points without knowledge of the final states. They conclude that minimum-mode-following methods are to be preferred if a more complete mapping of all existing (low-energy) saddle points is desired for large systems, which is typically the case for AKMC simulations. The method requires the lowest-eigenvalue eigenvector of the Hessian which can be calculated exactly or approximated by a method such as the dimer method,(173) Lagrange multipliers,(174)or the Lanczos method.(175)

As soon as a transition state is found, the rate of the corresponding process is calculated. Until now, harmonic transition state theory has been used for the rate estimation, but any other method may also be used to account for effects such as quantum tunneling and zero point motion. By performing a series of transition state searches, the table of events of the initial state is filled up until there is sufficient confidence that all low-lying saddle points have been found.(167) Once this is the case, time is advanced and the system proceeds to a new state according to the n-fold way algorithm.

The need to perform transition state searches every time a state is entered for the first time makes the AKMC method more demanding from a computational resource point of view than traditional lattice KMC. However, since the individual searches are fully independent, they can be distributed over a large number of computational nodes in a straightforward manner.(176) This feature is currently implemented in the EON code.(177) Naturally, when a state is entered which has been visited before, no new searches need to be performed as the table of events is already known for that particular state. For systems with only a limited number of thermally accessible states, this means that eventually the relevant part of the PES will be fully explored and then the KMC simulation can progress just as quickly as a lattice-gas simulation with a predefined event table. This situation was for example encountered in the simulation of CO diffusion of hexagonal ice, which will be discussed below. In this case, due to the low temperatures, only the weakly bound CO molecule was free to move, and when the surface binding sites were all explored, long time kinetic Monte Carlo trajectories could be run to extract the surface diffusion constant. Notice that, in this case, the AKMC simulation essentially reduces to a normal KMC simulation, with a predefined event table which has been constructed in an unbiased way, i.e., without any possible bias following from human intuition. Moreover, the rates are determined from a realistic interaction potential and are not estimated based on experimental information.

The aforementioned study of CO diffusion on hexagonal water ice, the first application of AKMC in an astrochemical context,(83) also showed the importance of an unbiased way to construct the table of events. Despite the near-crystalline nature of the water ice substrate, the dangling bond pattern on the surface leads to a large variation in both the CO orientations and the binding energies. This is illustrated in Figure 19, where both the configurations of the binding sites and a sketch of the potential energy landscape are shown. Since the typical binding energy of the CO molecule is much lower than the cohesive energy of the ice substrate, all the surface binding sites could be explored without the substrate itself evolving. Long-time kinetic Monte Carlo trajectories could therefore be obtained which enabled the extraction of the surface diffusion constant.

Figure 19

Figure 19. Binding sites and potential energy surface of CO on the basal plane of hexagonal ice, investigated with the AKMC method.(83) The top panel shows the binding sites on the substrate. The large variation in CO orientation on the surface is an example of the importance of having an unbiased method for finding potential energy minima and transition states. The bottom panel shows the potential energy surface. Again, due to the proton disorder in the substrate, a large variation in binding energies is observed. Reproduced from ref 83 by permission of the PCCP Owner Societies.

The required resources for AKMC simulations heavily depend on, besides the obvious system size, the method by which interactions are calculated and the number of states which have to be explored. For small systems with only a limited number of states, the method permits the use of density functional theory as was demonstrated by Xu et al., who calculated the dynamics of methanol decomposition on Cu(100)(178) and Pd island formation on MgO(100).(167) For larger systems or when many states need to be explored, simpler interaction models, such as effective pair potentials, can be used.(83, 166)

Several techniques have been developed to improve the efficiency of the transition state searches and thereby reduce the computational cost of AKMC simulation. One of these techniques is the recycling of previously found saddle points, a method proposed by Xu and Henkelman.(167) In this case, the knowledge of the transition state configurations from previous states is used as an initial guess for transition states in unexplored states. This method is particularly effective when the system is still rather ordered and only a few atoms are significantly involved in the considered process. Another method which can speed up the simulations is the cancellation of transition state searches, if they reach a point where there is a high certainty that the saddle point which will be found was already identified during a previous search.(179)

The power of the AKMC method lies in its ability to perform long-time-scale simulations of thermally driven rare-event systems while retaining information of all atomic coordinates as the system evolves in time. This makes for a powerful tool when studying kinetic processes in grain mantles. Of course, the high amount of structural detail limits the system size and simulation time compared to traditional KMC. However, individual processes can be systematically investigated at high accuracy. Dynamics of important processes such as diffusion, segregation, and even chemical reactions (given an interaction method capable of treating reactive processes) can be studied using AKMC, and the results may then be fed to larger scale models such as on-lattice KMC or rate equation methods.

7.4 Coarse Graining

A common problem in all KMC simulations occurs when the rates in the table of events differ from each other by several orders of magnitude. This situation can for example be encountered when energy barriers are very different or when the temperature is low. When simulating interstellar ices, both cases apply and the KMC simulation may become stuck in a certain subset of states, wasting many iterations without interesting dynamics. The problem is illustrated in Figure 20. Due to the low barrier (high rate) between states 1 and 2, the KMC simulation will repeatedly pick the process between these states without evolving to the other states, whereby wasting valuable computational resources.

Figure 20

Figure 20. Schematic potential energy surface with KMC states. Without coarse graining, the simulation will be stuck in states 1 and 2 for a long time, due to the low barrier. Coarse graining can solve the problem by merging states into composite states.

Novotny(180) provided a solution to this problem by using the theory of absorbing Markov chains. Using this theory, the states which are joined by the high-rate process are merged together to form a composite state. This composite state can in principle contain any number of “real” states. When the simulation visits a composite state, dynamics are not treated by the usual KMC theory that dictates the choice of the leaving state and the time evolution but rather by the use of absorbing Markov chains. This gives, in one step, the process through which the simulation will leave the composite state and the corresponding time it spends in the composite state. Hence, as Figure 20 illustrates, the processes between states inside the composite state are effectively removed from the event table, allowing the simulation to progress much more efficiently. It is important to note that the dynamics generated in a coarse-grained simulation follow exactly the same statistics as in a normal KMC simulation. However, the exact evolution inside the composite states is lost.

The decision to merge states into a composite can be based on several criteria. One such algorithm was proposed by Pedersen et al.(181) Here, each state is assigned a fictitious energy level, initially equal to the total energy of the state. This level is then raised by each subsequent visit to the state, and once it becomes higher than the transition state energy to a neighboring state, the states are merged into a composite. The application of coarse graining can accelerate KMC simulations by many orders of magnitude in terms of both CPU time and simulation time per iteration. Figure 21 shows the speedup factor of AKMC simulations of the aforementioned CO diffusion study on top of crystalline water ice as a function of temperature. Speedup is here defined as the factor by which the simulated time increases with respect to a simulation without coarse graining using the same number of iterations. Since the difference in rates increases with decreasing temperature, for thermally activated processes, the coarse graining algorithm is clearly more advantageous at lower temperature where speedups of 13 orders of magnitude can be achieved. Figure 21 further shows that the obtained speedup depends on the number of states inside the composites: the larger the composite states grow, the higher the speedup is. For most physical applications, however, one should not let the composite state size grow indefinitely because the exact dynamics inside the composites is not retained in the absorbing Markov chain theory. It is therefore advisable to put a limit on the number of states inside each composite. This limit should be chosen such that there is a good balance between the obtained speedup and the chemical detail in the simulations.(83)

Figure 21

Figure 21. Speedup of CO on water crystalline water ice AKMC simulations through coarse graining as a function of temperature. Reproduced from ref 83 by permission of the PCCP Owner Societies.

The application of coarse graining can significantly extend the time scales which can be reached by KMC simulations. It is particularly at low temperatures, when small differences in barrier heights lead to huge differences in process rates, that the method is at its most useful. This makes the method particularly applicable in KMC simulations of astrochemical systems. Since the underlying dynamics remain unaltered, the physical quantities extracted from the simulation remain the same; only the time over which statistics can be gathered becomes longer and the temperatures may be lower.

8 Conclusions and Outlook

ARTICLE SECTIONS

Jump To

Kinetic Monte Carlo was initially introduced in surface astrochemistry to circumvent the accretion limit, in which classical rate equations break down and obtain abundances which are off by many orders of magnitude. The technique, although still not generally applied, has come a long way since then. Since the method has no real restriction on how molecules are represented—in terms of number densities or in full atomistic detail—it has proven its additional strength by allowing layering and more detail in the simulations. We have discussed several examples in the previous sections. Mostly lattice KMC models are applied with the aim of (i) gaining more insight in the physicochemical mechanisms behind astrochemically relevant experiments or (ii) modeling grain surface chemistry with a relatively high level of detail. As discussed in the earlier sections, several of the interstellar ice studies have revealed that bulk diffusion is still an ill-understood phenomenon that cannot be well modeled using the lattice-gas KMC models. It appears however to be the limiting factor for obtaining better agreement between models and experiment. This agreement is a requirement for realistic models that can be applied to predict the chemical evolution of different astrophysical environments/objects. We therefore see the field moving in two directions: on one hand, gas–grain models include more microscopic detail and layering in order to introduce some distinction between surface and bulk diffusion—bulk diffusion is usually switched off entirely in these models; on the other hand, we see the move to a more atomistic description of the system, which allows for detailed investigation of the underlying mechanisms behind bulk diffusion. Unfortunately, the application of these atomistic models remains limited due to their high computational load. The real progress in the field will be made if these two approaches can be combined by feeding the results of the atomistic models to the gas–grain models.

The authors declare no competing financial interest.

Acknowledgment

ARTICLE SECTIONS

Jump To

H.M.C. is grateful for support from the VIDI research program 700.10.427, which is financed by The Netherlands Organisation for Scientific Research (NWO). H.M.C. and L.J.K. acknowledge the European Research Council (ERC-2010-StG, Grant Agreement No. 259510-KISMOL) for financial support. T.L. is supported by the Dutch Astrochemistry Network financed by The Netherlands Organisation for Scientific Research (NWO). The authors would like to thank Stefan Andersson and Eduardo Monfardini Penteado for general discussions concerning this review and Eric Herbst, Stephanie Cazaux, Qiang Chang, Kinsuk Acharyya, and Jochiam Krug for fruitful discussions regarding microscopic reversibility in particular.

References

ARTICLE SECTIONS

Jump To

This article references 181 other publications.

Viti, S.; Collings, M. P.; Dever, J. W.; McCoustra, M. R. S.; Williams, D. A. Mon. Not. R. Astron. Soc. 2004, 354, 1141

Translate

22 Aug 2020

The Kinetic Monte Carlo Method as a Way To Solve the Master Equation for Interstellar Grain Chemistry

Labels