Gas-grain chemistry in molecular clouds

A unified Monte Carlo treatment of gas-grain chemistry for large reaction networks. I. Testing validity of rate equations in molecular clouds

A.I. Vasyunin Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany D.A. Semenov Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany D.S. Wiebe Institute of Astronomy of the Russian Academy of Sciences, Pyatnitskaya str. 48, 119017 Moscow, Russia Th. Henning Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

In this study we demonstrate for the first time that the unified Monte Carlo approach can be applied to model gas-grain chemistry in large reaction networks. Specifically, we build a time-dependent gas-grain chemical model of the interstellar medium, involving about 6000 gas-phase and 200 grain surface reactions. This model is used to test the validity of the standard and modified rate equation methods in models of dense and translucent molecular clouds and to specify under which conditions the use of the stochastic approach is desirable.

Two cases are considered: (1) the surface mobility of all species is due to thermal hopping, (2) in addition to thermal hopping, temperature-independent quantum tunneling for H and H is allowed. The physical conditions characteristic for the core and the outer region of TMC1 cloud are adopted. The gas-phase rate file RATE 06 together with an extended set of gas-grain and surface reactions is utilized.

We found that at temperatures 25–30 K gas-phase abundances of HO, NH, CO and many other gas-phase and surface species in the stochastic model differ from those in the deterministic models by more than an order of magnitude, at least, when tunneling is accounted for and/or diffusion energies are 3x lower than the binding energies. In this case, surface reactions, involving light species, proceed faster than accretion of the same species. In contrast, in the model without tunneling and with high binding energies, when the typical timescale of a surface recombination is greater than the timescale of accretion onto the grain, we obtain almost perfect agreement between results of Monte Carlo and deterministic calculations in the same temperature range. At lower temperatures ( K) gaseous and, in particular, surface abundances of most important molecules are not much affected by stochastic processes.

astrochemistry ISM: abundances, clouds, molecules — molecular processes – radio lines: ISM — stars: formation

1 Introduction

Chemical processes are discrete in nature, and it has been realized long time ago that their microphysically correct theoretical treatment should rest upon stochastic methods (Gillespie, 1976). While in cellular biology Monte Carlo simulations of chemical processes in cells are widespread, most models in astrochemistry are based on deterministic rate equation (RE) approach, which has been proved to be inaccurate for chemical processes on “cell-like” objects like dust grains and PAHs (Tielens & Hagen, 1982; Caselli et al., 1998; Herbst & Shematovich, 2003; Lipshtat & Biham, 2003; Stantcheva & Herbst, 2004). The concept of mass-action kinetics, which lies behind rate equations, being adequate for gas-phase reactions, is not appropriate in the case of small surface population (number of particles per grain). This situation is often the case in grain surface catalysis. Rates of surface reactions obtained through this approach may be strongly incorrect, leading to an improper estimation of timescales of fundamental processes, like H formation, and significant errors in abundances of other species.

Up to now, a number of attempts have been made to assess the significance of stochastic effects for astrochemical simulations and to develop suitable numerical methods. In the context of astrochemistry first attempt to treat grain surface chemistry stochastically has been made by Tielens & Hagen (1982). In this study, equilibrium abundances of gas-phase species were used to obtain accretion rates and to calculate time-dependent populations of surface species through Monte Carlo approach. All surface species but H were assumed to remain on the grain surfaces.

Subsequent attempts to account for the discrete nature of surface reactions can be divided into two categories. The first category encompasses different modifications to the rate equations for grain surface chemistry (Caselli et al., 1998, 2002; Stantcheva et al., 2001). The main idea behind these modifications is to restrict rates of diffusive surface reactions by accretion or desorption rates of reactants. Thus, the so–called accretion-limited regime is taken into account, when the rate of a particular surface reaction is determined by the flux of reactants accreting on a grain and not by their diffusion over the surface. For simple reacting systems, results obtained with this semi-empirical approach were found, at least, to be in closer agreement with Monte Carlo (MC) computations than the treatment by rate equations alone. However, the applicability of the modified rate equations (MRE) to large chemical networks has not been tested so far. Moreover, it is impossible to investigate purely stochastic effects, like bistability, with this technique, still purely deterministic in its nature (Le Bourlot et al., 1993; Shalabiea & Greenberg, 1995; Lee et al., 1998; Boger & Sternberg, 2006).

In the second category of studies the master equation is solved with different methods. Because the chemical master equation describes probabilities of a chemically reacting system to be in all possible states of its phase space and the number of these states grows exponentially with the number of species, their direct integration is only possible for extremely simple systems, consisting of a very few different components. Biham et al. (2001), Green et al. (2001) and Lipshtat et al. (2004) have performed a direct integration of the master equation to investigate H formation. These studies demonstrated that in some cases this fundamental process cannot be adequately described with the rate equation approach. In Stantcheva et al. (2002), Caselli et al. (2002) and Stantcheva & Herbst (2003) direct integration of the master equation was applied to study the evolution of a simple H—O—CO surface chemical network and deuterium fractionation. Direct integration of the master equation in more complex networks of grain surface reactions is hampered by the lack of appropriate computing power. To the best of our knowledge, the only successful simulation of a more extended surface chemical network in astrochemistry, consisting of 19 reactions, with direct integration of the master equation was accomplished by Stantcheva & Herbst (2004).

A Monte Carlo approach based on continuous-time random-walk method has been developed by Chang et al. (2005) to model recombination of hydrogen atoms on interstellar grains. The advantage of this method is that it allows to model formation of H and other molecules on surfaces of arbitrary roughness, which has later been demonstrated by Cuppen & Herbst (2007).

To relax computational requirements needed to solve the master equation, the moment equation approximation was suggested by Lipshtat & Biham (2003) and Barzel & Biham (2007). Currently, this is the most promising approach to efficient simulations of gas-grain systems with complex grain surface chemistry, because the moment equations can be easily combined with rate equations used to simulate gas-phase chemical processes. But, strictly speaking, the moment equations are an approximation and their validity should be verified by a comparison with exact methods.

