Incorporation of stochastic chemistry on dust grains in the Meudon PDR code using moment equations
I. Application to the formation of H and HD
F. Le Petit^{1}  B. Barzel^{2}  O. Biham^{2}  E. Roueff^{1}  J. Le Bourlot^{1}
1  Observatoire de Paris, LUTH and Université Denis Diderot,
Place J. Janssen, 92190 Meudon, France
2  Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel
Received 3 April 2009 / Accepted 29 June 2009
Abstract
Context. Unlike gasphase reactions, chemical reactions taking place on interstellar dust grain surfaces cannot always be modeled by rate equations. Because of the small grain sizes and low flux, these reactions may exhibit large fluctuations and thus require stochastic methods such as the moment equations.
Aims. We evaluate the formation rates of H,
HD, and Dmolecules on dust grain surfaces and their abundances in the gas phase under interstellar conditions.
Methods. We incorporate the moment equations into the Meudon PDR code and compare the results with those obtained from the rate equations.
Results. We find that within the experimental constraints on the energy barriers for both diffusion and desorption and the density of adsorption sites on the grain surface, H,
HD and Dmolecules can be formed efficiently on dust grains.
Conclusions. In a wide range of conditions, the moment equation results agree with those obtained from the rate equations. However, for a range of relatively high grain temperatures, there are significant deviations: the rate equations fail, while the moment equations provide accurate results. The incorporation of the moment equations into the PDR code can be extended to other reactions taking place on grain surfaces.
Key words: astrochemistry  ISM: general  ISM: clouds  ISM: molecules
1 Introduction
Most of the detected interstellar molecules are in the gas phase.
The chemical balance between gas phase species in interstellar clouds
is commonly described by a set of rate equations, using
the known rate coefficients of gasphase reactions.
Many studies have been devoted to
the resolution of the chemical rate equations
to derive the compositions of interstellar clouds within
specific physical conditions.
Two main approaches are used by modellers:
the timedependent approach and the steadystate approach.
In the timedependent approach, one
follows the time evolution of the abundances of various atoms and
molecules by numerically integrating the corresponding set of
first order coupled differential
equations, from specific initial conditions.
In the steady state approach,
one directly solves the set of algebraic equations
that are obtained when the time derivatives of the
abundances are set to be zero.
Comparison with observations then allows them to
test the relevance of the various hypotheses.
However, some important reactions that give rise
to the formation of molecular hydrogen, ice mantles,
and certain organic molecules do not take place in
the gas phase, but on the surfaces of dust grains.
To account for these surface processes, one needs to
incorporate the grainsurface reactions into the gas
phase models of interstellar chemistry.
It is thus attractive to use rate equations for the surface reactions,
which enable us to couple surface reactions to gasphase
reactions in a straightforward way.
This approach is indeed valid in the case of sufficiently
large grains, when the number of reactive atoms and
molecules of each species on a grain is sizeable.
However, in the limit of small grains,
when the number of reactive atoms and molecules on a grain is small,
the rate equations are not always suitable.
Rate equations simply ignore the fluctuations
in the number of reactive species on the grains.
This problem has been
discussed by several authors (Stantcheva et al. 2002; Charnley et al. 1997; Stantcheva et al. 2001; Shalabiea et al. 1998; Tielens & Hagen 1982; Caselli et al. 1998).
To overcome these difficulties,
a master equation approach was proposed
and applied specifically to the formation of
molecular hydrogen.
This approach is suitable for the simulation of diffusive
chemical reactions on interstellar grains
in the accretion limit
(Biham et al. 2001; Green et al. 2001; Biham & Lipshtat 2002).
The master equation approach takes into account
the discrete nature of reactive species on the surfaces
of grains as well as the fluctuations in their populations.
This approach was later used for a larger chemical network
on grains, leading to the formation of methanol
and its deuterated versions
(Stantcheva et al. 2002; Stantcheva & Herbst 2003).
However, the number of equations included in the master
equation increases exponentially with the number of species
that are reactive on grain surfaces.
Therefore, it is not suitable for incorporation in codes
that include complex networks of grainsurface reactions.
In addition to the master equation method,
a Monte Carlo method was proposed by Tielens & Hagen (1982)
and applied to interstellar modelling by Charnley (2001),
Caselli et al. (2002), and Stantcheva & Herbst (2003).
Cuppen & Herbst (2005)
developed random walk models for molecular hydrogen formation
on grains that take into account the effect of surface roughness
on the diffusion and reaction rates.
They simulated these models using Monte Carlo methods
and showed that surface roughness tends to broaden the
range of temperatures in which molecular hydrogen formation
is efficient.
The advantage of the random walk models is that they enable them
to account in more detail for the microphysics of the grain
surface. In particular, they enable them to take into account the
entire distribution of binding energies and diffusion barriers
for H atoms. Detailed models of this type provide useful insight
but it is not feasible to include them in large models of
interstellar chemistry because of their complexity.
In contrast, the rate equation and moment
equation methods used in this paper, include a single energy
barrier for each process.
Finally, a semiempirical approach,
known as the modified rates method,
was proposed by
Caselli et al. (1998).
This method is easy to employ and was applied with mixed success by
Shalabiea et al. (1998), Stantcheva et al. (2002),
and Caselli et al. (2002).
Garrod (2008) studied a different version of the modified rate method.
The proposed moment equations method provides efficient stochastic simulations of grainsurface chemistry (Barzel & Biham 2007a,b). The method consists of only one equation for the average population size of each reactive species on a grain, and one equation for the rate of each reaction. It consists of a set of coupled ordinary differential equations, which resemble the rate equations. Therefore, they can be easily coupled to the rate equations of gasphase chemistry. Unlike the rate equations, the moment equations are linear, and thus easier to handle both in the timedependent and the steady state approaches. The moment equations were tested for the reaction network that gives rise to ice mantles on grains, which consist of water ice, carbon dioxide, and methanol (Barzel & Biham 2007a,b). The stability properties of the moment equations as well as the accuracy of their steadystate solution under astrophysically relevant conditions were examined by extensive computer simulations and by comparison with the master equation results.
In this paper, we incorporate the moment equations into the Meudon PDR code (Le Petit et al. 2006,2002), a stationary model of Photon Dominated Regions (PDRs).
The model considers a stationary planeparallel slab of gas and dust illuminated by an ultraviolet radiation field coming from one or both sides of the cloud. It solves, at each point in the cloud and in an iterative way, the radiative transfer in the UV, taking into account the absorption in both the continuum by dust and discrete transitions of H and H_{2}. Explicit treatment is performed for C and S photoionization, H_{2} and HD photodissociation, as well as CO (and its isotopomeres) predissociating lines. The model also computes the thermal balance by taking into account heating processes, such as the photoelectric effect on dust, cosmic rays, and chemistry. It also accounts for the cooling due to infrared and millimeter emission of the abundant ions, atoms, and molecules. The chemistry is solved under steady state conditions and the abundance of each species is computed at each point. The excitation states of a few important species are then computed. The column densities of these chemical species and their emissivities/intensities are then calculated. To examine the applicability and relevance of the moment equations within the PDR code, we consider a simple network of grain surface chemistry, involving only H, D, H, HD, and D_{2}. We compare the results obtained from the PDR code with the moment equations, with those obtained when the rate equations are incorporated into the same code. Unlike previous studies of molecular hydrogen formation that assumed a single grain size, the moment equations enable us to take into account the full distribution of grain sizes. We show that in the case of relatively high grain temperatures, the rate equations are not suitable and the moment equations should be used. We also demonstrate the importance of the Langmuir rejection effect at very low grain temperatures, where the grain surfaces are saturated by hydrogen atoms and reaction rates are low.
The paper is organized as follows. In Sect. 2, we briefly review the gasphase and surface reactions involved in the formation of H, HD and Dmolecules. We present in Sect. 3 the moment equations for the HD system and describe their incorporation into the PDR model. Simulation results are displayed in Sect. 4. The results are summarized and discussed in Sect. 5.
2 Formation of molecular hydrogen and its deuterated versions
2.1 Astrophysical context
No efficient gasphase mechanism is at work to
form molecular hydrogen in the gas phase.
Therefore, reactions on grain surfaces
are necessary to account for the observed abundance
of H.
Hollenbach & Salpeter (1971) were among the first to
quantitatively describe the formation of molecular
hydrogen on the surfaces of spherical dust grains.
Here we assume that
the grain size distribution follows a power law of the form
n(a) = c a^{q},  (1) 
as suggested by Mathis et al. (1977). The prefactor c is given by
(2) 
where (g cm^{3}) stands for the volumic mass of the grains, while and are the upper and lower cutoffs of the grain radii distribution, respectively. We denote the proton density in the gas phase by (cm^{3}) and the dust to gas mass ratio by G=0.01. In this case, one can formulate the formation rate of H_{2} as a function of various grain parameters, such as the dust to gas mass ratio, the minimum and maximum values of the grain radii, and the volumic mass of the grains. If the exponent of the power law is q=3.5, analytical formulae are obtained as derived by Le Bourlot et al. (1995). Le Petit et al. (2002) specifically studied the D/HD transition occurring at the edge of a translucent cloud and in dense PDRS. They also studied the formation of HD molecules on the surfaces of dust grains, following the procedure of Le Bourlot et al. (1995). The resulting formulae are very similar to the pure molecular hydrogen case. Following previous considerations by both Watson (1974) and Dalgarno, Black, & Weisheit (1973), Le Petit et al. (2002) found that HD is far more efficiently formed in the gas phase, in a succession of reactions initiated by cosmic ray ionization of atomic hydrogen. When protons are formed, the sequence of reactions involved is the following:
(3) 
The forward reaction is exothermic with an energy of 3.7 meV, corresponding to 43 K, so that the reverse reaction may take place at moderate temperatures, reducing the D^{+}formation efficiency. D^{+}reacts efficiently with H to form HD:
(4) 
The reverse reaction may also occur, but the endothermic barrier is now 40 meV, equivalent to 464 K, so that the corresponding probability is negligible at low temperatures. The main destruction channel of HD in diffuse environments is photodissociation, which was computed by Le Petit et al. (2002). This simple gasphase scheme accounts for the observed HD/H_{2} abundance ratio in diffuse clouds and translucent regions, as found from Copernicus and FUSE observations (Watson 1974; Dalgarno et al. 1973; Lacour et al. 2005). In dark and cold regions, HD is the major repository of deuterium. It reacts with H_{3}^{+}to form H_{2}D^{+}, which is the starting point of a rich chemistry of deuterium fractionation (Roueff et al. 2005,2000; Millar et al. 1989; Caselli et al. 2002; Roberts et al. 2003; Stantcheva & Herbst 2003). Although the formation of HD in diffuse clouds occurs primarily in the gas, it has proposed that under certain conditions, grainsurface reactions may also contribute to an enhanced production of HD and Dmolecules. The proposed mechanism is based on the assumption that D atoms stick more strongly than H atoms so that their desorption rate is lower. This isotope effect was observed in various experimental situations (Amiaud et al. 2007; Hoogers 1995; Koehler et al. 1988). As a result, the residence time of D atoms on grains is expected to be longer than that of H atoms. The D/H abundance ratio on the grains is then enhanced compared to the gasphase ratio, and newly adsorbed H (or D) atoms are more likely to find D atoms already residing on the grains. This may give rise to an enhanced production of HD molecules.
2.2 The interaction between hydrogen/deuterium atoms and dust grains
To quantitatively describe the formation of H,
HD, and Dmolecules on grain surfaces, one has to introduce several hypotheses and parameters,
which we recall below. The typical velocities of H and D
atoms in the gas phase are given by
and
,
respectively.
Based on the simplifying assumption that all
grains are of the same radius, a,
one can compute the numerical density of spherical dust grains
(cm^{3}) as a function of
to be
(5) 
For a power law distribution of grain sizes with the exponent q=3.5 (Mathis et al. 1977), the numerical density of the grains becomes
(6) 
The number of adsorption sites on a grain is denoted by S. Their density on the grain surface, s (sites cm^{2}), is given by . The density of adsorption sites, s, and the distance between adjacent sites, d, are related by d^{2} = 1/s. The fluxes and (atoms s^{1}) of H and D atoms respectively onto the surface of a single grain are given by , where I = H, D, and is the cross section of a grain, namely .
The atoms stick to the surface and hop as random walkers until they either desorb or recombine into molecules. The desorption rates of H and D atoms on the surface are given by
(7) 
where is the attempt rate, is the energy needed for desorption of an atom of isotope I, and is the grain temperature. The hopping rates of the atoms on the surface are
(8) 
where is the activation energy barrier for diffusion of the isotope I. The rate is the inverse of the residence time of an atom of isotope I in a single adsorption site. The sweeping rate is approximately the inverse of the time required for an atom of isotope I to visit nearly all the adsorption sites on the grain surface. Since the D atom is twice as heavy as the H atom, its ground state energy within an adsorption site on the surface is lower. Because of this isotope effect, the desorption barrier for D atoms is assumed to be higher by 5 meV than for H atoms. For the diffusion barriers, we assume that because the diffusion barrier balances the zero point energy of the potential well with that of the saddle point (or transition state), while desorption does not possess such a saddle point. Note that in the analysis of the experimental results, presented in Katz et al. (1999) and Perets et al. (2005), no distinction between H and D atoms was made to ensure that the number of fitting parameters was small.
Table 1: Parameters for the interaction of hydrogen and deuterium atoms with dust grains of different compositions and surface morphologies.
The possibility of the tunneling of adsorbed atoms between adsorption sites was studied by Cazaux & Tielens (2004).
The tunneling probability depends on both the distance between adjacent sites and the energy barrier, which is assumed to be the diffusion barrier.
For the sake of completeness, we present the corresponding tunneling rate,
which is given by (Eisbgerg 1961)
(9) 
where corresponds to the De Broglie wavelength and is given by
(10) 
and E is the kinetic energy of the adsorbed atom, determined by the grain temperature, namely . The efficiency of the tunnelling process is sizeable only when the ratio is of the order 1 or less. There is no specific reason why we should assume that the adsorption sites are uniformly distributed. However, for simplicity we assume that the distance d between adjacent sites is fixed and given by d^{2} = 1/s, as explained above. We have evaluated the tunneling rate for the parameters used in this paper, and it was found to be negligible. Therefore, in the simulation results presented in this paper, the mobility of H and D atoms on grains is only due to thermal hopping, while tunneling is ignored.
In Table 1 we display the parameters that describe the interaction of H and D atoms with grains of different compositions and surface morphologies, as reported in the literature. These parameters are based on a series of experiments and subsequent analysis, reported in Katz et al. (1999), Perets et al. (2005) and Perets et al. (2007). The attempt rate is often assumed to be equal to 10^{12} s^{1}. It may also be derived from an harmonic oscillator model, where .
The reported values of the desorption and diffusion barriers put severe constraints on the grain temperature range over which molecular hydrogen formation may occur. If one considers a distribution of grains of various sizes, such as the power law derived by Mathis et al. (1977), one has to calculate the average of the various quantities over the size distribution, keeping in mind that different temperatures may pertain to grains of different sizes. In the present paper, we assume that the grain temperature depends on the surrounding radiation field but not on the grain size. Based on this assumption, in a given location in the cloud, grains of all sizes have the same temperature. In the simulations reported in this paper, we use the parameters of amorphous carbon. Amorphous carbon is assumed to be a primary component of interstellar dust. Its surface properties are thus suitable for the simulation of molecular hydrogen formation in relatively warm regions in which the grain surfaces are not covered by ice mantles. The parameter values of amorphous carbon are also found to be close to those of low density amorphous ice, which were reported by Perets et al. (2005). Therefore, we use the same set of parameters for the interaction of hydrogen atoms with grain surfaces, independently of whether these grains are expected to be bare or covered by ice mantles.
3 Rate equations, master equation, and moment equations
3.1 Rate equations
After being adsorbed onto the surface, the H and D atoms hop between adsorption sites until they either desorb or form new molecules. The number of atoms of isotope I on the grain surface is denoted by . The rate equations account for the expectation values , of the population size of isotope I on a grain of a given radius, where I = H or D.
3.1.1 Basic processes
In the present version of the PDR code, we introduce three main physical processes to describe the formation of H_{2}, HD, and D_{2} on grain surfaces: adsorption, diffusionmediated reaction, and desorption:
 Adsorption
