A hybrid moment equation approach to gasgrain chemical modeling
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
email: fjdu@mpifr.de
Received: 5 December 2010
Accepted: 1 April 2011
Context. In addition to gas phase reactions, the chemical processes on the surfaces of interstellar dust grains are important for the energy and material budget of the interstellar medium. The stochasticity of these processes requires special care in modeling. Previously methods based on the modified rate equation, the master equation, the moment equation, and Monte Carlo simulations have been used.
Aims. We attempt to develop a systematic and efficient way to model the gasgrain chemistry with a large reaction network as accurately as possible.
Methods. We present a hybrid moment equation approach, which is a general and automatic method where the generating function is used to generate the moment equations. For large reaction networks, the moment equation is cut off at the second order, and a switch scheme is used when the average population of certain species reaches 1. For small networks, the third order moments can also be utilized to achieve a higher accuracy.
Results. For physical conditions in which the surface reactions are important, our method provides a major improvement over the rate equation approach, when benchmarked against the rigorous Monte Carlo results. For either very low or very high temperatures, or large grain radii, results from the rate equation are similar to those from our new approach. Our method is faster than the Monte Carlo approach, but slower than the rate equation approach.
Conclusions. The hybrid moment equation approach with a cutoff and switch scheme is a very powerful way to solve gasgrain chemistry. It is applicable to large gasgrain networks, and is demonstrated to have a degree of accuracy high enough to be used for astrochemistry studies. Further work should be done to investigate how to cut off the hybrid moment equation selectively to make the computation faster, more accurate, and more stable, how to make the switch to rate equation more mathematically sound, and how to make the errors controllable. The layered structure of the grain mantle could also be incorporated into this approach, although a full implementation of the grain microphysics appears to be difficult.
Key words: astrochemistry / ISM: abundances / ISM: clouds / ISM: molecules / radio lines: ISM / stars: formation
© 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 H_{2}, 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 gasgrain 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 populations^{1} of some reactants on a single grain is small. This failure of the rate equation is related to the treatment of twobody 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 grainsurface 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 gasgrain 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 x_{i} being the number of the ith species in the system. The evolution equation of P(x,t) is the socalled master equation, whose form is determined by the reaction network.
Many sophisticated methods have been proposed (mainly outside the astrochemical 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 H_{2} 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 literature^{2}. 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 socalled 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 gasgrain chemistry, named the hybrid moment equation (HME) approach. The goal is to find a systematic, automatic, and fast way to modeling gasgrain 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 semisteadystate 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 gasgrain 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 x_{j} 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 a_{i}(x) is called the propensity function, a_{i}(x)Δt 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 ⟨x_{j}⟩ , which is simply the average number of species j, ⟨x_{j}⟩ ≡ ∑ _{x}P(x,t)x_{j}, 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 singlebody reactions, a_{i}(x) is a linear function of x. In this case the ME is closed and can easily be solved. However, when twobody reactions are present, this is no longer true, as ⟨a_{i}(x)⟩ might be of a form ⟨x_{k}(x_{k} − 1)⟩ or ⟨x_{k}x_{l}⟩ , 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 finitedimensional 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 k_{ad}a, k_{evap}A, k_{AB}AB, and k_{AA}A(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. (3–6) ), 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. (7–11) 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. (7–11)) 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 k_{AB} ⟨AB⟩ in Eq. (10) and the term k_{AA} ⟨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 ^{3} ⟨A⟩ ⟨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 k_{AB} ⟨AB⟩ might be much larger than the retained terms such as ⟨A⟩ ^{2} ⟨B⟩ or ⟨A⟩ ⟨B⟩ ^{2}, and the omitted term k_{AA} ⟨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. (7–11) we do not write terms such as ⟨A(A − 1)⟩ in the split form ⟨A^{2}⟩ − ⟨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
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. 

