Free Access
Volume 530, June 2011
Article Number A131
Number of page(s) 11
Section Numerical methods and codes
Published online 25 May 2011

© ESO, 2011

1. Introduction

The chemistry of the interstellar medium can be roughly divided into two types: gas phase chemistry and grain surface chemistry. The two types of chemistry are coupled by the adsorption and desorption processes. Species adsorbed on the grain surface migrate in a random walk manner, and they may react with each other upon encounter at the same site (a local potential minimum). The products can be released back to the gas phase through certain desorption mechanisms. In addition to the gas phase chemistry, grain chemistry is important for the material and energy budget of the interstellar medium. For example, besides H2, molecules such as methanol are believed to be formed on the grain surfaces (Garrod et al. 2007), because its relatively high abundance (see, e.g., Menten et al. 1988) cannot be reproduced by gas phase chemistry.

Several methods have been used to model the gas-grain chemistry. In the rate equation (RE) approach (see, e.g., Hasegawa et al. 1992), the surface processes are treated the same way as the gas phase processes. This works fine when the number of reactants on a single grain is large (under the assumption that the system is well stirred; see Gillespie 2007), but might not be accurate enough when the average populations1 of some reactants on a single grain is small. This failure of the rate equation is related to the treatment of two-body reaction. For the REs to be applicable, the probability of one reactant being present should be independent of that of another being present. This is not always true, especially when the average populations of both reactants are low, in which case they might be highly correlated. The flaws in employing the RE for grain-surface chemistry were pointed out by Charnley et al. (1997) and Tielens & Charnley (1997).

To remedy this problem, modification schemes based on some empirical, heuristic, and/or physical reasoning have been applied to the RE approach (Caselli et al. 1998; Stantcheva et al. 2001), and are called modified rate equation (MRE) approach. The validity of this method has been questioned (Rae et al. 2003). A modification scheme developed by Garrod (2008) uses different functional forms for different surface populations, taking various competition processes and refinements into account. It has been shown to work very well, even for very large reaction networks (Garrod et al. 2009).

Mathematically, the gas-grain system should be viewed as a stochastic chemical system (see, e.g., McQuarrie 1967; Gillespie 1976; Charnley 1998), being described by a probability distribution P(x,t), which is the probability that the system has a population vector x at time t, with xi being the number of the ith species in the system. The evolution equation of P(x,t) is the so-called master equation, whose form is determined by the reaction network.

Many sophisticated methods have been proposed (mainly outside the astro-chemical community; see, e.g., the operator method described in Mattis & Glasser 1998; or the variational approach used by Ohkubo 2008) to solve the master equation. However, these methods work fine only when either the chemical network is small or some special assumptions are made in the derivation, thus their validity in the general case should be questioned. It is unclear whether these methods can be generalized to large complex networks.

The numerical solution of the master equation has also been performed (Biham et al. 2001; Stantcheva et al. 2002; Stantcheva & Herbst 2004). To limit the number of variables in the set of differential equations and to separate the deterministic and stochastic species, usually a priori knowledge of the system is required in these studies. The steady state solution of the master equation can also be obtained analytically in some very simple cases, such as the formation of H2 molecules on the grain surface (Green et al. 2001; Biham & Lipshtat 2002).

On the other hand, the master equation prescription can be “realized” through a stochastic simulation algorithm (SSA), proposed by Gillespie 1976 (see also Gillespie 1977, 2007). In this approach, the waiting time for the next reaction to occur, as well as which specific reaction will occur are random variables that are completely determined by the master equation, so this approach should be considered the most accurate. In principle, multiple runs are needed to average out the random fluctuations, but in practice this is unnecessary if one only cares about the abundant species. This approach has been applied successfully to astrochemical problems (Charnley 1998, 2001; Vasyunin et al. 2009), even in the case of very large networks (Vasyunin et al. 2009). Besides providing results that are accurate, this approach is very easy to implement. However, it requires a very long run time for large networks if a long evolution track is to be followed, although some approximate accelerated methods do exist (e.g. Gillespie 2000).

The SSA described above is somewhat different from a Monte Carlo (MC) approach which has also been applied to astrochemistry (e.g. Tielens & Hagen 1982); however, this approach is not rigorously consistent with the master equation (see the comment by Charnley 2005), and can lead to a reaction probability higher than 1 (Awad et al. 2005) in certain cases. The nomenclature of these two approaches is not always consistent in the astrochemical literature2. For example, the SSA used by Vasyunin et al. (2009) is called the Monte Carlo approach in their paper. Hereafter, we use the term “Monte Carlo” when referring to the rigorous stochastic simulation approach of Gillespie.

By taking various moments of the master equation, the so-called moment equation (ME) is obtained (Lipshtat & Biham 2003; Barzel & Biham 2007a,b). This set of equations describes the evolution of both the average population of each species and the average value of the products of the population of a group of species, usually cut off at the second order moments. Its formulation is similar to that of the RE, so it is relatively easy to implement. Furthermore, in this approach the gas phase chemistry and grain surface chemistry can be coupled together naturally. It has been tested on small surface networks.

