{'fraction_non_alphanumeric': 0.027225512528473803, 'fraction_numerical': 0.015908883826879272, 'mean_word_length': 4.4969448061704895, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gas-grain chemistry in cold interstellar cloud cores with a microscopic Monte Carlo approach to surface chemistry

AIM: We have recently developed a microscopic Monte Carlo approach to study surface chemistry on interstellar grains and the morphology of ice mantles. The method is designed to eliminate the problems inherent in the rate-equation formalism to surface chemistry. Here we report the first use of this method in a chemical model of cold interstellar cloud cores that includes both gas-phase and surface chemistry. The surface chemical network consists of a small number of diffusive reactions that can produce molecular oxygen, water, carbon dioxide, formaldehyde, methanol and assorted radicals. METHOD: The simulation is started by running a gas-phase model including accretion onto grains but no surface chemistry or evaporation. The starting surface consists of either flat or rough olivine. We introduce the surface chemistry of the three species H, O and CO in an iterative manner using our stochastic technique. Under the conditions of the simulation, only atomic hydrogen can evaporate to a significant extent. Although it has little effect on other gas-phase species, the evaporation of atomic hydrogen changes its gas-phase abundance, which in turn changes the flux of atomic hydrogen onto grains. The effect on the surface chemistry is treated until convergence occurs. We neglect all non-thermal desorptive processes. RESULTS: We determine the mantle abundances of assorted molecules as a function of time through 2x10^5 yr. Our method also allows determination of the abundance of each molecule in specific monolayers. The mantle results can be compared with observations of water, carbon dioxide, carbon monoxide, and methanol ices in the sources W33A and Elias 16. Other than a slight underproduction of mantle CO, our results are in very good agreement with observations.


Introduction
Not only are interstellar grain surfaces believed to be the sites where molecular hydrogen is formed in diffuse clouds, but in cold cores of dense clouds interstellar grain chemistry is an important mechanism for the formation of quite a few additional species, most of which are produced through gradual hydrogenation of gas-phase molecules (Tielens & Hagen 1982;Hasegawa et al. 1992;Herbst 1995;Ehrenfreund & Schutte 2000). Methanol is formed in such a manner, and its surface formation has received much attention over the years, both from theoreticians and modelers (Tielens & Whittet 1996;Stantcheva & Herbst 2004;Awad et al. 2005;) and experimentalists (Hiraoka et al. 2002;Watanabe & Kouchi 2002;Watanabe et al. 2003).
In order to simulate all of the chemistry that occurs in a cold region, it is necessary to include both the chemical kinetics in the gas and on the surfaces of dust grains. The basic method for simulating chemistry is to use coupled rate equations, one for each species in the simulation. Rate equations have been proven to be an accurate and efficient approach to simulate gas phase chemical kinetics (Herbst 1995;Roberts et al. 2004;Wakelam et al. 2006). In these Send offprint requests to: herbst@mps.ohio-state.edu equations, the average concentration of a species is progagated forward in time through a non-linear differential equation that contains both source and sink terms corresponding to chemical reactions. Their success comes from the fact that average concentrations are normally well defined in the gas phase as functions of time. In other words, the so-called mean-field approximation is valid. However, a straightforward and efficient generalization to simulate diffusive surface chemistry on interstellar grains in a similar manner (Pickles & Williams 1977;Hasegawa et al. 1992) can be in error (Tielens & Hagen 1982). In particular, when the number of surface species on an average grain is small, the rate equation method breaks down since its mean-field approach does not accurately describe the surface kinetics, which must be described by a method in which both the discrete number of species and fluctuations in this number are considered. Several alternative methods to rate equations have been proposed to overcome this difficulty. The simplest of these is the modified-rate-equation approach, which uses the rate equation method but instead of using the actual diffusion rate on the surface, scales this rate downward so that the rate of reaction does not on average exceed the rate of accretion of a reactant (Caselli et al. 1998). This semi-empirical approach is efficient computationally, but it is not reliable under all conditions (Stantcheva & Herbst 2004;. There are two more rigorous approaches, both of which are stochastic and macroscopic in nature. In other words, they allow a determination of the probability that a certain number of molecules of a given species exists on a grain at any time without distinguishing where on the grain the adsorbate might lie. In the Monte Carlo approach (Charnley 1998), random numbers are used to determine what processes occur in a time interval based on their relative rates (Caselli et al. 2002). It is difficult to develop a combined gas-grain model, in which the Monte Carlo method is used for the surface chemistry and the rate-equation approach for the gas-phase chemistry. The master equation approach (Biham et al. 2001;Green at al. 2001;Charnley 2001), on the other hand, is based on integrating a set of ordinary differential equations for the probabilities of surface populations, and can be coupled with rate equations for gas phase chemistry. Indeed, a coupling of the master equation approach to surface chemistry and a Monte Carlo approach to gas-phase chemistry has also been suggested (Charnley 1998). Although the master equation approach is believed to be the better candidate to solve the difficulties with the surface rate-equation approach, it has the drawback that there are many simultaneous equations to be solved, since the probabilities of surface populations for different species can be correlated. Nevertheless, with suitable approximations, two groups have successfully run gas-grain networks of cold interstellar cloud cores with the master equation approach to model surface chemistry (Stantcheva & Herbst 2004;Lipshtat & Biham 2004). In these simulations, only a small subset of the surface reactions used in full-size models (Hasegawa et al. 1992) were included, with methanol the most complex species synthesized.
Recently, we adopted a microscopic stochastic method for surface chemistry, known as continuous-time random walk (CTRW) Monte Carlo approach (Montroll & Weiss 1965), which tracks the specific trajectory and rate of diffusion of each adsorbate and the binding (desorption or evaporation) energy and barrier to diffusion at each local site (Chang et al. 2005;. In this method, grain surfaces are modeled as square lattices while hopping, evaporation and deposition are regarded as Poisson processes. We developed the CTRW approach for several reasons. First, for a surface with a continuous distribution of diffusion barriers and evaporation energies, neither the rate equation approach nor the macroscopic stochastic approaches can be easily implemented. Even for a surface with discrete values of the energy parameters, known as a rough surface, the earlier approaches are not readily useable. Moreover, it is facile to consider time-dependent barriers and desorption energies with the CTRW approach. Secondly, all previous methods neglected the "back diffusion" of hopping species on a granular surface. This neglect, however, turned out to be relatively unimportant for astronomical purposes (Lohmar & Krug 2006;Chang et al. 2006). Finally, the use of a microscopic Monte Carlo method allows us to follow the morphology of a mantle as it develops, so that only reactions in nearby monolayers can react with one another. For instance, a hydrogen atom on the surface is not likely to react with an oxygen atom deeply buried in the ice mantle. This constraint has been neglected in almost all previous gas-grain network calculations.
Although the CTRW Monte Carlo approach allows us in principle to study both the detailed microscopic chem-istry occurring in a grain mantle and the actual morphology of the mantle, use of the method presents some formidable computational challenges. Up to the present, we have used the method only to study the details of molecular hydrogen formation on surfaces relevant to diffuse clouds (Chang et al. 2005;Chang et al. 2006). In a work in preparation, we are expanding our use of the method to study the formation of mantles of water ice and related species when hydrogen and oxygen atoms accrete onto a grain surface (Cuppen & Herbst 2007). This extension still does not represent a true gas-grain model because the gas-phase chemistry is not followed. Indeed, there is no straightforward way to implement varying gasphase abundances at the same time the surface chemistry progresses, which makes it difficult to use the CTRW approach simultaneously with a gas phase kinetics simulation program.
The difficulty with time-varying gas-phase abundances is that these lead to time-varying fluxes of accreting species onto the grains. The CTRW method calculates deposition as a Poisson process with a time interval ∆t between events calculated from a random number X between 0 and 1 determined via the relation ∆t = −t ln X, wheret is the average time interval between deposition. If abundances change, this equation would be invalid sincet would no longer be constant and the deposition would not be a strictly defined Poisson process. This problem has been solved by Jansen (1995) for the simulation of a laboratory procedure known as temperature-programmed desorption (TPD), in which the temperature changes, resulting in a change of residence time for surface species. This technique, however, cannot be directly applied to a gas-grain chemical model because, unlike the TPD simulation, where the temperature is linearly dependent on time, no prior knowledge of the abundances of gas-phase species is available.
In the present paper, we present a first attempt to combine a treatment of surface chemistry by the CTRW approach with a treatment of gas-phase chemistry by the standard rate equation approach. The iterative method is based on rather weak coupling between the two chemistries, a condition that occurs primarily for cold cores, where thermal evaporation is negligible for most surface species. The present paper uses a variant of the so-called H, O and CO system as the surface reaction network (Charnley 2001;Stantcheva & Herbst 2004). In this system, the gasphase species atomic hydrogen, atomic oxygen, and carbon monoxide accrete onto grains and react to form species such as molecular oxygen, water, carbon dioxide, formaldehyde, methanol, and assorted radicals. The surface network, which is complex enough to mimic a full surface chemical reaction network, has been extensively studied by the modified rate equation approach (Caselli et al. 2002), the macroscopic Monte Carlo method (Caselli et al. 2002), and the master equation approach  with fixed gas-phase abundances of the accreting species. The master equation has also been used in a calculation, like the present, where gas-phase chemistry also occurs (Stantcheva & Herbst 2004;Lipshtat & Biham 2004). By carefully studying the influence of the desorption of species from the grain surfaces on the overall gas phase abundances, we are led to our weak coupling hypothesis that the gasphase kinetics is almost independent of the surface chemical kinetics. The weak coupling allows us to handle the chemistry in an iterative manner.
The paper is organized as follows. Section 2 reviews the H, O and CO system, while Section 3 discusses the initial conditions and basic assumptions used in our treatment. The details of our method of simulation and surface model are discussed in Section 4. We present and discuss our results in Section 5, including a comparison with observations of molecular ices in cold regions. Our conclusions are stated in Section 6.

The H, O, and CO system
In the stochastic simulations to be discussed later, we follow only the accretion of H, O and CO from the gas-phase and the surface reactions that result from this accretion. As will be discussed in detail, the accretion of all species is considered in a calculation prior to the Monte Carlo simulation. Table 1 shows the limited surface reaction network used in our stochastic simulations. With two exceptions, the reactions can occur through the diffusive (Langmuir-Hinshelwood) mechanism until the reactants find one another in the same site. The Eley-Rideal mechanism, in which a gas-phase species lands atop an adsorbed one and reacts with it, is also considered. Although the formation of molecular hydrogen via recombination of hydrogen atoms is listed, as we discuss below, this process need not be considered if one starts with hydrogen in its molecular form. Most of the reactions in Table 1 occur without activation energy E a barriers; for those that do possess activation energy, we use values taken from Caselli et al. (2002). The H, C, and CO network in Table 1 differs somewhat from that used by Stantcheva & Herbst (2004). These authors included two reactions to form carbon dioxide -the association reaction between CO and O and the reaction between O and the radical HCO -which are relatively unimportant in our calculations mainly because the large activation energy cannot be efficiently tunneled under by the heavy reactants. We also added a few reactions to make the system more realistic. The most important reactions in this network concern the gradual hydrogenation of CO to produce methanol, which has been studied in the laboratory (Watanabe & Kouchi 2002;Fuchs et al. 2007).
In our simulations, which pertain to a source at a temperature of 10 K or 15 K, we assume that only H and O can hop, and ignore the movement of all other species. The rate of hopping depends exponentially on the energy barrier between adsorption sites; this barrier is significantly lower for these physisorbed atoms than for the other surface species followed. At higher temperatures, the hopping of heavier species begins to become important . Likewise, thermal evaporation of almost all species is unimportant because it depends exponentially on the desorption energy, and this energy is much too large for species heavier than hydrogen at 10-15 K. In fact, we need only follow the evaporation of hydrogen atoms. Both the hopping barrier and the evaporation energy for H and O atoms are not constant during our simulation.
Since we follow the morphology of the growing mantle, we find that among monolayers that lie next to one another, products of a reaction in one of the monolayers may lie atop other reactive species in the next lower monolayer and be able to react with them. The OH + H 2 CO → HCO + H 2 O reaction in Table 1  H 2 CO molecule lies below the original oxygen atom, the newly-formed OH product will directly react with H 2 CO. What is interesting is that OH and H 2 CO are assumed to be stationary under these conditions, but they can react by the diffusion of H atoms.

Initial conditions and some basic assumptions
It is common in pseudo-time-dependent models of dense cloud chemistry to assume that the initial form of hydrogen is molecular, having been formed on grain surfaces during the diffuse cloud stage or somewhat later as the dense core was formed. In dense cloud cores, H 2 is consumed to some extent by the gas phase chemistry. Indeed, previous model calculations at 10 K show that cosmic ray ionization of H 2 leads to a residual atomic abundance of atomic hydrogen of ≈ 1 cm −3 while that for H 2 is maintained at 10 4 cm −3 . But does this maintenance depend on the continuing surface formation of H 2 under dense cloud conditions? Assuming a standard grain size of 0.1 µm and a standard gas-to-dust ratio, a sticking coefficient of unity for H atoms, and unit efficiency for the conversion of H to H 2 on grain surfaces, we obtain that at most 10% of the initial H 2 abundance can be produced within 10 8 yr. For a more reasonable time scale of 10 6 yr, the amount of H 2 produced is 100-fold less. Despite this minimal production, no depletion of H 2 has ever been reported even at higher temperatures when the formation efficiency drops significantly, which indicates that gas-phase chemistry during the dense cloud stage does not change the abundance of H 2 significantly. So, by this argument, molecular hydrogen formation on dust grains in dense clouds can be neglected, because it does not affect the gas. In order to confirm the unimportance of H 2 formation and to understand the strength of the coupling between gas-phase chemistry and surface chemistry, we performed two calculations at 10 K using the rate equation approach. In the first calculation, designated Model 1, we used the osu.2005 gas-grain code , which includes over 4000 gas-phase reactions (see http://www.physics.ohio-state.edu/∼eric/research.html) and hundreds of surface reactions, both of which are treated by the rate-equation method. Although all gasphase species were allowed to accrete with unit sticking efficiency, we included only those surface reactions in Table  1. In addition, atomic hydrogen was allowed to evaporate in competition with reaction. In the second calculation, designated Model 2, we removed all the surface reactions and ran the network again with the following additional changes: atomic hydrogen was not allowed to evaporate, and the accretion of H 2 and He onto grains was halted, because we assume that these two species simply land on grains and quickly evaporate again. The initial fractional abundances for both models are given in Table 2. These are based on elemental abundances considered by  based on unpublished work of Wakelam and Herbst. As compared with more standard low-metal elemental abundances, there is in general less depletion of the heavier elements. The constant density used is n H =2 × 10 4 cm −3 and the cosmic ray ionization rate ζ is set to 1.3 × 10 −17 s −1 . We will use the same initial abundances and total density throughout the paper. Table 3 shows the comparison of fractional abundances of selected species for these two models at a time of 10 6 yr. We can see that the difference is very limited; only the atomic hydrogen abundance is affected significantly; its abundance is greater in Model 1 because evaporation is allowed. But the extra atomic hydrogen does not make a significant difference for other species. This agreement shows that the gas-phase abundances remain mostly the same regardless of whether or not surface chemistry (Stantcheva & Herbst 2004) is considered as long as the temperature is sufficiently low that evaporation is not efficient for heavy species and non-thermal desorption is not included. This latter constraint is an important one because non-thermal desorption of surface ices may indeed be non-negligible for a few gasphase species such as methanol. Nevertheless, in this first attempt to wed the Monte-Carlo microscopic approach with gas-phase chemistry, we ignore the process and focus attention on the surface abundances. Under these conditions, we can infer that the details of surface kinetics are relatively insignificant for gas-phase abundances, whether or not they are handled by the rate-equation approach or by our more quantitative stochastic technique. The only important influence of grains on gas-phase chemistry occurs by the accretion of species onto grains, which eventually reduces the abundance of gas-phase species.

Calculation of Gas-Phase Chemistry
Since the two different approaches to surface chemistry give very similar gas-phase abundances except for atomic hydrogen, we can start with a rate-equation approach to the gas- 3.7 × 10 −10 4.4 × 10 −10 HCN 6.1 × 10 −10 7.5 × 10 −10 NH3 4.2 × 10 −9 4.3 × 10 −9 H2S 1.8 × 10 −12 2.0 × 10 −12 OH 1.6 × 10 −9 1.8 × 10 −9 C3H 4.1 × 10 −11 4.5 × 10 −11 phase abundances coupled with accretion of all species onto grains with unit sticking efficiency and no subsequent evaporation of atomic hydrogen, as in Model 2. The assumption that evaporation of H does not occur at all derives from the work of Chang et al. (2005) and , who showed that at a grain temperature of 10 K the recombination efficiency for molecular hydrogen is quite high for realistic surfaces. Nevertheless, it is not obviously true for all surfaces (see Table 3), and iteration with the actual Monte Carlo algorithm will be used to achieve better results.
The error in the assumption that hydrogen atoms do not return to the gas phase after accretion can be calculated as follows. Suppose that F evap (t) is the true time-dependent rate of atomic H evaporated from a grain. The evolution of the gas-phase abundance of atomic hydrogen n g (H) can then be written as where k acc is the accretion rate coefficient and the first two terms on the right-hand-side of the equation refer to gasphase formation and destruction processes. Formal integration of this equation yields represents the ratio between the integrated evaporation rate and accretion rate of atomic hydrogen. This ratio can be easily determined during a Monte Carlo simulation, as discussed later, by counting the number of evaporating and accreting atoms over the period of time. If ρ is much less than one, the error is small and the approximation of assuming F evap (t) to be zero is justified. If, however, ρ is large, the abundance of gas-phase atomic hydrogen can have a large error. Although this error will not affect other gas-phase abundances substantially, the incorrect hydrogen abundance will not be suitable for the subsequent surface chemistry simulation, which depends on the accretion flux of H atoms. In this case, an iteration has to be performed to obtain the correct gas-phase H abundance as a function of time. As a first-order approximation, we assume the ratio of evaporating H to incoming H to be constant throughout the simulation and equal to ρ determined over the entire time scale of the simulation. This approach corresponds to taking the average evaporating effect into account. It is undertaken because the Gear routine used for the gas-grain calculation cannot solve ordinary differential equations with variables explicit in time. In other words, a variable 1− ρ(t) factor cannot be inserted into the current gas-grain code. After ρ is obtained from the Monte Carlo simulation, we reduce the deposition rate coefficient of atomic hydrogen, k acc , to k acc (1 − ρ) and run the simulation again. If the new value of ρ is close to the previous one, the system is converged, otherwise the cycle is repeated with the new ρ until convergence is achieved. If ρ has no time dependence, the system converges to the true solution. However, ρ is typically time dependent and this time dependence cannot be included in the gasgrain network calculation. In order to gauge the resulting error, we define a time-dependent ρ 2 (t), which is the ratio of the number of evaporating H atoms to incoming H atoms counted during a "short" time period. We choose this "short" time to be 2 × 10 3 yr, which is 1/100 of the total simulation time of 2 × 10 5 yr. The time is short enough to determine the time dependence of ρ 2 (t) and long enough to suppress most of the statistical noise. A shorter time scale (1×10 3 yr) shows some more noise while a longer time scale (5 × 10 3 yr) shows less. With the rate equation approach, this method is equivalent to the integration If ρ 2 (t) is weakly dependent on time, t, then ρ has little time dependence, thus the resulting error is small. Otherwise, the error is large even if convergence has been achieved. With the current Gear algorithm used in the gas-grain model, we cannot use the parameter ρ 2 (t) in the actual simulation. A new gas-grain algorithm would have to be designed to follow the gas-phase abundance of H with the time-dependent ρ 2 .

Deposition of Gas-Phase Species
Assuming the sticking coefficients for each species to be 1, the deposition (accretion) flux of species i onto each absorption (lattice) site is equal to the product f i N , where f i is the incoming flux of species i in units of ML s −1 , while N , the number of sites on the grain, is given by where r is the radius of the grain and s is the surface site density. In our simulations, we start with bare olivine as the granular surface, and use a site density s of 1.5 × 10 15 cm −2 , which is also used in our standard gas-grain network. The flux f i can be calculated by where T g is the gas phase temperature and m i is the mass of species i. A factor of four, applicable to the case of spherical grains, is often used in the denominator. However, this factor is not used in the accretion term in our standard gasgrain network, and for consistency we do not use it here. If f i is independent of time, the deposition is an homogeneous Poisson process, and the deposition interval ∆t can then be calculated using where X is a random number between 0 and 1. In our simulations the deposition rate is time-dependent due to the time dependence of the gas-phase abundances of the impinging species. These gas-phase abundances are meanfield values, but the flux is still stochastic in arrival time.
In reality, there is already a distribution in the flux due to a distribution in velocity. The next deposition time, t next , satisfies the equation, with t prev the previous deposition time, and the process is referred to as an inhomogeneous Poisson one. It can easily be shown that Eq. (8) is also valid when f i is independent of time. With an analytical form of f i (t) known, t next can be determined numerically. However, only a numerical solution is obtained from the gas-grain network calculation, thus interpolation is needed to approximate f i . The current gasgrain network outputs gas phase abundance of species i, at different time t j . We employ linear interpolation to get the abundance of species i at time t, which is between t j and t j+1 . The flux f i (t) can then be calculated as where ∆t j = t j+1 − t j . We can formally write f i (t) for any t as, where S(t) is the step function and b ij = fi(tj+1)−fi(tj ) ∆tj . Combining Eqs. (8) and (10) gives Let us call Then where r = − ln X N . Eq. (12) can be solved numerically in the following manner. Suppose t k ≤ t prev < t k+1 . First, we calculate I(t prev , t k+1 ). If this is larger than r, we know that t next must be less than t k+1 . We can determine t next by expressing f i (t) in terms of t prev rather than in terms of t k . Integration then yields where If I(t prev , t k+1 ) < r, then deposition must happen later than t k+1 . We define a new quantity, r ′ , or If integration from t k+1 to t k+2 leads to a larger value than r ′ , t next < t k+2 and an expression for t next similar to eq. (14) can be obtained. If integration to t k+2 leads to a smaller value than r ′ , t next is larger than t k+2 , and the sequence must be continued until t next is bounded by two successive times t m and t m+1 .
To check for convergence, the original results of the gasgrain network are rerun twice with smaller time intervals, first 1/4 of the original ∆t j and then 1/4 of the second one. Convergence is achieved in all cases.

Surface Models
A detailed description of flat and rough surface models has been given in . Here we start with either a flat surface or their surface D, which is the roughest one and is generated by another Monte Carlo program. The diffusion barrier, E b , and the desorption energy, E D , for the two diffusing species -atomic hydrogen and atomic oxygen -are in our models dependent on the topological roughness. Starting from either a flat or rough bare olivine surface, a mantle gradually develops. The formation of surface species, especially water ice, on the surface has two consequences. First, the chemical nature of the surface changes, resulting in different energetics and therefore different values for E b and E D . Secondly, the surface topology changes, as more and more species are formed and roughness develops even from the initially flat surface due to fluctuations. Thus, E b and E D are inherently stochastic in the simulation in that they change for the site at the top of a given column as the stochastic calculation proceeds. For simplicity, we assume that only the olivine substrate and the H 2 O molecules (the main mantle constituent) influence E b and E D , either by lateral interactions or by a vertical interaction with the substrate.
The parameters E b and E D are given by respectively. Here E 0 b and E 0 D are the diffusion barrier and desorption energy on a surface without lateral bonds, while E 1 L and E 2 L are the lateral bonds imposed by horizontal olivine and H 2 O neighbors respectively. In this simulation, E L has a value of 0.1E D . The indices i and j are the number of olivine and H 2 O neighbors respectively. Since species other than H and O have large diffusion barriers even without lateral bonds (Ruffle & Herbst 2000), we ignore their movement in the temperature range studied here (10-15 K). The evaporation of O is ignored. Table 4 summarizes the parameters used for H and O on the ice and olivine surfaces. The diffusion barrier and desorption energy of H on olivine are from Katz et al. (1999). The desorption energy for H on water ice is not well constrained. Hollenbach & Salpeter (1970) found a value of 450 K, Buch & Czerminski (1991) found a value of 500 K, while Al-Halabi et al. (2002) found a value of 400 K. We take the average of these measurements. Recently, Perets et al. (2005) deduced a much higher value of 720 K. The oxygen values are from Tielens & Hagen (1982). For comparison, we also simulate a second system by ignoring the lateral bonds and only considering the change of the chemical nature of the surface, i.e. the vertical bond. This system starts from a flat surface. Although it gains some roughness from fluctuations as the ice layers grow, the system is kinetically equivalent to a flat system since there is no penalty for going uphill and no gain in falling into a valley.

Surface reactions with activation energy barriers
In gas-grain codes based on a rate equation or master equation approach to the surface chemistry, the common method to calculate diffusive rate coefficients for reactions with an activation energy barrier is to multiply the diffusion rate by a tunneling factor: where a is the width of the rectangular barrier, typically chosen to be 1Å, µ is the reduced mass, and E a is the reaction activation energy (Hasegawa et al. 1992). Note that the word "diffusion" is normally used to refer to motion over a whole grain, whereas the word "hopping" refers to motion from one potential minimum to the next. A physical interpretation of the rate coefficient is that the two species have a probability equal to the tunneling expression to react with each other if they diffuse into the same potential minimum. Thus, the rate coefficients are linearly dependent on the diffusion rate. Unlike the typical situation in the gas phase, two species that meet each other on a surface can stay in each other's vicinity for quite a while depending on the hopping rate to nearby lattice sites. During this time they can have many chances to react. It is therefore more physical to model reactions with barriers as processes competitive with hopping out of the potential minimum. Such an approach has been done by Awad et al. (2005) using rate equations. Here we incorporate this approach by considering the competition among Poisson processes.
In particular, the competition among reaction, hopping or evaporation can be implemented as follows. The reaction rate with tunneling can be written as where ν 1 is the attempt frequency for reaction. The competing rates for hopping and evaporation of the two species are given by and respectively, where T is the surface temperature, ν 2 is the attempt frequency for hopping, and ν 3 is the attempt frequency for evaporation, and the diffusion barrier and evaporation energy pertain to either of the two species. The stochastic algorithm for diffusion treats all directions with equal probability. Note that hopping is treated classically rather than by quantum mechanical tunneling because, unlike reaction, hopping is thought to occur over a broad and shallow potential. Quantum chemical calculations are being undertaken to verify this point (Woon, private communication). Assuming that only one species is moving, the probability for reaction to occur instead of diffusion or evaporation is We assume ν 1 = ν 2 = ν 3 = 10 12 s −1 , but since there three different types of processes the attempt frequencies can be different. The competition is implemented in the Monte Carlo algorithm as follows: a random number between 0 and 1 is generated. If this number is smaller than p the reaction occurs. Otherwise evaporating or hopping takes place. A competition between these two processes must also be evaluated.

Monte Carlo Simulation
The Monte Carlo algorithm starts with determining the first deposition time for H, O, and CO, as explained in Section 4.2. Then the first deposition is executed by placing the particle on a randomly picked site of the lattice, which is chosen to have a lattice of size 100 × 100. This size corresponds to a rather small grain of 0.015 µm, and a mean distance between atoms in adjacent lattice sites of ≈ 5Å. Typical activation energy barriers are shorter than this in width; the widths of the diffusion barriers are currently under investigation (Woon, private communication).
We have already shown that the lattice size used is large enough that there is little size dependence to the results for both flat and rough olivine surfaces (Chang et al. 2005;. Nevertheless, we redid our simulations for lattice sizes of 50 × 50 and 200 × 200 lattices at 15 K to confirm this point. For much smaller grains, the socalled accretion limit is reached, where the average mantle abundances of reactive species is less than unity, leading to a drop in reaction efficiency. For everything but the surface stochastic simulations, the dust-to-gas number density is taken to be 1.32 × 10 −12 and the grain size the nominal one of radius 0.1µm. As the clock moves forward, other species are deposited on the lattice, and species already on the lattice can hop from site to site or evaporate. By hopping on a rough surface, we mean motion on the top level of any occupied site, whether is is occupied by the original silicate or by any other species. Thus hopping can include climbing or falling in a vertical sense. For very rough surfaces, the climbing or falling length can be up to 20-30 molecules. Evaporation and hopping are assumed to be Poisson distributed, with time intervals between events ∆t = − ln(X) b , where X is a random number between 0 and 1, and b is the occurring rate of that process as given by Eqs. (21) and (22). For H atoms, which are allowed to both hop and evaporate, these two processes are treated first as one combined process with a rate b 1 + b 2 and then another random number X ′ is called to decide which one will occur first. If X ′ is smaller than b1 b1+b2 , the hydrogen atom hops and it evaporates in the other case. The total number of H atoms that evaporate is counted to calculate ρ. Oxygen atoms are only allowed to hop. The time interval between hopping is solely determined by the hopping rate.
If a hydrogen or oxygen atom hops to a site where there is already one species that can react with it without a barrier, the reaction happens immediately. If either atom hops into the same site as a species with which it can react with an activation energy barrier, the competition among reaction and other processes occurs as discussed in the previous section. Reaction products other than H 2 remain on the surface and can cover species on the same sites. Some reaction products can react further with species they cover. For instance, if OH is produced by a reaction in which H hops into a site occupied by an O atom, and there is an H 2 CO molecule lying under the O atom, then another reaction is allowed to take place, producing HCO and H 2 O (see Table 1).
If a gas-phase species lands atop a bound surface species, reaction can occur according to the Eley-Rideal mechanism. No distinction is made between reaction by hopping (Langmuir-Hinshelwood mechanism) and by landing. If reaction does not occur upon landing, the incoming particle covers the one already present. No reaction is allowed between species covered by at least one other species. For instance, assume that an oxygen atom is buried by a CO molecule. After being buried, the oxygen atom cannot react with H or O even if they occupy the site above the CO.
The Monte Carlo simulation is performed by executing the following sequence of steps: 1) The gas-grain code is run for Model 2 through 10 6 yr, 2) The gas-phase results for H, O, CO are used to run the Monte Carlo simulation of deposition, surface chemistry, and evaporation of atomic hydrogen for 2 × 10 5 yr, 3) An average value of ρ over the full time of the stochastic calculation is obtained, 4) If ρ is close to zero or is converged, the calculation is stopped, otherwise: 5) The gas-grain code is rerun with a deposition rate reduced by a factor of 1 -ρ, and the simulation is also run with a reduced deposition rate for H,  Fig. 1. Fractional abundances of gas-phase species that accrete onto grains for Model 2. Solid lines are for 10 K while dotted lines are for 15 K.
6) The procedure then returns to step 3.
Note that the gas-grain code is used to obtain gasphase results; surface results are subsequently obtained by stochastic simulation. The Monte Carlo code is run for a shorter time than the gas-grain code for two main reasons: (i) it is very time-consuming computationally, and (ii) the parameter ρ 2 (t) becomes significantly time-dependent after 2 × 10 5 yr. Finally, the iterative procedure used is forced on us by the current structure of the gas-grain code, which does not allow parameters with explicit time dependence.

