© ESO, 2012

1. Introduction

Despite the general agreement that H₂ formation occurs on grain surfaces, the actual formation mechanisms at work in interstellar conditions remain unclear. In this paper we present a new chemical model of H₂ formation that encompasses recent laboratory experimental and theoretical results on this topic and show some consequences of the model results. A follow-up paper will extend these results to other molecules on grains. We recall in Sect. 2 the previous assumptions concerning H₂ formation, in relation to observational issues. Section 3 describes two new implementations concerning i) the contribution of different dust grain sizes to the formation process within the Langmuir-Hinshelwood mechanism (labeled LH in the following), ii) the introduction of the Eley-Rideal mechanism (labeled ER) for H₂ formation. In Sect. 4 we display and compare the results obtained for different pressures and radiation fields in pertaining to both diffuse gas and dense photon dominated regions (PDRs) conditions. Particular focus is placed on the values of the H₂ formation rates themselves and on observational tests such as the emissivities of fine structure transitions of the abundant atoms and ions and infrared H₂ transitions. We present our conclusions in Sect. 5.

2. Historical status

2.1. H₂ formation efficiency

Hollenbach & Salpeter (1971) estimated the binding energies of H and H₂ on various surfaces as well as the corresponding sticking efficiencies and subsequently developed the first theory of the formation mechanism of H₂ on dust grain surfaces in interstellar conditions by diffusion of H atoms on the surface of the grains. The formation rate was subsequently derived thanks to Copernicus far ultraviolet observations of H₂ in absorption in diffuse interstellar clouds (Jura 1974). This mean formation rate, confirmed by FUSE observations (Gry et al. 2002), is 3 × 10^-17   n_H   n(H)   cm³   s^-1. The launch of the ISO satellite and the development of infrared astronomy using ground-based telescopes has allowed detections of warm excited H₂ via its electric quadrupolar transitions in the infrared in PDRs. At long wavelengths (~7−28   μm) pure rotational transitions take place and are detected from space, whereas rovibrational transitions at shorter wavelengths can be observed from the ground in the I, J, H, and K photometric bands. Le Bourlot et al. (1995b) showed that this reaction rate coefficient can be readily expressed in terms of the dust grain properties, i.e. density ρ_gr, dust-to-gas mass ratio G, and the minimum (a_min) and maximum (a_max) values of the grain radii, for which the size distribution is a power-law dn_gr    = A_gr   a^-3.5   da (Mathis et al. 1977), subsequently quoted as the MRN distribution. Le Bourlot et al. (1995b) derived the normalization factor and the formation rate of H₂This expression includes a sticking coefficient s for impinging hydrogen atoms, the mass of atomic hydrogen m_H, the Boltzmann constant k, and the gas temperature T. For typical dust properties (Table 1), this expression becomes . Assuming that the sticking factor is 1 up to 10   K and then decreases with temperature as , the formation rate becomes equal to 4.4 × 10^-17   n_H   n(H)   cm³   s^-1. We note that if we take a single radius of grains with a value of 0.1   μm, the formation rate is now 1.3 × 10^-17   n_H   n(H)   cm³   s^-1. The present prescription also allows us to derive the mean-cross section per H atom For the standard parameters, the corresponding value is 1.95 × 10^-21   n_H.

While this simple physical picture reflects the coherence between the observations and the basic properties of dust, a new perspective arose when the quantitative values of the desorption and diffusion barriers were obtained from temperature desorption experiments reported by Katz et al. (1999); Biham et al. (2001); Biham & Lipshtat (2002); Lipshtat & Biham (2003) and Lipshtat et al. (2004) with various types of surfaces. The corresponding values are reported in Table 1 and compared to older values from Hasegawa et al. (1992). These values were included in theoretical models of H₂ formation by Biham & Lipshtat (2002), who also discussed the relevance of using the rate equations formalism within the LH diffusion mechanism. In the case of small grains, Biham & Lipshtat (2002) argue that only a few hydrogen atoms can stick to the surface so that master equations describing the various probabilities of having N adsorbed H atoms are required. This formalism leads to a large number of coupled equations. Biham & Lipshtat (2002) further introduced the use of moment equations, reducing the number of equations and leading, in the case of H₂ formation, to only two coupled equations. This procedure compares very satisfactorily with the master equation solution (Biham et al. 2005). During the diffusion process, physisorbed H atoms explore the available adsorption sites with a sweeping rate , where ν₀ is a typical vibration frequency and T_d = E_d/k the diffusion energy threshold in Kelvin. The number of available sites, 4πa²n_s, depends on the radius of the grain. The sweeping rate probability was subsequently revisited by Lohmar & Krug (2006) and Lohmar et al. (2009) for various surface geometries.The desorption probability is assumed to follow a similar exponential decrease with the same ν₀ frequency but the binding energy in Kelvin T_b = E_b/k.

The frequency ν₀ is often assumed to be equal to 10¹²   s^-1 (Biham et al. 2001). Hasegawa et al. (1992) introduced the vibrational frequency corresponding to an harmonic oscillator within the binding well potential of the form . The usual expression is where the force constant at the minimum of the interaction potential U(r). The last expression defines a typical size of the bottom of the well d₀ by identification with the expression of Hasegawa et al. (1992). Kim et al. (2011) and Sakong & Kratzer (2010) showed that the stretching vibration wavenumber of H on graphene is about 2600   cm^-1, and the bending about 1200   cm^-1. This translates to vibrational frequencies of a few 10¹³   s^-1, which are rather higher than usual estimates. We deduced a typical size of d₀ ~ 1    Å . The corresponding values of are given in Table 1 for various dust compositions.We verified that model results are almost insensitive to the value of ν₀.

Le Petit et al. (2009) included the possibility of the rejection of the impinging hydrogen atoms, when they land on an already occupied site, and introduced these mechanisms into the Meudon PDR code together with the size distribution of grains, while assuming a uniform dust temperature for all grain sizes.

