{'fraction_non_alphanumeric': 0.03774699665200552, 'fraction_numerical': 0.02126961202652137, 'mean_word_length': 4.37733851041299, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Simulations of Interstellar Gas-Grain Chemistry Using Moment Equations

Networks of reactions on dust grain surfaces play a crucial role in the chemistry of interstellar clouds, leading to the formation of molecular hydrogen in diffuse clouds as well as various organic molecules in dense molecular clouds. Due to the sub-micron size of the grains and the low flux, the population of reactive species per grain may be very small and strongly fluctuating. Under these conditions rate equations fail and the simulation of surface-reaction networks requires stochastic methods such as the master equation. However, the master equation becomes infeasible for complex networks because the number of equations proliferates exponentially. Here we introduce a method based on moment equations for the simulation of reaction networks on small grains. The number of equations is reduced to just one equation per reactive specie and one equation per reaction. Nevertheless, the method provides accurate results, which are in excellent agreement with the master equation. The method is demonstrated for the methanol network which has been recently shown to be of crucial importance.


Introduction
Chemical networks in interstellar clouds consist of gas-phase and grain-surface reactions (Hartquist & Williams 1995;Tielens 2005). Reactions that take place on dust grains include the formation of molecular hydrogen (Gould & Salpeter 1963;Hollenbach et al. 1971) as well as reaction networks producing ice mantles and various organic molecules. Unlike gas phase reactions in cold clouds that mainly produce unsaturated molecules, surface processes are dominated by hydrogen-addition reactions that result in saturated, hydrogen-rich molecules, such as H 2 CO, CH 3 OH, NH 3 and CH 4 . In particular, recent experiments show that methanol cannot be efficiently produced by gas phase reactions (Geppert et al. 2006).
On the other hand, there are indications that it can be efficiently produced on ice-coated grains (Watanabe 2005). Therefore, the ability to perform simulations of the production of methanol and other complex molecules on grains is of great importance (Garrod et al. 2006). Unlike gas-phase reactions, simulated using rate equation models (Pickles & Williams 1977;Hasegawa et al. 1992), grain-surface reactions require stochastic methods such as the master equation (Biham et al. 2001;Green et al. 2001), or Monte Carlo (MC) simulations (Charnley 2001). This is due to the fact that under interstellar conditions, of extremely low gas density and sub-micron grain sizes, surface reaction rates are dominated by fluctuations which cannot be accounted for by rate equations (Tielens & Hagen 1982;Charnley et al. 1997;Caselli et al. 1998;Shalabiea et al. 1998). A significant advantage of the master equation over MC simulations is that it consists of differential equations, which can be easily coupled to the rate equations of gas-phase chemistry. Furthermore, unlike MC simulations that require the accumulation of statistical information over long times, the master equation provides the probability distribution from which the reaction rates can be obtained directly. However, the number of equations increases exponentially with the number of reactive species, making the simulation of complex networks infeasible (Stantcheva et al. 2002;Stantcheva & Herbst 2003). The recently proposed multi-plane method dramatically reduces the number of equations, by breaking the network into a set of fully connected sub-networks (Lipshtat & Biham 2004), enabling the simulation of more complex networks. However, the construction of the multi-plane equations for large networks turns out to be difficult.
In this Letter we introduce a method based on moment equations which exhibits crucial advantages over the multi-plane method. The number of equations is further reduced to the smallest possible set of stochastic equations, including one equation for the population size of each reactive specie (represented by a first moment) and one equation for each reaction rate (represented by a second moment). Thus, for typical sparse networks the complexity of the stochastic simulation becomes comparable to that of the rate equations. Unlike the master equation (and the multi-plane method) there is no need to adjust the cutoffs -the same set of equations applies under all physical conditions. Unlike the multi-plane equations, the moment equations are linear and for steady state conditions can be easily solved using algebraic methods. Moreover, for any given network the moment equations can be easily constructed using a diagrammatic approach, which can be automated (Barzel & Biham 2007).

The Method
To demonstrate the method we consider a simple network, shown in Fig. 1, that involves three reactive species: H and O atoms and OH molecules Shalabiea et al. 1998;Stantcheva et al. 2002). For simplicity we denote the reactive species by X 1 = H, X 2 = O, X 3 = OH, and the non-reactive product species by X 4 = H 2 , Consider a spherical grain of diameter d, exposed to fluxes of H and O atoms and OH molecules. The cross-section of the grain is σ = πd 2 /4 and its surface area is πd 2 . The density of adsorption sites on the surface is denoted by s (sites cm −2 ). Thus, the number of adsorption sites on the grain is S = πd 2 s. The desorption rates of atomic and molecular species from the grain are given by where ν is the attempt rate (standardly taken to be 10 12 s −1 ), E 1 (i) is the activation energy for desorption of specie X i and T (K) is the grain temperature. The hopping rate of adsorbed atoms between adjacent sites on the surface is a is the activation energy for hopping of X i atoms (or molecules). Here we assume that diffusion occurs only by thermal hopping, in agreement with experimental results (Katz et al. 1999;Perets et al. 2005). For small grains it is convenient to define the scanning rate, A i = a i /S, which is approximately the inverse of the time it takes an X i atom to scan the surface of the entire grain.
The master equation provides the time derivatives of the probabilities P (N 1 , N 2 , N 3 ) that on a random grain there will be N i adsorbed atoms/molecules of the reactive specie X i . It takes the forṁ (1) The terms in the first sum describe the incoming flux, where F i (atoms s −1 ) is the flux per grain of the specie X i . The second sum describes the effect of desorption. The third sum describes the effect of diffusion mediated reactions between two atoms of the same specie and the last two terms account for reactions between different species. The rate of each reaction is proportional to the number of pairs of atoms/molecules of the two species involved, and to the sum of their scanning rates. The moments of P (N 1 , N 2 , N 3 ) are given by N a where a, b, c are integers. In particular, N i is the average population size of the specie X i on a grain. The production rate per grain, R(X k ) (molecules s −1 ), of X k molecules produced by the reaction In numerical simulations the master equation must be truncated in order to keep the number of equations finite. This can be done by setting upper cutoffs N max i , i = 1, . . . , J on the population sizes, where J is the number of reactive species. However, the number of coupled equations, , grows exponentially with the number of reactive species. This severely limits the applicability of the master equation to interstellar chemistry (Stantcheva & Herbst 2003). To reduce the number of equations one tries to use the lowest possible cutoffs under the given conditions. In any case, to enable all reaction processes to take place, the cutoffs must satisfy N max i ≥ 2 for species that form homonuclear diatomic molecules (H 2 ,O 2 , etc.) and N max i ≥ 1 for other species.
The average population sizes of the reactive species and the reaction rates are completely determined by all the first moments and selected second moments of the distribution P (N 1 , N 2 , N 3 ). Therefore, a closed set of equations for the time derivatives of these first and second moments could provide complete information on the population sizes and reaction rates. For the simple network considered here one needs equations for the time derivatives of the first moments N 1 , N 2 and N 3 and of the second moments N 2 1 , N 2 2 , N 1 N 2 and N 1 N 3 [nodes and edges, respectively, in the graph shown in Fig. 1]. Such equations are obtained by taking the time derivative of each moment and using Eq. (1) to express the time derivatives of the probabilities (Lipshtat & Biham 2003). Here we show two of the resulting moment equations: In these equations, the time derivative of each moment is expressed as a linear combination of several other moments. However, the right hand sides of these equations include third order moments for which we have no equations. In order to close the set of moment equations we must express the third order moments in terms of first and second order moments (Lipshtat & Biham 2003). This can be done by imposing the following constraint on the master equation: at any given time, at most two atoms or molecules can be adsorbed simul-taneously on the surface. Furthermore, these two atoms or molecules must be from species that react with each other. The resulting cutoffs allow only eight non-vanishing probabilities, namely, P (0, 0, 0), P (0, 0, 1), P (0, 1, 0), P (1, 0, 0), P (2, 0, 0), P (0, 2, 0), P (1, 1, 0) and P (1, 0, 1). The third moments in Eq.
(2) can now be expressed in terms of these non-vanishing probabilities, giving rise to the following rules: Using these rules, which are general and apply to any network of binary reactions, one can modify Eqs. (2), and obtain a closed set of the form where i = 1, 2, j = 2, 3 and δ i,j = 1 if i = j and 0 otherwise. This set includes one equation that accounts for the population size of each reactive species and one equation that accounts for the rate of each reaction. Although these equations were derived using strict cutoffs, that are expected to apply only in the limit of very small grains and low flux, they provide accurate results for a very broad range of conditions. The point is that once the set of moment equations is derived, the probabilities do not appear anymore, so the constraint is not explicitly enforced. In fact, the equations maintain their accuracy even when the populations sizes of the reactive species are well beyond the constraints imposed above.
In Fig. 1(a) we present the population sizes of H (+) and O (×) atoms on a grain, vs. grain diameter, obtained from the moment equations. In Fig. 1(b) we present the production rates of H 2 (squares), O 2 (triangles) and H 2 O (circles) vs. grain diameter, obtained from the moment equations. The results are in excellent agreement with the master equation (solid lines). In the limit of large grains they also coincide with the rate equations (dashed lines). Note that the moment equations apply even when there are as many as 10 hydrogen atoms on a grain. The parameters used in the simulations are s = 5 × 10 13 (sites cm −2 ), F 1 = 5.0 × 10 −10 S (atoms s −1 ), F 2 = 0.1F 1 and F 3 = 0. The activation energies for diffusion and desorption were taken as E 0 (1) = 44, E 1 (1) = 52, E 0 (2) = 47, E 1 (2) = 54, E 0 (3) = 47 and E 1 (3) = 54 meV. The grain temperature was T = 15K. The parameters used for hydrogen are the experimental results for low density amorphous ice (Perets et al. 2005).
For the other species there are no concrete experimental results, and the values reflect the tendency of heavier species to bind more strongly. The fluxes and grain temperatures are suitable for dense molecular clouds.