In the present paper, we propose yet another approach to modeling gas-grain chemistry, named the hybrid moment equation (HME) approach. The goal is to find a systematic, automatic, and fast way to modeling gas-grain chemistry as accurately as possible. Our method is based on the ME approach. Different approximations are applied to the MEs at different time depending on the overall populations at that specific time. It is hybrid in the sense that the RE and the ME are combined together. The basic modification and competition scheme presented in Garrod (2008) can be viewed as a semi-steady-state approximation to our approach (by assuming that the time derivatives of certain second order moments are equal to zero), while our approach can also be viewed as a combination of the ME approach of Barzel & Biham (2007a) and the RE. In our approach, the MEs are generated automatically with the generating function technique, and in principle MEs up to any order can be obtained this way. We benchmark our approach against the exact MC approach (i.e. the SSA of Gillespie).

The remaining part of this paper is organized as follows. In Sect. 2, we review the chemical master equation and ME, then describe the main steps of the HME approach. In Sect. 3, we benchmark the HME approach with a cutoff at the second order and the RE approach against the MC approach with a large gas-grain network; we also tested the HME approach with a cutoff at the third order with a small network. In Sect. 4, we discuss the performance of the HME, and its relation with previous approaches, as well as possibilities for additional improvements. Our way of generating the MEs is described in Appendix A. A surface chemical network we used for benchmark is listed in Appendix B.

2. Description of the hybrid moment equation (HME) approach

In this section, we first review both the chemical master and moment equations. Although this content can be found in many other papers (e.g., Charnley 1998; Gillespie 2007), we present them here as they are the basis of our HME approach. We then describe the MEs and REs for a simple set of reactions as an example, to demonstrate how the HME approach naturally arise as a combination of ME and RE. Finally we show the main steps of the HME approach.

2.1. The chemical master equation and the moment equation (ME)