H + grain H: in s D + grain D: in s ,
where H: and D: represent the hydrogen and deuterium atoms that are adsorbed on the grain surface. The corresponding rate equations are
(11)
(12)
The adsorption coefficients and are directly proportional to the grain cross section and to the sticking probability , namely
(13)
For a distribution of grain sizes, one has to calculate the average of for the distribution function, which leads to the following expression:
(14)
Similar equations hold for deuterium adsorption.
 Reaction on the grain surface due to diffusion
The diffusionmediated formation of molecular hydrogen and its deuterated versions is described by
+ ( in cm^{3}s^{1}) + ( in cm^{3}s^{1}) + ( in cm^{3}s^{1}). Assuming that all formed molecules are directly released into the gas phase, the formation rates of H_{2}, HD, and D_{2} are given by
(15)
(16)
(17)
where ( ) is the number of adsorbed hydrogen (deuterium) atoms on a single grain. The rate coefficients corresponding to these surface reactions can be obtained from the relations and .
 These rate coefficients take the form
= = =
These formulae are valid only if the grains have the same temperature because the diffusion coefficient is temperature dependent.
 Desorption
The desorption processes
H: H + grain ( in s^{1}) D: D + grain ( in s^{1}),
(18) (19)
where
= (20) = (21)
The overall variations in the number of adsorbed H and D atoms per cm^{3} in the rate equation formalism are
3.1.2 Rejection effects
Since the number of adsorption sites on a grain is finite,
one must consider that
hydrogen atoms incident on already occupied sites
may be rejected rather than adsorbed onto the grain.
This mechanism is often referred to as Langmuir rejection.
When the rejection is taken into account,
the flux terms are modified according to
where I stands for H or D. Since the abundance of D atoms is very low ( ) relative to H atoms, it is sensible to consider first, for simplicity, the case in which the rejection is caused only by already adsorbed H atoms.
 Rejection caused only by adsorbed H atoms