The net result, for H₂ formation, is that the reaction is only efficient for a narrow window of dust temperatures, ranging typically from 11   K to 19   K, depending on the values of desorption and diffusion barriers, as also pointed out by Chang et al. (2006). This restricted range of grain temperatures is difficult to reconcile with H₂ observations towards dense PDRs (Habart et al. 2004, 2005; Cazaux & Tielens 2004, 2010; Habart et al. 2011), where the dust temperatures reach values close to and even higher than 30   K. As a result, new scenarios have been proposed such as the reaction between a physisorbed and a chemisorbed atom (Habart et al. 2004) or reactions between two chemisorbed atoms (Cazaux & Tielens 2004, 2010).

Other attempts to improve the physics describing the actual process have considered morphological aspects of the dust grain surfaces, such as roughness effects probed through continuous-time random-walk Monte Carlo simulation studies, resulting in effective formation-rate coefficients that depend on the surface temperature and impinging atomic flux (Chang et al. 2006). In summary, despite the large number of laboratory and theoretical studies that have been completed on this subject, a detailed description of H₂ formation is still missing. As a consequence, several authors (Hollenbach et al. 2009; Islam et al. 2010; Kaufman et al. 2006; Sheffer et al. 2011) prefer to use an ad hoc “standard” value of 3 × 10^-17   n_H   n(H)   cm³   s^-1, determined by Copernicus and FUSE observations for diffuse clouds, and/or introduce an additional multiplication factor of two for dense and bright PDRs.

2.2. H₂ excitation in the formation process

The excitation state of H₂ in the formation process is another concern. Experimental studies (Hornekær et al. 2003; Islam et al. 2007; Latimer et al. 2008; Lemaire et al. 2010) demonstrated that H₂ can be rovibrationally excited in the formation process, where the experiments were conducted with two beams of H and D and various types of surfaces. These experiments were performed on cold surfaces, indicating that the formation mechanism involves mainly physisorbed atoms, so a LH mechanism. Alternatively, Farebrother et al. (2000); Bachellerie et al. (2009); Sizun et al. (2010) performed theoretical calculations for the ER formation mechanism of H₂ on graphenic surfaces and derived the rovibrational excitation of the nascent H₂ molecule for different energies of the impinging hydrogen atoms. The most recent study of Sizun et al. (2010) included explicitly the role of the zero point energy, which leads to a slight increase in the rotational excitation within v = 5 and 6. In the astrophysical community, Duley & Williams (1993) proposed that infrared fluorescence of H₂ could be detected towards dark clouds, as a signature of its formation process. Le Bourlot et al. (1995a) showed that the excitation due to the secondary electrons generated by cosmic rays was much more efficient, even for the standard value of the cosmic ionization rate. Similar conclusions were reached by Islam et al. (2010) in a detailed modeling study including the most recent experimental information on the excitation of H₂ and assuming however the “standard” formation rate of 3 × 10^-17   n_H   n(H)   cm³   s^-1. Congiu et al. (2009) also found that internal excitation is quenched before desorption from experiments involving amorphous water ice. However, formation pumping signatures can also be searched for in PDRs and were reported as such by Burton et al. (2002) in the Messier 17 PDR through the mapping of the 6–4 O(3) infrared transition of H₂ at 1.733   μm. We investigate the possible consequences of this result in our present study.

3. Model calculations

Le Petit et al. (2009) studied the efficiency of H₂ formation for different dust temperatures when H₂ formation takes place on amorphous carbon, assuming that the grains have the same temperature, independently of their size. This hypothesis is inadequate as discussed by different authors (Cuppen et al. 2006). Compiègne et al. (2011) explicitly computed the dust temperatures resulting from the balance between photoelectric effects and radiative infrared emission for the different grain sizes. The results are very close to the analytic expressions given by Hollenbach et al. (1991) and used in the Meudon PDR code (Le Petit et al. 2006). We describe now two major extensions, i.e. inclusion of the different dust grain temperatures for surface reactions involving physisorbed hydrogen atoms and the introduction of the ER mechanism for H₂ formation. Other recent updates of the PDR code are described in Appendix A.

3.1. Langmuir-Hinshelwood mechanism

In a pure diffusion process (LH mechanism), physisorbed hydrogen atoms, labelled as “H:”, can either diffuse on the surface, encounter another physisorbed hydrogen atom to form H₂ or leave the surface through thermal desorption or another desorption process such as photodesorption. We derived the detailed rate equations in Appendix B, taking into account size-dependent dust temperatures and rejection effects. The formation rate of H₂ is no longer analytic and involves a numerical integration over the dust grain sizes. In this paper, we assume a MRN distribution law, but the expressions are given for any distribution. The number of particles X physisorbed on the surface of a grain of size a is N_X:(a) and is one of the unknowns we have to compute.

Considering only the balance between H and H₂ (without deuterium here), we can derive, as shown in Appendix B, Eqs. (B.2) and (B.4) In this expression, k_rej, k_ev, k_ph, and k_CR are the chemical reaction rate constants holding respectively for rejection, thermal evaporation, photo-desorption, and cosmic-ray-induced desorption. The parameter k_rej is proportional to both the thermal velocity of the impinging hydrogen atoms and the sticking coefficient, which is also a function of the gas temperature. The parameter S_gr, the total grain surface per cubic centimeter, is equal to , for the MRN distribution law, and d_s is the mean distance between adsorption sites. Among the ejection processes, only thermal evaporation (k_ev(a)) depends on the grain size. The formation of gas phase H₂ is assumed to occur directly after the encounter of adsorbed hydrogen atoms. The corresponding rate is given by Eq. (B.6) The hoping time t_H depends on both diffusion and possible tunneling. The diffusion time, t_D, equal to , now depends on the size of the grain. The tunneling time is derived for a rectangular barrier of height E_b and width d_s¹, where and λ_D is the de Broglie length, given by . The total hoping frequency is a function of the grain size.