The only feasible technique allowing us to obtain an exact solution to the master equation for complex chemically reacting systems is Monte Carlo method developed by Gillespie (1976). This approach is widely used in molecular biology for simulation of chemical processes in living cells and can also be applied to astrochemical problems. Gillespie’s Stochastic Simulation Algorithm (SSA) was first used in complex astrochemical networks by Charnley (1998), where SSA was applied to gas-phase chemistry only. Later a simple grain surface H—O—CO chemical model, similar to the one considered in Caselli et al. (1998), was developed by Charnley (2001). In this model, gas-phase abundances of species were assumed to be constant and desorption was neglected. No coupled gas-grain chemical model was studied.

The chemical master equation is usually considered as a method to model surface reactions. However, due to its universality, it can be used to model the gas-phase chemistry as well. The reason why the rate equations are so popular as a tool to study gas-phase chemistry in astronomical environments is the fact that they are easier to implement and much more computationally effective. However, if one uses the rate equations to model gas-phase chemistry and some stochastic approach to model surface chemistry one has to find a way to match the two very different computational techniques. The first study in which a Monte Carlo treatment of grain surface chemistry was coupled to time-dependent gas-phase chemistry is Chang et al. (2007). In this paper iterative technique is utilized to combine very different Monte Carlo and rate equations approaches. Calculations are performed for the time interval of years and a relatively simple H—O—CO surface network consisting of 12 reactions is used.

Even though remarkable progress has been made to develop stochastic approaches to astrochemistry, up to now there are no studies in which a stochastic treatment of a complex grain surface chemical network is fully coupled to time-dependent gas-phase chemistry. The situation gets less complicated when the same method is used for both gas-phase and surface chemistry. In this paper, for the first time, we present a “complete” time-dependent chemical gas-grain model, calculated with Monte Carlo approach. We solve the master equation with a Monte Carlo technique, using it as a single method to treat gas-phase reactions, surface reactions, and gas-grain exchange processes simultaneously. The gas-phase chemical network includes more than 6000 reactions while surface network consists of more than 200 reactions.This Monte Carlo model is used to test the validity of rate equations and modified rate equations for a range of physical conditions, typical of diffuse and dense clouds. So far, no comparison of RE and MRE techniques against rigorous stochastic approach over a wide range of physical conditions, typical for astrophysical objects, has been made.

There is clear observational evidence that interstellar grains have complex mineral composition and may have either amorphous or crystalline structure. In general, it is believed that interstellar dust particles are made of olivine-like silicates and some form of amorphous carbon or even graphite (Draine & Lee, 1984). The location on a surface which may be occupied by an accreted molecule is commonly referred to as a surface site. The density of sites and binding energies for specific species strongly depend on the dust material and its structure. Binding energies , that is, potential barriers between adjacent sites, have been assumed by Hasegawa et al. (1992) to be roughly 0.3 of corresponding desorption energies . In Pirronello et al. (1997a, b, 1999), Katz et al. (1999) and Perets et al. (2005) experimental results on the H formation, the ratio of diffusion energy to desorption energy, , was found to be close to 0.77 for H atoms. This ratio determines the surface mobility of species and is of utter importance for estimating surface reaction rates.

Another factor defining rates of reactions involving H and H (and their isotopologues) is related to the possibility of their tunneling. Non-thermal quantum tunneling of these species allows them to scan grain surface and find a reacting partner quickly. This mechanism has been suggested as an explanation of efficient H formation in the ISM. Later analysis of the cited experimental results has lead Katz et al. (1999) to the conclusion that tunneling does not happen. The absence of tunneling implies lower mobility, so that surface processes are essentially reduced to accretion and desorption. However, this conclusion turned out to be not the final one. In other studies (e.g. Cazaux & Tielens, 2004) alternative explanations of these experiments have been proposed and some shortcomings of the theoretical analysis performed in Katz et al. (1999) have been found. So, up to now the question of presence or absence of tunneling effects on grain surfaces is not settled, and the ratio is not well constrained. With this in mind, two different models of grain surface reactions are considered: (1) the surface mobility of all species is caused by thermal hopping only with high ratio equal to 0.77, (2) temperature-independent quantum tunneling for H and H is included in addition to thermal hopping with low ratio equal to 0.3.

The organization of the paper is as follows. In Section 2 we describe physical conditions adopted in simulations and chemical model used. Section 3 contains basics of stochastic reaction kinetics and detailed description of the Monte Carlo code used for simulations. In Section 4 the validity of rate equations and modified rate equations is checked against Monte Carlo method. First, global agreement between methods is investigated. Then, some interesting species are discussed separately. In Section 5 a discussion is presented. Section 6 contains the conclusions.

2 Modeling

2.1 Physical conditions

In the present study, we consider grain surface reactions under physical conditions typical of irradiated translucent clouds, cold dark cores and infrared dark clouds. We consider temperatures between 10 and 50 K, densities between 200 and  cm (three values of and five values of ), and visual extinctions of 0.2 mag at  cm, 2.0 mag at  cm and 15 mag at  cm (e.g., Snow & McCall, 2006; Hassel et al., 2008). This roughly corresponds to the distance of order of 0.4–0.5 pc from the cloud boundary (at constant density). The cloud is illuminated by the mean interstellar diffuse UV field. The dust temperature is assumed to be equal to gas temperature. We do not study the earliest stage of the molecular cloud evolution, which is essentially the stage when atomic hydrogen is converted into hydrogen molecules. The chemical evolution is simulated assuming static physical conditions. While it is still a matter of debate how long does it take for a typical isolated cloud core to become gravitationally unstable and to start collapsing, we assume a short timescale of 1 Myr which seems to be appropriate for the TMC1 cloud (see e.g. Roberts et al., 2004; Semenov et al., 2004).

2.2 Chemical model