A chemical system at a given time t can be described by a state vector x which changes with time, with its jth component xj being the number of the jth species in this system. As a chemical system is usually stochastic, x should be viewed as a random variable, whose probability distribution function P(x,t) evolves with time according to the master equation (Gillespie 2007) (1)where ai(x) is called the propensity function, ai(xt is the probability that given a current state vector x an ith reaction will happen in the next infinitesimal time interval Δt, and νi is the stoichiometry vector of the ith reaction. The sum is over all the reactions, and M is the total number of reactions.

The ME is derived by taking moments of the master equation. For example, for the first order moment  ⟨xj⟩ , which is simply the average number of species j,  ⟨xj⟩  ≡  ∑ xP(x,t)xj, its evolution is determined by (Gillespie 2007) (2)where νij is the jth component of the stoichiometry vector of the ith reaction, i.e. the number of jth species produced (negative when being consumed) by the ith reaction. For higher order moments, their corresponding evolution equations can be similarly derived, although the final form will be more complex. In Appendix A, we present another method based on the generating function technique to derive the MEs, which is more suitable for programming.

For the simplest network, in which all the reactions are single-body reactions, ai(x) is a linear function of x. In this case the ME is closed and can easily be solved. However, when two-body reactions are present, this is no longer true, as  ⟨ai(x)⟩  might be of a form  ⟨xk(xk − 1)⟩  or  ⟨xkxl⟩ , which is of order two and cannot be determined in general by the lower order moments. Hence additional equations governing their evolution should be included, i.e., they should be taken to be independent variables. The evolution equation of these second order moments may also involve moments of order three, and this process continues without an end, thus the ME is actually an infinite set of coupled equations (although in principle they are not completely independent if the chemical system being considered is finite, which leads to a finite-dimensional space of state vectors). The equation cannot be solved without a compromise, e.g., a cutoff procedure, except for the simplest cases in which an analytical solution is obtainable in the steady state.

2.2. The MEs and REs for a set of reactions

We take the following symbolic reactions as an illustrative example where the ks are the reaction rates of each reaction, A–E are assumed to be surface species that are distinct from each other, and “a” is the gas phase counterpart of A.

In the following we first write down the MEs and REs for this system, then discuss the relations and differences between them, as well as the relation between a cutoff of MEs and a cutoff of master equations in previous studies. These discussions will be essential to developing our HME approach.

2.2.1. The MEs for this system

The propensity functions for the above four reactions are kada, kevapA, kABAB, and kAAA(A − 1), respectively. Here for convenience we use the letter “A” to represent both the name of a species and the population of the corresponding species.

For the first order moments, we have Other similar equations are omitted. The symbol  ⟨ ∗ ⟩  is used to represent the average population of “*” in the system; the average should be understood as an ensemble average. The second order moments  ⟨AB⟩  and  ⟨A(A − 1)⟩  have their own evolution equations, which are (10)(11)For this simple example set of reactions (Eqs. (36) ), the above equations can be easily obtained from the master equation (see, e.g., Lipshtat & Biham 2003, p. 8). In the general case (e.g., when A–E are not completely distinct from each other), an automatic way of obtaining the MEs is described in Appendix A. The method described there is also applicable to moments with any order, and to all the common reaction types in astrochemistry.

In general, the third order moments in the above equations cannot be expressed as a function of the lower order moments, so they need their own differential equations. In the case of a cutoff at the second order, the chain of equations, however, stops here. We describe the method required to evaluate them in Sect. 2.3.

2.2.2. The REs for this system

When using REs, Eqs. (711) are replaced by The equations for  ⟨A⟩  ⟨B⟩  and  ⟨A⟩ 2 are of course not needed in the RE approach but are simply derived from Eq. (7′) (and an omitted similar equation for  ⟨B⟩ ) using the chain rule of calculus.

2.2.3. The relation between MEs and REs

The differences between the MEs (Eqs. (711)) and the REs (Eqs. (7′–11′)) in the present case are as follows: all the  ⟨AB⟩  are replaced by  ⟨A⟩  ⟨B⟩ , all the  ⟨A(A − 1)⟩  are replaced by  ⟨A⟩ 2, all the  ⟨A(A − 1)B⟩  are replaced by  ⟨A⟩ 2 ⟨B⟩ , the  ⟨AB(B − 1)⟩  is replaced by  ⟨A⟩  ⟨B⟩ 2, and the  ⟨A(A − 1)(A − 2)⟩  is replaced by  ⟨A⟩ 3. Furthermore, the term kAB ⟨AB⟩  in Eq. (10) and the term kAA ⟨A(A − 1)⟩  in Eq. (11) disappear in the RE (10′) and (11′).

These differences make clear why the REs are accurate when the involved species are abundant (namely when  ⟨A⟩     ≫    1 and  ⟨B⟩     ≫    1). This is because, in this case,  ⟨AB⟩  can be approximated well by 3A⟩  ⟨B⟩ , and  ⟨A(A − 1)⟩  can be approximated well by  ⟨A⟩ 2.

The RE approach will be erroneous when  ⟨A⟩  or  ⟨B⟩  are smaller than 1 because, in this case, the correlation between A and B might cause  ⟨AB⟩  to differ considerably from  ⟨A⟩  ⟨B⟩ , and the fluctuation in A might cause  ⟨A(A − 1)⟩  to differ considerably from  ⟨A⟩ 2. It can also be viewed like this: in Eqs. (10′) and (11′) that govern the evolution of second order moments, the omitted term kAB ⟨AB⟩  might be much larger than the retained terms such as  ⟨A⟩ 2 ⟨B⟩  or  ⟨A⟩  ⟨B⟩ 2, and the omitted term kAA ⟨A(A − 1)⟩  might be much larger than the retained term  ⟨A⟩ 3.

2.2.4. The relation between a cutoff of MEs and a cutoff of possible states in previous master equation approaches

In Eqs. (711) we do not write terms such as  ⟨A(A − 1)⟩  in the split form  ⟨A2⟩  −  ⟨A⟩ . We keep terms such as  ⟨A(A − 1)B⟩  and  ⟨A(A − 1)(A − 2)⟩  in their present forms intentionally. One reason for this is that terms such as  ⟨A(A − 1)⟩  look more succinct and follow naturally from our way of deriving them (see Appendix A). When  ⟨A⟩  ≫ 1,  ⟨A(A − 1)⟩  and  ⟨AB⟩  can be directly replaced by  ⟨A⟩ 2 and  ⟨A⟩  ⟨B⟩ , respectively, to obtain the RE formulation.

More importantly, this formulation can be directly connected to the cutoff schemes in the previous master equation approaches (e.g., Biham et al. 2001; Stantcheva et al. 2002). For example, in a scheme in which no more than two particles of A are expected to be present on a single grain at the same time, we have P(A    >    2) = 0. In this case, . Thus we see that a cutoff at a population of two in the master equation approach corresponds naturally to assigning a zero value to moments containing A more than twice, as far as the moments are defined in the form presented above.

2.3. The HME approach

thumbnail Fig. 1

A flow chart of the main components of our HME code. Steps 2, 5, 6, and 8 are described in detail in the text.

The HME approach is a combination of the ME and RE approaches. The basic idea is that, for deterministic (average population  > 1) species, the REs are used, while for stochastic (average population  < 1) species, the stochastic effects are taken into account by including higher order moments in the equations. Since a deterministic species may become stochastic as time goes by, and vice versa, the set of differential equations governing the evolution of the system also changes with time, and is determined dynamically. A flow chart of our HME code is shown in Fig. 1.

We first set up all potentially needed MEs (using the procedure described in Appendix A), with a cutoff of moments at a prescribed order (usually two). After this and some other initialization work, the program enters the main loop.

The main loop contains an ODE solver because the system of MEs is a set of ordinary differential equations (ODEs). We use the solver from the ODEPACK package4.

Not all MEs and moments are used at all times; which ones are used is determined dynamically. In each iteration of the main loop, we verify whether some surface species have changed from stochastic to deterministic, or from deterministic to stochastic. The gas phase species are always treated as deterministic, regardless of how small their average populations are. In either of these two cases, we re-examine all the moments, and determine the way to treat them. There are four cases:

  • 1.

    all the first order moments are treated as independent variables;

  • 2.

    if a moment consists of only stochastic species, and its order is no larger than the prescribed highest allowed order, it will be treated as an independent variable, and the corresponding moment equation will be included and solved. For the sake of numerical stability, its value should be no larger than its deterministic counterpart. For example, if the ODE solver yields a value of  ⟨AB⟩  >  ⟨A⟩  ⟨B⟩ , then the latter value will be assigned to  ⟨AB⟩ ;

  • 3.

    if a moment consists of only stochastic species, and its order is larger than the prescribed highest allowed order, its value will be set to zero, and of course, its moment equation will not be solved. For example, if  ⟨A⟩     <    1 and  ⟨B⟩     <    1, then, with a highest allowed order set to two, moments such as  ⟨A(A − 1)(A − 2)⟩  and  ⟨A(A − 1)B⟩  will be set to zero. This follows from the discussion in Sect. 2.2.4;

  • 4.

    if a moment contains at least one deterministic species, it will not be treated as an independent variable, and its moment equation will not be solved. It can be evaluated in the following way: assuming that the moment under consideration has a form  ⟨AB(B − 1)⟩ , and that A is deterministic (i.e.  ⟨A   ⟩  >    1), then the value of  ⟨AB(B − 1)⟩  is set to be  ⟨A⟩  ⟨B(B − 1)⟩ . If B is also deterministic, then it will be evaluated as  ⟨A⟩  ⟨B⟩ 2. This follows from the discussion in Sect. 2.2.3.

From these procedures, we see that the number of equations, as well as the form of these equations will change when a transition between stochastic and deterministic state of certain species occurs. Each time the ODE system is updated, the ODE solver must therefore be re-initialized.

It seems possible to replace the sharp transition between the stochastic and deterministic state of a species (based on whether its average population is smaller than 1) with a smooth transition, e.g., using a weight function similar to that in Garrod (2008). However, it is not mathematically clear which weight function we should choose, and an arbitrary one might cause some artificial effects, so we prefer not to use this formulation.

3. Benchmark with the Monte Carlo approach

We compare the results of our HME approach with those from the exact stochastic simulation (Gillespie 2007; Charnley 1998, 2001; Vasyunin et al. 2009). The RE results are also compared for reference. As in previous studies (Charnley 2001; Vasyunin et al. 2009), we consider a closed chemical system in a volume containing exactly one grain particle. The number of each species in this volume is called a “population”, which can be translated into an abundance relative to H nuclei by multiplying it by the dust-to-gas ratio, which is 2.8    ×    10-12(0.1  μm/r)3, where r is the grain radius, assuming an average molecular weight of 1.4, a dust-to-gas mass ratio 0.01, and a density of grain material of 2 g cm-3.

In the MC approach, the number of each species in this volume at any time is an integer. Owing to the large amount of time steps (> 109), it is impractical to store all intermediate steps, so we average the population of each species in time, weighted by the time intervals (remember that the lengths of time intervals between reactions are also random in MC). Because of this weighted average (rather than merely saving the state vector at certain instants), the MC approach can resolve average populations much smaller than one, although the fluctuations that are intrinsic to the MC approach can be larger than the average populations when the latter is small.

We first demonstrate how our method works for a large gas-grain network. We then show that for a small surface network, third order moments can also be included to improve the accuracy.

3.1. Test of the HME approach truncated at the second order on a large gas-grain network

Table 1

Percentage of agreement between the results from MC and those from HME/RE.

Table 2

Same as Table 1 except a smaller grain radius of 0.02 μm is taken.

We use the “dipole-enhanced” version of the RATE06 gas phase reaction network5 (Woodall et al. 2007), coupled with a surface network of Keane (1997, see Appendix B). The surface network contains 44 reactions between 43 species, which is basically a reduced and slightly revised version of the network of Tielens & Hagen (1982), containing the formation routes of the most common grain mantle species, such as H2O, CH3OH, CH4, NH3, etc. This surface network is not really large in comparison with some of the previous works, such as that used by Garrod et al. (2009). However, it is already essential for the most important species. The energy barriers for thermal desorption and diffusion are taken from Stantcheva & Herbst (2004). Diffusion of H atoms on the surface through quantum tunneling is included. Desorption by cosmic rays is taken into account following the approach of Hasegawa & Herbst (1993). The rate coefficients of the gas phase reactions are calculated according to Woodall et al. (2007), while the rate coefficients of the surface reactions are calculated following Hasegawa et al. (1992). The initial condition is the same as in Stantcheva & Herbst (2004).

We assumed a dust-to-gas mass ratio of 0.01. The grain mass density is taken to be 2 g cm-3, with a site density 5 × 1013 cm-2. Two grain sizes have been used: 0.1 μm and 0.02 μm. A cosmic ray ionization rate of 1.3 × 10-17 s-1 is adopted. Four different temperatures (10, 20, 30, 50 K) and three different densities (103, 104, 105 cm-3) have been used. In total, the comparison has been made for 24 different sets of physical parameters. These conditions are commonly seen in translucent clouds and cold dark clouds.

As in Garrod et al. (2009), we make a global comparison between the results of MC, HME, and RE. For each set of physical parameters, the comparisons are made at a time of 103, 104, and 105 years. We calculate the percentage of species for which the agreement between MC and HME/RE is within a factor of 2 or 10. Only species with a population (either from MC or from HME/RE) larger than 10 are included for comparison. This is because for species with smaller populations, the intrinsic fluctuation in the MC results can be significant. For several different sets of physical parameters, we repeated the MC several times to get a feeling for how large the fluctuation magnitude would be, although this is impractical for all the cases.

thumbnail Fig. 2

Typical time evolution of the average populations of certain species from MC (solid lines), HME to the 2nd order (dotted lines), and RE (dashed lines). Note that the Monte Carlo has been repeated twice. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 2.8 × 10-12. Physical parameters used: T = 20 K, n = 2 × 105 cm-3, grain radius = 0.1 μm.

thumbnail Fig. 3

Cases in which the agreement between the results of MC and those of HME are not so good, especially at t = 103 year. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 3.5 × 10-10. Physical parameters used: T = 20 K, n = 2 × 105 cm-3, grain radius = 0.02 μm.

The comparison results are shown in Table 1 (grain radius = 0.1 μm) and Table 2 (grain radius = 0.02 μm). The HME approach always has a better global agreement (or the same for several cases) than the RE approach in the cases we tested. The typical time evolution of certain species is shown in Fig. 2. In each panel of the figure, the species with a name preceded with a “g” means it is a surface species.

The poorest agreement of HME (Fig. 3) is at t = 103 year for T = 20 K, nH = 105 cm-3, and grain radius = 0.02 μm. This is mostly because at the time of comparison the populations of certain species were changing very rapidly, so a slight mismatch in time leads to a large discrepancy. This mismatch is probably caused by the truncation of higher-order terms in the HME (see Sect. 3.2). For gN2 in Fig. 3, its population seems to be systematically smaller in HME than in MC during the early period, although the HME result matches the one from MC at a later stage (after 3 × 103 years).

The RE is as effective as the HME in several cases, when the temperature is either relatively low (~10 K) or high (~50 K) (see also Vasyunin et al. 2009), and generally works better for a grain radius of 0.1 μm than of 0.02 μm. When the temperature is very low, many surface reactions with barriers cannot happen (at least in the considered timescales). On the other hand, when the temperature is high, the surface species evaporate very quickly and the surface reactions are also unable to occur. In these two extreme cases, the surface processes are inactive, and the RE works fine.

thumbnail Fig. 4

Comparison of the results from MC, HME to the 2nd order, HME to the 3rd order, and RE. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 3.5 × 10-10. Physical parameters used when running these models include T = 10 K, nH = 2 × 105 cm-3, grain radius = 0.02 μm.

The RE becomes problematic in the intermediate cases, when the temperature is high enough for many surface reactions to occur, but not too high to evaporate all the surface species; in these cases the HME represents a major improvement over the RE. For a smaller grain radius, the population of each species in a volume containing one grain will be smaller, thus the stochastic effect will play a more important role, and the RE will tend to fail.

We note that, in the HME approach, there is no elemental leakage except those caused by the finite precision of the computer. In all the models that we have run, all the elements (including electric charge) are conserved with a relative error smaller than 5    ×    10-14. The reason why elemental conservation is always guaranteed is that either the rate equations or the moment equations for the first order moments conserve the elements.

When comparing the results from the HME approach with those from MC simulation, it is important to see how the intrinsic fluctuation in MC behaves. If we assume the probability distribution of the population of a species, say A, is Poissonian, then the variance of A is σ2(A) =  ⟨A⟩ . Hence if  ⟨A⟩  is small, the relative fluctuation of the MC result can be quite large. This fluctuation might be smoothed out by means of a weighted average in time, but this procedure does not always work. This is why we choose to only compare species with a population higher than 10, corresponding to an abundance relative to H nuclei of 2.8 × 10-11 (for grain radius = 0.1 μm) or 3.5 × 10-9 (for grain radius = 0.02 μm). For a real reaction network, it is usually difficult to predict the intrinsic fluctuation in a MC simulation, unless it is repeated many times. These fluctuations will not have any observational effects, because along a line of sight there are always a large number of a certain species (as far as it is detectable) and the fluctuations are averaged out.

We note that the gas phase processes are not treated identically in our HME approach and MC simulation. In the MC approach, the gas phase processes are always treated as being stochastic (see, e.g., Charnley 1998; Vasyunin et al. 2009), in the same way as the surface processes. However, in our HME approach, the gas phase species are treated in a deterministic way, i.e., REs are always applied to them. This means that even if two reacting gas phase species A and B both have average populations much smaller than one, we still assume that  ⟨AB⟩  =  ⟨A⟩  ⟨B⟩ . This is physically quite reasonable, because the presence of large amounts of reacting partners in the gas phase (if not limited to a volume containing only one dust grain; see, e.g., Charnley 1998) ensures that the RE is applicable. However, although it might sound a bit pedantic, mathematically this is not equivalent to the MC approach, and some discrepancies caused by this are expected. For a large network, it is impractical to treat the gas phase processes in the same way as the surface processes in the HME approach, because in that case the number of independent variables in the ODE system will be quite large (at least no less than the number of two-body reactions), and the performance of the ODE solver will be degraded.

3.2. Test of the HME approach truncated at the third order on a small surface network

To test the improvement in accuracy when the cutoff is made at a higher order, we compare the results of the HME approach with a cutoff at the second order to those obtained from the same approach with a cutoff at the third order. We use a small surface reaction network of Stantcheva & Herbst (2004), containing 17 surface reactions between 21 species, producing H2O, CH3OH, CH4, NH3, and CO2. No gas phase reactions are included, except adsorption and desorption processes. The initial gas phase abundances of the relevant species are obtained from the steady state solution of the RATE06 network under the corresponding physical conditions.

As before, we run the HME, RE, and the MC code for different sets of physical parameters. Although by transferring from the RE to the second order HME a major improvement in accuracy can be obtained, the inclusion of the third order moments to the HME usually only improves slightly over the second order case. In Fig. 4, we show an example (T = 10 K, nH = 2 × 105 cm-3, grain radius =  0.02 μm), in which the distinctions between the results from the second and third order HME are relatively large.

For several species, we note that the third order HME is still unable to match the MC results perfectly, and for gHCO (Fig. 4) the third order HME even produces an artificial spike in the time evolution curve. The results from the third order HME are otherwise of greater accuracy than the second order one, the abundances of gH2CO and gCH3OH in particular being in almost perfect agreement with those from the MC approach. In the case of gHCO, the timescale mismatch between HME and MC is alleviated by including the third order moments.

It might be useful to see the difference between the second order HME and the third order HME in a computational sense. For the current reaction network with physical parameters described above, the number of variables (same as the number of equations, which changes with time) is 145 initially in the second order case, and this number becomes 705 for the third order case. To reach a time span of  ~106 year, the second order HME takes about 3 s, while the third order one takes about 220 s on a standard desktop computer (a CPU @ 3.00 GHz with double cores, 4 GBytes memory). The number of variables depends on the network structure, and it is not straightforward to derive a formula to calculate it. Qualitatively, this number (as a function of the number of reactions or the number of species) seems to increase with the cutoff order less quickly than exponential growth. However, such a “mild” increase affects the behavior of the ODE solver quite significantly. This is partly because the solver contains operations (such as matrix inversion) that become slower as the number of variables become larger. An increase in the number of variables might also increase the stiffness of the problem, let alone the memory limitation of the computer. For the larger reaction network described in the previous section, the third order HME would involve about 5000 variables and has not been tested successfully.

4. Discussion

In our HME approach, we have used a general and automatic way to derive the MEs. For large gas-grain networks, the MEs are cut off at the second order. For small networks, a cutoff at the third order is possible and higher accuracy can be obtained. We incorporate a switching scheme between the ME and RE when the average population of a species reaches 1.

The results from HME are more accurate than those from the RE in the cases we have tested, when benchmarked against the exact MC results. The abundances of almost all the abundant species (≳ 2.8    ×    10-11 for a grain radius of 0.1 μm and  ≳ 3.5    ×    10-9 for a grain radius of 0.02 μm) from HME are accurate to within a factor of two, especially at later stages of the chemical evolution, while in some cases nearly 40% of the results from RE are incorrect by a factor of at least ten.

In terms of computation time, our approach usually takes several tens of minutes to reach a evolution time span of 106 years, so it is slower than the RE, but faster than the MC approach (which usually takes from several hours to days). Our approach may also be slower than the MRE approach of Garrod (2008), because more variables (namely the moments with orders higher than one) are present in our method, and the ODE system in HME is usually stiffer. For example, for a moderate temperature many surface reactions can be much faster than any gas phase reactions, and yield a very large coefficient in some of the MEs. However, this is not the case in the stochastic regime of the MRE approach, because when a competition scheme is used in MRE, such a large coefficient does not appear. In this sense, we also advocate the MRE approach of Garrod (2008).

Mathematically, our approach is partially equivalent to the master equation approach of Stantcheva & Herbst (2004) in two respects. 1) They separated stochastic and deterministic species, which is similar to our adopting RE for the abundant species. 2) They set a cutoff for the possible states of the stochastic species. This is in essence equivalent to letting moments containing these species with order higher than a certain number equal to zero.