The rate equation for adsorbed H on a single grain becomes:
The net result is an apparent increase in the desorption term by . When the incoming flux of deuterium atoms is similarly modified, the rate equation for the adsorbed deuterium becomes
The correction caused by rejection introduces a term proportional to . The corresponding equations for the quantities expressed in volumic density becomes the following
 Rejection by adsorbed H and D atoms
We now include the adsorbed deuterium atoms in the rejection term according to Eq. (24). Since the equations are very similar, we display only the resulting rate equations
The values of the effective desorption coefficients are
=  (33)  
=  (34)  
=  (35) 
This is a generalized form of the equations given in Lipshtat et al. (2004).
3.2 Master and moment equations when only the rejection term for H atoms is included
In the two following sections, we provide the moment and rate equations when rejection is only caused by H atoms.
3.2.1 Master equation
When the number of adsorbed atoms is small,
the rate equations may become inaccurate.
In this case, one may consider using the master equation
or the moment equations derived from it.
The master equation describes the temporal
evolution in the probabilities
that
hydrogen atoms
and
deuterium atoms reside on
the surface of a given grain.
Here we write the master equation
only when the Langmuir rejection due to hydrogen atoms
is taken into account. It takes the form
=  
+  
+  
+  
+  
+  
+  (36) 
The set of equations is written for the various integer values of and . These equations are almost identical to those derived previously by Lipshtat et al. (2004), except that we have explicitly introduced the possibility of the Langmuir rejections viathe terms. Using the relation , we rewrite the master equation in the following way:
=  
(37) 
This set of equations is identical to the standard master equation system (Barzel & Biham 2007b) with two additional terms proportional to and . By a suitable summation over the probabilities in the master equation, one can obtain the moments of the distribution of adsorbed species populations, defined by
(38) 
The order of each moment is defined by the sum l+k. The firstorder moments and represent the mean number of adsorbed H and D atoms, respectively.
3.2.2 Moment equations
Lipshtat & Biham (2003) demonstrated that the moment equation
formalism can adequately describe the evolution of the
system and allows a significant reduction in the number
of coupled equations that need to be solved.
We apply the same technique as in
Barzel & Biham (2007b) to derive the corresponding
moment equations when the additional terms introduced
by the Langmuir rejection are included.
The time derivatives of the moments introduced
by the term proportional to
are
=  
=  
=  
=  
=  (39) 
where the dependent terms are given by,
=  
+  (40) 
Adding these terms to the standard moment equations, one obtains the following system
=  
=  
=  
  