The chemical model is the same as described in Vasyunin et al. (2008). The gas-phase reactions and their rates are taken from the RATE 06 database, in which the effects of dipole-enhanced ion-neutral rates are taken into account (Woodall et al., 2007). All reactions with large negative activation barriers are excluded. The rates of photodissociation and photoionization of molecular species by interstellar UV photons are taken from van Dishoeck et al. (2006). The self- and mutual shielding of CO and H against UV photodissociation are computed as described in van Zadelhoff et al. (2003) using pre-calculated factors from Tables 10 and 11 from Lee et al. (1996). Ionization and dissociation by cosmic ray particles are also considered, with a cosmic-ray ionization rate of  s (Spitzer & Tomasko, 1968).

The gas-grain interactions include accretion of neutral species onto dust grains, thermal and photodesorption of mantle materials, dissociative recombination of ions on charged grains and grain re-charging processes. The dust grains are uniform m spherical particles made of amorphous silicates of olivine stoichiometry (Semenov et al., 2003), with a dust-to-gas mass ratio of . The sticking probability is 100%. Desorption energies and surface reaction list are taken from Garrod & Herbst (2006).

The grain surface is assumed to be compact, with surface density of  sites cm, which gives sites per grain (Biham et al., 2001). We employ two models to calculate the rates of surface reactions (Table 1). In Model T (T for tunneling) the tunneling timescale for a light atom to overcome the potential barrier and migrate to another potential well is computed using Eq. (10) from Hasegawa et al. (1992), with the barrier thickness of 1Å. In Model H (H for hopping) we do not allow H and H to scan surface sites by tunneling. The diffusion timescale for a molecule is calculated as the timescale of thermal hopping multiplied by the total number of surface site and is given by Eq. (2) and Eq. (4) from Hasegawa et al. (1992). The hopping rates are sensitive to the adopted values of diffusion energy, which are not well constrained. Thus, we consider low and high diffusion energies which are calculated from adopted desorption energies by multiplying them by factors of 0.3 (like in Hasegawa et al., 1992) (Model T) and 0.77 (Katz et al., 1999) (Model H), respectively. For all considered models, the total rate of a surface reaction is calculated as a sum of diffusion or tunneling rates divided by the grain number density and multiplied by the probability of reactions (100% for processes without activation energy).

Overall, our network consists of 422 gas-phase and 157 surface species made of 13 elements, and 6002 reactions including 216 surface reactions. As initial abundances, we utilize the “low metallicity” set of Lee et al. (1998), where abundances of heavy elements in the gas are assumed to be severely depleted. All hydrogen is molecular initially. The chemical evolution for 1 Myr in the classical deterministic approaches is computed with the fast “ALCHEMIC” code111Available upon request: in which the modified rate equations are implemented according to Caselli et al. (1998). No further modification for reactions with activation energy barriers is used (see Caselli et al., 2002).

3 Stochastic reaction kinetics

3.1 Theoretical foundations

We consider a chemically reacting system which consists of different types of species {} and chemical reactions {} (Gillespie, 1976). All these species are contained in a constant volume , in which local thermal equilibrium is reached, so that the system is well mixed. We denote the number of species of type at time as . The ultimate goal is to determine the state vector of the system at any given time , assuming certain initial conditions .

Let us assume that each reaction retains properties of a Markov chain and can be considered as a set of independent instantaneous events. Each chemical reaction is described by two quantities, the discrete state vector and the propensity function . The components of the state vector represent the net change in populations of the th species due to the th reaction. By definition, the propensity function


is the probability for a reaction to occur in the volume over the time interval . In analogy to the reaction rate in the rate equation approach, for a one-body reaction the propensity function can be expressed through rate constants and numbers of reactants as


Here, is the rate constant of the th reaction and is the absolute population of reagent. For a two-body heterogeneous reaction this expression changes as follows:


where is the reaction rate constant, and are the abundances of first and second reactants. If we deal with a two-body homogeneous reaction (like H + H  H), the expression for its propensity function becomes somewhat more complicated:


This expression reflects the fact that the rate of a homogeneous two-body reaction is proportional to the number of all possible pairs of its reactants. The term () is important for a proper calculation of the reaction rate in the stochastic regime when the average reagent population is close to unity and cannot be properly reproduced by the rate equation approach. Note that the abundances of species are always integer numbers in these expressions.

These definitions allow to characterize the microphysical nature of chemical processes and to establish a basis for stochastic chemical kinetics. A vector equation that describes the temporal evolution of a chemically reacting system by stochastic chemical kinetics is the chemical master equation. To derive this equation one has to introduce the probability of a system to be in the state at a time :


Here is the conditional probability density function of the time-dependent value . Temporal changes of are caused by chemical reactions with rates defined by Eq. (1). Therefore, the change of over the time interval is the sum of probabilities of all possible transitions from the state to the state :

The first term in this equation is derived from the fact that the system is already in the state and describes the probability for the system to leave the state during the time interval . The second term defines the total probability for the system to reach the state from states during the time interval due to a reaction . Finally, we formulate the chemical master equation:


This vector equation looks similar to the conventional set of balance equations. However, it properly takes into account the discrete nature of chemical processes and, thus, has a much wider range of applicability. Unfortunately, an analytic integration of the chemical master equation is only possible in a very limited number of cases (several species and reactions, see e.g. Stantcheva & Herbst, 2004). Thus, one has to rely on various numerical techniques to solve this equation.

3.2 Implementation of the Monte Carlo algorithm

In principle, any normalization of abundances can be used in Eq. (6). In particular, if we take to be a unit volume, we end up with the usual number densities, which are widely used in astrochemical models. When dealing with mixed gas-dust chemistry, one has to be more careful with the normalization. As surface reactions can only proceed when both reactants reside on a surface of the same grain, we define as the volume of the interstellar medium that contains exactly one dust grain,


where is the grain radius, is the mass density of grain material, is the mass gas density, and is the dust-to-gas (0.01) mass ratio. Here grains are assumed to be spheres of equal size. As we also assume that the interstellar medium is well-mixed, so that the volume is representative for any volume of the real interstellar medium with the same physical conditions.