Open with DEXTER 
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 package^{4}.
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 reexamine 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 reinitialized.
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 dusttogas 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 dusttogas 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 (> 10^{9}), 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 gasgrain 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 gasgrain network
Percentage of agreement between the results from MC and those from HME/RE.
We use the “dipoleenhanced” version of the RATE06 gas phase reaction network^{5} (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 H_{2}O, CH_{3}OH, CH_{4}, NH_{3}, 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 dusttogas mass ratio of 0.01. The grain mass density is taken to be 2 g cm^{3}, with a site density 5 × 10^{13} 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 (10^{3}, 10^{4}, 10^{5} 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 10^{3}, 10^{4}, and 10^{5} 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.
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 yaxis 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 × 10^{5} cm^{3}, grain radius = 0.1 μm. 

Open with DEXTER 
Fig. 3
Cases in which the agreement between the results of MC and those of HME are not so good, especially at t = 10^{3} year. The yaxis 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 × 10^{5} cm^{3}, grain radius = 0.02 μm. 

Open with DEXTER 
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 = 10^{3} year for T = 20 K, n_{H} = 10^{5} 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 higherorder terms in the HME (see Sect. 3.2). For gN_{2} 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 × 10^{3} 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.
Fig. 4
Comparison of the results from MC, HME to the 2nd order, HME to the 3rd order, and RE. The yaxis 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, n_{H} = 2 × 10^{5} cm^{3}, grain radius = 0.02 μm. 

Open with DEXTER 
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 twobody 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 H_{2}O, CH_{3}OH, CH_{4}, NH_{3}, and CO_{2}. 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, n_{H} = 2 × 10^{5} 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 gH_{2}CO and gCH_{3}OH 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 ~10^{6} 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 gasgrain 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 10^{6} 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 semisteady 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 nonnegative integer values, while “average population” is a nonnegative real number.
For a discussion about the relations and differences between “stochastic simulation” and “Monte Carlo”, see Kalos & Whitlock (2008) and Ripley (2008).
Downloaded from www.netlib.org
Acknowledgments
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/11. The first author acknowledges travel funding from the IMPRS for astronomy and astrophysics at the Universities of Bonn and Cologne.
References
 Awad, Z., Chigai, T., Kimura, Y., Shalabiea, O. M., & Yamamoto, T. 2005, ApJ, 626, 262 [NASA ADS] [CrossRef] (In the text)
 Barzel, B., & Biham, O. 2007a, ApJ, 658, L37 [NASA ADS] [CrossRef] (In the text)
 Barzel, B., & Biham, O. 2007b, J. Chem. Phys., 127, 144703 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Barzel, B., & Biham, O. 2011, Phys. Rev. Lett., 106, 150602 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Biham, O., & Lipshtat, A. 2002, Phys. Rev. E, 66, 056103 [NASA ADS] [CrossRef] (In the text)
 Biham, O., Furman, I., Pirronello, V., & Vidali, G. 2001, ApJ, 553, 595 [NASA ADS] [CrossRef] (In the text)
 Caselli, P., Hasegawa, T. I., & Herbst, E. 1998, ApJ, 495, 309 [NASA ADS] [CrossRef] (In the text)
 Chang, Q., Cuppen, H., & Herbst, E. 2005, A&A, 434, 599 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Charnley, S. B. 1998, ApJ, 509, L121 [NASA ADS] [CrossRef] (In the text)
 Charnley, S. B. 2001, ApJ, 562, L99 [NASA ADS] [CrossRef] (In the text)
 Charnley, S. B. 2005, Adv. Space Res., 36, 132 [NASA ADS] [CrossRef] (In the text)
 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 (In the text)
 Charnley, S. B., Tielens, A. G. G. M., & Rodgers, S. D. 1997, ApJ, 482, L203 [NASA ADS] [CrossRef] (In the text)
 Cuppen, H., & Herbst, E. 2007, ApJ, 668, 294 [NASA ADS] [CrossRef] (In the text)
 Garrod, R. 2008, A&A, 491, 239 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Garrod, R., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Garrod, R., Vasyunin, A., Semenov, D., Wiebe, D., & Henning, T. 2009, ApJ, 700, L43 [NASA ADS] [CrossRef] (In the text)
 Gillespie, D. 1976, J. Comp. Phys., 22, 403 [NASA ADS] [CrossRef] (In the text)
 Gillespie, D. 1977, J. Phys. Chem., 81, 25 [CrossRef] (In the text)
 Gillespie, D. T. 2000, J. Chem. Phys., 113, 297 [NASA ADS] [CrossRef] (In the text)
 Gillespie, D. 2007, Ann. Rev. Phys. Chem., 58, 35 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Green, N., Toniazzo, T., Pilling, M., et al. 2001, A&A, 375, 1111 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] (In the text)
 Hasegawa, T., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] (In the text)
 Kalos, M. H., & Whitlock, P. A. 2008, Monte Carlo Methods: Second Revised and Enlarged Edition (WileyVCH Verlag) (In the text)
 Keane, J. 1997, Ph.D. Thesis, Rijks Univ., Leiden (In the text)
 Lipshtat, A., & Biham, O. 2003, A&A, 400, 585 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Mattis, D., & Glasser, M. 1998, Rev. Mod. Phys., 70, 3 [CrossRef] (In the text)
 McQuarrie, D. A. 1967, J. Appl. Probab., 4, 413 [CrossRef] (In the text)
 Menten, K. M., Walmsley, C. M., Henkel, C., & Wilson, T. L. 1988, A&A, 198, 253 [NASA ADS] (In the text)
 Ohkubo, J. 2008, J. Chem. Phys., 129, 044108 [NASA ADS] [CrossRef] [PubMed] (In the text)
 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] (In the text)
 Ripley, B. D. 2008, Stochastic Simulation (John Wiley & Sons, Inc.) (In the text)
 Stantcheva, T., & Herbst, E. 2004, A&A, 423, 241 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Stantcheva, T., Caselli, P., & Herbst, E. 2001, A&A, 375, 673 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Stantcheva, T., Shematovich, V., & Herbst, E. 2002, A&A, 391, 1069 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tielens, A., & Hagen, W. 1982, A&A, 114, 245 [NASA ADS] (In the text)
 Tielens, A. G. G. M., & Charnley, S. B. 1997, OLEB, 27, 23 [NASA ADS] [CrossRef] [PubMed] (In the text)
 van Kampen, N. 2007, Stochastic processes in physics and chemistry, 3rd edn. (North Holland) (In the text)
 Vasyunin, A., Semenov, D., Henning, T., et al. 2008, ApJ, 672, 629 [NASA ADS] [CrossRef] (In the text)
 Vasyunin, A., Semenov, D., Wiebe, D., & Henning, T. 2009, ApJ, 691, 1459 [NASA ADS] [CrossRef] (In the text)
 Woodall, J., Agúndez, M., MarkwickKemper, A., & Millar, T. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
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 z_{i}s should be thought of as merely symbols without any physical meaning, and they have a onetoone correspondence with the x_{i}s.
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 z_{i}s.
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 (x_{j} − 1), while the third should be understood as (x_{k} − 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 a_{i}(x) (see Eq. (1)) are allowed to take any functional form, then this is not straightforward. Fortunately, in practice a_{i}(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 x_{i}s and y_{i}s 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 x_{1}, x_{2}, ..., x_{n}, the probability that the above reaction will happen in a unit time is kx_{1}x_{2}...x_{n} (namely, the propensity function a(x) = kx_{1}x_{2}...x_{n}; the product should be understood as explained in the sentence following Eq. (A.5)), then the generating function^{6} 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
All Tables
All Figures
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. 

Open with DEXTER  
In the text 
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 yaxis 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 × 10^{5} cm^{3}, grain radius = 0.1 μm. 

Open with DEXTER  
In the text 
Fig. 3
Cases in which the agreement between the results of MC and those of HME are not so good, especially at t = 10^{3} year. The yaxis 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 × 10^{5} cm^{3}, grain radius = 0.02 μm. 

Open with DEXTER  
In the text 
Fig. 4
Comparison of the results from MC, HME to the 2nd order, HME to the 3rd order, and RE. The yaxis 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, n_{H} = 2 × 10^{5} cm^{3}, grain radius = 0.02 μm. 

Open with DEXTER  
In the text 