=  
  
=  (41) 
This is the set of equations that is actually used in the interstellar chemistry model when only rejection by H atoms is included. This system of coupled differential equations accounts for the populations of H and D atoms on a grain (given by the two first order moments) and for the reaction rates (given by the three second order moments). One then obtains the formation rates , , and , of H, HD, and D, respectively, on a single grain (Lipshtat et al. 2004)
The formation rate of H, HD, and Don grains is also obtained by summing over the number of grains.
3.3 Master and moment equations with rejection terms for both H and D atoms
When the Langmuir rejection terms caused by both H and D atoms adsorbed onto the surface are included, the master equation takes the form
=  
+  
+  
+  
+  
+  
+  (43) 
Additional terms proportional to and are introduced into the master equations to account for the Langmuir rejection caused by D atoms. The derivation of the moment equations is similar to the previously described procedure, leading to the following set of coupled equations for the first and second moments:
This is the set of equations that is used in the interstellar chemistry model when rejection due to both H and D atoms is included. The production rates of H, HD, and Don a single grain are given by the same equations as those reported previously (Eq. (42)). Integrating over the grain size distribution, we obtain the formation rates (in cm^{3}s^{1}), as well as the adsorption and desorption rates, and the number of adsorbed hydrogen and deuterium atoms (in cm^{3}) as described in the next section.
4 Calculations and results
4.1 Fixed dust temperatures
Table 2: The parameters used in the model. Numbers in parentheses refer to powers of ten. is the total column density of protons (N(H) + 2 N(H_{2})) and is the visual extinction.
Table 3: The three models simulated using the rate equations and the moment equations and their legends.
We first consider a diffuse cloud model whose parameters are given in Table 2. The chemical network comprises of 134 chemical species containing H, D, C, N, O, S and a typical metal ``Fe'' undergoing charge exchange reactions only in addition to photoionisation and electronic recombination.
We compare the results obtained by means of several approximations within the rate and moment equations on the formation of H, HD, and Don amorphous carbon grains (see Table 1 for their properties) for different fixed relevant dust temperatures (between 8 and 20 K) and their link to the gas phase chemistry. We use the Meudon PDR code (Le Petit et al. 2006,2002) in which we have introduced the appropriate changes. The model corresponds to a slab of gas of 6.06 pc (corresponding to a total visual magnitude of 1) irradiated from both sides by the standard interstellar radiation field as given by Draine (1978).
4.1.1 Incorporation of the moment equations into the Meudon PDR code
We couple the full gas phase network and the surface chemistry of H, HD, and Din the PDR code by solving explicitly the steady state of Eq. (44), describing the evolution in the first and second moments of the distribution of hydrogen and deuterium atoms adsorbed on a single grain. The values of the first order moment infer the average numbers of H and D atoms adsorbed on a single grain, whereas no obvious physical meaning is associated with the second order moments. We perform the integration on the size distribution to derive the appropriate quantities such as the number densities of adsorbed H, D, and gas phase atomic and molecular number densities. We first checked the relevance of the treatment by considering only H formation by comparing with the analytical solution obtained from the resolution of the two coupled equations (Lipshtat & Biham 2003).
4.1.2 Formation of H and HD on grains
By using both the moment equations and the rate equations, we examined the assumptions described above. More specifically, we compared three models: a model that includes no Langmuir rejection, a model that includes rejection only caused by adsorbed H atoms, and a model that includes rejection due to both H and D atoms adsorbed onto grains. The three models are listed in Table 3.
We present the results in the form of the standard
rate law for the production of molecular hydrogen
in interstellar clouds. In this rate law, the production rate
(cm^{3} s^{1})
is expressed by
where (cm^{3}) is the total density of H nuclei, in atomic and molecular forms. It can be approximated by . The parameter (cm^{3} s^{1}) is the effective rate coefficient. In the case of diffuse clouds, this rate coefficient is often crudely approximated by . The formation rate of HD molecules can be expressed in a similar way, by
The parameter (cm^{3} s^{1}) is the effective rate coefficient for HD formation on dust grains.
In Fig. 1, we present the rate coefficient
for the formation of H_{2} molecules on grains,
and in Fig. 2 we present the rate coefficient
for the formation of HD molecules.
These rate coefficients are plotted
versus the assumed fixed temperature of the grains
for the three models specified above,
simulated by both the rate equations and
the moment equations.
The conditions used in the simulations
are those obtained at the edge of the modelled cloud,
i.e. A
.
When the Langmuir rejection effect is taken into account,
the formation of molecular hydrogen on grains is efficient
within a narrow window of the grain temperatures,
in agreement with previous studies
(Lipshtat et al. 2004; Katz et al. 1999).
At higher temperatures, the hydrogen atoms do not
remain long enough on the grain to encounter each
other and form molecules.
At lower temperatures, the grain surface is saturated by
nearly immobile hydrogen atoms, and the reaction rate
is very low.
These adsorbed atoms reject new atoms that are incident on
the surface and prevent the formation of additional
layers of hydrogen atoms on the surface.
When the Lan
gmuir rejection term is removed from the
equations, the range of high efficiency recombination is
extended towards lower temperatures.
However, in this case the model enables the accumulation
of many layers of adsorbed hydrogen atoms on the surface.
This accumulation is unphysical. Furthermore,
under these circumstances, the formation of molecular
hydrogen is not mediated by diffusion, and thus the
model itself is not valid.
Laboratory experiments provide strong evidence of
Langmuir rejection
(Katz et al. 1999).
Therefore, the standard forms of both the rate equations and
of moment equations that we incorporate in the PDR
code are those that include the Langmuir rejection term.
In the presence of this term, we observe that
the high efficiency window is somewhat larger for HD
formation than for H_{2} formation.
This is caused by the somewhat higher desorption barrier for
D atoms than for H atoms. The maximum value of the formation rate of H shown in Fig. 1 is in excellent agreement with
the effective formation rate
(cm^{3}s^{1}), expressed by (Le Bourlot et al. 1995) to be
(46) 
for the chosen parameters. The corresponding value for HD is (cm^{3}s^{1}), as seen in Fig. 2. The rate and moment equations are almost identical for the maximum value when desorption effects are insignificant. However, significant discrepancies appear at the edges. Rejection effects play a significant role at the lowest temperatures of the grains, which may not be very relevant for astrophysical purposes. We display in Table 4 the resulting column densities of gas phase H, H, D, and HD when integration of the abundance densities is performed over the width of the cloud.
Figure 1: The formation rate (H_{2}) (cm^{3}s^{1}) of H molecules on grains, obtained from the rate equations (solid line) and from the moment equations (dashed lines), where the grain sizes follow the MRN size distribution. Results obtained from the three models of Table 3 are presented: Model A which includes no rejection (open circles), model B which includes rejection only due to adsorbed H atoms (open squares), and model C that includes rejection caused by both adsorbed H and adsorbed D (open diamonds). 