Our approach can also be viewed as a combination of some of the ideas of Garrod (2008) and Barzel & Biham (2007a). The basic modification and competition scheme in Garrod (2008) can be derived from the MEs, with a semi-steady state assumption for the second order moments. Barzel & Biham (2007a) used the MEs, but they did not include a switch scheme, and their way of deriving the MEs is different from the one in the present paper.

There are still many possibilities for improvement. Although in principle moments with any order can be included, the number of equations grows quite quickly with the cutoff order, which makes the system of equations intractable with a normal desktop computer. It is unclear whether it is possible to include the moments selectively. It is unclear whether there are better, and more mathematically well founded strategies than the switch at an average population of 1. The present approach is usually stable numerically. However, this is not always guaranteed, especially if higher order moments are to be included. The behavior of the numerical solution also depends on other factors, such as the ODE solver being used and the tunable parameters for it, while the MC approach does not have such issues. In this sense, the MC approach is the most robust.

Even in the accurate MC approach described above, the detailed morphology of the grain surface and the detailed reaction mechanism is not taken into account. One step in this direction would be to take into account the layered structure of the grain mantle. This was done by Charnley (2001; see also Charnley & Rodgers 2009) by means of stochastic simulation. It could also be included in the HME approach, as far as the underlying physical mechanism could be described by a master equation.