Results and Discussion
The simulation starts in earnest with runs at 10 K and 15 K in which the fluxes of hydrogen atoms, oxygen atoms, and CO striking the grain as a function of time are taken from Model 2. Figure 1 shows the fractional abundances of gaseous H, O and CO as functions of time at 10 K and 15 K obtained from the Model 2 run. As can be seen, the atomic hydrogen abundances do not change much over the first 10 5 yr, while the decline in the atomic oxygen abundance turns on after 10 4 yr. A significant amount of gas-phase CO is produced after 10 3 yr. These abundances are those used to calculate fluxes for the Monte Carlo simulation of the surface chemistry. In this run, both surface chemistry and evaporation of atomic hydrogen occur. As discussed in Section 4.1, the evaporation of atomic hydrogen is followed stochastically and averaged over time to determine the parameter ρ. If this parameter is significantly different from zero, then further runs are necessary with a reduced hydrogen flux to achieve convergence. The closer the value of ρ to unity, the slower the convergence.
We first simulate the surface chemistry starting with a flat olivine surface and using a surface-adsorbate system with no lateral bonds as the ice develops. We then redo the simulation starting from a rough olivine surface and containing lateral bonds. The latter surface, which is rough at all stages, will just be referred to as "rough".

Surface without lateral bonds
At 10 K, ρ is initially calculated to be 3 × 10 −4 , which is significantly less than 1. Thus, the gas-phase fractional abundances shown in Fig. 1 are also reasonably good solutions of the gas-grain network with evaporation. However, at 15 K, ρ is 0.37, which is comparable to 1, because evaporation is more important. Following the prescription in Section 4.1, we reduce the hydrogen-atom accretion rate to 0.63 of the previous value, and perform the gas-grain simulation a second time with the Monte Carlo approach to the surface chemistry. In this simulation, ρ now is calculated to be 0.40, which indicates that reasonable average convergence has been reached. Moreover, we can look at the time-dependent ρ 2 (t) to check that it is reasonably independent of time, another test of convergence. Figure 2 shows ρ 2 (t) as a function of time for 15 K. From this figure, we see that the error of our approach is small since ρ 2 (t) converges and does not change much with time, although there is an increase after 10 5 yr. The increase of ρ 2 (t) derives from the decrease of O in the gas phase, which decreases its flux onto grain surfaces, so that more H atoms can evaporate before combining with surface O. The gaseous fractional abundance of atomic hydrogen after the initial Monte Carlo simulation and its iteration at 15 K is shown in Fig. 3. As expected, the H fractional abundance increases after the initial surface simulation since evaporation is allowed. For example, at 10 5 yr, the H fractional abundance exceeds 10 −4 whereas it is closer to 7 × 10 −5 for Model 2. After the iteration, however, the atomic H abundance has become closer to its original abundance at all times. Even with the increased gas-phase H abundance after the initial simulation, the abundance for each surface species changes less than 15% at 15 K. This indicates that the change in the gas-phase atomic hydrogen abundance due to the change in H evaporation from the grain surface has little effect on the surface abundances. The other gas-phase abundances do not change much either because atomic hydrogen only plays a minor role in the gas-phase chemistry. Figure 4 shows the surface abundances of CO, H 2 CO, CH 3 OH and CO 2 on grain surfaces as functions of time at 10 K and 15 K. After some initial fluctuations, especially  Fig. 4. Mantle abundances for selected species in monolayers as a function of time on a surface-adsorbate system with no lateral bonds. Circles are for CO, triangles pointed upward are for CO 2 , squares are for H 2 CO, while diamonds are for CH 3 OH. The left panel refers to results at 10 K while the right panel refers to results at 15 K. at 10 K, the population of each carbon-bearing species (in monolayers per grain) increases monotonically with time. At 10 K, methanol exceeds CO after 500 yr and eventually becomes the most abundant carbon-bearing surface species, while CO 2 is also abundantly produced. At 15 K, on the other hand, methanol is always the least abundant of the four carbon-bearing species shown. Table 5 shows the populations of all the mantle species at a time of 2 × 10 5 yr. It is easily seen that the mantles are dominated by water ice, with molecular oxygen second in abundance. The total number of species at 15 K is slightly larger because of the larger accretion rate at higher temperatures and the still minimal evaporation rate. As already seen in Fig. 4, the abundance of methanol is highly dependent on temperature; its abundance drops by a factor of seven as the temperature increases from 10 K to 15 K. The higher temperature makes the formation of methanol less efficient because diffusion rates increase quickly at higher temperatures while the tunneling reaction rates for the critical reactions with activation energy -H + CO and H + H 2 CO -are independent of temperature. Thus, competition favors hopping of H over reaction as the temperature  increases. This decrease in the methanol abundance with increasing temperature is not exactly what is measured in the laboratory, although the conditions there are somewhat different from those studied here (Fuchs et al. 2007). Our analysis for the calculated temperature dependence of methanol also leads to more HCO and H 3 CO at 10 K than at 15 K. In addition, the increasing temperature makes hydrogen atoms hop faster, thus the activation energy-less reactions with HCO and H 3 CO become more efficient. For HCO and H 3 CO, there is another mechanism for destruction: O atoms become mobile at 15 K, and can react with both radicals. Unlike other simulations, much CO 2 is produced in our calculation, mainly because products from one reaction can react further with species under them; thus, species that are not mobile can react with each other. In this case, the reaction of importance is that between CO and OH, the barrier of which is low enough that the probability of reaction is near unity. For the case of O 2 , its formation by surface recombination of oxygen atoms is helped by an increase in temperature since the surface atoms move around significantly more quickly at 15 K. At 10 K, there is a significant amount of atomic oxygen in the grain mantles, most of which is buried and not able to react further. Figure 5 shows the abundances of the major carbonbearing species in individual monolayers, starting from the surface and moving upward. These differential abundances of each carbon-containing species typically increase initially as a function of height from the original surface but peak before the top layers are reached at both 10 K and 15 K. The bottom layers hardly have any molecules with carbon, because gas-phase CO takes some time to produce. The top layers contain fewer species overall because these layers are not yet fully occupied.