Open with DEXTER 
Figure 2: The formation rate (HD) (cm^{3}s^{1}) of HD molecules on grains, obtained from the rate equations and the moment equations, where the grain sizes follow the MRN size distribution. Same conventions as in Fig. 1. 

Open with DEXTER 
The column densities of H exhibit significant variations with the assumed dust temperature as reflected by the values of the molecular fraction defined by . The resulting column densities of HD become less sensitive as soon as some molecular hydrogen is present. This is because in presence of H_{2}, the formation of HD occurs mainly via gas phase processes. The formation of HD on grains is more efficient than in the gas phase only at the edge where no H has yet formed. We discuss this in more detail in the next section where we introduce a variable dust temperature as a function of the visual magnitude.
4.2 dependent dust temperatures
4.2.1 Homogeneous temperature distribution
It is well recognised that the grain temperatures decrease when the strength of the incident radiation field decreases for increasing visual magnitude. We introduce this dependence by using the simple analytic formula presented by Hollenbach et al. (1991), which infers the grain temperature as a function of the radiation field strength and the visual magnitude. As our model involves a slab of gas, irradiated from both sides, we extend the formula of Hollenbach et al. (1991) by introducing for the dependence, where and are the strengths of the radiation field on both sides, expressed in Draine units, and is the total visual magnitude of the cloud. In this approximation, grains of all sizes exhibit the same average temperature. We consider possible fluctuations effects in Sect. 4.2.2.
Table 4: Column densities (in cm^{2}) of H, H_{2}, D and HD and the molecular fraction f for model C, obtained from the rate equations (RC) and from the moment equations (MC) for different grain temperatures. Numbers in parentheses refer to powers of ten.
To study a significant range of dust temperatures, we consider a dense photodissociation region (PDR) with a radiation field of intensity of 100 impinging on the left side ( ), and a standard radiation field ( ) impinging on the right side. The parameters of the model are displayed in Table 5. As we want to discuss the different effects of the rate versus moment equations systems, we keep all the physical parameters identical.
We display in Figs. 3 and 4 the dust temperature variation and the formation rates through gas phase and grain surface reactions of H and HD obtained with the RC and MC treatments, because these are expected to be the most physical within the rate equation and moment equation formalisms. The dust temperature profile spans a range of values between about 27 K on the left side, reaching a minimum of 10 K in the shielded region at in the range 57 , and reaching a value of about 10 K on the right edge. This range corresponds to very low to high efficiencies of formation of H2 and HD via grain surface reactions. As far as H formation is concerned, we see that the difference between rate and moment equation formalisms is significant at the left hand side of the cloud, where the dust temperature reaches a value of 27 K, corresponding to a situation where H atoms desorb efficiently from the surface. The results obtained via the rate equations are higher by a factor of about 2. In this particular case, this leads to competition between gas phase and dust processes. However, as soon as the dust temperature reaches a value of about 25, the formation efficiency on dust becomes preponderant. Both treatments are equivalent for a range of temperature between 20 and 11K as already shown in Fig. 1. Significant discrepancies occur again for the range of visual magnitude 58 where the dust temperature reaches a value of 10 K.
These behaviours have a direct impact on the HD formation mechanisms, as displayed in Fig. 4. As long as no significant formation of H occurs, HD is formed mainly via dust processes. However, as soon as H is present, formation of HD in gas phase is much more efficient. In the shielded region where reaches a value of a few and the dust temperature is about 10 K, significant variations occur due to rejection effects. Table 6 displays the resulting column densities. Whereas significant differences exist for the column densities of H and H, the values relevant to the other species are quite similar.
4.2.2 Temperature fluctuation effects
Table 5: Physical conditions relevant to a dense PDR.
Figure 3: Comparison between gas phase formation (circles) and grainsurface formation (squares) rates of H_{2} (in cm^{3}s^{1}), obtained from the rate equations (solid lines and full symbols) and from the moment equations (dottedlines and open symbols). The thin solid line shows the grain temperature (the temperature scale is on the right hand side). 