However, a microscopic MC approach has also been used to study the grain chemistry (see, e.g., Chang et al. 2005; Cuppen & Herbst 2007). In this approach, the morphology of the grain mantle and the interaction between species are modeled in detail. As far as we know, this approach is only practical when the network is small. It remains unclear whether it is possible to incorporate these details into the current HME approach.

In some cases, errors caused by uncertainties in the reaction mechanism and rate parameters might be larger than those introduced by the modeling method (Vasyunin et al. 2008). Hence, further experimental study and a more sophisticated way of interpreting those results would be indispensable.


Here by “population” we mean the number of a species in a volume of interest, and by “average” we mean an ensemble average (i.e. average over many different realizations of the same system setup). Hence “population” can only take non-negative integer values, while “average population” is a non-negative real number.


For a discussion about the relations and differences between “stochastic simulation” and “Monte Carlo”, see Kalos & Whitlock (2008) and Ripley (2008).


Assuming Poisson statistics, we have


Downloaded from


Instead of using zis as symbols for the independent variables of the generating function f, we use xis and yis instead. This will not cause any confusion.


We thank Rob Garrod, Julia Roberts, Tom Millar, Guillaume Pineau des Forêts, and Malcolm Walmsley for answering questions of the first author during the early stage of this work. We also thank A. G. G. M. Tielens, Anton Vasyunin, and Eric Herbst for discussions related to the present work, and Antoine Gusdorf for useful comments. This work is financially supported by the Deutsche Forschungsgemeinschaft Emmy Noether program under grant PA1692/1-1. The first author acknowledges travel funding from the IMPRS for astronomy and astrophysics at the Universities of Bonn and Cologne.


  1. Awad, Z., Chigai, T., Kimura, Y., Shalabiea, O. M., & Yamamoto, T. 2005, ApJ, 626, 262 [NASA ADS] [CrossRef] [Google Scholar]
  2. Barzel, B., & Biham, O. 2007a, ApJ, 658, L37 [NASA ADS] [CrossRef] [Google Scholar]
  3. Barzel, B., & Biham, O. 2007b, J. Chem. Phys., 127, 144703 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  4. Barzel, B., & Biham, O. 2011, Phys. Rev. Lett., 106, 150602 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  5. Biham, O., & Lipshtat, A. 2002, Phys. Rev. E, 66, 056103 [NASA ADS] [CrossRef] [Google Scholar]
  6. Biham, O., Furman, I., Pirronello, V., & Vidali, G. 2001, ApJ, 553, 595 [NASA ADS] [CrossRef] [Google Scholar]
  7. Caselli, P., Hasegawa, T. I., & Herbst, E. 1998, ApJ, 495, 309 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chang, Q., Cuppen, H., & Herbst, E. 2005, A&A, 434, 599 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Charnley, S. B. 1998, ApJ, 509, L121 [NASA ADS] [CrossRef] [Google Scholar]
  10. Charnley, S. B. 2001, ApJ, 562, L99 [NASA ADS] [CrossRef] [Google Scholar]
  11. Charnley, S. B. 2005, Adv. Space Res., 36, 132 [NASA ADS] [CrossRef] [Google Scholar]
  12. Charnley, S. B., & Rodgers, S. B. 2009, in ASP Conf. Ser. 420, ed. K. J. Meech, J. V. Keane, M. J. Mumma, J. L. Siefert, & D. J. Werthimer, 29 [Google Scholar]
  13. Charnley, S. B., Tielens, A. G. G. M., & Rodgers, S. D. 1997, ApJ, 482, L203 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cuppen, H., & Herbst, E. 2007, ApJ, 668, 294 [Google Scholar]
  15. Garrod, R. 2008, A&A, 491, 239 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Garrod, R., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Garrod, R., Vasyunin, A., Semenov, D., Wiebe, D., & Henning, T. 2009, ApJ, 700, L43 [Google Scholar]
  18. Gillespie, D. 1976, J. Comp. Phys., 22, 403 [Google Scholar]
  19. Gillespie, D. 1977, J. Phys. Chem., 81, 25 [Google Scholar]
  20. Gillespie, D. T. 2000, J. Chem. Phys., 113, 297 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gillespie, D. 2007, Ann. Rev. Phys. Chem., 58, 35 [Google Scholar]
  22. Green, N., Toniazzo, T., Pilling, M., et al. 2001, A&A, 375, 1111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hasegawa, T., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kalos, M. H., & Whitlock, P. A. 2008, Monte Carlo Methods: Second Revised and Enlarged Edition (Wiley-VCH Verlag) [Google Scholar]
  26. Keane, J. 1997, Ph.D. Thesis, Rijks Univ., Leiden [Google Scholar]
  27. Lipshtat, A., & Biham, O. 2003, A&A, 400, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Mattis, D., & Glasser, M. 1998, Rev. Mod. Phys., 70, 3 [Google Scholar]
  29. McQuarrie, D. A. 1967, J. Appl. Probab., 4, 413 [CrossRef] [Google Scholar]
  30. Menten, K. M., Walmsley, C. M., Henkel, C., & Wilson, T. L. 1988, A&A, 198, 253 [NASA ADS] [Google Scholar]
  31. Ohkubo, J. 2008, J. Chem. Phys., 129, 044108 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  32. Rae, J. G. L., Green, N. J. B., Hartquist, T. W., Pilling, M. J., & Toniazzo, T. 2003, A&A, 405, 387 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Ripley, B. D. 2008, Stochastic Simulation (John Wiley & Sons, Inc.) [Google Scholar]
  34. Stantcheva, T., & Herbst, E. 2004, A&A, 423, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Stantcheva, T., Caselli, P., & Herbst, E. 2001, A&A, 375, 673 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Stantcheva, T., Shematovich, V., & Herbst, E. 2002, A&A, 391, 1069 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Tielens, A., & Hagen, W. 1982, A&A, 114, 245 [NASA ADS] [Google Scholar]
  38. Tielens, A. G. G. M., & Charnley, S. B. 1997, OLEB, 27, 23 [Google Scholar]
  39. van Kampen, N. 2007, Stochastic processes in physics and chemistry, 3rd edn. (North Holland) [Google Scholar]
  40. Vasyunin, A., Semenov, D., Henning, T., et al. 2008, ApJ, 672, 629 [NASA ADS] [CrossRef] [Google Scholar]
  41. Vasyunin, A., Semenov, D., Wiebe, D., & Henning, T. 2009, ApJ, 691, 1459 [NASA ADS] [CrossRef] [Google Scholar]
  42. Woodall, J., Agúndez, M., Markwick-Kemper, A., & Millar, T. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: A method to generate the moment equations based on the generating function