Rough surface
The calculated ρ at 10 K is about 1 × 10 −4 , which is much smaller than unity, so that no iteration is required. At 15 K, ρ is 0.14, which is still smaller than 1. But we perform an iteration anyway. After one iteration, ρ is still 0.14, so the calculation is converged.
The time-dependent surface abundances of CO, H 2 CO, CH 3 OH and CO 2 are shown in Fig. 6. At 10 K, the figure is similar to Fig. 4, which represents the analogous results on a surface-adsorbate system without lateral bonds, while at 15 K, on the other hand, more methanol is produced than in the system without lateral bonds, because there are stronger binding sites due to lateral bond interactions and these stronger binding sites help to add hydrogen atoms to molecules when there is a reaction barrier, since the competitive diffusion rate is lowered.
The surface abundances of all studied species at 2 × 10 5 yr are shown in Table 6. Compared with Table 5, there is not much difference at 10 K except that more hydrogen atoms are trapped on the surface. Indeed, for any reasonably-sized grain, the mantle population of 0.019 corresponds to far greater than 1 atom per grain. At 15 K, however, as in Fig. 6, one can see that more methanol is produced than without lateral bonds.
The abundance per monolayer of the carbon-bearing species at different layers is shown in Fig. 7 for a time of 2 × 10 5 yr. For 10 K, this gives a very similar result to the system without lateral bonds in Fig. 5. At 15 K, more methanol is formed at each layer because of the slower hopping rates caused by the lateral bond interaction.  Fig. 7. Abundance per monolayer in fractions of a monolayer for selected species at each layer on a rough grain surface at 2 × 10 5 yr. Monolayers are numbered from zero at the initial surface. The symbols representing the species are the same as in Fig. 5. The left and right panels represent results at 10 K and 15 K, respectively.
At the rough surface at 15 K, the populations of species as a function of monolayer have a distinct feature, in that they drop to zero abruptly, unlike the other cases. This effect happens because the oxygen atoms are able to hop efficiently out of sites with small lateral bond interaction, located above surrounding sites, to sites with larger lateral bond interaction, located in valleys below surrounding sites, where they are trapped. So, sites that are lower in vertical height above the original surface tend to have a larger probability to be occupied by oxygen atoms. These atoms can remain or react with other atoms to fill the up the valleys. Thus, the diffusion-reaction mechanism automatically "smooths" the surface strongly at 15 K.
In addition to the rapid dropping in monolayer abundance of the four major carbon-bearing species, there is a second distinct feature at 15 K. Both the monolayer abundances of CO and H 2 CO increase strongly towards the highest monolayers before suddenly dropping. This effect occurs because it is relatively difficult to bury CO and H 2 CO under higher layers due to their reactivity with OH radicals.