To include these equations in the chemical scheme, the mean number of H or H₂ physisorbed on a grain of size a, namely N_H:(a) and N_H2:(a) have to be computed. We obtain from Eqs. (B.1), (B.3) and (B.5): where k_ad is the chemical rate constant for adsorption. After discretization of the grain size distribution, the equations can be solved as ordinary rate equations and are thus easy to incorporate within the other chemical rate equations. However, the number of variables grows significantly.

3.2. Eley-Rideal mechanism

We explicitly introduce adsorption of H atoms onto chemisorbed sites in the PDR model. They are denoted H::. We introduce also the subsequent reactions with impinging gas-phase atomic hydrogen leading to H₂ formation. This process has been studied theoretically by different authors, who have derived the interaction potential with the surface and the reaction probabilities. The static properties, potential well, and possible barrier height depend significantly on the nature of the surface, which allows some liberty in the choice of the actual value to be used in the modeling.

There is a significant difference from the adsorption process described previously, as trapping a gas-phase hydrogen atom into a chemisorption site involves the crossing of a barrier where the exponential term now involves the gas temperature. Given the variety of grain types and surfaces properties found in the ISM, it is probable that this barrier is not unique. Those properties are poorly known and trying to derive a detailed model would, arbitrarily, increase the number of free parameters. Therefore we limit ourselves to a simple approximation that takes into account the major aspects of this process.

With the assumption that H is removed from a chemisorption site only by the formation of a molecule and using results from Appendix C, Eq. (C.2), the H₂ formation rate becomes (1)where α(T) is an efficiency factor that depends on the gas temperature T and  ⟨ nσ_gr ⟩  is the total surface of grains available per cubic centimeter (see the Appendix for a discussion about the parameters). If a more complicated set of species and reactions is used (e.g., formation of HD), then the equations given allow the computations of all necessary abundances, but the rate is no longer analytic.

The corresponding rate is with the standard values displayed in Table 1. We note that the numerical factor is twice the value reported in Sect. 2.1 as H₂ is now produced from the collision between two different partners, a gas-phase H atom and a chemisorbed H::. The order of magnitude of this chemical rate constant is again within the “standard” value derived by Jura (1974). The actual value of k_ER depends on both the chemisorption barrier T₁ and the sticking efficiency α(T) but not the distance between adsorption sites. The influence of T₁ and α(T) on the formation efficiency and line intensities is discussed in Appendix C.

We see below that this reaction has a major effect on the thermal balance at the edge of the cloud as each H₂ formation process as well as each destruction contributes to the gas heating. This is thus a very efficient way to couple the strong ultraviolet radiation field to the gas.

3.3. Excitation during formation

The H₂ formation process releases about 4.5   eV. Most current prescriptions assume an equipartition of this energy between the translational energy, the H₂ internal energy, and the heating of the grain. In this way, about 17   000   K is spread in H₂ rovibrational states. Alternative scenarios were proposed by Black & van Dishoeck (1987). Recent experimental results obtained for cold graphitic surfaces (Islam et al. 2010) and silicates (Lemaire et al. 2010) show that high vibrational levels are preferentially populated in the formation process.

The Meudon PDR code permits us to choose between several scenarios. In this paper, for the LH mechanism, we adopt the default option in which one third of the formation enthalpy is spread in internal energy with a Boltzmann distribution among rovibrational levels (see the description in Le Petit et al. 2006). We note that this is physically relevant only if a sufficient number of levels of H₂ are included in the computations, even if their steady-state populations are negligible.

For the ER process, we follow the energy distribution obtained by Sizun et al. (2010, their Table 1) and consider that the uncertainties are larger than the differences that they computed at incident energies of 15   meV and 50   meV for the impinging H atom. The energy available after formation of a H₂ molecule is thus split into:

• 0.3   eV   (~3500   K): depth of the chemical bound potential well;

• 2.7   eV   (~31   300   K): internal energy, spread following a Boltzmann distribution (mean level: v ~ 5, J ~ 6);

• 0.6   eV   (~7000   K): translational kinetic energy (heating term for the gas);

• 1   eV   (~11   600   K): heating of the grain.

We note that this last term is not yet included in the thermal balance of the grains.

3.4. Heating

This enhanced formation rate at the edge of PDRs provides a new contribution to the gas heating. At the edge of clouds, since the H₂ photodissociation rate is important and the ER mechanism is efficient, the formation-destruction cycle of H₂ has a short timescale. Each destruction by a UV photon leads to about 0.6   eV of kinetic energy, but each subsequent formation event also provides a way to tap the formation enthalpy reservoir. As seen above (Sect. 3.3), most of this energy is in high lying rovibrational levels. Hence, one H₂ formation event is usually followed by the emission of an IR photon, which may be detected. However, for high enough densities, collisional deexcitation with atomic H comes into competition with spontaneous radiative decay. Above a critical density, the collisional process dominates² and the released kinetic energy provides an efficient heating mechanism. This may lead to gas temperatures in the range of a few thousands, which further increases the formation rate (which is proportional to v_th, see Eq. (1)), decreases the cycling timescales, and enhances yet further the heating rate.

This positive loop is quenched, our sticking coefficient decreasing to zero with increasing temperature, which allows an equilibrium temperature to be reached.

For a high density and a high radiation field, this process is by far the dominant one at the edge of the cloud. Figure 1 shows the fraction of heating provided by the ER H₂ formation versus the photo-electric effect in a cloud where the temperature profile is fixed. Close to the edge, the high gas temperature allows the ER process to be fast and most of the heating comes from collisional de-excitation. Once the temperature drops, that process is quenched and photo-electric heating dominates as long as the optical depth is not too high.

The transition region is chemically very active (H and H₂ both being abundant, and the temperature still being quite high). The formation rates of molecules such as OH or CH ⁺  may be enhanced in particular if H₂ internal energy can contribute to overcome the energy thresholds of the chemical reaction. Heat released by exothermic reactions also provides most of the heating missing from the two processes of Fig. 1.