Open with DEXTER 
Figure 4: Comparison between gas phase formation (circles) and grainsurface formation (squares) rates of HD (in cm^{3}s^{1}) obtained from the rate equations (solid lines and full symbols) and from the moment equations (dotted lines and open symbols). The thin solid line shows the grain temperature (the temperature scale is on the right hand side). 

Open with DEXTER 
The consideration of a homogeneous distribution of dust temperatures as a function of A is reasonable for low incident radiation fields and grains of radii larger than cm (Horn et al. 2007). For grains of radii smaller than cm, the absorption of a single UV photon may cause a temperature spike that heats up the grains by over 20 K (Horn et al. 2007). This is likely to cause the immediate desorption of all the H and D atoms adsorbed on the grain surface.
This conclusion is also derived from detailed kinetic Monte Carlo simulations of H_{2} formation on stochastically heated olivine grains (Cuppen et al. 2006; Herbst & Cuppen 2006). To quantify the effect of temperature fluctuations on small grains, we ran models similar to those presented previously, except that we cut the size distribution of grains at the minimum value of cm. In Table 7, we present the corresponding column densities. In this way, we do not take into account production of any molecular hydrogen on grains of radii smaller than cm both in rate equations and systems of moments approach.
In these models, computed column densities of H_{2} and HD are smaller than those obtained with the MRN distribution extended to
cm (cf. Table 7) because of a smaller integrated grain crosssection value. We find that rate and moment equations infer significantly different values of the total column densities of H_{2}
even for these larger grains. These differences arise essentially from the regions when
is higher than about 18 K, i.e., at the edges of the cloud.
The present results infer an upper limit to the effect of temperature fluctuations because these can only produce a decrease in molecular formation at the surface of grains. In addition, these effects will not occur at
values higher than 1 and lower than 9, because the radiation field is then considerably reduced.
5 Summary
We have incorporated the moment equations for the formation of H_{2} and its deuterated versions into the Meudon PDR code, examined their applicability and relevance, and compared our results with those obtained from rate equations incorporated into the same code under identical conditions. Previous analysis has shown that as long as the temperatures of the dust grains are not too high, the reaction rates obtained from the moment equations coincide with those obtained from the rate equations. At grain temperatures of around 18 K or higher, on the hightemperature end of the efficiency window there are significant differences, where the rate equations over estimate the formation rates of H_{2} and HD. These deviations are caused mainly by the very small grains, on which the population sizes of adsorbed H and D atoms exhibit large fluctuations. In these conditions, it is important to use the moment equations rather than the rate equations. For low grain temperatures, below 12 K or so, we find that the Langmuir rejection term makes a crucial difference. In this range, the rejection term is required to prevent the freezing of multilayers of hydrogen atoms onto the grains. This freezing is unphysical. Here, for the first time, we have derived a set of moment equations in which the Langmuir rejection term is incorporated. Comparison was made with rate equations, in which a suitable rejection term was also incorporated. For very low temperatures, where the grains are covered by a layer of H and D atoms and the rejection term plays a major role, we also find a discrepancy between the formation rates of H_{2} and HD obtained from the rate equations and the moment equations. In these conditions, deep inside a molecular cloud, the formation rate and dissociation rate of H_{2}molecules are low, while HD molecules form mainly by gas phase reactions. These are conditions in which both the rate equations and the moment equations are accurate. The observed discrepancy is caused by the different physical consitions that emerge as a result of the discrepancies in warmer regions, where the moment equations are accurate and the rate equations are not. Therefore, the moment equation results are those that one can rely on in the cold regions also.
Table 6: Column densities in cm^{2} and molecular fraction f for the dense PDR models. Numbers in parentheses refer to power of ten.
Table 7: Column densities in cm^{2} and molecular fraction f for the dense PDR models in which surface chemistry is cut for very small grains ( cm). Numbers in parentheses refer to power of ten.
We conclude that the moment equations provide an efficient method for determining molecular hydrogen formation in interstellar chemistry codes. The equations provide accurate results for the reaction rates on both large and small grains. The equations are easy to construct and are efficient in terms of computational resources. They can be easily coupled to the rate equations of gasphase chemistry. Since the moment equations are linear, their stability and convergence properties are often superior even in compared to the rate equations.
The present study has been achieved by considering that molecular hydrogen can be formed only by means of diffusion of adsorbed atoms, within the strong constraints of present experimental knowledge. Other formation mechanisms may be at work, such as those suggested by Cazaux & Tielens (2004) and Habart et al. (2004) in their interpretation of observations of warm molecular hydrogen.
The method can now be extended to more complex reaction networks on grains. The network that involves H, O, and CO molecules, from which ice mantles containing H_{2}O, CO_{2} and CH_{3}OH are formed by successive hydrogenation and oxydation reactions, is promising (Barzel & Biham 2007b). Unlike the H and D atoms, these heavier atoms and molecules are more strongly bound to both the surface and each other and do not exhibit the Langmuir rejection property. Therefore, in the construction of moment equations for more complex networks, one only needs to maintain the rejection terms for H and D, presented above and there is no need to incorporate these terms for any other atomic or molecular species. Recent studies have shown that CH_{3}OH molecules do not form in the gas phase as efficiently as previously expected. The approach presented here will enable us to evaluate the formation rate of CH_{3}OH on dust grains in various physical conditions, taking into account the contributions of grains of different sizes.
Acknowledgements
This work was supported in part by the FranceIsrael High Council for Science and Technology Research via contract 07 AST F. Evelyne Roueff, J. Le Bourlot and F. Le Petit thank the french national program Physics and Chemistry of Interstellar Medium (PCMI). This work is related the ANR project SCHISM.
References
 Amiaud, L., Dulieu, F., Fillion, J.H., Momeni, A., & Lemaire, J. L. 2007, J. Chem. Phys., 127, 144709 [NASA ADS] [CrossRef]
 Barzel, B., & Biham, O. 2007a, ApJ, 658, L37 [NASA ADS] [CrossRef]
 Barzel, B., & Biham, O. 2007b, J. Chem. Phys., 127, 144703 [NASA ADS] [CrossRef]
 Biham, O., & Lipshtat, A. 2002, Phys. Rev. E, 66, 056103 [NASA ADS] [CrossRef]
 Biham, O., Furman, I., Pirronello, V., et al. 2001, ApJ, 553, 595 [NASA ADS] [CrossRef]
 Caselli, P., Hasegawa, T. I., & Herbst, E. 1998, ApJ, 495, 309 [NASA ADS] [CrossRef]
 Caselli, P., Stantcheva, T., Shalabiea, O., Shematovich, V. I., & Herbst, E. 2002, Planet. Space Sci., 50, 1257 [NASA ADS] [CrossRef] (In the text)
 Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222 [NASA ADS] [CrossRef] (In the text)
 Charnley, S. B. 2001, ApJ, 562, L99 [NASA ADS] [CrossRef] (In the text)
 Charnley, S. B., Tielens, A. G. G. M., & Rodgers, S. D. 1997, ApJ, 482, L203 [NASA ADS] [CrossRef]
 Cuppen, H. M., & Herbst, E. 2005, MNRAS, 361, 565 [NASA ADS] [CrossRef] (In the text)
 Cuppen, H. M., Morata, O., & Herbst, E. 2006, MNRAS, 367, 1757 [NASA ADS] [CrossRef]
 Dalgarno, A., Black, J. H., & Weisheit, J. C. 1973, Astrophys. Lett., 14, 77 [NASA ADS] (In the text)
 Draine, B. T. 1978, ApJS, 36, 595 [NASA ADS] [CrossRef] (In the text)
 Eisbgerg, R. 1961, Fundamentals of Modern Physics (New York: John Wiley & sons) (In the text)
 Garrod, R. T. 2008, A&A, 491, 239 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Green, N. J. B., Toniazzo, T., Pilling, M. J., et al. 2001, A&A, 375, 1111 [NASA ADS] [CrossRef] [EDP Sciences]
 Habart, E., Boulanger, F., Verstraete, L., Walmsley, C. M., & Pineau des Forêts, G. 2004, A&A, 414, 531 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Herbst, E., & Cuppen, H. M. 2006, Proceedings of the National Academy of Science, 103, 12257 [NASA ADS] [CrossRef]
 Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155 [NASA ADS] [CrossRef] (In the text)
 Hollenbach, D. J., Takahashi, T., & Tielens, A. G. G. M. 1991, ApJ, 377, 192 [NASA ADS] [CrossRef] (In the text)
 Hoogers, G. 1995, Surface Science, 327, 47 [NASA ADS] [CrossRef]
 Horn, K., Perets, H. B., & Biham, O. 2007, ArXiv eprints (In the text)
 Katz, N., Furman, I., Biham, O., et al. 1999, ApJ, 522, 305 [NASA ADS] [CrossRef] (In the text)
 Koehler, B. G., Mak, C. H., Arthur, D. A., Coon, P. A., & George, S. M. 1988, J. Chem. Phys., 89, 1709 [NASA ADS] [CrossRef]
 Lacour, S., André, M. K., Sonnentrucker, P., et al. 2005, A&A, 430, 967 [NASA ADS] [CrossRef] [EDP Sciences]
 Le Bourlot, J., Pineau des Forets, G., Roueff, E., et al. 1995, A&A, 302, 870 [NASA ADS] (In the text)
 Le Petit, F., Nehmé, C., Le Bourlot, J., et al. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef]
 Le Petit, F., Roueff, E., & Le Bourlot, J. 2002, A&A, 390, 369 [NASA ADS] [CrossRef] [EDP Sciences]
 Lipshtat, A., & Biham, O. 2003, A&A, 400, 585 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Lipshtat, A., Biham, O., & Herbst, E. 2004, MNRAS, 348, 1055 [NASA ADS] [CrossRef] (In the text)
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] (In the text)
 Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906 [NASA ADS] [CrossRef]
 Perets, H. B., Biham, O., Manicó, G., et al. 2005, ApJ, 627, 850 [NASA ADS] [CrossRef] (In the text)
 Perets, H. B., Lederhendler, A., Biham, O., et al. 2007, ApJ, 661, L163 [NASA ADS] [CrossRef] (In the text)
 Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41 [NASA ADS] [CrossRef]
 Roueff, E., Lis, D. C., van der Tak, F. F. S., Gerin, M., & Goldsmith, P. F. 2005, A&A, 438, 585 [NASA ADS] [CrossRef] [EDP Sciences]
 Roueff, E., Tiné, S., Coudert, L. H., et al. 2000, A&A, 354, L63 [NASA ADS]
 Shalabiea, O. M., Caselli, P., & Herbst, E. 1998, ApJ, 502, 652 [NASA ADS] [CrossRef]
 Stantcheva, T., Caselli, P., & Herbst, E. 2001, A&A, 375, 673 [NASA ADS] [CrossRef] [EDP Sciences]
 Stantcheva, T., & Herbst, E. 2003, MNRAS, 340, 983 [NASA ADS] [CrossRef]
 Stantcheva, T., Shematovich, V. I., & Herbst, E. 2002, A&A, 391, 1069 [NASA ADS] [CrossRef] [EDP Sciences]
 Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 [NASA ADS]
 Watson, W. D. 1974, ApJ, 188, 35 [NASA ADS] [CrossRef] (In the text)
