EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 505 / No 3 (October III 2009) A&A, 505 3 (2009) 1153-1165 Full HTML
Free Access
Issue
A&A
Volume 505, Number 3, October III 2009
Page(s) 1153 - 1165
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/200912269
Published online 18 August 2009

Incorporation of stochastic chemistry on dust grains in the Meudon PDR code using moment equations

I. Application to the formation of H$_{\sf 2}$ and HD

F. Le Petit1 - B. Barzel2 - O. Biham2 - E. Roueff1 - J. Le Bourlot1

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 gas-phase 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$_{\rm 2}$, HD, and D$_{\rm 2}$molecules 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$_{\rm 2}$, HD and D$_{\rm 2}$molecules 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 gas-phase 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 time-dependent approach and the steady-state approach. In the time-dependent 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 grain-surface 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 gas-phase 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 grain-surface 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 semi-empirical 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 grain-surface 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 gas-phase chemistry. Unlike the rate equations, the moment equations are linear, and thus easier to handle both in the time-dependent 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 steady-state 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 plane-parallel 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 H2. Explicit treatment is performed for C and S photoionization, H2 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$_{\rm 2}$, HD, and D2. 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 gas-phase and surface reactions involved in the formation of H$_{\rm 2}$, HD and D$_{\rm 2}$molecules. We present in Sect. 3 the moment equations for the H-D 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 gas-phase 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$_{\rm 2}$. 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