We describe our means of getting the MEs. Our method is automatic, and can be easily coded into a computer program. It is applicable to moments of any order and all the common astrochemical reactions. It makes use of the probability generating function. While preparing the present paper, we noted that Barzel & Biham (2011) also proposed a binomial formulation of ME, which in essence is partly equivalent to our approach presented here, although our way of deriving the MEs is quite different from theirs.

For a probability distribution P(x,t), the corresponding generating function is defined as (van Kampen 2007) (A.1)here all the zis should be thought of as merely symbols without any physical meaning, and they have a one-to-one correspondence with the xis.

It is obvious that f(z = 1,t) ≡ 1, which is just the normalization condition for probability. It is also easy to see that the average population of the ith species is (A.2)The right hand side of the above equation means taking the partial derivative first, then assigning a value one to all the zis.

For the second order moment between two distinct species i and j, we have (A.3)If i equals j in the above equation, then what we actually get is (A.4)In general, we have (A.5)If several of the subscripts are the same in the left hand side of the above equation, say, i = j = k, then the second should be understood as (xj − 1), while the third should be understood as (xk − 2), and so on.

From the master equation in Eq. (1), it seems possible to get an equation for the evolution of f(z,t) in the general case. However, if the propensity functions ai(x) (see Eq. (1)) are allowed to take any functional form, then this is not straightforward. Fortunately, in practice ai(x) usually has a very simple form. On the other hand, we note that in the right hand side of the master equation in Eq. (1) the contributions from all the reactions are added linearly. Hence the contribution of each reaction to the evolution of generating function can be considered independently of each other.