This setup allows to construct the unified master equation for all three kinds of processes listed above (gas phase reactions, accretion/desorption processes and grain surface reactions). This equation is solved with the stochastic Monte Carlo algorithm described in Gillespie (1976). We implemented this technique in a FORTRAN77 code which allows to simulate all chemical reactions in the network in a self-consistent manner. Such a “brute-force” approach requires substantial CPU power and cannot be utilized in massive calculations. In the present study, it is used as a benchmark method to simulate chemical evolution in astrophysical objects. A typical run on a single Xeon 3.0GHz CPU takes between 10 hours and several days of computational time, and involves several billions of time steps. If one is only interested in a single point then the model can be used not only for benchmarking, but also for practical purposes. Unfortunately, it becomes impractical if abundances are computed for a number of spatial locations. High density also slows down the computation significantly.

Due to the fact than the Monte Carlo technique operates with integers, the smallest abundance of a molecule is 1. Given the dust properties in our model (radius cm, density g cm), a dust-to-gas ratio of 1/100 means that volume per one grain is about cm, where is the number density of hydrogen nuclei. Therefore, an absolute population of 1 corresponds to a relative abundance of with respect to the total number of hydrogen atoms (see Eq. 7). This is the lowest abundance directly resolvable in our model. However, the huge amount of tiny time steps taken during the calculations allows to average stochastic abundances over wider time spans of 10–1 000 years and push the smallest abundance resolvable by the Monte Carlo code below :


where is the abundance of species after the th time step, the is the time step in seconds and is the amount of time steps over which the averaging is performed. In such an averaging procedure, the noise in the random variable decreases with as and enables us to resolve abundances below . In our simulations we averaged abundances over time steps. Therefore, the resolution of our calculations is about 10 cm. The average value of the Monte Carlo time step in Gillespie’s algorithm is equal to an inverse sum of rates of all chemical reactions in the model. For example, in the dense cloud model many molecules form rapidly at initial times and less actively toward the end of evolution (Fig. 1). Time intervals, corresponding to time steps in a medium with density cm vary from 10 to years for early and late time moments, correspondingly.

4 Comparison of the Monte Carlo method with rate equations

While the Monte Carlo (MC) technique seems to be the most adequate method to solve the chemical master equation, it is computationally extremely demanding. Rate equations (RE) do not account for the stochastic nature of surface processes and thus may produce spurious results in some circumstances. On the other hand, they do require much less computer power and are easier to handle. Therefore rate equations will remain an important tool in theoretical astrochemistry.

Given the importance of surface processes, it is thus necessary to isolate regions in the parameter space where rate equations should not be used at all or, at least, should be used with caution. Similar comparisons, having been made so far, are based on a very limited number of species and reactions (e.g. Charnley, 1998, 2001; Ruffle & Herbst, 2000; Green et al., 2001; Caselli et al., 2002; Stantcheva et al., 2002; Stantcheva & Herbst, 2003, 2004; Lipshtat & Biham, 2003; Barzel & Biham, 2007; Chang et al., 2007). As we implemented a unified Monte Carlo approach to the gas-phase and surface chemistry, which is capable to treat a “regular size” chemical network, we are able to perform a comprehensive comparison between predictions of the MC and the (M)RE approaches.

For both the Model T and the Model H, the rate equations (RE), the modified rate approach (MRE) (for surface chemistry) and the Monte Carlo method are used to simulate gas-phase and surface chemistry. The differences in abundances discussed are related to surface chemistry and/or dust-gas interactions. For the gas-phase chemistry, all models produce identical results. Each run represents a combination of the tunneling treatment (T or H) and the surface chemistry treatment (RE, MRE, or MC). In total, we considered six models denoted as T-RE, T-MRE, etc.

4.1 Global Agreement

To perform an overall comparison, in Figures 2 and 3 we show diagrams of global agreement between MC, RE, and MRE calculations as a function of and for evolutionary times , , and years. For a given species, the two methods, being compared, are assumed to agree at time if they produce abundances, that differ by no more than an order of magnitude. As we have mentioned above, in the MC method we are able to resolve abundances down to , which would lead to a formal disagreement between abundances, say, in the MC method and (or smaller) in the RE method (all abundances are relative to the total number of H atoms). This kind of disagreement is not meaningful from the observational point of view because such small relative abundances are extremely difficult to observe with a satisfactory accuracy. In addition, even though lower abundances can be calculated by the RE method, their actual accuracy is limited by the numerical interpolation. In the following discussion, species, for which both methods predict relative abundances less than , are excluded from further consideration. In Figures 2 and  3, contours labeled with 0.9 mean that the two methods give abundances that differ by less than an order of magnitude for 90% of species with abundances higher than in both methods.

An order-of-magnitude agreement criterion may seem to be too coarse. However, when stochastic effects do not play an important role, all three method produce results, which are nearly identical. To illustrate this, we show in Figure 4 the evolution of the HCCCC abundance at  cm and  K, simulated with the MC and RE methods. The species is chosen because its abundance strongly varies with time. It can be seen that the difference between the two methods is very small.

4.1.1 Model T

In Model T the species are highly mobile—atomic and molecular hydrogen because of tunneling, and all the other species because of the low energy barrier for diffusion. Mobility drives rapid surface reactions with rates which sometimes exceed the accretion/desorption rate. Surface reactions, which are on average faster than accretion, are not permitted in the MC approach when surface populations of their reagents are low ( 1), but they are allowed in the RE approach, where the surface abundance of a species can be much less than one. Because of that, differences between Monte Carlo and rate equation methods are quite noticeable in Model T.

At temperatures below 20 K and above  K the agreement between the stochastic method and the rate equation method (both RE and MRE) is about 85% or better at all times and densities (Figure 2). Around  K the residual discrepancy is mainly caused by complex surface species, which have zero abundances in Model T-MC and are overproduced in Model T-RE and Model T-MRE, so that their abundances are just above the adopted cutoff of . If we would raise the cutoff to , the agreement would be almost 100%. At  K the agreement is also nearly perfect, as the chemistry is almost a purely gas-phase chemistry under these conditions (dust temperature being equal to gas temperature).