Comparison with observations
Our primary interest in this work has been to introduce a new method to combine rate equations and microscopic Monte Carlo simulations in a gas-grain model of the chemistry of cold interstellar cores. Toward this goal, a small grain-surface chemical network was chosen for simplicity. However, we would still like to compare our results with astronomical observations. Since, again for simplicity, we do not consider non-themal desorption mechanisms here, as we have done recently with our latest gas-grain rate-equation code , it is unnecessary to compare our gas-phase results with observations since the results are very similar to those of Ruffle & Herbst (2000). Rather, it makes more sense at this stage to compare our mantle results with infrared absorption measurements. These results are insensitive to non-thermal desorption, which only affects mantle abundances slightly. We are particularly interested in determining whether or not a rough surface model can help to obtain abundances closer to the observations. In Table 7 we list calculated and observed fractional abundances for four ices: H 2 O, CO 2 , CO, and CH 3 OH. The water abundance is tabulated as the fractional abundance with respect to n H while the abundances of the other ices are listed as percentages of the water ice abundance. Our theoretical results are calculated at 2 × 10 5 yr. For comparison, master equation results at 10 K for a time of 3.2 × 10 5 yr are also listed (Stantcheva & Herbst 2004). The observational data for W33 A are from Dartois et al. (1999) and Gibb et al. (2000), while the data for Elias 16 are from Gibb et al. (2000).
As compared with the master equation results, which are based on a macroscopic stochastic approach in which no distinction is made for local site differences and differences according to monolayer, the Monte Carlo results are a distinct improvement. The master equation abundance for CO is very small at 10 K, although better results are obtained at 20 K. As for CO 2 , the results are low at both 10 K and 20 K. For the Monte Carlo models, agreement is reasonably good regardless of the type of surface model except perhaps for CO. The agreement for CO 2 is especially important given the ubiquity and high abundance of this ice component. The W33A results are closer to our 10 K model predictions whereas the Elias 16 results are closer to the 15 K predictions. The higher temperature leads to more CO and less methanol. Still, CO is somewhat underproduced in the models. For the rough surface, increasing the lateral bond strength to 0.4E D from its standard value of 0.1E D can increase CO from 0.75 percent to 3.0 percent of the H 2 O ice while keeping other species in good agreement with observation at 10 K.