The Methanol Network
Consider the case in which a flux of CO molecules is added to the network. This gives rise to the network shown in Fig. 2, which includes the following sequence of hydrogen addition reactions (Stantcheva et al. 2002): H + CO → HCO, H + HCO → H 2 CO, H + H 2 CO → H 3 CO and H + H 3 CO → CH 3 OH. Two other reactions that involve oxygen atoms also take place: O + CO → CO 2 and O + HCO → CO 2 + H. This network was studied before using the multiplane method, which required about a thousand equations compared to about a million equations in the master equation with similar cutoffs (Lipshtat & Biham 2004). The moment equations include one equation for each node and one for each edge, namely, the network shown in Fig. 2 requires only 17 equations. We have performed extensive simulations of this network using the moment equations and found that they are in excellent agreement with the master equation. In Fig. 2(a) we present the moment-equation results for the population sizes of H, O and CO on a grain vs. grain diameter for the methanol network. In Fig. 2(b) we present the moment-equation results for the production rates per grain of some of the final products of the network. The results are in excellent agreement with the master equation (solid lines), and coincide with the rate equations (dashed lines) for large grains. The activation energies for diffusion and desorption of CO (X 7 ), HCO (X 8 ), H 2 CO (X 9 ) and H 3 CO (X 10 ) were taken as E 0 (7) = 50, E 1 (7) = 55, E 0 (8) = 52, E 1 (8) = 58, E 0 (9) = 53, E 1 (9) = 59, E 0 (10) = 55 and E 1 (10) = 62 meV. The flux of CO molecules was taken as F 7 = 0.2F 1 . Note that experiments on CO desorption from ice surfaces indicate that E 1 (7) should be slightly higher than the value used here (Collings et al. 2004). However, it turns out that using a higher value would compromise the feasibility of the master equation simulations, which are required here in order to establish the validity of the moment equation method.
When simulating highly complex networks one encounters the problem of obtaining the equations themselves. The ability to automate the construction of the equations becomes important. A crucial advantage of the moment equations is that they can be easily obtained using a diagrammatic approach. A detailed presentation of the diagrammatic method will be given in Barzel & Biham (2007).

Summary and Discussion
In summary, we have introduced a method, based on moment equations, for the simulation of chemical networks taking place on dust-grain surfaces in interstellar clouds. The method provides highly efficient simulations of complex reaction networks under the extreme conditions of low gas density and sub-micron grain sizes, in which the reaction rates are dominated by fluctuations and stochastic simulations are required. The number of equations is reduced to one equation for each reactive specie and one equation for each reaction, which is the lowest possible number for such networks. This method enables us to efficiently simulate networks of any required complexity without compromising the accuracy. It thus becomes possible to incorporate the complete network of surface reactions into gas-grain models of interstellar chemistry. To fully utilize the potential of this method, further laboratory experiments are needed, that will provide the activation energy barriers for diffusion, desorption and reaction processes not only for hydrogen but for all the molecules involved in these networks.