All the major discrepancies are concentrated in the temperature range between 25–30 K. In this range the accretion rate (which depends on gas temperature) is high enough to allow accumulation and some processing of surface species, while the correspondingly high desorption rate (which depends on dust temperature and ) precludes surface production of complex molecules. The latter process is adequately described in the MC runs only.

At earlier times ( years) the largest differences in this temperature range are observed at lower densities. They are caused by an overproduction of ‘terminal’ surface species, like water, ammonia, hydrogen peroxide, carbon dioxide, etc., in the RE and MRE calculations. At later times ( years) situation changes drastically. The discrepancies shift toward higher density, and the overall agreement falls off to about 60% for both RE and MRE calculations. Even though the stochastic chemistry is only important on dust surfaces, by this time it also influences many gas-phase abundances due to effective gas-grain interactions.

At  years, in the RE calculation at  K, among 234 species with abundances higher than , 95 species (including 81 gas-phase species) have an order of magnitude disagreement with the MC run. In the MRE run at  K, the number of species with abundances above is 221. Among them 83 species disagree with the MC results by more than an order of magnitude. In both calculations these are mostly carbon-bearing species, like CO, CO, CS, to name a few, and, in particular, carbon chains (like some cyanopolyynes). It is noteworthy that not only the abundances of neutral species disagree, but also the abundances of some ions, including C and S.

4.1.2 Model H

In Model H the mobility of all species on the surface (including H and H) is caused by thermal hopping only, which is slow because of the high ratio. One of the consequences of these slow surface reaction rates is that results produced by the modified rate equations almost do not differ from results of conventional rate equations. The essence of the modification is to artificially slow down surface reactions to make their rates consistent with accretion/desorption rates. In Model H all surface reactions are slow anyway, so that modifications never occur. Thus, we only compare the results of the MC and RE calculations.

At years the percentage agreement between the two runs never falls below 80% and most often is actually much better than the adopted order-of-magnitude criterion. At later times, the discrepancy appears in the same temperature range as in Model T, i.e., between 25 and 30 K. Unlike Model T, both at  years and  years the disagreement shows up at low density. The set of discrepant species is nearly identical at both times and consist mostly of ices, overproduced by rate equations. These ices include some key observed species. For example, NH, HO, CO ices are all abundant in the Model H-RE despite the low density, with a surface HO abundance reaching by years. On the other hand, the Model H-MC results in very small or zero abundances for the same species. More chemically rich ices in Model H-RE bring about enhanced gas-phase abundances for the same molecules. The only species which is underproduced by the rate equation approach is surface CO.

4.2 Selected Species

The diagrams presented above indicate areas in the parameter space, relevant to molecular clouds, where deterministic methods fail to describe stochastic surface processes. However, a global agreement does not necessarily imply that the abundances of key species are also correctly calculated. In the following we discuss a few important gas-phase species and ices, for which an order of magnitude (or more) disagreement has been found.

4.2.1 Model T

Because surface reaction rates are high in this model, there are many species for which stochastic and deterministic methods give quite different results. The most abundant trace molecule, CO, is among these species. In Figure 5, we show agreement diagrams for Models T-RE and T-MRE (top row) and the CO abundance evolution for combinations corresponding to the worst agreement. Both the RE and MRE methods fail to reproduce the late CO evolution at high density and temperatures of 25–30 K (where the global agreement is also worst in this model), but each in a different way. While in Model T-RE CO is underabundant by an order of magnitude, in Model T-MRE it is overabundant by the similar amount. This difference is related to the treatment of CO CO conversion on dust surfaces, where underabundance of the CO ice in Model T-RE and overabundance in Model T-MRE are observed at almost all times. At later time, these surface abundances just start to propagate to gas-phase abundances.

Overproduction of the CO ice in Model T-RE also consumes carbon atoms which would otherwise be available not only for CO, but also for more complex molecules, in particular, carbon chains, starting from C, observed in diffuse clouds, and CS, used as a diagnostics in prestellar cores. In the RE and MRE calculations these molecules, like CO, show trends opposite to that of CO ice: when surface CO is overabundant with respect to the MC model, carbon chains are underabundant and vice versa.

When CO molecule sticks to a grain, it either desorbs back to the gas-phase, where it may participate in further processing, or is converted into CO. As the desorption energy is quite high for CO, it acts as a sink for carbon atoms at these temperatures. In the accretion-limited regime, rate equations overproduce CO with respect to the MC method. Thus, we see less CO and other carbon-bearing species in the gas-phase. The MRE method helps, probably, too much to account for accretion-limited CO processing and essentially quenches surface CO CO conversion. It is worth noting that at high density in this particular temperature range (and for adopted CO desorption energy) CO balances between complete freeze-out (at  K) and near absence on dust surface (at  K). Thus, even relatively minor changes in treatment may lead to noticeable consequences.

Another carbon-bearing molecule that shows an interesting behavior is methanol in the gas phase. Diagrams, comparing its abundance computed with the MC approach and with the rate equations, are presented in Figure 6. The agreement between results of the MC calculations and the RE calculations is almost perfect, at least, within the scatter produced by the MC simulations. However, the modified rate equations, which tend to improve agreement between stochastic and deterministic methods for many other species, in this particular case underpredict methanol abundance by more than two orders of magnitude in comparison with the MC run. The same is true for surface methanol as well. It is interesting to note that the region in the parameter space, where MC and MRE methanol abundances disagree, extends down to lowest temperatures considered.

An interesting pattern is presented by the abundance of water, ammonia and their ices (Figure 7). These species are severely overproduced by the RE methods at  years, but after this time the modified rate equations are able to restore the agreement with the MC method, so that by years difference between MC and MRE runs is not significant. This is an example of a situation where the MRE help to produce realistic results, at least, at later times. We show plots only for gas-phase abundances because these are what is really observed. However, these abundances are only a reflection of surface processes, and surface water and ammonia abundances behave in exactly the same way. It is hard to name a single reaction which is responsible for the difference between results of RE and MRE methods in this case. Let’s consider water as its chemistry is somewhat simpler. The primary formation reaction for water is H + OH, with only a minor contribution from HO (see below). However, hydroxyl is not produced in H + O reaction, as one might have expected. Because of the paucity of H atoms on the grain surface, hydroxyl formation is dominated by reaction