Conclusions
For the first time, a gas-grain chemical simulation of a dense interstellar cloud core has been performed with a combined rate equation-microscopic Monte Carlo approach. In this approach, the local characteristics of the surface and all monolayers accreted onto it can be followed. The surface network is chosen to be the H, O and CO system, which has been widely studied and is the key network to form methanol on grain surfaces. Our results for the total mantle abundances of the carbon-bearing species CO, CO 2 , and CH 3 OH as well as for water are in reasonable agreement with observation in two well-studied protostellar sources. In fact, the agreement is superior to that achieved by use of a macroscopic Monte Carlo approach known as the master equation treatment (Stantcheva & Herbst 2004).
The approach utilized here is based on a number of simplifications made possible mainly by the low temperatures (10 K, 15 K), the starting chemical abundances, and an unusual aspect of the gas-phase chemistry. First, if we start with hydrogen mainly in molecular form, we need not consider subsequent H 2 formation since it is a small effect. Secondly, at the low temperatures used here, the only reactive species that can evaporate from grain surfaces is atomic hydrogen. Thirdly, the residual atomic hydrogen abundance in the gas does not strongly affect the abundances of all other gas-phase species. These simplifications allow an iterative solution in which gas-phase abundances computed initially do not change significantly except for the case of atomic hydrogen. Our iterative procedure involves starting with a gas-grain calculation in which no surface chemistry or H evaporation is allowed. The time-dependent abundances of H, O, and CO are then allowed to deplete onto grain surfaces and the surface chemistry and evaporation of H are treated by the continuous-time random walk Monte Carlo procedure. The evaporation rate of hydrogen is used to refine the abundance of atomic hydrogen in the gas by an average correction over the time period of the calculation, an approximation made possible by the neglect of H 2 formation. This refined abundance is then used, if necessary, to redo the Monte Carlo calculation of the surface chemistry until convergence is achieved. If the evaporation of atomic hydrogen from grain surfaces is strongly time dependent, a more complex bootstrap procedure would be called for. Were other species to evaporate, as occurs at higher temperatures, and were non-thermal desorption to be included, a process that occurs even at low temperatures, we would probably have to devise a more intricate, non-iterative, method of integrating the gas-phase abundances simultaneously with the Monte Carlo determination of surface abundances. Such improvements are under active consideration, including an approach similar to the effective rate coefficient method Chang et al. 2006). Finally, only a small number of surface processes are currently used; inclusion of a much larger number of processes by the Monte Carlo procedure would be exceedingly time-consuming whether or not changes in gas-phase abundances are also calculated. Whether a suitable procedure that will allow us to remove all of these simplifications can be achieved in the near future remains uncertain at his point. a(-b) refers to a × 10 −b 1 water abundance with respect to nH; other ices listed as percentage of the water ice abundance