\begin{displaymath}c = \frac{3}{4\pi} \frac{1.4 \cdot m_{\rm H} \cdot G \cdot (4... ... \left(a_{\rm max}^{4-q} - a_{\rm min}^{4-q}\right)}n_{\rm H}, \end{displaymath} (2)

where $\rho$ (g cm-3) stands for the volumic mass of the grains, while $a_{\rm max}$ and $a_{\rm min}$ are the upper and lower cutoffs of the grain radii distribution, respectively. We denote the proton density in the gas phase by $n_{\rm H}$ (cm-3) and the dust to gas mass ratio by G=0.01. In this case, one can formulate the formation rate of H2 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:

\begin{displaymath}{\rm H^{+} + D \rightleftharpoons D^{+} + H.} \end{displaymath} (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$_{\rm 2}$ to form HD:

\begin{displaymath}{\rm D^{+} + H_2 \rightarrow HD + H^{+}.} \end{displaymath} (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 gas-phase scheme accounts for the observed HD/H2 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 H3+to form H2D+, 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, grain-surface reactions may also contribute to an enhanced production of HD and D$_{\rm 2}$molecules. 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 gas-phase 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$_{\rm 2}$, HD, and D$_{\rm 2}$molecules 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 $\upsilon_{\rm H}$ and $\upsilon_{\rm D}$, 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 $n_{\rm gr}$ (cm-3) as a function of $n_{\rm H}$ to be

\begin{displaymath}n_{\rm gr} =1.4 \cdot \frac{3 \cdot m_{\rm H} \cdot G} {4\pi \cdot a^3 \cdot \rho } \cdot n_{\rm H}. \end{displaymath} (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

\begin{displaymath}n_{\rm gr} = 1.4 \cdot \frac{3 \cdot m_{\rm H} \cdot G \cdo... ...t (\sqrt{a_{\rm max}} - \sqrt{a_{\rm min}}) } \cdot n_{\rm H}. \end{displaymath} (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 $s = S / 4 \pi a^2$. The density of adsorption sites, s, and the distance between adjacent sites, d, are related by d2 = 1/s. The fluxes $F_{\rm H}$ and $F_{\rm D}$ (atoms s-1) of H and D atoms respectively onto the surface of a single grain are given by $F_{\rm I} = n_{\rm I} \upsilon_{\rm I} \sigma$, where I = H, D, and $\sigma$ is the cross section of a grain, namely $\sigma = \pi a^2$.

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

\begin{displaymath}W_{\rm I} = \nu_{\rm I} \cdot \exp\left( - \frac{E_{\rm I}^{\rm des}}{k T_{\rm gr}} \right), \end{displaymath} (7)

where $\nu_{\rm I}$ is the attempt rate, $E_{\rm I}^{\rm des}$is the energy needed for desorption of an atom of isotope I, and $T_{\rm gr}$ is the grain temperature. The hopping rates of the atoms on the surface are

\begin{displaymath}a_{\rm I} = \nu_{\rm I} \cdot \exp\left( - \frac{E_{\rm I}^{\rm diff}}{k T_{\rm gr}} \right), \end{displaymath} (8)

where $E_{\rm I}^{\rm diff}$ is the activation energy barrier for diffusion of the isotope I. The rate $a_{\rm I}$ is the inverse of the residence time $\tau_{\rm I}$ of an atom of isotope I in a single adsorption site. The sweeping rate $A_{\rm I} = a_{\rm I }/ S$ is approximately the inverse of the time $S \tau_{\rm I}$ 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 $E_{\rm D}^{\rm diff} = E_{\rm H}^{\rm diff}$ 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)

\begin{displaymath}k_{\rm tun} = {\nu_{\rm I}} \left[ {1+\frac{(E_{\rm I}^{\rm d... ...-{E})} \cdot {\rm sinh^2}(d/ \lambda_{\rm DB})} \right]^{-1}, \end{displaymath} (9)

where $\lambda_{\rm DB}$ corresponds to the De Broglie wavelength and is given by

\begin{displaymath}\lambda_{\rm DB} = \frac{\hbar}{\sqrt{2m(E_{\rm I}^{\rm diff}-E)}}, \end{displaymath} (10)

and E is the kinetic energy of the adsorbed atom, determined by the grain temperature, namely $E = k_{\rm B} T_{\rm gr}$. The efficiency of the tunnelling process is sizeable only when the ratio $d/\lambda_{\rm DB}$ 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 d2 = 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 $\nu_{\rm I}$ is often assumed to be equal to 1012 s-1. It may also be derived from an harmonic oscillator model, where $\nu_{\rm I} = \sqrt{2E^{\rm des}_{\rm I}/m} ~ / ({\pi ~d}) $.

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 $N_{\rm I:}$. The rate equations account for the expectation values $\langle N_{\rm I:} \rangle$, 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 H2, HD, and D2 on grain surfaces: adsorption, diffusion-mediated reaction, and desorption:

  • Adsorption
    H + grain $\rightarrow$ H: $\Big(k^{\rm H}_{\rm ads}$ in s $^{-1}\Big)$
    D + grain $\rightarrow$ D: $\Big(k^{\rm D}_{\rm ads}$ in s $^{-1}\Big)$,


    where H: and D: represent the hydrogen and deuterium atoms that are adsorbed on the grain surface. The corresponding rate equations are

    \begin{displaymath}{\frac{{\rm d}n({\rm H:})}{{\rm d}t} = +k^{\rm H}_{\rm ads} n({\rm H})} \end{displaymath} (11)


    \begin{displaymath}{\frac{{\rm d}n({\rm D:})}{{\rm d}t} = +k^{\rm D}_{\rm ads} n({\rm D})}. \end{displaymath} (12)

    The adsorption coefficients $k^H_{\rm ads}$ and $k^D_{\rm ads}$ are directly proportional to the grain cross section $\sigma = \pi a^2$ and to the sticking probability $\gamma$, namely

    \begin{displaymath}k^{\rm H}_{\rm ads} = \frac{F_{\rm H}}{n({\rm H})} n_{\rm gr} = \gamma \upsilon_{\rm H} \pi a^2 n_{\rm gr}. \end{displaymath} (13)

    For a distribution of grain sizes, one has to calculate the average of $\sigma n_{\rm gr}$ for the distribution function, which leads to the following expression:

    \begin{displaymath}k^{\rm H}_{\rm ads} = \left\langle \frac{F_{\rm H}}{n({\rm H}... ...\pi \cdot \rho \sqrt{a_{\rm min} a_{\max}} } \cdot n_{\rm H} \end{displaymath} (14)

    Similar equations hold for deuterium adsorption.

  • Reaction on the grain surface due to diffusion

    The diffusion-mediated formation of molecular hydrogen and its deuterated versions is described by

    $\rm H:$ + $\rm H:$ $\rightarrow$ $\rm H_2$       ( $k^{\rm H_2}_{\rm surf}$ in cm-3s-1)
    $\rm H:$ + $\rm D: $ $\rightarrow$ $\rm HD$       ( $k^{\rm HD}_{\rm surf}$ in cm-3s-1)
    $\rm D: $ + $\rm D: $ $\rightarrow$ $\rm D_2$       ( $k^{\rm D_2}_{\rm surf}$ in cm-3s-1).

    Assuming that all formed molecules are directly released into the gas phase, the formation rates of H2, HD, and D2 are given by

    \begin{displaymath}{\frac{{\rm d}n({\rm H_2})}{{\rm d}t} = k^{\rm H_2}_{\rm sur... ... H}:)^2} = A_{\rm H} \langle N_{\rm H:} \rangle^2 n_{\rm gr}, \end{displaymath} (15)


    \begin{displaymath}{\frac{{\rm d}n({\rm HD})}{{\rm d}t} = k^{\rm HD}_{\rm surf} ... ...gle N_{\rm H:} \rangle \langle N_{\rm D:} \rangle n_{\rm gr}, \end{displaymath} (16)


    \begin{displaymath}{\frac{{\rm d}n({\rm D_2})}{{\rm d}t} = k^{\rm D_2}_{\rm sur... ... D}:)^2} = A_{\rm D} \langle N_{\rm D:} \rangle^2 n_{\rm gr}, \end{displaymath} (17)

    where $\langle N_{\rm H:} \rangle$ ( $\langle N_{\rm D:} \rangle$) 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 $n({\rm H}:) = \langle N_{\rm H:} \rangle \cdot n_{\rm gr}$ and $n({\rm D}:) = \langle N_{\rm D} \rangle \cdot n_{\rm gr}$.
These rate coefficients take the form
                           $\displaystyle k^{\rm H_2}_{\rm surf}$ = $\displaystyle \frac{A_H \langle N_{\rm H:} \rangle^2}{n({\rm H}:)^2}n_{\rm gr}= \frac{A_{\rm H}}{n_{\rm gr}}$  
$\displaystyle k^{\rm HD}_{\rm surf}$ = $\displaystyle \frac{{(A_{\rm H}+A_{\rm D})} \langle N_{ \rm H:} \rangle \langle... ...gle}{n({\rm H}:)n({\rm D}:)}n_{\rm gr} = \frac{A_{\rm H}+A_{\rm D}}{n_{\rm gr}}$  
$\displaystyle k^{\rm D_2}_{\rm surf}$ = $\displaystyle \frac{A_{\rm D} \langle N_{\rm D:}\rangle^2}{n(\rm D:)^2}n_{\rm gr}= \frac{A_{\rm D}}{n_{\rm gr}}\cdot$  

These formulae are valid only if the grains have the same temperature because the diffusion coefficient is temperature dependent.
  $\mbox{$\bullet$}$
Desorption The desorption processes
H:   $\rightarrow$ H + grain     ( $k^{\rm H}_{\rm des}$ in s-1)
D:   $\rightarrow$ D + grain     ( $k^{\rm D}_{\rm des}$ in s-1),
are described by the equations
$\displaystyle \frac{{\rm d}n(\rm H:)}{{\rm d}t} = -k^{\rm H}_{\rm des}~n(\rm H:)$     (18)
$\displaystyle \frac{{\rm d}n(\rm D:)}{{\rm d}t} = -k^{\rm D}_{\rm des}~n(\rm D:),$     (19)

where
$\displaystyle k^{\rm H}_{\rm des}$ = $\displaystyle W_{\rm H}$ (20)
$\displaystyle k^{\rm D}_{\rm des}$ = $\displaystyle W_{\rm D}.$ (21)

The overall variations in the number of adsorbed H and D atoms per cm3  in the rate equation formalism are
 
$\displaystyle \frac{{\rm d} n({\rm H:}) }{{\rm d}t}$ = $\displaystyle k^{\rm H}_{\rm ads}n({\rm H}) -W_{\rm H} n({\rm H:}) - 2\frac{A_{... ...\rm H:}) ^2 - \frac{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H:}) n({\rm D:}),$ (22)
$\displaystyle \frac{{\rm d} n({\rm D}:)} {{\rm d}t}$ = $\displaystyle k^{\rm D}_{\rm ads}n({\rm D}) -W_{\rm D} n({\rm D}:) - 2\frac{A_{... ...\rm D:}) ^2 - \frac{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H}:) n({\rm D}:),$ (23)

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

\begin{displaymath}F^{\rm eff}_{\rm I}= F_{\rm I}\left(1 - \frac{ N_{\rm H} + N_{\rm D} }{S}\right), \end{displaymath} (24)

where I stands for H or D. Since the abundance of D atoms is very low ( $1.5 \times 10^{-5}$) 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:

    \begin{displaymath} \frac{{\rm d} \langle N_{\rm H:} \rangle}{{\rm d}t} = F_{\r... ...m D}) \langle N_{\rm H:} \rangle \langle N_{\rm D:} \rangle. \end{displaymath} (25)

    The net result is an apparent increase in the desorption term by $f_{\rm H} = {F_{\rm H}}/{S}$. When the incoming flux of deuterium atoms is similarly modified, the rate equation for the adsorbed deuterium becomes

    \begin{displaymath} \frac{{\rm d} \langle N_{\rm D:} \rangle}{{\rm d}t} = F_{\r... ...m D}) \langle N_{\rm H:} \rangle \langle N_{\rm D:} \rangle. \end{displaymath} (26)

    The correction caused by rejection introduces a term proportional to ${\langle N_{\rm H:} \rangle}/{\langle N_{\rm D:} \rangle}$. The corresponding equations for the quantities expressed in volumic density becomes the following

    \begin{displaymath} \frac{{\rm d} n({\rm H:}) }{{\rm d}t} = k^H_{\rm ads}n({\rm... ...c{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H:}) n({\rm D:}), \end{displaymath} (27)


    \begin{displaymath} \frac{{\rm d} n({\rm D:}) }{{\rm d}t} = k^D_{\rm ads}n({\rm ... ...c{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H:}) n({\rm D:}). \end{displaymath} (28)

  • 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

    \begin{displaymath} \frac{{\rm d} n({\rm H:}) }{{\rm d}t} = k^{\rm H}_{\rm ads}... ...c{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H:}) n({\rm D:}), \end{displaymath} (29)


    \begin{displaymath} \frac{{\rm d} n({\rm D:}) }{{\rm d}t} = k^{\rm D}_{\rm ads}... ...c{(A_{\rm H}+A_{\rm D})}{n_{\rm gr}} n({\rm H:}) n({\rm D:}). \end{displaymath} (30)

    The values of the effective desorption coefficients are

    \begin{displaymath} W_{\rm H}^{\rm eff} = {W_{\rm H}+\frac{\gamma \upsilon_{\rm ... ...psilon_{\rm H} n({\rm H})}{4s}\frac{n({\rm D}:)}{n({\rm H}:)}, \end{displaymath} (31)


    \begin{displaymath} W_{\rm D}^{\rm eff} = {W_{\rm D}+\frac{\gamma \upsilon_{\rm ... ...on_{\rm D} n({\rm D})}{4s}\frac{n({\rm H}:)}{n({\rm D}:)}\cdot \end{displaymath} (32)

The production rates of H$_{\rm 2}$, HD, and D$_{\rm 2}$ (in units of cm-3s-1) are
                             $\displaystyle R_{\rm H_2}$ = $\displaystyle {(A_{\rm H }/n_{\rm gr})} n({\rm H:})^2$ (33)
$\displaystyle R_{\rm HD}$ = $\displaystyle {((A_{\rm H }+A_{\rm D})/n_{\rm gr})}n({\rm H:})n({\rm D:})$ (34)
$\displaystyle R_{\rm D_2}$ = $\displaystyle {(A_{\rm D }/n_{\rm gr})} n({\rm D:})^2.$ (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 $P(N_{\rm H:}, N_{\rm D:})$ that $N_{\rm H:}$ hydrogen atoms and $N_{\rm D:}$ 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

                          $\displaystyle \frac{{{\rm d} P}(N_{\rm H:}, N_{\rm D:})}{{\rm d}t}$ = $\displaystyle F_{\rm H}\left[\left(1-\frac{N_{\rm H:}-1}{S}\right) P(N_{\rm H:}... ...\rm D:}) - \left(1-\frac{N_{\rm H:}}{S}\right) P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle F_{\rm D}\left[\left(1-\frac{N_{\rm H:}}{S}\right) P(N_{\rm H:}, ... ...rm D:}-1) - \left(1-\frac{N_{\rm H:}}{S}\right) P(N_{\rm H:}, N_{\rm D:}\right]$  
  + $\displaystyle W_{\rm H} \left[(N_{\rm H:}+1)P(N_{\rm H:}+1, N_{\rm D:}) - N_{\rm H:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle W_{\rm D} \left[(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+1) - N_{\rm D:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle A_{\rm H} \left[(N_{\rm H:}+2)(N_{\rm H:}+1)P(N_{\rm H:}+2, N_{\rm D:}) - N_{\rm H:}(N_{\rm H:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle A_{\rm D} \left[(N_{\rm D:}+2)(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+2) - N_{\rm D:}(N_{\rm D:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle {(A_{\rm H}+A_{\rm D}) \left[(N_{\rm H:}+1)(N_{\rm D:}+1) P(N_{\rm H:}+1, N_{\rm D:}+1) - N_{\rm H:} N_{\rm D:} P(N_{\rm H:}, N_{\rm D:}) \right]}.$ (36)

The set of equations is written for the various integer values of $N_{\rm H:}$ and $N_{\rm D:}$. 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 ${-F \cdot N/S}$ terms. Using the relation $f_{\rm I} = F_{\rm I} /S$, we rewrite the master equation in the following way:
                          $\displaystyle \frac{{{\rm d} P}(N_{\rm H:}, N_{\rm D:})}{{\rm d}t}$ = $\displaystyle F_{\rm H}\left[P(N_{\rm H:}-1, N_{\rm D:})-P(N_{\rm H:}, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle F_{\rm D}\left[(P(N_{\rm H:}, N_{\rm D:}-1) -P(N_{\rm H:}, N_{\rm D:}\right]$  
  $\textstyle \quad +$ $\displaystyle f_{\rm H} \left[N_{\rm H:}P(N_{\rm H:}, N_{\rm D:}) - (N_{\rm H:}-1) P(N_{\rm H:}-1, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle f_{\rm D} \left[N_{\rm H:}P(N_{\rm H:}, N_{\rm D:}) - N_{\rm H:} P(N_{\rm H:}, N_{\rm D:}-1)\right]$  
  $\textstyle \quad +$ $\displaystyle W_{\rm H} \left[(N_{\rm H:}+1)P(N_{\rm H:}+1, N_{\rm D:}) - N_{\rm H:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle W_{\rm D} \left[(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+1) - N_{\rm D:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle A_{\rm H} \left[(N_{\rm H:}+2)(N_{\rm H:}+1)P(N_{\rm H:}+2, N_{\rm D:}) - N_{\rm H:}(N_{\rm H:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle A_{\rm D} \left[(N_{\rm D:}+2)(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+2) - N_{\rm D:}(N_{\rm D:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  $\textstyle \quad +$ $\displaystyle {(A_{\rm H}+A_{\rm D}) \left[(N_{\rm H:}+1)(N_{\rm D:}+1) P(N_{\rm H:}+1, N_{\rm D:}+1) - N_{\rm H:} N_{\rm D:} P(N_{\rm H:}, N_{\rm D:}) \right]}.$ (37)

This set of equations is identical to the standard master equation system (Barzel & Biham 2007b) with two additional terms proportional to $f_{\rm H}$ and $f_{\rm D}$. 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

\begin{displaymath}\langle N_{\rm H:}^kN_{\rm D:}^l \rangle = \sum_{N_{\rm H:},N_{\rm D:}} N_{\rm H:}^kN_{\rm D:}^l P(N_{\rm H:}, N_{\rm D:}). \end{displaymath} (38)

The order of each moment is defined by the sum l+k. The first-order moments $\langle N_{\rm H:} \rangle$ and $\langle N_{\rm D:} \rangle$ 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 $f_{\rm I}$ are

                            $\displaystyle \frac{{\rm d}\langle {N_{\rm H:}} \rangle_{\rm f}}{{\rm d}t}$ = $\displaystyle \sum_{N_H,N_D} N_{\rm H:} \dot{P_{\rm f}}(N_{\rm H:}, N_{\rm D:}) = - f_{\rm H} \langle N_{\rm H:} \rangle$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm D:}} \rangle_{\rm f}}{{\rm d}t}$ = $\displaystyle \sum_{N_H,N_D} N_{\rm D:} \dot{P_{\rm f}}(N_{\rm H:}, N_{\rm D:}) = - f_{\rm D} \langle N_{\rm H:} \rangle$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm H:}}^2 \rangle_{\rm f}}{{\rm d}t}$ = $\displaystyle \sum_{N_H,N_D} N_{\rm H:}^2 \dot{P_{\rm f}}(N_{\rm H:}, N_{\rm D:}) = - f_{\rm H} (\langle N_{\rm H:} \rangle + 2 \langle N_{\rm H:}^2 \rangle )$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm D:}}^2 \rangle_{\rm f}}{{\rm d}t}$ = $\displaystyle \sum_{N_H,N_D} N_{\rm D:}^2 \dot{P_f}(N_{\rm H:}, N_{\rm D:}) = - f_{\rm D} (\langle N_{\rm H:} \rangle + 2 \langle N_{\rm H:} N_{\rm D:} \rangle )$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm H:}N_{\rm D:}} \rangle_{\rm f}}{{\rm d}t}$ = $\displaystyle \sum_{N_H,N_D} N_{\rm H:} N_{\rm D:} \dot{P_f}(N_{\rm H:}, N_{\rm... ...} \langle N_{\rm H:}N_{\rm D:} \rangle -f_{\rm D} \langle N_{\rm H:}^2 \rangle,$ (39)

where the $f_{\rm I}$ dependent terms are given by,
                            $\displaystyle \dot{P_{\rm f}}(N_{\rm H:}, N_{\rm D:})$ = $\displaystyle f_{\rm H} \left[N_{\rm H:}P(N_{\rm H:}, N_{\rm D:}) - (N_{\rm H:}-1) P(N_{\rm H:}-1, N_{\rm D:})\right]$  
  + $\displaystyle f_{\rm D} \left[ N_{\rm H:}P(N_{\rm H:},N_{\rm D:}) -N_{\rm H:}P(N_{\rm H:},N_{\rm D:}-1) \right].$ (40)

Adding these terms to the standard moment equations, one obtains the following system
                            $\displaystyle \frac{{\rm d} \langle {N_{\rm H:}} \rangle}{{\rm d}t}$ = $\displaystyle F_{\rm H} + (2 A_{\rm H} - W_{\rm H}-f_{\rm H}) \langle N_{\rm H:... ...m H:} ^2 \rangle - (A_{\rm H}+ A_{\rm D}) \langle N_{\rm H:}N_{\rm D:} \rangle,$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm D:}} \rangle}{{\rm d}t}$ = $\displaystyle F_{\rm D} -f_{\rm D} \langle N_{\rm H:} \rangle + (2 A_{\rm D}- W... ...rm D:}^2 \rangle - (A_{\rm H}+ A_{\rm D}) \langle N_{\rm H:}N_{\rm D:} \rangle,$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm H:}}^2 \rangle}{{\rm d}t}$ = $\displaystyle F_{\rm H} + (2 F_{\rm H} + W_{\rm H} + 4 A_{\rm H} -f_{\rm H}) \l... ...\rangle - (2 W_{\rm H} + 4 A_{\rm H}+2f_{\rm H} ) \langle N_{\rm H:} ^2 \rangle$  
  - $\displaystyle (A_{\rm H} + A_{\rm D}) \langle N_{\rm H:}N_{\rm D:} \rangle,$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm D:}}^2 \rangle}{{\rm d}t}$ = $\displaystyle F_{\rm D} - f_{\rm D} \langle N_{\rm H:} \rangle + (2 F_{\rm D} +... ...e N_{\rm D:} \rangle - (2 W_{\rm D} + 4 A_{\rm D}) \langle N_{\rm D:}^2 \rangle$  
  - $\displaystyle (A_{\rm H} + A_{\rm D}+2f_{\rm D}) \langle N_{\rm H:}N_{\rm D:} \rangle,$  
$\displaystyle \frac{{\rm d} \langle {N_{\rm H:}N_{\rm D:}} \rangle}{{\rm d}t}$ = $\displaystyle F_{\rm D} \langle N_{\rm H:} \rangle + F_{\rm H} \langle N_{\rm D... ...rm D}+ A_{\rm H} + A_{\rm D} + f_{\rm H}) \langle N_{\rm H:}N_{\rm D:} \rangle.$ (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 $r_{\rm H_2}$, $r_{\rm HD}$, and $r_{\rm D_2}$, of H$_{\rm 2}$, HD, and D$_{\rm 2}$, respectively, on a single grain (Lipshtat et al. 2004)
                   $\displaystyle r_{\rm H_2}$ = $\displaystyle A_{\rm H} \left (\langle N_{\rm H}^2 \rangle - \langle N_{\rm H} \rangle \right)$  
$\displaystyle r_{\rm HD}$ = $\displaystyle (A_{\rm H}+A_{\rm D} ) \langle N_{\rm H}N_{\rm D} \rangle$  
$\displaystyle r_{\rm D_2}$ = $\displaystyle A_{\rm D} \left ( \langle N_{\rm D} ^2 \rangle - \langle N_{\rm D} \rangle \right).$ (42)

The formation rate of H$_{\rm 2}$, HD, and D$_{\rm 2}$on 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

                            $\displaystyle \frac{{{\rm d} P}(N_{\rm H:}, N_{\rm D:})}{{\rm d}t}$ = $\displaystyle F_{\rm H}\left[ \left(1-\frac{N_{\rm H:}-1}{S} -\frac{N_{\rm D:}}... ...rac{N_{\rm H:}}{S}-\frac{N_{\rm D:}}{S}\right) P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle F_{\rm D}\left[\left(1-\frac{N_{\rm H:}}{S} -\frac{N_{\rm D:}-1}{... ...rac{N_{\rm H:}}{S} -\frac{N_{\rm D:}}{S}\right)P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle W_{\rm H} \left[(N_{\rm H:}+1)P(N_{\rm H:}+1, N_{\rm D:}) - N_{\rm H:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle W_{\rm D} \left[(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+1) - N_{\rm D:}P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle A_{\rm H} \left[(N_{\rm H:}+2)(N_{\rm H:}+1)P(N_{\rm H:}+2, N_{\rm D:}) - N_{\rm H:}(N_{\rm H:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle A_{\rm D} \left[(N_{\rm D:}+2)(N_{\rm D:}+1)P(N_{\rm H:}, N_{\rm D:}+2) - N_{\rm D:}(N_{\rm D:}-1)P(N_{\rm H:}, N_{\rm D:})\right]$  
  + $\displaystyle {(A_{\rm H}+A_{\rm D}) \left[(N_{\rm H:}+1)(N_{\rm D:}+1) P(N_{\rm H:}+1, N_{\rm D:}+1) - N_{\rm H:} N_{\rm D:} P(N_{\rm H:}, N_{\rm D:}) \right]}.$ (43)

Additional terms proportional to $f_{\rm H}$ and $f_{\rm D}$ 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:

 \begin{eqnarray*}\frac{{\rm d} \langle {N_{\rm H:}}\rangle}{{\rm d}t} &=& F_{\r... ...D} + f_{\rm H}+f_{\rm D}) \langle N_{\rm H:}N_{\rm D:} \rangle. \end{eqnarray*}


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$_{\rm 2}$, HD, and D$_{\rm 2}$on 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-3s-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. $N_{\rm H}$ is the total column density of protons (N(H) + 2 N(H2)) and $A_{\rm v}$ 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$_{\rm 2}$, HD, and D$_{\rm 2}$on 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$_{\rm 2}$, HD, and D$_{\rm 2}$in 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$_{\rm 2}$ 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$_{\rm 2}$ 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 $R_{\rm H_2}$ (cm-3 s-1) is expressed by

\begin{displaymath}R_{\rm H_2} = \alpha({\rm H_2}) n({\rm H}) n_{\rm H}, \end{displaymath} (44)

where $n_{\rm H}$ (cm-3) is the total density of H nuclei, in atomic and molecular forms. It can be approximated by $n_{\rm H}= n({\rm H}) + 2 n({\rm H_2})$. The parameter $\alpha({\rm H_2})$(cm3 s-1) is the effective rate coefficient. In the case of diffuse clouds, this rate coefficient is often crudely approximated by $\alpha({\rm H_2}) = 5 \times 10^{-17} \sqrt{T_{\rm gas}/300}$. The formation rate of HD molecules can be expressed in a similar way, by

\begin{displaymath}R_{\rm HD} = \alpha({\rm HD}) n({\rm D}) n_{\rm H}. \end{displaymath} (45)

The parameter $\alpha({\rm HD})$ (cm3 s-1) is the effective rate coefficient for HD formation on dust grains.

In Fig. 1, we present the rate coefficient $\alpha({\rm H_2})$ for the formation of H2 molecules on grains, and in Fig. 2 we present the rate coefficient $\alpha({\rm HD})$ 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 $_{\rm V} = 0$. 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 H2 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$_{\rm 2}$ shown in Fig. 1 is in excellent agreement with the effective formation rate $\alpha_{H_2}$ (cm3s-1), expressed by (Le Bourlot et al. 1995) to be

\begin{displaymath}\alpha_{H_2} = 1.4 \cdot \frac{3 \cdot m_{\rm H} \cdot G } {4... ...a_{\max}}} \upsilon_{\rm H} = 1.62 \times 10^{-16} n_{\rm H}, \end{displaymath} (46)

for the chosen parameters. The corresponding value for HD is $\alpha_{\rm HD} = \sqrt{2} \alpha_{\rm H_2} = 2.3 \times 10^{-16} n_{\rm H}$ (cm3s-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$_{\rm 2}$, D, and HD when integration of the abundance densities is performed over the width of the cloud.
 \begin{figure} \par\includegraphics[width=8.8cm,clip]{12269fg1.eps} \end{figure} Figure 1:

The formation rate $\alpha $(H2) (cm3s-1) of H$_{\rm 2}$ 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
 \begin{figure} \par\includegraphics[width=8.8cm]{12269fg2.eps} \end{figure} Figure 2:

The formation rate $\alpha $(HD) (cm3s-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
We display only the results of model C, which are the most accurate when rejection effects are included both for H and D impinging atoms.

The column densities of H$_{\rm 2}$ exhibit significant variations with the assumed dust temperature as reflected by the values of the molecular fraction defined by $f = {{2 N({\rm H}_2)}/{[N({\rm H}) + 2 N({\rm H}_2)]}}$. The resulting column densities of HD become less sensitive as soon as some molecular hydrogen is present. This is because in presence of H2, 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$_{\rm 2}$ 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 $A_{\rm v}$ 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 $[\chi_{\rm left} \exp(-1.8 A_{\rm v}) + \chi_{\rm right} \exp(-1.8(A_{\rm v}^{\rm tot}-A_{\rm v})]$ for the $A_{\rm v}$ dependence, where $\chi_{\rm left}$ and $\chi_{\rm right}$ are the strengths of the radiation field on both sides, expressed in Draine units, and $A_{\rm v}^{\rm tot}$ 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, H2, 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 ( $\chi_{\rm left}$), and a standard radiation field ( $\chi_{\rm right}$) 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$_{\rm 2}$ 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 $A_{\rm v}$ in the range 5-7 , 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$_{\rm 2}$  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 5-8 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$_{\rm 2}$ occurs, HD is formed mainly via dust processes. However, as soon as H$_{\rm 2}$  is present, formation of HD in gas phase is much more efficient. In the shielded region where $A_{\rm v}$ 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$_{\rm 2}$, 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.

 \begin{figure} \par\includegraphics[width=8.8cm]{12269fg3.eps} \end{figure} Figure 3:

Comparison between gas phase formation (circles) and grain-surface formation (squares) rates of H2 (in cm-3s-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

 \begin{figure} \par\includegraphics[width=8.8cm]{12269fg4.eps} \end{figure} Figure 4:

Comparison between gas phase formation (circles) and grain-surface formation (squares) rates of HD (in cm-3s-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$_{\rm v}$ is reasonable for low incident radiation fields and grains of radii larger than $1.0 \times 10^{-6}$ cm (Horn et al. 2007). For grains of radii smaller than $\sim$ $1.0 \times 10^{-6}$ 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 H2 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 $1.0 \times 10^{-6}$ 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 $1.0 \times 10^{-6}$ cm both in rate equations and systems of moments approach.

In these models, computed column densities of H2 and HD are smaller than those obtained with the MRN distribution extended to $3.0 \times 10^{-7}$ cm (cf. Table 7) because of a smaller integrated grain cross-section value. We find that rate and moment equations infer significantly different values of the total column densities of H2 even for these larger grains. These differences arise essentially from the regions when $T_{\rm dust}$ 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 $A_{\rm v}$ values higher than $\sim$1 and lower than $\sim$9, because the radiation field is then considerably reduced.

5 Summary

We have incorporated the moment equations for the formation of H2 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 high-temperature end of the efficiency window there are significant differences, where the rate equations over estimate the formation rates of H2 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 multi-layers 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 H2 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 H2molecules 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 ( $r < 1.0 \times 10^{-6}$ 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 gas-phase 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 H2O, CO2 and CH3OH 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 CH3OH 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 CH3OH 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 France-Israel 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

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. $N_{\rm H}$ is the total column density of protons (N(H) + 2 N(H2)) and $A_{\rm v}$ 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, H2, 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 ( $r < 1.0 \times 10^{-6}$ cm). Numbers in parentheses refer to power of ten.

All Figures

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{12269fg1.eps} \end{figure} Figure 1:

The formation rate $\alpha $(H2) (cm3s-1) of H$_{\rm 2}$ 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

  \begin{figure} \par\includegraphics[width=8.8cm]{12269fg2.eps} \end{figure} Figure 2:

The formation rate $\alpha $(HD) (cm3s-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

  \begin{figure} \par\includegraphics[width=8.8cm]{12269fg3.eps} \end{figure} Figure 3:

Comparison between gas phase formation (circles) and grain-surface formation (squares) rates of H2 (in cm-3s-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

  \begin{figure} \par\includegraphics[width=8.8cm]{12269fg4.eps} \end{figure} Figure 4:

Comparison between gas phase formation (circles) and grain-surface formation (squares) rates of HD (in cm-3s-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

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

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

Initial download of the metrics may take a while.