Abundance of HNO is restored in reaction

These two reactions form a semi-closed loop for which the sole result is OH synthesis out of an H atom and an O atom. Evolution of surface O abundance seems to be the key to the difference between the RE and MRE models. In the RE model the number of O atoms on a surface is determined by the balance between their accretion from the gas phase and consumption in O + HNO reaction. The same processes define O abundance in the MRE model as well, but only at years. After this time O+HNO reaction slows down significantly due to modifying correction, and O abundance is controlled purely by accretion and desorption thereafter (as it is in the MC model at all times). As a result, OH and HO abundances decrease, getting closer to the MC model prediction. A similar mechanism is at work in the case of ammonia as well.

4.2.2 Model H

The only parameter region where a discrepancy occurs for Model H is located at  cm and  K. Because of the low density, the degree of molecular complexity at these conditions is also low, and by years only 73 species have abundances higher than at least in one of the calculations. Of these 73 species, 25 show disagreement between the MC calculations and the RE calculations.

In the top row of Figure 8 we show the low density evolution of gas-phase water and molecular oxygen abundances. Both plots show a similar trend. In the MC calculations the abundance of the molecule stays nearly constant. At earlier times this steady-state behavior is reproduced in the RE run as well, but later the abundance grows and exceeds the MC abundance by more than an order of magnitude. It is interesting to note that the average O abundance in the MC run seems to decrease somewhat at later times, and this trend is reproduced by the RE model. The behavior of gas-phase abundances is a reflection of grain-surface abundances. It must be kept in mind that surface water abundance is quite high in the H-RE model, reaching almost . Even though evaporation rate is not very high, it does increase gas-phase abundance up to a level of a few time . Also, it is a low density model, so desorption is primarily caused by photons.

In Figure 9 we show the evolution of HO surface abundance and O gas-phase abundance in Model T-RE at  cm and the gas-phase ammonia abundance in Model H-RE at  cm. These plots demonstrate that disagreement does not always mean simply higher or lower abundances. Different processes determine the HO abundance at different times, and while the RE method is able to capture some of these processes, others are obviously missed. Our analysis shows that this behavior is related to the treatment of reaction with H atoms. Specifically, surface abundance of HO is controlled by three reactions only. It is produced in a reaction H + OH and destroyed by atomic hydrogen in reactions producing either HO + OH or OH + H. (Note that in the considered model water ice is almost irreversible sink for O atoms as there are no surface reactions destroying water, and desorption is negligible.) The formation reaction is fast as it does not have an activation barrier. So, initially, as H and O atoms start to stick to a grain, surface abundance of HO steadily grows. At the same time, there is not enough H atoms for either of destruction reactions. Then, in the MC model at around 1300 years the first destruction reaction (H + HO HO + OH), which has a slightly lower activation barrier than the other one (H + HO OH + H), starts to “steal” some hydrogen atoms from the formation reaction, irreversibly removing one O atom per reaction from the HO (re)formation process. This is the reason for the sharp fall-off of HO abundance. On the other hand, in the RE model both destruction and formation reactions may occur simultaneously, which prevents the first destruction reaction from having such a dramatic effect. Note that the second reaction restores an OH molecule which may react again with an H atom to restore an HO molecule.

HO is not unique in this kind of behavior. There are some other surface species which show more or less similar trends, that is, almost perfect agreement at some times and quite noticeable disagreement at other times. These are some carbon chains (CH, CN) or simpler molecules (CS, NO). Another equally dramatic example is represented by CH. The mechanism is similar, that is, it is related to the sequence of hydrogen additions which is treated differently in the MC model and the RE model.

Evolution of the O abundance computed with the RE method follows the evolution computed with the MC model, but most of the time the RE curve is somewhat higher than the MC curve. Finally, the plot for ammonia shows the situation when the RE method correctly describes the average abundance evolution. However, the abundance predicted by the MC method fluctuates so wildly, that at each particular we can detect disagreement with high probability, using our formal order-of-magnitude criterion.

5 Discussion

In the previous sections we investigated the validity of the (modified) rate equation approach to grain surface chemistry under different physical conditions encompassing diffuse clouds, giant molecular clouds, infrared dark clouds etc. Evolution of the medium is studied with an extended gas-grain chemical model over a time period of 10 years. We utilize the unified stochastic Monte Carlo approach, applied simultaneously to gas-phase and grain-surface reactions. The results are used to test the validity of conventional deterministic approaches. In general, differences in results obtained with deterministic and stochastic methods strongly depend on the adopted microphysical model of surface chemistry. In Model T, where tunneling for atomic and molecular hydrogen is permitted and diffusion/desorption energies ratio is low (), discrepancy both for RE and MRE methods is very significant. At low temperatures (10 - 20 K) only abundances of surface species are discrepant by more than an order of magnitude. At moderate temperatures (25–40 K), due to active gas-grain interaction, incorrect treatment of grain surface chemistry becomes important for gas-phase abundances, too, and leads to dramatic decrease of overall agreement. At even higher temperatures agreement becomes better again due to limited surface chemistry. In Model H with no tunneling and high ratio, surface species are much less mobile, and the agreement between deterministic and stochastic methods is better. The only exception is low density ( cm) with moderate temperatures (25 – 30 K), where average residence time of species on grain surface is still lower than the average interval between reaction events. In general, the stochastic chemistry severely affects abundances and, thus, must be taken into account in chemical models of a warm and moderately dense medium.

Modifications in rate equations, aimed at taking into account the accretion-limited regime, do not provide a significant improvement over ‘canonical’ rate equations in Model T. Close inspection of Figure 2 shows that at earlier times ( years) modified rate equations mostly produce results, which are rather less consistent with the ‘exact’ Monte Carlo solution, than results of the RE method. Only at later times ( years) results of the MRE method become more consistent with the ‘exact’ solution than the RE results. In the previous section, we have already shown some examples of species, for which the MRE method actually worsens the agreement.