We assume there is only one reaction in the network, which has a form (A.6)where the xis and yis represent the reactants and products, which do not have to be different from each other. We also use these symbols to represent the populations of the corresponding species. If, given a population of x1, x2, ..., xn, the probability that the above reaction will happen in a unit time is kx1x2...xn (namely, the propensity function a(x) = kx1x2...xn; the product should be understood as explained in the sentence following Eq. (A.5)), then the generating function6 will evolve according to (A.7)It is not difficult to derive the above equation from the master equation and the definition of generating equation and our assumption about the propensity function.

We note that Eq. (A.7) has a very simple pattern that is easy to remember: (a) the constant coefficient is the rate coefficient; (b) in the parenthesis, the symbols of all the products are multiplied together, with a coefficient  + 1, while the symbols of all the reactants are multiplied together, with a coefficient  −1; (c) in the differential part, all the reactants are present as they are in the left hand side of Eq. (A.6), while none of the products appear.

We now attempt to derive the ME. We obtain the evolution equation of each moment (as defined in Eq. (A.5)) by simply differentiating both sides of Eq. (A.7) with respect to the relevant components in the moment, then setting all the symbols to a value of one.

For example, for the reaction , the evolution equation of the generating function f is For  ⟨AB⟩, we differentiate both sides of the above equation by A and B. We obtain Next we assign a value of one to all the symbols (A–D) appearing in the resulting expressions, and “translate” the remaining terms into moments (recalling the remark about Eq. (A.5)), obtaining Although the above derivation involves differentiations, these operations can be easily translated into some combinatorial rules and written as a computer program. A recursive procedure is needed to generate all the potentially needed moments up to a given order.