All Tables
Table 1: Parameters for the interaction of hydrogen and deuterium atoms with dust grains of different compositions and surface morphologies.
Table 2: The parameters used in the model. Numbers in parentheses refer to powers of ten. is the total column density of protons (N(H) + 2 N(H_{2})) and is the visual extinction.
Table 3: The three models simulated using the rate equations and the moment equations and their legends.
Table 4: Column densities (in cm^{2}) of H, H_{2}, D and HD and the molecular fraction f for model C, obtained from the rate equations (RC) and from the moment equations (MC) for different grain temperatures. Numbers in parentheses refer to powers of ten.
Table 5: Physical conditions relevant to a dense PDR.
Table 6: Column densities in cm^{2} and molecular fraction f for the dense PDR models. Numbers in parentheses refer to power of ten.
Table 7: Column densities in cm^{2} and molecular fraction f for the dense PDR models in which surface chemistry is cut for very small grains ( cm). Numbers in parentheses refer to power of ten.
All Figures
Figure 1: The formation rate (H_{2}) (cm^{3}s^{1}) of H molecules on grains, obtained from the rate equations (solid line) and from the moment equations (dashed lines), where the grain sizes follow the MRN size distribution. Results obtained from the three models of Table 3 are presented: Model A which includes no rejection (open circles), model B which includes rejection only due to adsorbed H atoms (open squares), and model C that includes rejection caused by both adsorbed H and adsorbed D (open diamonds). 

Open with DEXTER  
In the text 
Figure 2: The formation rate (HD) (cm^{3}s^{1}) of HD molecules on grains, obtained from the rate equations and the moment equations, where the grain sizes follow the MRN size distribution. Same conventions as in Fig. 1. 

Open with DEXTER  
In the text 
Figure 3: Comparison between gas phase formation (circles) and grainsurface formation (squares) rates of H_{2} (in cm^{3}s^{1}), obtained from the rate equations (solid lines and full symbols) and from the moment equations (dottedlines and open symbols). The thin solid line shows the grain temperature (the temperature scale is on the right hand side). 

Open with DEXTER  
In the text 
Figure 4: Comparison between gas phase formation (circles) and grainsurface formation (squares) rates of HD (in cm^{3}s^{1}) obtained from the rate equations (solid lines and full symbols) and from the moment equations (dotted lines and open symbols). The thin solid line shows the grain temperature (the temperature scale is on the right hand side). 

Open with DEXTER  
In the text 
Copyright ESO 2009