Because at many of considered physical conditions results obtained with Monte Carlo method deviate significantly from those obtained with (modified) rate equations, the natural question to ask is whether the inability of the RE method to treat these combinations of parameters might have caused mis-interpretation of observational data. To answer this question, one would need to consider an object, for which a large volume of observational data is available, so that simultaneous comparison for many molecules is possible. So far, there seems to be just one molecular cloud, which is studied in necessary details. This is a sub-region of the TMC1 cloud where cyanopolyynes reach high concentrations, the so-called TMC1-CP. The physical conditions in this object ( K,  cm) do fall into the range adopted in our model (Benson & Myers, 1989). However, the TMC1-CP is probably too cold to show significant dependence of gas-phase abundances from surface processes.

We may expect that the disagreement will be more significant in hotter objects. The range of physical conditions studied in this paper is quite narrow, as it has only been chosen to provide an initial view on the importance of stochastic treatment for surface chemistry. There are two directions which are worth to follow for further study. One of them is related to diffuse clouds. In this paper we most of the time assumed initial conditions in which all the hydrogen is molecular even at low density. Another simplification, which is, strictly speaking, not valid in low density medium, is equality of dust and gas temperature. So, it is interesting to check if stochastic effects work in predominantly atomic gas with different gas and dust temperatures. Another direction is chemistry at later stages of star formation, that is, in hot cores, hyper- and ultracompact HII regions, protoplanetary disks. In particular, in dark cold disk midplanes, which are poor in atomic hydrogen, stochastic effects may severely influence surface synthesis of hydrogenated species, like formaldehyde, which are later transferred by turbulent mixing in warmer disk areas and desorbed into gas phase (Aikawa, 2007). In these circumstances one would definitely need to use more accurate methods to model surface chemistry. Stochastic effects may be even more pronounced in the dynamically evolving medium, like during slow warm-up phase in hot cores or “corinos”. Complex (organic) molecules will be among the most affected species.

While the Monte Carlo method is the most direct approach, it is by no means practical. So, other methods are to be developed, in particular, the method of moment equations, if we want to understand deeper the chemical evolution in warm and dense astrophysical objects.

6 Conclusions

In this study, a gas-grain chemical model is presented, consisting of about 600 species and 6000 reactions with stochastic description of grain surface chemistry. For the first time both gas phase and grain surface reactions are simulated with a unified Monte Carlo approach. This unified model is used to test the validity of rate equations and modified rate equations for a set of physical conditions, relevant to translucent clouds, dark dense cores, and infrared dark clouds.

A comparison of results obtained with deterministic rate equations (RE) and modified rate equations (MRE) with results of Monte Carlo simulations was performed for this set of physical conditions and two different models of surface chemistry. We found that results obtained with both RE and MRE approaches sometimes significantly deviate from those obtained with Monte Carlo calculations in Model T, where surface species are highly mobile and most of the grain chemistry occurs in the accretion limited regime. While at low temperatures (10 K – 20 K) mainly RE and MRE abundances of surface species show significant discrepancies, at moderately high temperatures (25 K – 30 K) even abundant gas-phase species like CO considerably deviate from the Monte Carlo results. At these temperatures grain surface chemistry is still very active and at the same time processes of accretion and desorption are very efficient. The gas-phase abundances of many species are influenced by surface chemistry that is correctly described in Monte Carlo simulations only. In Model H surface species are much less mobile than in Model T and the agreement between the (modified) rate equations approach and Monte Carlo simulations is almost perfect. The only parameter region where this discrepancy is significant is at a temperature of  K and at a low density cm. On the whole, results of modified rate equations seem to be in somewhat closer agreement with Monte Carlo simulations than results of the conventional rate approach. But still, the agreement is far from being perfect, so modified rates should be used with care. We conclude that stochastic effects have a significant impact on chemical evolution of moderately warm medium and must be properly taken into account in corresponding theoretical models. In the studied parameter range, stochastic effects are, in general, not important at low ( K) and high ( K) temperatures.