Appendix B: A surface reaction network we used to test our code

Table B.1

The surface network used in Sect. 3.1 of this paper.

All Tables

Table 1

Percentage of agreement between the results from MC and those from HME/RE.

Table 2

Same as Table 1 except a smaller grain radius of 0.02 μm is taken.

Table B.1

The surface network used in Sect. 3.1 of this paper.

All Figures

thumbnail Fig. 1

A flow chart of the main components of our HME code. Steps 2, 5, 6, and 8 are described in detail in the text.

In the text
thumbnail Fig. 2

Typical time evolution of the average populations of certain species from MC (solid lines), HME to the 2nd order (dotted lines), and RE (dashed lines). Note that the Monte Carlo has been repeated twice. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 2.8 × 10-12. Physical parameters used: T = 20 K, n = 2 × 105 cm-3, grain radius = 0.1 μm.

In the text
thumbnail Fig. 3

Cases in which the agreement between the results of MC and those of HME are not so good, especially at t = 103 year. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 3.5 × 10-10. Physical parameters used: T = 20 K, n = 2 × 105 cm-3, grain radius = 0.02 μm.

In the text
thumbnail Fig. 4

Comparison of the results from MC, HME to the 2nd order, HME to the 3rd order, and RE. The y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 3.5 × 10-10. Physical parameters used when running these models include T = 10 K, nH = 2 × 105 cm-3, grain radius = 0.02 μm.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.