4. Model results

We study the effect of different prescriptions for molecular hydrogen formation on grains. We model isobaric clouds since this is the simplest equation of state for dilute gas in a steady state embedded in a stationary environment where gravitation is negligible: large pressure variations (such as found in a constant density cloud with large temperature variations) would lead to some kind of relaxation. Smaller entities (clumps, filaments, etc.) may be better described as isochoric clouds within these large structures. A grid of these isochoric models, with the same physics as the isobaric models presented below has also been produced. They are not discussed in this paper but can be found in the PDR Data Base (http://pdr.obspm.fr).

Models labelled “A” correspond to a fixed formation rate of molecular hydrogen of 3 × 10^-17   cm³   s^-1, following the prescriptions of Hollenbach et al. (2009) and Islam et al. (2010); in models “B”, H₂ is formed only by LH mechanism, and in models “C” both the LH and ER mechanisms are activated to form H₂. The model geometry is a plane parallel infinite slab with an isotropic radiation field impinging on both sides.

We explore a large domain of gas pressure and intensity of the incident UV radiation field. The gas pressure is given by P = n × T = [n(H) + n(H₂) + n(He)] × T. The radiation field is scaled versus the ISRF evaluated by Mathis et al. (1983) in the UV and visible part of the spectra. Longward of  ~6000  Å  the ISRF is corrected using the infrared observations of Hauser et al. (1998, Boulanger & Gonzalez, priv. comm.). The scaling factor χ is applied only to the VUV range, from 13.6   eV to 5.17   eV (corresponding to the wavelength range 912    Å  to 2400    Å ) to mimic the presence of nearby bright OB stars. It represents the scaling factor of the UV part of the radiation field at a point that receives the radiation field from 4π steradians.

Dust properties, elemental abundances (where {X} represents the total abundance of the X element), and other parameters are given in Tables 1 and 2. The adopted dust extinction curve corresponds to the parameters of Fitzpatrick & Massa (1990) for the Galaxy. The photodissociation probabilities of H₂ and CO are explicitly computed from the different discrete photodissociating transitions. We include a turbulent component to the thermal line width and take a common value of 3   km   s^-1 for all transitions. We choose to display results of isobaric models with values relevant to both diffuse and dense gas.

4.1. Diffuse gas

The chemical network includes only gas phase processes, except for H and H₂ as well as the neutralization of atomic ions onto grain surfaces. It involves 131 chemical species and about 2650 chemical reactions. The analysis of the UV absorption of neutral carbon lines towards a sample of translucent and diffuse line of sights (Jenkins & Shaya 1979; Jenkins & Tripp 2001, 2007) allows us to derive the range of thermal pressures relevant to the Galactic interstellar gas. Table 3 displays results for two isobaric chemical models corresponding to a total visual magnitude and χ = 1 on both sides of the cloud. The impinging flux at the edge is about 0.56 times the Mathis field as some photons are coming from the opposite side of the cloud or are back-scattered by dust. The corresponding radiation pressure at the edge is u/3, where u is the energy density per unit volume, i.e. 1.7 × 10^-14   dynes   cm^-2, always smaller than the considered thermal pressures (respectively 4.14 × 10^-13 and 1.38 × 10^-12   dynes   cm^-2). The size in pc is noted l (see conversion in Appendix D). The grain temperatures depend only on the radiation field intensity. Their common values at the edge of the clouds span a range between 10.7   K and 12.5   K.

We note that the results are typical of cold neutral medium (CNM) properties, following Wolfire et al. (1995). We also find that the dependence of the displayed quantities depend marginally on the assumed formation scenario of molecular hydrogen. This results directly from the range of gas and dust temperatures and the assumed dust size distribution. Thus, in diffuse clouds, the role of the ER mechanism is negligible except at the edge of the low pressure models (P = 3 × 10³   (cm^-3   K) = 4.14 × 10^-13   dynes   cm^-2). We also report T₀₁, the excitation temperature derived from the computed column densities of H₂ in levels J = 0 and J = 1: T₀₁    =    170.5/log (9   N₀/N₁). The values agree with the observations (T ≃ 67   K from FUSE observations, Burgh et al. 2007; Rachford et al. 2009; and Sheffer et al. 2008), especially for models with P = 3 × 10³   cm^-3   K. For the low pressure models, we note that H₂ formation is enhanced in models B and C compared to model A leading to a significantly higher value of the molecular fraction. The FUSE survey of diffuse clouds (Rachford et al. 2009) has not detected clouds with a molecular fraction over 0.6. This is one of the difficulties in reconciling models of diffuse clouds and observations. First, molecular fractions determined by observations also include the atomic hydrogen content over the full line of sight. Second, the molecular fraction depends strongly on the size of the clouds and the UV illumination. The detailed study is beyond the scope of this paper.

4.2. Dense PDRs

Dense and bright PDRs offer the opportunity to detect infrared transitions of molecular hydrogen and fine structure transitions of neutral and ionized atoms (O, C, C ⁺  etc.), which are then used to derive the relevant physical conditions by the observers. We again consider isobaric models irradiated on the observer side by three different radiation fields, respectively 10², 10³, and 10⁴ times the standard interstellar radiation field of Mathis et al. (1983). The backside side is illuminated by a standard radiation field (χ = 1). A total visual magnitude of 10 is assumed.

4.2.1. H₂ formation

We first check the role of the formation process in the atomic to molecular transition properties. Tables 4–7 display results for four different pressures and three radiation field enhancement factors.

The range of dust temperatures at the edge of the clouds has a dramatic effect on the formation rate of H₂, when only diffusion mechanisms are involved (B models), as expected. In this latter case, there is not even a range of visual magnitudes where the local density of H₂ exceeds that of atomic H, when the pressure is 10⁶   cm^-3   K or lower (shown as an ellipsis in the corresponding cells in Tables 4 and 5). The different trends can be understood from analytic developments obtained in a single grain size approximation, as given in Appendices B and C. The effect of pressure is best seen with the χ = 10³ results:

• Grains are warm (from 26   K to 44   K), so that thermal evaporation ofphysisorbed hydrogen is efficient. Therefore, the ER processdominates over the LH one.

• All things equal, increasing the pressure in case B (where only LH is active) gives results in accordance with Eq. (B.8). The temperature at the edge does not vary much, so that [H]  ∝ P and R_H2 is indeed proportional to [H]  as expected.

• Including the ER process, R_H2 first increases with increasing P owing to a (at first) slow rise in temperature. For the highest pressure, the heating mechanism described in Sect. 3.4 is fully operational, and R_H2 decreases as the temperature stabilizes on the right side of the formation rate curve of Fig. C.1.

Increasing χ leads to similar results, with still higher grain temperatures. However, complex UV pumping processes have an increasing role and Eq. (B.8) is less accurate. For a lower radiation field, the range of H densities spanned by increasing pressure is within the critical densities defined in B and the full Eq. (B.7) applies. In both cases, the ER rate closely follows the temperature behavior with a positive feedback being most efficient at high pressure.

We see that allowing for the possibility of chemisorption compensates, in a natural way, the deficiencies of the LH mechanism when dust grains become too warm. The values of the formation rate of molecular hydrogen attained in such bright PDR conditions then become very comparable to the values deduced from the observations (Habart et al. 2004) with an increase by a factor of two to five in the standard 3 × 10^-17   n_H   n(H)   cm³   s^-1 value depending on the impinging incident radiation field. Even if the chemisorption properties are somewhat empirical (the dependence is discussed in Appendix C), we feel that such a mechanism allows us to solve a major difficulty in the theory of H₂ formation.

4.2.2. Emission lines

We now investigate the implications of the new formation scenario on the various emissivities. We display in Tables 8 and 9 the emissivities derived from the previous models. We focus on transitions emitted in the PDR region where photo-chemical and -physical effects are taking place. In line with the results already displayed, we find that the fine-structure emission lines are not very dependent on the formation model of H₂ except in a few cases, where temperatures become of the order of thousands K, owing to the enhanced heating mechanisms described in Sect. 3.4. The molecular hydrogen emission transitions are much more sensitive. We display the emissivities of the pure rotational lines detected routinely by Spitzer, as well as the 1–0 S(1) transition at 2.12 μm detected from the ground. Ratios of these emissivities with respect to those of other vibrational lines detected in the K band are given to compare with past (Burton et al. 2002) and future observations. We have also run 190 type C models with different pressures (10⁵ − 10⁸   cm^-3   K) and intensity scaling factors of the incident radiation field (1 to 10⁶ Mathis ISRF) on the observer side to display the corresponding contours in different figures. The back-side radiation field is set to χ = 1 and the total visual magnitude, in order to neglect radiative effects from the back-side. The other model parameters are given in Table 2 and the dust properties correspond to the standard values displayed in Table 1. We display results when both LH and ER mechanisms are activated for the formation of H₂ (C type). Figures 2 to 6 display H₂, C ⁺ , and O line intensities contours in erg   cm^-2   s^-1   sr^-1 for a face-on geometry.

Intensity maps for other line intensities and other observation angles can be found on our VO-Theory database, PDR Database (PDRDB), under the project “H2 formation LH+ER – 2011”³. The complete set of output data of the models (giving access to abundances profiles, temperature profile, and formation and destruction rates, ...), input parameters, and source code used to produce these models are also available through PDRDB.

Figures 2–4 show that the intensities (excepted for the 0–0 S(0) line) depend mostly on P for high values of χ/P and on χ for low values of χ/P. Habart et al. (2004) explained this behavior in the case of isochoric models. At the position of the H/H₂ transition, and for high χ/n_H, photons are mostly absorbed by grains, whereas for low χ/n_H they are absorbed in H₂ lines (self-shielding). It can be shown that in the first case, the H₂ line intensities are proportional to the H₂ formation rate and so to n_H, whereas, in the second case, they are proportional to the intensity of the incident radiation field. In the case of isobaric models, an analytical derivation is less obvious, but, as we can see in the intensity maps, this trend is valid over three orders of magnitudes of P. In this P range, T varies only by a factor smaller than 10 (Tables 4 to 7). Thus, the variation in P is essentially dominated by the variation in n_H, and the same argument applies.

The orders of magnitude of the H₂ line intensities are in good agreement with those of the observations by Spitzer presented in Habart et al. (2011). For example, comparisons of our intensity maps to observations towards the Horsehead for H₂ lines (Habart et al. 2011) and O at 63   μm (Goicoechea et al. 2009) gives χ of a few hundreds in Mathis units, which is consistent with previous models. The gas pressure is more difficult to constrain in this way since, for moderate χ, the intensities of those lines depend weakly on P. Of all H₂ lines, the 0–0 S(3) line is the more difficult to reproduce. We also note that the intensity of the 1–0 S(1) line, as well as the ratios of vibrational transitions displayed in Table 9 as potential tests of formation pumping, are within the values reported for R1 and R2 towards M17 by Burton et al. (2002) for a pressure of 10⁶   cm^-3   K and radiation field enhancement factor of 10⁴. However, the corresponding modeled R3 ratio is significantly smaller. The models in this paper were developed using typical properties of the interstellar medium (such as grain properties) that may not be applicable to a detailed interpretation of the observations towards specific line of sights. As an example, the star illuminating the Horsehead is a O 9.5 type star with an effective temperature of 33   000   K. Scaling a Mathis ISRF as done here does not spread the UV energy in the same way as in a stellar spectrum, which can affect level excitation. We note that most of the PDR observations presented in Habart et al. (2011) are seen edge on, whereas with our plane-parallel model we provide in this paper line intensities for PDRs seen face on and in PDRDBs of up to 60 degrees. Applications to specific lines of sight will be subsequently performed by including all possible constraints derived from the observations.

5. Conclusions

We have presented implementations of the Langmuir-Hinshelwood and Eley-Rideal mechanisms to describe H₂ formation on grain surfaces within the Meudon PDR code. This is the first of a series of intended studies to investigate the formation of other molecules on grains the light of recent experimental and theoretical progresses. The Eley-Rideal mechanism involving chemisorbed sites is an efficient process for forming H₂ in regions where the gas is hot and dust grains are warm. In previous versions of the code, the Langmuir-Hinshelwood mechanism alone had allowed us to recover the value of the mean formation rate of H₂ derived from VUV observations of diffuse and translucent clouds. However, it led to significantly lower values of H₂ formation rates in denser and strongly illuminated clouds due to the increasing importance of H evaporation relative to diffusion on the surface of the dust grains. Observations of these high excitation PDRs by, e.g., Habart et al. (2011) and references therein inferred mean formation rates typically a factor of two to five higher (with large uncertainties). This is now recovered by our treatment of the Eley-Rideal process, whose efficiency rises with increasing gas temperature (up to a threshold). Thus, considering that observed rates are mean values along the line of sight and across a large variety of physical conditions, regardless of whether our theoretical rates are local values computed from micro-physics considerations alone, the agreement between model and observations is now also quite satisfactory in high excitation regions.

The Eley-Rideal mechanism also provides an efficient means of heating the gas. Following the Sizun et al. (2010) prescription for spreading the energy released by the formation of a H₂ molecule in terms of kinetic energy and internal energy, we have shown that this formation process can be the major heating mechanism under some specific conditions. In particular, in regions where H₂ self-shielding is not yet fully efficient and at densities above the critical density of collisions with H, a significant amount of energy can be transferred into the kinetic energy of the gas. The signatures of formation pumping are not straightforward in the light of present knowledge about the internal state of newly formed H₂.

The efficiency of the Langmuir-Hinshelwood process depends on the temperature of the grains. As shown by many authors, this formation mechanism can only take place in a small window of grain temperatures (Biham & Lipshtat 2002; Biham et al. 2005; Le Petit et al. 2009), namely 10–20 K for amorphous carbons. As a consequence, the LH mechanism becomes the major H₂ formation process at moderate visual extinctions. The strong dependence of the rate on T_gr implies that a proper computing of grain temperature has to be done in numerical models, including the temperature fluctuations of the smallest grains.

We have compared the results of PDR models in which H₂ formation is treated by these detailed mechanisms to models using a constant formation rate with R_H2 = 3 × 10^-17   cm³   s^-1. For diffuse clouds parameters, the main effect is to shift the H/H₂ transition and to increase the molecular fraction. For stronger PDRs, we have shown that for strong radiation fields, the line intensities depend strongly on gas pressure. The difference between the detailed treatment of H₂ formation and a constant formation rate can lead to differences of up to a factor of three in the most commonly observed H₂ line intensities. We provided maps of these line intensities for a set of 190 models for different pressures and incident radiation field intensities.

---

Acknowledgments

This work was partially funded by grant ANR-09-BLAN-0231-01 from the French Agence Nationale de la Recherche as part of the SCHISM project and by the French CNRS national program PCMI. The authors acknowledge computing support in term of resources and services of France Grilles and the EGI e-infrastructure as well as MESOPSL, the computing center for Paris Sciences et Lettres. Some kinetic data have been downloaded from the online KIDA (KInetic Database for Astrochemistry, http://kida.obs.u-bordeaux1.fr) database. We thank our colleagues at the “Laboratoire des Collisions Atomiques et Moléculaires” in Orsay, V. Sidis, N. Rougeaud, and D. Bachellerie for fruitful discussions on the ER process.We thank David Languignon, Nicolas Moreau, Benjamin Ooghe, and Emeric Bron for their help in running some of the models and setting them up on the on-line database.

References

Online material

Appendix A: Other recent updates to the PDR code

Progress have been achieved in computing fine structure and rotational excitation due to collisions with H and/or H₂ as reported in the BASECOL (Dubernet et al. 2006) and LAMDA (Schöier et al. 2005) databases. We have updated and/or implemented the collisional excitation rates that play an active role in the cooling processes and computed explicitly their emission spectrum by solving the statistical equilibrium equations, including radiative pumping by the cosmic background radiation field and dust infrared emission (Gonzalez Garcia et al. 2008). We point out that the new fine-structure excitation collision rates of atomic oxygen computed by Abrahamsson et al. (2007) have a significant impact on the temperature at the edge of PDRs. As an example, for a typical proton density of 10⁴   cm^-3 and a radiation scaling factor of 10, the temperature at the edge is 90   K with the old Launay & Roueff (1977) atomic oxygen collision rates and only 67   K with the values displayed in Abrahamsson et al. (2007). At the present stage, the emission spectra of millimeter and submillimeter transitions of HCN, OH, CH ⁺ , and O₂ are readily computed and the implementation of other molecules is straightforward and depends only on the availability of the collision rates by the relevant perturbers. Another significant issue is the inclusion of the thermal and charge balance of the grains in the overall ionization fraction. The thermal balance is obtained through a coupling with the DustEM program Compiègne et al. (2011), which can be switched on through the F_Dustem parameter (F_Dustem=1). If not (F_Dustem=0), the temperature of the various grain bin sizes is obtained from the formula given in Eq. (5) of Hollenbach et al. (1991), where the actual value of the radiation field is introduced. The determination of the grain charge is obtained from the balance between photoelectric effect and recombination on dust particles, expanding on the treatment of Draine & Sutin (1987) and Bakes & Tielens (1994).

The code also computes the photodissociation rates from the integration of the photodissociation cross-sections, when available, with the interstellar radiation field. The attenuation by dust particles is then directly obtained from the dust properties considered in the model, i.e. their absorption and extinction coefficients, which depend on the size and the nature of the dust particles. Different options are proposed depending on the treatment of the grain temperatures. If F_Dustem=1, we use the absorption and scattering coefficients computed by the DustEM code (Compiègne et al. 2011). If F_Dustem=0, we derive the albedo and dust properties from the extinction curve given by the Fitzpatrick and Massa analytic expansion (Fitzpatrick & Massa 2007), extended towards longer wavelengths by the data from Weingartner & Draine (2001)⁴. It is remarkable that we recover the dust-free photodissociation rates displayed in van Dishoeck (1988) for the Mathis or Draine incident radiation field. The A_V dependence of the photodissociation rates then directly reflects the appropriate dust environment.

Appendix B: Langmuir-Hinshelwood mechanism

Upon landing on a grain, most heavy species may build an ice mantle. In this case, it is possible to account for the total number of physisorbed molecules by integrating over the grain size distribution. This is the subject of a forthcoming paper on surface chemistry (Le Petit et al., in prep.). However, this is most probably not the case for the lightest species (H, H₂, D, HD, ...). We assume here that they only build a single monolayer above either the grain surface or the ice mantle.

In that case, two effects must be taken into account:

• upon landing on a site already occupied by a light species, theimpinging species is rejected to the gas phase;

• binding depends on the (size-dependent) temperature of the grain and thus the steady state depends on the size (and characteristics) of the grain.

Hence, we must compute the number of physisorbed particles of type X on a grain of size a. This is N_X:(a) in the following (in particles per grain, and not in particles per cubic centimeter). The total amount of X on all grains, that follows by integration, is where the last expression is for a MRN size distribution with dn_g = A_gr   n_H   a ^− α   da. A_gr is a normalization factor, n_H is the gas density (in cm^-3), and a the grain radius (in cm). In the following, we use that case as an example, but it is easy to generalize to any distribution. Numerical integration is performed by discretizing the size. When needed, we use where the weights w_i and abscissae a_i are chosen according to the distribution. The number of abscissae is npg. Table B.1 displays the abscissae and weights computed for the parameters of the MRN distribution given in Table 1. These coefficients must be computed anew whenever one changes the range of sizes⁵. Here Remark that the power a^-3.5 does not appear in the discrete sum.

Three types of reactions must be considered:

• Adsorption;

• Ejection;

• Reaction.

We do not consider reactions with heavy atoms nor molecules (including ices) which will be the subject of a follow-up paper.

B.1. Adsorption

If the outer layer of the grain is populated by light species (H, H₂), then any one of them may lead to the rejection of an impinging atom. We consider n_j of these species. The number of accretion events of a species X per unit time interval on a single grain of size a is then where s(X) is the sticking coefficient of species X, d_s the mean distance between adsorption sites (assumed to be identical for all grains), and the term in parentheses takes into account rejection by any species Y_j that is already on the grain. The term is the total number of adsorption sites on a grain of size a. Formally, this equation may be divided into a first-order formation reaction of rate , and n_j different second-order destruction reactions of rate . The relevant creation and destruction equations are thus (B.1)or (B.2)where S_gr is the total surface of grains per unit volume. One can see that, although the form of the reaction terms is preserved (first or second order polynomials in the variables) the total number of individual contributions becomes large (one accretion leads to (n_j × npg + 1) reactions)⁶.

B.2. Desorption processes

Ejection can occur spontaneously (thermal evaporation) or by either photo-desorption or cosmic-ray ejection. All processes are similar in the sense that they only involve a single variable N_X(a) for a grain of size a.

If the vibration frequency of the adsorbed particle is ν₀, the temperature of the grain is T_gr(a), and the binding energy is T_b(X), then the number of evaporation per unit time is If the flux of photons (respectively cosmic rays) is F_ph (respectively F_CR) and the number of particles desorbed by impact is η_ph (respectively η_CR), then the number of desorption is (for a photon) Writing , and , we have (B.3)(B.4)

B.3. Surface reactions

We seek to compute the number of encounters per grain and per unit time. We adopt first the point of view of X:. On a single grain of size a, the number of encounter per s is proportional to 1/t_X the inverse hoping time of X:, the probability of finding a Y: upon landing and the number of X: During the same time, from the point of view of Y:, the number of encounters made is So the total number of encounters is where the factor of takes care of each encounter having been counted twice. Thus, for two surface species, we can write , and (B.5)The production rate of Z in the gas phase, occurring directly after the encounter of two adsorbed atoms, is obtained after integration on the grain size distribution (B.6)

B.4. Approximate H₂ formation rate

For a single grain size, and negligible photodesorption and cosmic-ray desorption, we can derive an analytic approximation to the H₂ formation rate in the spirit of the discussion of Biham & Lipshtat (2002). If H is the only atom sticking to a grain of size a with a sticking probability of 1, then In a steady state, this leads to (B.7)where is the maximum number of H on the grain, and the two critical densities [H] _a and[H] _b are defined as Table B.2 gives the values of the critical densities for different grain temperatures.

The H₂ formation rate per grain is then (with n_g the number of grains per cubic centimeter of gas = ) and can be given analytically from the previous formulae. Both [H] _a and[H] _b vary slowly with the gas temperature, but very strongly with the grain temperature. Thus, we can define two limiting regimes for the gas phase atomic hydrogen density, namely

• Case 1: [H]  ≫ [H] _a > [H] _b. Then

This is possible only for cold grains (typically below 15   K). Hence, it requires both a high density (or pressure) and a weak radiation field.

• Case 2: [H]  ≪ [H] _b < [H] _a. Then (B.8)

This is the case for all grain sizes as soon as T_g is higher than about 25   K. Hence, it applies to all strong radiation field models. In that case, the formation rate increases as the square of the density of H.

The usual expression for the formation rate R_H2 in cm³   s^-1 follows from Since n_g is proportional to n_H, we see that in the first case R_H2 ∝ 1/[H] , whereas in the second R_H2 ∝ [H] . These relations apply only for the approximations made here.

Appendix C: Eley-Rideal mechanism

C.1. Formalism

We consider the impact of a fast atom (hot gas) on a grain. Since a fully detailed description (taking into account all possible kinds of surfaces) is far beyond the capabilities of our model, we search for an approximate mechanism that takes into account the following requirements:

• it is efficient in “hot” gas and on “hot” grains, hence the impinging Hatom must eventually reach a chemisorbed site on the grain;

• it leads to H₂ formation rates that are consistent with observational constraints;

• the number of free parameters remains at the lowest possible number.

Since this process takes place at the edge of the cloud, we assume that the grains are essentially bare (without ice coating) and that the process does not depend on the grain temperature. This approximation is justified since, at the edge of PDRs, grain temperatures (at most 100   K) are much lower than either the gas temperature or chemical binding energies on grain surfaces.

On impact, the gas phase H can find either a free chemisorpsion site or an already chemisorbed H. In the second case, since the formation of H₂ releases 4.5 eV, an energy far higher than the chemisorbed well, we assume that a newly formed H₂ is immediately released in the gas phase. Hence where [H::] _max is the maximum number of chemisorbed H atoms (saturated grains) and [H::]  the corresponding abundance. In addition v_th is the thermal velocity of the gas phase H and we consider the geometrical cross-section to compute the total amount of grain surface per unit volume  ⟨ nσ_gr ⟩ . If the mean distance between chemisorption sites is d_s and is the same on all types of grains, one can see (from purely geometric considerations) that This is true for any grain size distribution. If the gas phase atom impacts a free chemisorption site, we have to estimate the probability that it sticks to the grain. The simplest hypothesis requires that it be proportional to the number of collisions of H with grains per unit of time (v_th    ⟨ nσ_gr ⟩ ), possibly with a barrier to cross (where T is the gas temperature and T₁ the threshold), with a temperature dependent sticking coefficient α(T) and proportional to the “free room” . Hence we write This equation is split in the code into two: a direct formation reaction and a “pseudo” rejection reaction. The corresponding rates are

C.2. Sticking coefficient and choice of T₁

There is not much information on how to define the sticking function α(T), but we expect that it goes to 0 for very high temperatures (the atom just bounces on the grain without there being any time to evacuate the excess kinetic energy). We introduce the empirical form (C.1)In this expression, the index β controls the steepness of the decrease in α(T), and T₂ defines the temperature such that .

We may constrain the value of β with the following considerations:

• an estimate of the velocity v₂ above which theatom bounces back to the gas is given by

• the sticking coefficient is approximated as the fraction of gas phase atom with velocity lower than v₂. Using a Maxwell distribution at temperature T, we have

At low temperatures, this expression tends to 1 as expected. At high temperatures, we may expand this expression as a function of . Defining T₂ by , we get This shows that in Eq. (C.1) is exact. The prefactor is 1.08, which we can take as 1 given the approximations involved; thus Eq. (C.1) is a good approximation over the whole temperature range. The estimate of v₂ gives T₂ in the range 400 − 500   K⁷.

Given α(T), we may investigate which barrier T₁ gives a “standard” formation rate of 3 × 10^-17   cm³   s^-1 at a given temperature T. We find that Postulating that where grains are warm, the gas is warm too, we may require that this standard rate is reached for a temperature in the range [150:450]    K. This translate into a range [100:800]    K for T₁. Our choice of T₁ = 300   K favors an efficient formation, and reflects the idea that chemisorption is easy (but not instantaneous) on grains with plenty of surface defects. This leads to a higher formation rate at high gas temperature as found observationally by Habart et al. (2004).

C.3. Analytical approximation

If hydrogen is the only chemisorbed species, the abundance of H:: can be analytically derived at a steady state which leads to an H₂ formation rate of (C.2)with In this expression, v_th    ⟨ nσ_gr ⟩    [H]  refers to a purely geometric collisional process. Figure C.1 displays the variation in the chemisorption rate as a function of gas temperature (relative to the one at 100   K). The rate is negligible at low temperature owing to the exponential barrier. It grows as the square root of T where the barrier is negligible, then is quenched by the sticking cutoff. The chemisorption efficiency, κ(T), is displayed in Fig. C.1 (right axis). It peaks at a few hundred Kelvin and remains significantly high up to a few thousands of Kelvin.

This behavior is qualitatively very similar to results found by Cuppen et al. (2010, their Fig. 2) from Monte Carlo simulations of H₂ formation including both physisorption and chemisorption. They also found that the formation efficiency increases with temperature for gas temperatures around a thousand K.

In Sect. 4.2.2, we present several line intensities computed from models in which H₂ is formed by the ER and LH mechanisms. In these models, we adopted β = 1.5, T₁ = 300   K, and T₂ = 464   K. Figure C.2 presents the effect of variations in T₁ and T₂ on the intensity of one line of H₂. We note that T₁ is the most important parameter and that the line intensity can be reduced by

a factor  ≃ 5 if this parameter is increased from 300   K to 1000   K. For other lines, such as 0–0 S(0) this decrease can reach a factor of ten. As mentioned above, for real interstellar grains, we can expect to have a large range of T₁ depending on the nature and structure of the grains surfaces. Even if most chemisorbed sites have high thresholds, it only requires a few low thresholds sites for the ER mechanism to be efficient.

Appendix D: A_V to size conversion

For constant dust properties along the line of sight, it is possible to convert optical depth in the visible to a distance (in pc) analytically. We define C_d to be the total proton column-density to color index ratio, and R_V the usual extinction to color index ratio: Then, using A_V = 2.5   log ₁₀(e)   τ_V, we have hence, with , we have We note that , where κ_V + σ_V is the extinction (sum of absorption plus scattering) per H atom by grains at the wavelength of the photometric band V.

All Tables

All Figures