Part of this work was supported by the Russian Dynasty foundation. DW is supported by the RFBR grant 07-02-01031. Authors are grateful to Sergey Koposov for valuable and fruitful discussions and assistance in data analysis software development. Authors are thankful to the anonymous referee for valuable comments and suggestions. This research has made use of NASA’s Astrophysics Data System.


  • Aikawa (2007) Aikawa, Y. 2007, ApJ, 656, L93
  • Barzel & Biham (2007) Barzel, B. & Biham, O. 2007, ApJ, 658, L37
  • Benson & Myers (1989) Benson, P. J. & Myers, P. C. 1989, ApJS, 71, 89
  • Biham et al. (2001) Biham, O., Furman, I., Pirronello, V., & Vidali, G. 2001, ApJ, 553, 595
  • Boger & Sternberg (2006) Boger, G. I. & Sternberg, A. 2006, ApJ, 645, 314
  • Caselli et al. (1998) Caselli, P., Hasegawa, T. I., & Herbst, E. 1998, ApJ, 495, 309
  • Caselli et al. (2002) Caselli, P., Stantcheva, T., Shalabiea, O., Shematovich, V. I., & Herbst, E. 2002, Planet. Space Sci., 50, 1257
  • Cazaux & Tielens (2004) Cazaux, S. & Tielens, A. G. G. M. 2004, ApJ, 604, 222
  • Chang et al. (2005) Chang, Q., Cuppen, H. M., & Herbst, E. 2005, A&A, 434, 599
  • Chang et al. (2007) —. 2007, A&A, 469, 973
  • Charnley (1998) Charnley, S. B. 1998, ApJ, 509, L121
  • Charnley (2001) —. 2001, ApJ, 562, L99
  • Cuppen & Herbst (2007) Cuppen, H. M. & Herbst, E. 2007, ApJ, 668, 294
  • Draine & Lee (1984) Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89
  • Garrod & Herbst (2006) Garrod, R. T. & Herbst, E. 2006, A&A, 457, 927
  • Gillespie (1976) Gillespie, D. T. 1976, Journal of Computational Physics, 22, 403
  • Green et al. (2001) Green, N. J. B., Toniazzo, T., Pilling, M. J., Ruffle, D. P., Bell, N., & Hartquist, T. W. 2001, A&A, 375, 1111
  • Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167
  • Hassel et al. (2008) Hassel, G. E., Herbst, E., & Garrod, R. T. 2008, ApJ, 681, 1385
  • Herbst & Shematovich (2003) Herbst, E. & Shematovich, V. I. 2003, Ap&SS, 285, 725
  • Katz et al. (1999) Katz, N., Furman, I., Biham, O., Pirronello, V., & Vidali, G. 1999, ApJ, 522, 305
  • Le Bourlot et al. (1993) Le Bourlot, J., Pineau des Forets, G., Roueff, E., & Schilke, P. 1993, ApJ, 416, L87+
  • Lee et al. (1996) Lee, H.-H., Herbst, E., Pineau des Forets, G., Roueff, E., & Le Bourlot, J. 1996, A&A, 311, 690
  • Lee et al. (1998) Lee, H.-H., Roueff, E., Pineau des Forets, G., Shalabiea, O. M., Terzieva, R., & Herbst, E. 1998, A&A, 334, 1047
  • Lipshtat & Biham (2003) Lipshtat, A. & Biham, O. 2003, A&A, 400, 585
  • Lipshtat et al. (2004) Lipshtat, A., Biham, O., & Herbst, E. 2004, MNRAS, 348, 1055
  • Perets et al. (2005) Perets, H. B., Biham, O., Manicó, G., Pirronello, V., Roser, J., Swords, S., & Vidali, G. 2005, ApJ, 627, 850
  • Pirronello et al. (1997a) Pirronello, V., Biham, O., Liu, C., Shen, L., & Vidali, G. 1997a, ApJ, 483, L131+
  • Pirronello et al. (1999) Pirronello, V., Liu, C., Roser, J. E., & Vidali, G. 1999, A&A, 344, 681
  • Pirronello et al. (1997b) Pirronello, V., Liu, C., Shen, L., & Vidali, G. 1997b, ApJ, 475, L69+
  • Roberts et al. (2004) Roberts, H., Herbst, E., & Millar, T. J. 2004, A&A, 424, 905
  • Ruffle & Herbst (2000) Ruffle, D. P. & Herbst, E. 2000, MNRAS, 319, 837
  • Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
  • Semenov et al. (2004) Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93
  • Shalabiea & Greenberg (1995) Shalabiea, O. M. & Greenberg, J. M. 1995, A&A, 296, 779
  • Snow & McCall (2006) Snow, T. P. & McCall, B. J. 2006, ARA&A, 44, 367
  • Spitzer & Tomasko (1968) Spitzer, L. J. & Tomasko, M. G. 1968, ApJ, 152, 971
  • Stantcheva et al. (2001) Stantcheva, T., Caselli, P., & Herbst, E. 2001, A&A, 375, 673
  • Stantcheva & Herbst (2003) Stantcheva, T. & Herbst, E. 2003, MNRAS, 340, 983
  • Stantcheva & Herbst (2004) —. 2004, A&A, 423, 241
  • Stantcheva et al. (2002) Stantcheva, T., Shematovich, V. I., & Herbst, E. 2002, A&A, 391, 1069
  • Tielens & Hagen (1982) Tielens, A. G. G. M. & Hagen, W. 1982, A&A, 114, 245
  • van Dishoeck et al. (2006) van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, in Faraday Discussions 133: Chemical Evolution of the Universe, ed. I. R. Sims & D. A. Williams, 231–243
  • van Zadelhoff et al. (2003) van Zadelhoff, G.-J., Aikawa, Y., Hogerheijde, M. R., & van Dishoeck, E. F. 2003, A&A, 397, 789
  • Vasyunin et al. (2008) Vasyunin, A. I., Semenov, D., Henning, T., Wakelam, V., Herbst, E., & Sobolev, A. M. 2008, ApJ, 672, 629
  • Woodall et al. (2007) Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197
Model Based on Tunneling
T Hasegawa et al. (1992) Yes 0.3
H Katz et al. (1999) No 0.77
Table 1: Adopted surface chemistry models
Figure 1: The Monte Carlo time step as a function of time in the model of the chemical evolution of a dense cloud.
Figure 2: Global agreement diagrams for Model T at three time moments. Percentage of species is shown in contour labels, for which abundances computed with the rate equations and with the Monte Carlo method differ by no more than an order of magnitude. Grayscale map with arbitrary contours is added for clarity. Darker color corresponds to worse agreement. Results are presented for conventional rate equations (top row) and for modified rate equations (bottom row).
Figure 3: Same as in Figure 2, but for Model H.
Figure 4: Evolution of HCCCC gas-phase abundance at  cm and  K, simulated with the MC and RE methods (Model T). Dotted line indicates raw MC results, solid line shows MC results after averaging. Results of RE integration are shown with dashed line.
Figure 5: (Top row) Agreement diagrams for CO in Model T. Contour labels correspond to the logarithm of the ratio of gas-phase CO abundance in the MC run to that in the RE (top left) and MRE (top right) runs. (Bottom row) Gas-phase CO abundance evolution in the RE run ( cm and  K; bottom left) and in the MRE run ( cm and  K; bottom right).
Figure 6: (Top row) Agreement diagrams for methanol. Contour labels correspond to the logarithm of the ratio of gas-phase CHOH abundance in the MC run to that in the RE (top left) and MRE (top right) runs. (Bottom row) Gas-phase methanol abundance evolution in the RE run ( cm and  K; bottom left) and in the MRE run ( cm and  K; bottom right).
Figure 7: Water and ammonia gas-phase abundance evolution in the RE run in Model T ( cm and  K; left) and in the MRE run ( cm and  K; right).
Figure 8: Water and molecular oxygen abundances at  cm and  K in Model H.
Figure 9: Abundance evolution for selected species. Models and physical conditions are indicated on the diagrams.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description