© ESO, 2014

1. Introduction

H₂ is the most abundant molecule in the interstellar medium. It is found in a wide variety of astrophysical environments in which it often plays a leading role in the physics and evolution of objects (see the book by Combes & Pineau des Forets 2001, for a review). As an efficient coolant of hot gas, it controls for instance through thermal balance the collapse of interstellar clouds leading ultimately to star formation. Furthermore, its formation is the first step of a long sequence of reactions leading to the great chemical complexity found in dense clouds. In addition to its physical importance, H₂ is also of great observational usefulness as a diagnostic probe for many different processes in various environments.

The detailed mechanism of its formation is thus a key part of the understanding and modeling of the interstellar medium. Direct gas phase formation is very inefficient and cannot explain the observed abundances. The ion-neutral reaction H + H⁻ is more efficient, but the low H⁻ abundance does not allow sufficient H₂ formation this way. The H₂ molecule is thus thought to form mainly on the surface of dust grains (Gould & Salpeter 1963; Hollenbach & Salpeter 1971). Grains act as catalysts. They provide a surface on which adsorbed H atoms can meet and react, and they absorb that excess energy released by the formation that prevented the gas phase formation of a stable molecule.

This process is usually modeled using the Langmuir-Hinshelwood (LH) mechanism in which weakly bound adsorbed atoms migrate randomly on the surface, and one H₂ molecule is formed each time two atoms meet. This mechanism has been studied in detail using the master equation approach (Biham & Lipshtat 2002), the moment equation approach (Lipshtat & Biham 2003; Le Petit et al. 2009) and Monte Carlo simulations (Chang et al. 2005; Cuppen & Herbst 2005; Chang et al. 2006). Laboratory experiments have also been conduced and their results modeled (Pirronello et al. 1997a,b, 1999; Katz et al. 1999). It is found to be efficient over a limited range of grain temperatures of about 5−15 K for flat surfaces. However, observations of dense photon dominated regions (PDRs) found efficient formation despite higher grain temperatures (Habart et al. 2004, 2011). Various modifications of the mechanism have been proposed to extend the range of efficient formation. Some authors introduced sites of higher binding energies due to surface irregularities (Chang et al. 2005; Cuppen & Herbst 2005). Some introduced a reaction between chemisorbed atoms and atoms migrating from a physisorption site (Cazaux & Tielens 2004; Iqbal et al. 2012). Others proposed that the Eley-Rideal (ER) mechanism for chemisorbed atoms is a relevant formation process at the edge of PDRs (Duley 1996; Habart et al. 2004; Bachellerie et al. 2009; Le Bourlot et al. 2012). In this mechanism, chemisorbed H atoms, which are attached to the surface by a chemical bond, stay fixed on the surface until another atom of the gas falls into the same site, forming a H₂ molecule.

The grain temperatures in illuminated environments are not only higher but also strongly fluctuating. Small grains have small heat capacities, and each UV-photon absorption makes their temperature fluctuate widely. This effect has received much attention because some components of dust emission are produced during these transient temperature spikes (Desert et al. 1986; Draine & Li 2001; Li & Draine 2001; Compiègne et al. 2011). Those fluctuations are stronger in strong radiation field environments like PDRs. Moreover, usual dust size distributions favor small grains so that they contribute the most to the total dust surface, making the small grains dominant for surface reactions. Those fluctuations can thus have an important effect on H₂ formation, especially in strongly illuminated environments, as temperature spikes are likely to cause massive desorption. This effect has so far only been studied by Cuppen et al. (2006) using Monte Carlo simulations for the LH mechanism, computing the H₂ formation rate in one single specific diffuse cloud condition under a standard interstellar radiation field. The dependance on the astrophysical parameters was not investigated. Very recently, similar models using the same Monte Carlo method were presented in Iqbal et al. (2014) for silicate surfaces. Those models also included the ER mechanism. The authors remain limited by the prohibitive computation time (several weeks) required by the kinetic Monte Carlo method when including temperature fluctuations. We study here this problem for both ER and LH in a wider range of environments, including PDRs and their high radiation fields, and study the effect induced on cloud structure and observable line intensities by introducing the computation of the fluctuation effect into full cloud simulations.

We propose an analytical approach of the joint fluctuations problem (temperature and surface coverage) based on the master equation, which takes the form of an integral equation. Its numerical resolution incurs a high computational cost, which grows with grain size. In the case of the ER mechanism, the absence of non-linear term in the chemistry allows a decomposition of the two-dimensional full master equation into two one-dimensional equations : a marginal master equation for the temperature fluctuations and an equation on the conditional average population (conditional on T). This makes the numerical computation much more tractable. In the case of the LH mechanism, such simplification is not possible. The computationally heavy full resolution is used to build and verify a fast and accurate approximation of the solution. The final method proposed here allows computation of the total formation rate with a computation time of the order of one minute.

In Sect. 2, we introduce the simple physical model that we use and the exact resolution method. In Sect. 3, approximations of the solution are constructed for both mechanisms. Section 4 presents the results of the numerical computation of the H₂ formation rate using the full method for ER and the approximation for LH. Section 5 then shows how this computed effect affects the structure of the cloud in PDR simulations. Section 6 is our conclusion.

2. The model

2.1. Physical description

We use a simple physical model of H₂ formation on grains to be able to solve the problem of the coupled fluctuations and to obtain an estimate of the importance of dust temperature fluctuations for H₂ formation.

Our system is a spherical grain of radius a, density ρ_gr, mass $\mathit{m}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{a}}^{\mathrm{3}}\hspace{0.17em}\mathit{\rho }$, and heat capacity C(T). Its thermal state can be equivalently described by its temperature T or by its thermal energy $\mathit{E}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{T}}\mathrm{d}{\mathit{T}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{C}\mathrm{\left(}{\mathit{T}}^{\mathrm{\prime }}\mathrm{\right)}$ ¹. For polycyclic aromatic hydrocarbons (PAHs), we define the size a as the equivalent size of a sphere of the same mass (see Compiègne et al. 2011).

We consider two groups of phenomena:

• Interaction with the radiation field through photon absorption and emission, which are both described as discrete events.

• Adsorption of H atoms of the gas onto the surface of the grain and H₂ formation on the surface (physisorption and chemisorption are considered separately, as are the two corresponding H₂ formation mechanisms).

2.1.1. Photon absorption and emission

Under an ambient isotropic radiation field of specific intensity I_U(U) in W m^-2 J^-1 sr^-1 (in this article, U always denotes a photon energy), the power received by the grain at a certain photon energy U is P_abs(U) = 4π²a²Q_abs(U) I_U(U), where Q_abs(U) is the absorption efficiency coefficient of the grain at energy U. In later parts of the paper, we are interested in transition rates between states of the grain. The rate of photon absorptions at this energy U is simply $\begin{array}{ccc}{\mathit{R}}_{\mathrm{abs}}\mathrm{\left(}\mathit{U}\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{abs}}\mathrm{\left(}\mathit{U}\mathrm{\right)}}{\mathit{U}}\mathrm{·}& & \end{array}$We postulate that the grain emits according to a modified black body law with a specific intensity Q_abs(U) B_U(U,T), where B_U(U,T) is the usual black body specific intensity. The power emitted at energy U is P_em(U,T) = 4π²a²Q_abs(U) B_U(U,T) and the photon emission rate is $\begin{array}{ccc}{\mathit{R}}_{\mathrm{em}}\mathrm{\left(}\mathit{U,T}\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{em}}\mathrm{\left(}\mathit{U,T}\mathrm{\right)}}{\mathit{U}}\mathrm{·}& & \end{array}$These events, which are supposed to occur as Poisson processes, cause fluctuations of the temperature. We treat these fluctuations exactly in Sect. 2.2.

When neglecting these fluctuations, the equilibrium temperature T_eq of the grain is defined as balancing the instantaneous emitted and absorbed powers ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{U}\hspace{0.17em}{\mathit{P}}_{\mathrm{abs}}\mathrm{\left(}\mathit{U}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{E}}_{\mathrm{grain}}}\mathrm{d}\mathit{U}\hspace{0.17em}{\mathit{P}}_{\mathrm{em}}\mathrm{\left(}\mathit{U,}{\mathit{T}}_{\mathrm{eq}}\mathrm{\right)}\mathit{,}$(1)where the upper bound on the right hand side accounts for the finite total energy of a grain. This upper cutoff of emission frequencies can become very important for small cold grains, effectively stopping their cooling.

In the following, we use a standard interstellar radiation field (Mathis et al. 1983) and apply a scaling factor χ to the UV component of the field. We measure the UV intensity of those fields using the usual ${\mathit{G}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{u}}_{\mathrm{Habing}}}\hspace{0.17em}{{}^{\mathrm{\int }}}_{\hspace{0.17em}\mathrm{912}\AA }^{\mathrm{2400}\AA }\mathrm{d}\mathit{\lambda }\hspace{0.17em}{\mathit{u}}_{\mathit{\lambda }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}$, where u_Habing = 5.3 × 10^-15 J m^-3.

The dust properties (C(T), Q_abs(U) and ρ) are taken from Compiègne et al. (2011). We consider PAHs, amorphous carbon grains, and silicates dust populations and use the properties used in this reference (see references therein, in their Appendix A) for each of those populations.

2.1.2. Surface chemistry via Langmuir-Hinshelwood mechanism

We first consider the physisorption-based LH mechanism. In this mechanism, atoms bind to the grain surface due to Van der Waals interactions, leading to the so-called physisorption. These weakly bound atoms are able to migrate from site to site. Formation can then occur when two such migrating atoms meet.

Real grain surfaces are irregular, and binding properties vary from site to site due to surface defects. To keep the model simple enough to allow a detailed modeling of temperature fluctuations, we assume identical binding sites uniformly distributed on the grain surface. The number of sites is ${\mathit{N}}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{a}}^{\mathrm{2}}}{{\mathit{d}}_{\mathrm{s}}^{\mathrm{2}}}$, where d_s is the typical distance between sites. We take d_s = 0.26 nm as in Le Bourlot et al. (2012). We assume the physisorption sites to be simple potential wells of depth E_phys without barrier.

Atoms from a gas at temperature T_gas collide with the grain at a rate k_coll = πa²n(H) v_th, where ${\mathit{v}}_{\mathrm{th}}\mathrm{=}\sqrt{\frac{\mathrm{8}\mathit{k}\hspace{0.17em}{\mathit{T}}_{\mathrm{gas}}}{\mathit{\pi }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}}}$ is the thermal velocity. We call s(T_gas) the sticking probability and discuss later our choice for this function. The probability to land on an empty site is $\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}$ and the physisorption rate is thus $\begin{array}{ccc}{\mathit{R}}_{\mathrm{phys}}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\hspace{0.17em}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\right)\mathrm{·}& & \end{array}$We assume that gas atoms falling on an occupied site are never rejected and react with the adsorbed atom to form a H₂ molecule, which is immediately desorbed by the formation energy. This assumption of Eley-Rideal-like reaction for physisorbed atoms is similarly made in the model of Cuppen et al. (2006) with the difference that they do not consider desorption at the formation of the new molecule. This direct Eley-Rideal process is expected to have a large cross section at temperatures relevant to the interstellar medium (Martinazzo & Tantardini 2006), and it is thus a reasonable assumption to assume that it dominates over rejection. To estimate the influence of this assumption on our results, we computed the contribution of this direct ER-like reaction process to the total mean formation rate through physisorption and found it to be a negligible fraction (always less than 1% of the total average formation rate). Despite the very small contribution of this Eley-Rideal-like process to the physisorption-based formation rate in our model, we keep calling the physisorption-based mechanism LH mechanism to distinguish it from the chemisorption-based ER mechanism.

Surface atoms can then evaporate. The adsorbed atoms are assumed to be thermalized at the grain temperature and must overcome the well energy E_phys. We take this desorption rate for one atom to be of the form ${\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}\hspace{0.17em}\exp \left(\mathrm{-}\frac{{\mathit{T}}_{\mathrm{phys}}}{\mathit{T}}\right)$, where T is the grain temperature, T_phys = E_phys/k, and ν₀ is a typical vibration frequency taken as ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}\hspace{0.17em}\sqrt{\frac{\mathrm{2}\hspace{0.17em}{\mathit{E}}_{\mathrm{phys}}}{{\mathit{d}}_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}}}$ with a typical size d₀ = 0.1 nm (Hasegawa et al. 1992). The total desorption rate on the grain is then $\begin{array}{ccc}{\mathit{R}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathit{n}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{.}& & \end{array}$To migrate from site to site, an atom must cross a barrier of height E_migr. It can do so by thermal hopping or by tunneling for which we use the formula from Le Bourlot et al. (2012). The migration rate for one physisorbed atom is thus $\begin{array}{ccc}{\mathit{k}}_{\mathrm{migr}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}\hspace{0.17em}\exp \left(\mathrm{-}\frac{{\mathit{T}}_{\mathrm{migr}}}{\mathit{T}}\right)\mathrm{+}\frac{{\mathit{\nu }}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{migr}}^{\mathrm{2}}{\sinh }^{\mathrm{2}}\left({\mathit{d}}_{\mathrm{s}}\sqrt{\mathrm{2}{\mathit{m}}_{\mathrm{H}}\mathit{k}\mathrm{(}{\mathit{T}}_{\mathrm{migr}}\mathrm{-}\mathit{T}\mathrm{)}}\mathit{/}ħ\right)}{\mathrm{4}\mathit{T}\mathrm{(}{\mathit{T}}_{\mathrm{migr}}\mathrm{-}\mathit{T}\mathrm{)}}}\mathit{,}& & \end{array}$where T_migr = E_migr/k.

The main formation process is the meeting of two physisorbed atoms during a migration event. In a simple approximation (see Lohmar & Krug 2006; Lohmar et al. 2009, for a more detailed treatment of the sweeping rate), we take this formation rate to be $\begin{array}{ccc}{\mathit{R}}_{\mathrm{form}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{migr}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\frac{{\mathit{n}}^{\mathrm{2}}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}& & \end{array}$We also take direct reaction of a physisorbed atom with a gas atom in an ER-like mechanism into account, with a rate $\begin{array}{ccc}{\mathit{R}}_{\mathrm{form}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}& & \end{array}$This term is only significant for very low dust temperatures (a few K, depending on the collision rate with gas atoms).

We assume immediate desorption of the newly formed molecule because of the high formation energy compared to the binding energy, thus making a similar assumption to recent models (Iqbal et al. 2012; Le Bourlot et al. 2012). Experimental results indicate that a small fraction of the formed molecules may stay on the surface at formation (Katz et al. 1999; Perets et al. 2007). Theoretical studies (Morisset et al. 2004a, 2005) suggest that even when the formed molecule does not desorb immediately at formation, it desorbs after a few picoseconds by interacting with the surface due to its high internal energy. It is thus equivalent to immediate desorption compared to the timescales of the other events.

If the temperature fluctuations are neglected, an equilibrium surface population can be computed using the equilibrium rate equation ${\mathit{R}}_{\mathrm{phys}}\mathrm{-}{\mathit{R}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{R}}_{\mathrm{form}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{-}{\mathit{R}}_{\mathrm{form}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathrm{0}$ (the factor of 2 comes from the fact that one formation through migration event reduces the population by two). We find $\begin{array}{ccc}{\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{(}{\mathit{T}}_{\mathrm{gas}}{\mathrm{)}}^{\mathrm{)}}{\mathit{k}}_{\mathrm{coll}}}{\mathrm{4}\hspace{0.17em}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\left)}\mathrm{\right)}}\right)\\ & & \end{array}$(2)and the corresponding equilibrium formation rate on one grain at temperature T is ${\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{migr}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{\mathrm{\left[}{\mathit{n}}_{\mathrm{eq}}\mathrm{\right(}\mathit{T}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{{\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}$(3)Once again this expression is not valid for a fluctuating temperature and the full calculation is described in the next Sect. 2.2.

For the binding and migration energies, we use the values given in Table 1 associated to different substrates.

Many different expressions of the sticking coefficient have been given in the literature. For simplicity, we use the same sticking probability for the physisorption-based LH mechanism and the chemisorption-based ER mechanism, except for the additional presence of a barrier in the ER case. We use the expression from Le Bourlot et al. (2012): $\begin{array}{ccc}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\left(\frac{{\mathit{T}}_{\mathrm{gas}}}{{\mathit{T}}_{\mathrm{2}}}\right)}^{\mathit{\beta }}}& & \end{array}$using T₂ = 464 K and β = 1.5, which gives results close to the expression of Sternberg & Dalgarno (1995).

We also simply assume an equipartition of the formation energy among the kinetic energy of the newly formed molecule, its excitation energy, and the heating of the grain.

The total volumic H₂ formation rate is often parametrized as R_H2 = R_fn_Hn(H) to factor out the dependency to the atomic H abundance n(H) and to the dust abundance, which is proportional to the total proton density n_H, coming from the collision rate between gas H atoms and dust grains. The standard value of R_f is 3 × 10^-17cm³ s^-1. We also use the formation efficiency defined for one grain as $\frac{\mathrm{2}\hspace{0.17em}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{k}}_{\mathrm{coll}}}$, which represents the fraction of the H atoms colliding with the grain that are turned into H₂.

2.1.3. Surface chemistry via Eley-Rideal mechanism

We also consider the chemisorption-based ER mechanism. The H atoms of the gas phase that hit the grain can bind to the surface by a covalent bond after overcoming a repulsive barrier. This process is called chemisorption (see for instance Jeloaica & Sidis 1999). Once chemisorbed, the atoms either evaporate and return to the gas phase, or react if another gas H atom lands on the same adsorption site.

We assume the chemisorption sites to be similarly distributed with the same inter-site distance d_s = 0.26 nm. Each site is a potential well of depth E_chem and is separated from the free state by a potential barrier of height E_barr.

The collision rate k_coll is the same as in the LH case. However, we have to take the presence of the chemisorption barrier into acount in the sticking coefficient. Following Le Bourlot et al. (2012), we use a sticking probability $\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\mathrm{=}\exp \left(\mathrm{-}\frac{{\mathit{E}}_{\mathrm{barr}}}{\mathit{k}\hspace{0.17em}{\mathit{T}}_{\mathrm{gas}}}\right)\hspace{0.17em}{\left(\mathrm{1}\mathrm{+}{\left(\frac{{\mathit{T}}_{\mathrm{gas}}}{{\mathit{T}}_{\mathrm{2}}}\right)}^{\mathit{\beta }}\right)}^{-1}$, where the second term accounts for rebound of high energy atoms without sticking with T_barr = E_barr/k = 300 K, T₂ = 464 K, and β = 1.5. The energy barrier for chemisorption on perfect graphitic surfaces has been found to be of the order of 0.15−0.2 eV (but adsorption in para-dimer configuration has a much lower barrier of 0.04 eV) (Sha & Jackson 2002; Kerwin & Jackson 2008). As discussed in Le Bourlot et al. (2012), the choice made here of a lower barrier for chemisorption is based on the idea that real grain surfaces are not perfect flat surfaces like graphite, but present a large number of defects. Surface defects are expected to induce a strongly reduced barrier for chemisorption (Ivanovskaya et al. 2010; Casolo et al. 2013). Edge sites on PAHs behave in a similar way (Bonfanti et al. 2011). Moreover, Mennella (2008) found that chemisorption in aliphatic CH_2,3 groups could also lead to efficient Eley-Rideal formation with a very low chemisorption barrier.

Finally, only atoms arriving on an empty adsorption site stick, corresponding to a probability $\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}$, where n is the number of chemisorbed atoms.

The chemisorption rate is thus $\begin{array}{ccc}{\mathit{R}}_{\mathrm{chem}}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\hspace{0.17em}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\right)\mathrm{·}& & \end{array}$Atoms hitting an occupied site react with the adsorbed atom to form a H₂ molecule. Some theoretical studies show a very small activation barrier of 9.2 meV (Morisset et al. 2004b) for this reaction, but the existence of this barrier is under debate (Rougeau et al. 2011). To stay consistent with the model of Le Bourlot et al. (2012) to which we compare our results in Sect. 5, we neglect this barrier compared to the stronger adsorption barrier. We found that including it in the model reduces the formation efficiency by less than 10% in the range of gas temperatures where the ER mechanism is important compared to the physisorption-based LH mechanism.

As the formation reaction releases 4.5 eV, which is much more than the binding energy, we assume that the new molecule is immediately released in the gas. The formation rate is then $\begin{array}{ccc}{\mathit{R}}_{\mathrm{form}}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}& & \end{array}$Last, chemisorbed atoms can evaporate. This effect is usually negligible at typical grain temperatures, but the fluctuating temperatures of small grains can easily make excursions into the domain where evaporation is significant. The adsorbed atoms are supposed to be thermalized at the grain temperature, and must overcome the well energy E_chem. The desorption process is similar to the LH case, and we take again ${\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}\hspace{0.17em}\exp {}^{\mathrm{(}}\mathrm{-}{\frac{{\mathit{T}}_{\mathrm{chem}}}{\mathit{T}}}^{\mathrm{)}}$, where T_chem = E_chem/k, and ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}\hspace{0.17em}\sqrt{\frac{\mathrm{2}\hspace{0.17em}{\mathit{E}}_{\mathrm{chem}}}{{\mathit{d}}_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}}}$ (with d₀ = 0.1 nm). The total desorption rate on the grain is then $\begin{array}{ccc}{\mathit{R}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathit{n}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{.}& & \end{array}$Theoretical studies find the chemisorption energy on graphene to be in the range 0.67−0.9 eV (Sha & Jackson 2002; Lehtinen et al. 2004; Casolo et al. 2009). On small PAHs, Bonfanti et al. (2011) find binding energies in the range 1−3 eV for edge sites and in the range 0.5−1 eV for non edge sites. In this overall range 5800 −35 000 K, we choose a low value of 7000 K to estimate the maximum effect of temperature fluctuations on the ER mechanism.

We sum up the model parameters in Table 2.

When neglecting the fluctuations of the grain temperature, we can calculate the equilibrium surface population n_eq using the equilibrium rate equation R_chem − R_form − R_des(T) = 0, which gives ${\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}\frac{\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}}}\mathit{,}$(4)and the corresponding equilibrium formation rate on one grain at temperature T is ${\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\frac{{\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}$(5)However, when temperature fluctuations are taken into account, this equilibrium situation is not valid. The method we use to calculate this effect is described in Sect. 2.2.

The energy liberated from the formation reaction is split between the excitation and kinetic energy of the molecule and the heating of the grain. As in Le Bourlot et al. (2012), we take the results of Sizun et al. (2010) for the ER mechanism : 2.7 eV goes into the excitation energy of the molecule, 0.6 eV into its kinetic energy, and 1 eV is given to the grain.

We must note that the part of the formation energy given to the grain creates another coupling between surface population and grain temperature. This fraction of the formation energy (1 eV) is not completely negligible compared to photon energies. This additional retro-coupling between surface chemistry and temperature is ignored in the statistical calculation presented in this article. However, its effect on the equilibrium situation when neglecting fluctuations is evaluated in Appendix A. It is found to be negligible as the power given to the grain is usually much smaller that the radiative power it receives.

We can again define the formation parameter R_f as R_H2 = R_fn_Hn(H) and the formation efficiency for one grain as $\frac{\mathrm{2}\hspace{0.17em}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{k}}_{\mathrm{coll}}}$.

2.2. Method

To take into account temperature fluctuations, we adopt a statistical point of view, and we describe the statistical properties of the fluctuating variables.

The necessary statistical information on our system is contained in the probability density function (PDF) f(X), giving the probability to find the system in the state X that is defined by the two variables E, the thermal energy of the grain (equivalent to its temperature T), and n, the number of adsorbed H atoms on its surface. The function f is thus the joint PDF in the two variables. As we treat our two formation mechanisms independently, n counts chemisorbed or physisorbed atoms depending on which mechanism we are discussing.

The time evolution of this PDF is governed by the master equation: $\begin{array}{ccc}\frac{\mathrm{d}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathrm{\int }}_{\mathrm{states}}\mathrm{d}\mathit{Y}\hspace{0.17em}{\mathit{p}}_{\mathit{Y}\mathrm{\to }\mathit{X}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{Y}\mathrm{\right)}\mathrm{-}{\mathrm{\int }}_{\mathrm{states}}\mathrm{d}\mathit{Y}\hspace{0.17em}{\mathit{p}}_{\mathit{X}\mathrm{\to }\mathit{Y}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathit{,}& & \end{array}$where p_{X → Y} is the transition rate from state X to Y. This equation expresses that the probability of being in a given state varies in time due to the imbalance between the rate of arrivals and the rate of departures. In a stationary situation, the two rates, which we call as the gain and loss rates respectively, must compensate each other.

We are thus looking for a solution to the equation $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{states}}\mathrm{d}\mathit{Y}\hspace{0.17em}{\mathit{p}}_{\mathit{Y}\mathrm{\to }\mathit{X}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{Y}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{states}}\mathrm{d}\mathit{Y}\hspace{0.17em}{\mathit{p}}_{\mathit{X}\mathrm{\to }\mathit{Y}}\mathrm{=}\mathrm{0.}& & \end{array}$Two different kinds of transitions occur: the emission or absorption of a photon, which changes only the thermal energy E of the grain, and the adsorption, desorption, or reaction of hydrogen atoms, which changes only the surface population n. The two variables are only coupled by the transition rates for the chemical events, which depend on the grain temperature. The population fluctuations are thus entirely driven by the temperature fluctuations with no retroactions (as mentioned before we neglect grain heating by the formation process).

As the transition mechanisms modify only one variable at a time, we split each of the gain and loss terms into two parts, respectively for photon processes (affecting only E) and for surface chemical processes (affecting only n), ${\mathit{G}}_{\mathit{E}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{+}{\mathit{G}}_{\mathit{n}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{-}{\mathit{L}}_{\mathit{E}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{-}{\mathit{L}}_{\mathit{n}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(6)with $\begin{array}{ccc}{\mathit{G}}_{\mathit{E}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathrm{+}\mathrm{\infty }\\ \mathrm{0}\end{array}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{p}}_{{\mathit{E}}^{\mathrm{\prime }}\mathrm{\to }\mathit{E}}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\\ {\mathit{G}}_{\mathit{n}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \sum _{{\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{{\mathit{N}}_{\mathit{s}}}{\mathit{p}}_{{\mathit{n}}^{\mathrm{\prime }}\mathrm{\to }\mathit{n}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{E,}{\mathit{n}}^{\mathrm{\prime }}\mathrm{\right)}\\ {\mathit{L}}_{\mathit{E}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{p}}_{\mathit{E}\mathrm{\to }{\mathit{E}}^{\mathrm{\prime }}}\\ {\mathit{L}}_{\mathit{n}}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\sum _{{\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{{\mathit{N}}_{\mathit{s}}}{\mathit{p}}_{\mathit{n}\mathrm{\to }{\mathit{n}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathit{,}\end{array}$where we have simplified the notations for the transition rates as transitions only affect one variable at a time and the transitions affecting E (photon absorptions or emissions) have rates that are independent of n.

We now first show how the thermal energy marginal PDF f_E(E) (or the equivalent temperature PDF f_T(T)) can be derived from our formalism, which reduces to equations similar to previous works (Desert et al. 1986; Draine & Li 2001). We then show how the formation rates in the LH case and in the ER case can be derived.

2.2.1. Thermal energy PDF

Summing Eq. (6) over all n values, we obtain ${\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{E}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{p}}_{\mathit{E}\mathrm{\to }{\mathit{E}}^{\mathrm{\prime }}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{p}}_{{\mathit{E}}^{\mathrm{\prime }}\mathrm{\to }\mathit{E}}\hspace{0.17em}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(7)where f_E(E) = ∑ _nf(E,n) is the marginal PDF for the single variable E. The marginal PDF f_E can easily be converted into a grain temperature PDF f_T through the relationship $\begin{array}{ccc}{\mathit{f}}_{\mathit{T}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\hspace{0.17em}\mathit{C}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{.}& & \end{array}$We thus obtain an independent master equation for the temperature fluctuations.

Detailing the transition rates in Eq. (7), we can rewrite it as $\begin{array}{ccc}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{+}{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}{\mathit{R}}_{\mathrm{em}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E,T}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}}{{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathit{,T}\mathrm{(}\mathit{E}\mathrm{\left)}\mathrm{\right)}}\mathrm{·}& & \end{array}$(8)This equation is an eigenvector equation for the linear integral operator defined by the right-hand side: f_E(E) = ℒ [ f_E ] (E). A stationary PDF for the grain thermal energy is a positive and normalized eigenvector of this operator for the eigenvalue 1. The existence and unicity of this eigenvector is proven in Appendix B. Moreover, the other eigenvalues all have real parts that are lower than 1 (see Appendix B). We can thus converge toward the solution by building a sequence ℒⁿ [ f_E ] (the exponent refers to operator composition) from an initial guess. We solve this equation numerically by choosing an energy grid {E_i} and solving for piecewise linear functions on this grid.

2.2.2. Eley-Rideal mechanism

We first present the method for the ER mechanism as the linearity of the chemical rates in n simplifies the problem significantly. We can avoid directly solving the full master equation for the joint PDF in two variables (Eq. (6)).

We are interested in the average ER formation rate $\begin{array}{ccc}⟨{\mathit{R}}_{\mathrm{form}}⟩\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{E}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{N}}_{\mathrm{s}}}\mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\hspace{0.17em}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}& & \end{array}$Knowing f_E(E), we can rewrite it as: $\begin{array}{ccc}⟨{\mathit{R}}_{\mathrm{form}}⟩\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{E}\hspace{0.17em}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\hspace{0.17em}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\frac{⟨\mathit{n}\mathrm{|}\mathit{E}⟩}{{\mathit{N}}_{\mathrm{s}}}\mathit{,}& & \end{array}$where $⟨\mathit{n}\mathrm{|}\mathit{E}⟩\mathrm{=}{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{N}}_{\mathrm{s}}}\mathit{n}\frac{\mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}}{{\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{E}\mathrm{\right)}}$ is the conditional expectation of the surface population at thermal energy E. It is by definition the expected value (or ensemble average) of the surface population of a grain, knowing that this grain has thermal energy E. Multiplying Eq. (6) by n and summing over all n values gives the equation governing this quantity ⟨n | E⟩, $\frac{{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}}{\mathit{M}\mathrm{\left(}\mathit{E}\mathrm{\right)}}\mathrm{=}⟨\mathit{n}\mathrm{|}\mathit{E}⟩\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}⟨\mathit{n}\mathrm{|}{\mathit{E}}^{\mathrm{\prime }}⟩\hspace{0.17em}\mathit{K}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(9)where $\begin{array}{ccc}\mathit{M}\mathrm{\left(}\mathit{E}\mathrm{\right)}& \mathrm{=}{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}& \\ & \mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{)}}{{\mathit{N}}_{\mathit{s}}}\mathrm{+}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{E}\mathrm{\right)}& \end{array}$and $\begin{array}{ccc}\mathit{K}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ \\ & \end{array}& & \end{array}$After a similar discretization, this is a linear system of equations. However, while being nonsingular, it converges exponentially fast toward a singular system as the grain size a grows, making standard numerical procedures unusable.

To avoid this problem, we rewrite this equation. Multiplying by M(E)f_E(E) and integrating it over E yields ${\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\left[\frac{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{)}}{{\mathit{N}}_{\mathit{s}}}\mathrm{+}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\right]{\mathit{f}}_{\mathit{E}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}⟨\mathit{n}\mathrm{|}{\mathit{E}}^{\mathrm{\prime }}⟩\mathit{.}$(10)Dividing again by M(E) and subtracting it from the initial equation gives $\mathrm{0}\mathrm{=}⟨\mathit{n}\mathrm{|}\mathit{E}⟩\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}⟨\mathit{n}\mathrm{|}{\mathit{E}}^{\mathrm{\prime }}⟩\hspace{0.17em}{\mathit{K}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(11)where ${\mathit{K}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathit{K}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{)}\mathrm{+}{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}\mathit{M}\mathrm{\left(}\mathit{E}\mathrm{\right)}}\hspace{0.17em}{\mathit{f}}_{\mathit{E}}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}$. This is again an eigenvector equation for the eigenvalue 1 of a linear integral operator. After discretization, we numerically compute this eigenvector and use Eq. (10) as a normalization condition. This procedure proved to work on the entire range of relevant grain sizes.

2.2.3. Langmuir-Hinshelwood mechanism

For the LH mechanism, the formation rate contains a quadratic term. When trying to apply a similar method to the LH mechanism, the equivalent of Eq. (11) is then an infinite system of equations on the conditional moments of the population ⟨n | E⟩, ${}^{\mathrm{⟨}}\mathit{n}^{\mathrm{2}}\mathrm{|}{\mathit{E}}^{\mathrm{⟩}}$, ${}^{\mathrm{⟨}}\mathit{n}^{\mathrm{3}}\mathrm{|}{\mathit{E}}^{\mathrm{⟩}}\mathit{...}$ We thus directly solve the main master equation (Eq. (6)), despite the computational burden.

Writing explicitly the transition rates in Eq. (6), we get $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \int \begin{array}{c}\mathit{E}\\ \mathrm{0}\end{array}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\\ & & \mathrm{+}{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E,T}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\\ & & \mathrm{+}\mathit{f}\mathrm{(}\mathit{E,n}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}\mathrm{-}\mathrm{1}}{{\mathit{N}}_{\mathrm{s}}}\right)\\ & & \mathrm{+}\mathit{f}\mathrm{(}\mathit{E,n}\mathrm{+}\mathrm{1}\mathrm{)}\left[\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{\mathit{n}\mathrm{+}\mathrm{1}}{{\mathit{N}}_{\mathrm{s}}}\right]\\ & & \mathrm{+}\mathit{f}\mathrm{(}\mathit{E,n}\mathrm{+}\mathrm{2}\mathrm{)}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{2}\mathrm{)}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}{{\mathit{N}}_{\mathrm{s}}}\\ & & \mathrm{-}\left[{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}\right]\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\\ & & \mathrm{-}\mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}[{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\right)\mathrm{+}\mathit{n}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{\mathit{n}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}{{\mathit{N}}_{\mathrm{s}}}]\mathit{,}\end{array}$(12)where, as boundary conditions, all expressions n − 1, n + 1 and n + 2 are implicitly taken to be 0 if they become negative or greater than N_s.

We define the total loss rate as $\begin{array}{ccc}\mathit{Q}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathrm{+}\mathrm{\infty }\\ \mathit{E}\end{array}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\right)\mathrm{+}\mathit{n}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{\mathit{n}}{{\mathit{N}}_{\mathrm{s}}}\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{\mathit{n}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}{{\mathit{N}}_{\mathrm{s}}}\mathit{,}\end{array}$the integral operator (for photon induced transitions) as $\begin{array}{ccc}\mathrm{𝒢}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{E}\\ \mathrm{0}\end{array}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\\ & & \mathrm{+}{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E,T}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\mathit{,}\end{array}$and the discrete jump operator (for chemical transitions) as $\begin{array}{ccc}\mathrm{𝒥}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}& \mathrm{=}& \mathit{f}\mathrm{(}\mathit{E,n}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{\mathit{n}\mathrm{-}\mathrm{1}}{{\mathit{N}}_{\mathrm{s}}}\right)\\ & & \mathrm{+}\mathit{f}\mathrm{(}\mathit{E,n}\mathrm{+}\mathrm{1}\mathrm{)}\left[\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{\mathit{n}\mathrm{+}\mathrm{1}}{{\mathit{N}}_{\mathrm{s}}}\right]\\ & & \mathrm{+}\mathit{f}\mathrm{(}\mathit{E,n}\mathrm{+}\mathrm{2}\mathrm{)}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{2}\mathrm{)}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}{{\mathit{N}}_{\mathrm{s}}}\mathrm{·}\end{array}$We can then rewrite Eq. (12)in a simplified form as $\mathit{f}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{Q}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}}\left[\mathrm{𝒢}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\mathrm{+}\mathrm{𝒥}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E,n}\mathrm{\right)}\right]\mathit{.}$(13)This is an eigenvector equation for the eigenvalue 1, which is formally similar to Eq. (8). For the same reasons as previously, we expect the operator defined by the right hand side to have all its eigenvalues with real parts that are strictly lower than 1, except the eigenvalue 1, which is simple. We can thus find the solution by iterating the application of the operator from an initial guess.

The equation is solved numerically using this iterative procedure. The computation time is found to explode when increasing the grain size a. In addition to the number of possible values of n increasing as a², the number of iterations necessary to converge toward the stationary solution is also found to quickly increase with a. Grouping the values of n in bins to reduce the matrix size does not solve the problem, as the reduction of the computation time due to the smaller matrices is found to be balanced by a roughly equivalent increase in the number of iterations necessary to converge to a given stationarity threshold.

We can thus only use this method up to moderate grain sizes (~20 nm with 2 days of computation). In later sections, we see, however, that this is sufficient to observe the range of sizes for which the temperature fluctuations significantly impact the chemistry. In addition, a much quicker yet sufficiently accurate approximation is presented in Sect. 3.

3. Approximations

In this section, we construct fast approximations of the formation rate based on simple physical arguments and compare their results with those of the exact method presented in Sect. 2.2. An approximation is especially necessary for the LH mechanism to avoid the prohibitive computational cost of the exact method. A similar approximation is given for completeness for the ER mechanism. Those methods assume the temperature PDF f_T(T) to be known. It can be computed, for instance, by the method described in Sect. 2.2.1.

3.1. Approximation of the Langmuir-Hinshelwood formation rate

We build our method around an approximation of ⟨n|T⟩, the average population on the grains that are at a given temperature T. We first write a balance equation for this average population, taking into account both chemical processes (changing only the grain’s population), and grains leaving or reaching this temperature T.

We first consider the chemical processes. Those are described in Sect. 2.1.2:

• The grain can gain atoms through adsorption (average rate ${\mathit{k}}_{\mathrm{ads}}\mathrm{=}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}{}^{\mathrm{(}}\mathrm{1}\mathrm{-}{\frac{⟨\mathit{n}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}}^{\mathrm{)}}$).

• The grain can lose atoms through desorption (rate ⟨n|T⟩k_des(T)).

• The grain can lose atoms through LH reaction, a term we express later. For now, let us write it as ${\mathrm{⟨}}^{}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{LH}}|{\mathit{T}}^{\mathrm{⟩}}$. We discuss this term in depth below.

• The grain can lose atoms through direct ER reaction (rate ${\mathit{k}}_{\mathrm{coll}}\frac{⟨\mathit{n}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}$).

On the other hand, the grain can also leave the temperature T. At equilibrium, the rate at which grains leave the state T is exactly balanced by the rate of arrivals from other states. The net effect on the conditional average ⟨n|T⟩ can be a loss or a gain depending on whether the average population of the grains arriving from other states is higher or lower than ⟨n|T⟩. In general, noting k_leave(T), the rate at which grains leave the state T, and ⟨n_arrivals|T⟩, the average population on grains arriving in state T, the net gain rate is k_leave(T)(⟨n_arrivals|T⟩ − ⟨n|T⟩).

Up to now, no approximation has been made.

As we expected high temperatures states, where desorption dominates all other processes, to have negligible formation and extremely low average population, we focus our approximation on a low temperature regime.

We want to estimate what fraction of the grains leaving state T will come back bare (or with a negligible population compared to ⟨n|T⟩). We make the following approximation: grains leaving T to reach the regime where desorption dominates, come back with no population. Fluctuations that do not reach the desorption regime leave the surface population unchanged. We define the temperature T_lim that delimits these regimes as the temperature at which the desorption rate for one atom becomes equal to the adsorption rate for the grain (we call E_lim the corresponding thermal energy): $\begin{array}{ccc}{\mathit{T}}_{\lim }\mathrm{=}\frac{{\mathit{T}}_{\mathrm{phys}}}{\ln \left(\frac{{\mathit{\nu }}_{\mathrm{0}}}{{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}}\right)}\mathrm{·}& & \end{array}$We thus define ${\mathit{k}}_{\mathrm{phot}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\max \mathrm{(}{\mathit{E}}_{\lim }\mathrm{-}\mathit{E}\mathrm{(}\mathit{T}\mathrm{\left)}\mathit{,}\mathrm{0}\mathrm{\right)}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{U}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{\left(}\mathit{U}\mathrm{\right)}\mathit{,}$(14)the rate at which grains leave state T to reach a temperature above the limit. The net loss rate due to fluctuations is then k_phot(T)⟨n|T⟩.

We have thus obtained a rate equation for the conditional average population at temperature T: $\begin{array}{ccc}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{⟨\mathit{n}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}\right)& \mathrm{=}& \mathrm{(}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}{\mathrm{)}}^{\mathrm{)}}⟨\mathit{n}|\mathit{T}⟩\\ & & \end{array}$(15)This equation is in closed form if we can express ${\mathrm{⟨}}^{}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{LH}}|{\mathit{T}}^{\mathrm{⟩}}$ as a function of ⟨n|T⟩ only. We distinguish the two regimes ⟨n|T⟩ > 1 and ⟨n|T⟩ < 1. From now on, we simplify the notations by noting s = s(T_gas).

3.1.1. Regime with ⟨n|T⟩ > 1

In this regime, we assume that the discrete nature of the number of surface atoms is negligible and treat this variable as continuous. We then have ${}^{\mathrm{⟨}}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{LH}}|{\mathit{T}}^{\mathrm{⟩}}\mathrm{=}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\frac{⟨{\mathit{n}}^{\mathrm{2}}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}$. We also know that when ⟨n|T⟩ ≫ 1, direct ER reaction is likely to dominate so that a precise determination of the LH reaction rate is not necessary. We thus make the simple approximation ^⟨n²|T^⟩ = ⟨n|T⟩². Using this approximation in Eq. (15), we then get $\begin{array}{ccc}⟨\mathit{n}|\mathit{T}⟩& \mathrm{=}& \frac{\mathrm{1}\mathrm{+}\mathit{s}}{\mathrm{4}}\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{sites}}\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{)}}\right)\\ & & \mathrm{\times }\left[\sqrt{\mathrm{1}\mathrm{+}\frac{\mathrm{8}\mathit{s}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}{\mathrm{)}}^{\mathrm{2}}}\frac{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}}\frac{{\mathit{N}}_{\mathrm{sites}}}{{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{sites}}\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{s}\mathrm{)}}\right)}^{\mathrm{2}}}}\mathrm{-}\mathrm{1}\right]\mathrm{·}\end{array}$(16)If this results becomes smaller than 1, we switch to the other regime.

3.1.2. Regime with ⟨n|T⟩ < 1

In this regime, discretisation effects are very important. Formation is dominated by LH reaction, which can only happen if two atoms are present on the grain at the same time. We use reasoning that is similar to the modified rate equation approach developed by Garrod (2008).

As the average population is low, we can make the approximation ⟨n|T⟩ ≃ p(n = 1 | T), where the right hand side is the conditional probability of having one surface atom knowing that the temperature is T. The LH formation rate can then be computed the following way. The formation rate is the rate at which gas atoms fall on a grain that already had one adsorbed atom AND reaction occurs before any other process removes one of the atoms. The first part of the sentence gives a rate $\begin{array}{ccc}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{sites}}}\right)\hspace{0.17em}\mathit{p}\mathrm{(}\mathit{n}\mathrm{=}\mathrm{1}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{T}\mathrm{)}\mathrm{\simeq }{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{sites}}}\right)\hspace{0.17em}⟨\mathit{n}|\mathit{T}⟩& & \end{array}$using the previous approximation. Once we have two atoms on the grain, the probability that they react before anything else removes one atom is $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{\frac{\mathrm{2}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{sites}}}}{\frac{\mathrm{2}{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{sites}}}\mathrm{+}\frac{\mathrm{2}{\mathit{k}}_{\mathrm{coll}}}{{\mathit{N}}_{\mathrm{sites}}}\mathrm{+}\mathrm{2}\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}\mathrm{·}& & \end{array}$This gives $\begin{array}{ccc}\mathrm{⟨}{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{LH}}|{\mathit{T}}^{\mathrm{⟩}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{(}\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{sites}}}\mathrm{)}}{\mathrm{1}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{sites}}\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}}\hspace{0.17em}⟨\mathit{n}|\mathit{T}⟩\mathit{.}& & \end{array}$We inject this expression in Eq. (15)and finally find the solution for this regime: $⟨\mathit{n}|\mathit{T}⟩\mathrm{=}\frac{{\mathit{N}}_{\mathrm{sites}}}{{\mathit{N}}_{\mathrm{sites}}\frac{\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}\mathit{s}}\mathrm{+}\frac{\mathrm{1}\mathrm{+}\mathit{s}}{\mathit{s}}\mathrm{+}\frac{\mathrm{2}\mathrm{(}{\mathit{N}}_{\mathrm{sites}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{1}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{+}{\mathit{N}}_{\mathrm{sites}}\frac{{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{k}}_{\mathit{mig}\mathrm{\left(}\mathit{T}\mathrm{\right)}}}}}\mathrm{·}$(17)

3.1.3. Total formation rate

The total formation rate is then computed by integrating over the temperature PDF f(T) and distinguishing the two regimes: $\begin{array}{ccc}⟨{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}⟩& \mathrm{=}& \int \begin{array}{c}{\mathit{T}}_{\mathrm{switch}}\\ \mathrm{0}\end{array}\mathrm{d}\mathit{T}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\left(\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{N}}_{\mathrm{sites}}}⟨\mathit{n}|\mathit{T}⟩\mathrm{+}\frac{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{sites}}}{⟨\mathit{n}|\mathit{T}⟩}^{\mathrm{2}}\right)\\ \mathrm{+}{\mathrm{\int }}_{{\mathit{T}}_{\mathrm{switch}}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{Tf}\mathrm{\left(}\mathit{T}\mathrm{\right)}⟨\mathit{n}|\mathit{T}⟩& & \left(\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{N}}_{\mathrm{sites}}}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{(}\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{sites}}}\mathrm{)}}{\mathrm{1}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{coll}}\mathrm{+}{\mathit{N}}_{\mathrm{sites}}\mathrm{\left(}{\mathit{k}}_{\mathrm{des}}\mathrm{\right(}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}\mathrm{\left)}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{mig}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}}\right)\mathit{,}\end{array}$(18)where T_switch is the temperature at which ⟨n|T⟩ becomes <1 (⟨n|T⟩ > 1 for T<T_switch, and ⟨n|T⟩ < 1 for T>T_switch).

This approximation is compared to the exact method of Sect. 2.2.3 in Appendix C. It is found to give a very accurate estimate of the total average formation rate (within at most 6%). This approximation is used in all results shown in the following sections as the exact method is not practically usable.

3.2. Approximation of the Eley-Rideal formation rate

A similar approximation can be constructed in the case of the ER mechanism. We again write a rate equation for the conditional average ⟨n|T⟩, including a fluctuation loss term. We use the same approximation for this loss term with the limiting temperature being $\begin{array}{ccc}{\mathit{T}}_{\lim }\mathrm{=}\frac{{\mathit{T}}_{\mathrm{chim}}}{\ln \left(\frac{{\mathit{\nu }}_{\mathrm{0}}}{{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}}\right)}\mathit{,}& & \end{array}$and the loss rate being k_phot(T)⟨n|T⟩ with k_phot(T) as defined as previously by Eq. (14). The resulting rate equation is then directly obtained in closed form: $\begin{array}{ccc}{\mathit{k}}_{\mathrm{coll}}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{⟨\mathit{n}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}\right)\mathrm{=}\mathrm{(}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{(}\mathit{T}{\mathrm{)}}^{\mathrm{)}}⟨\mathit{n}|\mathit{T}⟩\mathrm{+}{\mathit{k}}_{\mathrm{coll}}\frac{⟨\mathit{n}|\mathit{T}⟩}{{\mathit{N}}_{\mathrm{sites}}}\mathit{,}& & \end{array}$and no further approximation is needed. The solution is $\begin{array}{ccc}⟨\mathit{n}|\mathit{T}⟩\mathrm{=}\frac{\mathit{s}}{\mathrm{1}\mathrm{+}\mathit{s}}{\mathit{N}}_{\mathrm{sites}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{sites}}}{\mathrm{1}\mathrm{+}\mathit{s}}\left(\frac{{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{phot}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}}\right)}\mathit{,}& & \end{array}$and the average H₂ formation rate can then be computed as $\begin{array}{ccc}⟨{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}⟩\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{T}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{coll}}}{{\mathit{N}}_{\mathrm{sites}}}⟨\mathit{n}|\mathit{T}⟩\mathrm{·}& & \end{array}$A comparison of this method with the exact result is also performed in Appendix C. However, as the exact method is easily tractable, the approximation is not used in the results presented in the rest of this article.

4. Results

The method described in the previous two sections was implemented as a stand alone code called Fredholm. The temperature PDF is computed using the exact method of Sect. 2.2.1 and the ER formation rate using the exact method of Sect. 2.2.2, while we use the approximation presented in Sect. 3 for the LH formation rate. This code is now used to study the effect of the temperature fluctuations on the formation rate and the influence of the external conditions (gas conditions and radiation field). We first present a few results on the temperature PDFs before presenting the results on H₂ formation by the ER and LH mechanisms.

4.1. Dust temperature probability density functions and grain emission

The thermal energy PDFs that we compute are shown in Fig. 1, which highlights the change of shape as we go from small grains to bigger grains. For small sizes, most of the grains are in the lowest energy states (the non-zero limit at zero energy can be easily derived from the master equation), while the high energy tail extends very far (the abrupt cut around 2 × 10^-18 J is due to the Lyman cutoff present in the radiation field). As we go to bigger grains, the fraction of grains in the low energy plateau decreases, and the PDF takes a more peaked form but keeps a clear asymmetry. For the bigger sizes, the PDF becomes very sharp around its maximum and its relative width decreases towards zero, but the shape remains asymmetrical and could be asymptotically log-normal.

As expected, the average temperature converges toward the usual equilibrium temperature at large grain sizes but is significantly lower at small sizes, as shown in Fig. 2. High temperature grains have a much higher contribution to the average energy balance through their high emissivity. Thus, a very low fraction of high temperature grains is sufficient to lower the average temperature of the PDF compared to the equilibrium temperature. The temperature increase at very low grain sizes, which we can see on both the equilibrium and average temperatures, is due to the property that the smallest sizes have decreased cooling efficiency (as they cannot emit photons with higher energies than their own internal energy, their emission spectrum is reduced).

This part of the computation has been done before (Desert et al. 1986; Draine & Li 2001) and especially in the DustEM code (see Compiègne et al. 2011) to compute the infrared emission of dust. We checked the results of our code, which is named Fredholm, by comparing both the temperature PDFs, and the derived grain emissivity to the public DustEM code. We found a good agreement between the results, and despite small differences in the temperature PDFs (due mainly to the continuous cooling approximation in DustEM), the final emissivities are extremely close. Figure 3 shows the final comparison on a dust population corresponding to Compiègne et al. (2011) with PAHs, small carbonaceous grains, large carbonaceous grains, and large silicate grains.

4.2. H₂ formation rate: LH mechanism

We apply the approximation method described in Sect. 3 to compute the H₂ formation rate on physisorption sites. The results presented here concern unshielded gas in which the equilibrium rate equation method gives very inefficient LH formation rate due to high dust temperatures (i.e., PDR edges). For results in a more shielded gas, see Sect. 5 describing the coupling with full cloud simulations.

The results shown in this section use binding and barrier energies corresponding to either ices or amorphous carbon surfaces (see Table 1) as specified for each figure. Silicates have intermediate values and are thus not shown here. On the other hand, the optical and thermal properties of the grains are those of amorphous carbon grains.

We first show results for ice surfaces. Figure 4 shows the formation efficiency on one carbonaceous grain as a function of the grain size under various radiation field intensities and compares the equilibrium rate equation method with our estimation of the fluctuation effects. The results for small grains are different by several orders of magnitude; the method presented here giving a much more efficient formation. For large grains, the two methods converge as expected.

To understand those differences, we show how the temperature of a 3 nm dust grain is distributed compared to the formation efficiency domain in Fig. 5. We show both the efficiency curve that would be obtained from an equilibrium computation without fluctuations (dashed line) and the true conditional mean efficiency $\frac{\mathrm{2}⟨{\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}|\mathit{T}⟩}{{\mathit{k}}_{\mathrm{coll}}}$(solid line), for fluctuating grains when they are at the given temperature T. We see two different effects caused by the fluctuations. First, the efficiency domain is reduced in both range and maximum value. This occurs because grains do not stay in the low temperature regime but undergo frequent fluctuations going to the high temperature domain where their surface population evaporates quickly before cooling down. Their average coverage when they are at low temperatures is thus decreased. The amplitude of this effect depends on the competition between the fluctuation timescale and the adsorption timescale.

The other effect comes from the spread of the temperature PDF. While the average temperature of the PDF falls in the low efficiency domain (and the equilibrium temperature, being higher, is even worse), a significant part of the temperature PDF actually falls in the high efficiency domain, resulting in a quite efficient formation as seen in Fig. 4. The grain spends most of its time in the cold states and makes quick excursions into the high temperatures that shift the average toward high temperatures. The resulting formation rate is thus determined by the fraction of the time spent in the high efficiency domain and not by the equilibrium or average temperature of the grain. As described in Sect. 4.1, the smaller the grain, the more asymmetric the temperature PDF with most of the PDF below its average and a long high temperature tail. The smallest grains are thus the most affected, and we see (Fig. 4) that formation on the very small grains is significantly improved even for the highest radiation field intensities. As the radiation field intensity is decreased, larger grains are affected and become efficient. Grains up to 20 nm are affected. We note here that Cuppen et al. (2006) do not consider grains smaller than 5 nm and thus do not include part of the affected range; they, however, find a similar effect from 5 nm to ~10 nm on their flat surface model.

To evaluate the importance of this effect on the global formation rate, we integrate our formation rate on a MRN-like power law distribution (Mathis et al. 1977) of carbonaceous grains with an extended size range that goes from 1 nm to 0.3 μm (see Table 3) and obtain a volumic formation rate R_H2.

As this distribution strongly favors small grains in term of surface, this integrated rate is strongly affected by the effect on small grains. Figure 6 shows this total formation rate as a function of the radiation field intensity for various gas densities (still using ice energy values). We obtain very efficient formation from the LH mechanism in unshielded PDR edge conditions for high gas densities. In all cases, the formation rate is increased by several orders of magnitude compared to the equilibrium computation. As the smallest grains remain efficient even for strong radiation fields, the total formation rate decreases more slowly when increasing the radiation field.

Moreover, we note a difference in the dependency to the gas density. The gas density dependency that comes from the collision rate is already factored out in the definition of the R_f parameter shown on this graph. The equilibrium rate equation result is in the low efficiency regime, where the efficiency is determined by the competition between the adsorption rate (proportional to n_H) and the desorption rate, and is thus proportional to the gas density. The result from the method presented here is mainly determined by the fraction of the time spend in the high efficiency regime. This fraction is fully determined by the photon absorption and emission processes and is independent of the gas density. The final result for high gas densities (n_H> 10³ cm^-3) is thus largely independent of the gas density. For low densities, the reduction of the conditional average efficiency curve (noted in Fig. 5) becomes dominant. As this effect results from the competition between the fluctuation rate (due to photons) and the adsorption rate (proportional to n_H), the resulting formation parameter R_f becomes proportional again to the gas density for low densities (n_H< 10³ cm^-3). This change of behavior corresponds to the point where the fluctuation timescale becomes significantly smaller than the adsorption timescale.

Figure 7 shows the same result when using binding and barrier energies corresponding to amorphous carbon surfaces. As a consequence of the higher binding energy for amorphous carbon surface, the equilibrium result yields higher formation rates. The increase caused by the fluctuations is thus less dramatic but remains very large for strong radiation fields.

We can already note here that the formation rates obtained in those unshielded PDR edge conditions are comparable to the typical value of 3 × 10^-17 cm³ s^-1. The LH and the ER mechanisms thus become of comparable importance in PDR edges for high densities, as seen in the results of full cloud simulations shown in Sect. 5.

4.3. H₂ formation rate: ER mechanism

The first step in our computation of the formation rate is the conditional mean of the surface population knowing the grain temperature. We show this quantity (as a fraction of the total number of sites) in Fig. 8, as compared to the coverage fraction for a grain at equilibrium at constant temperature (this equilibrium coverage fraction is independent of the grain size). As expected, the equilibrium coverage shows a plateau as long as the desorption rate is negligible compared to the adsorption rate. This plateau is lower than unity as the sticking is less than unity. Then the equilibrium coverage decreases as desorption becomes dominant. The dynamical effect of the temperature fluctuations has a consequence that the actual average populations are not at equilibrium. The grains undergoing a fluctuation to high temperature arrive in their new state with their previous population and then go through a transient desorption phase, which lasts a certain time, thus increasing the average population at the high temperatures compared to the equilibrium. When cooling down back to the low temperature states, re-accretion of their population takes some time, and this effect lowers the average population at low temperatures. Because this lowering effect is important, as seen in Fig. 8, it means that the characteristic time between two temperature spikes is sufficiently short compared to the adsorption timescale. The mechanism at play is thus that grains that go to high temperatures lose their population but then cool down faster than they re-adsorb their surface population and do not spend enough time in the low temperature state before the next fluctuation to have regained their equilibrium population. Of course, as the fluctuations are less dramatic for bigger grains (as a given photon energy makes a smaller temperature change: a Lyman photon at 912 Å brings a 1 nm grain to ~500 K but a 2 nm grain to only ~200 K), the effect disappears for larger grains, as seen in Fig. 8.

The formation rate on one grain is then simply ${\mathit{r}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{T}\hspace{0.17em}{\mathit{f}}_{\mathit{T}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}{\mathit{k}}_{\mathrm{coll}}\frac{⟨\mathit{n}\mathrm{|}\mathit{T}⟩}{{\mathit{N}}_{\mathrm{s}}}$. The resulting formation efficiency on one grain is shown in Fig. 9 as a function of grain size.

As explained above, the smaller grains have a lower coverage than equilibrium due to temperature fluctuations and, thus, a decreased formation rate. Only grains below 2 nm are affected, as the same photon absorption causes a smaller temperature fluctuation for bigger grains. For stronger radiation fields, fluctuations occur with a shorter timescale, leaving less time for the surface population to reach equilibrium. The effect is then stronger and extends to higher sizes.

The affected range seems very small, but we have to remember that usual dust size distributions are such that the small grains represent most of the surface (e.g., the MRN distribution). Thus, the contribution from the smallest grain sizes dominates the global formation rate. By integrating over the dust size distribution (same as in Sect. 4.2), we obtain the volumic formation rate in the gas, which we can express as R_H2 = R_fn(H) n_H defining R_f. In the following, we compare the integrated formation rate that we obtain to the one obtained with the rate equation at constant temperature as used in Le Bourlot et al. (2012) and plot our results as a fraction of this reference formation rate.

The two main parameters are the radiation field intensity and the H atom collision rate (and the sticking coefficient, which modifies the effective collision rate). We use standard interstellar radiation fields with a scaling factor as explained in Sect. 2.1.1.

Figure 10 shows the effect of these two parameters. We show the formation rate as a function of the radiation field intensity for different gas densities to probe the effect of the collision rate. Note that the reference rate equation result at constant temperature does not depend on either of these parameters: the density dependency (from the collision rate) is eliminated in the definition of R_f, and the radiation field has no effect when the fluctuations are neglected (the equilibrium temperature is too low to trigger desorption). The effect we show is thus fully the effect of grain temperature fluctuations. As expected, the formation rate decreases when either the radiation intensity is increased or the collision rate is decreased, as the effect results from a competition between the fluctuation rate and the collision rate. The observed reduction remains limited and reaches a ~40% reduction or lower in the range of parameters that we explore. The result that an effect affecting only a tiny fraction of the range of grain sizes can have a significant effect on the integrated formation rate illustrates how the smallest sizes dominate the chemistry in the usual MRN size distribution.

The choice of the lower limit of this size distribution is thus of crucial importance. The lower this minimal size, the larger the usual equilibrium rate equation formation rate will be, as the total dust surface is increased for a constant total dust mass. However, as we lower this limit, we extend the range affected by the fluctuations, and thus obtain a stronger decrease compared to the rate equation result, as shown in Fig. 11. We observe no effect if the size distribution starts above the affected size (e.g., for a_min = 3 nm), while a lower limit causes a stronger effect (up to ~60% with a_min = 0.5 nm).

Last, as noted in Sect. 2.1.3, the microphysical parameters of the adsorption process are not well known. The adsorption barrier affects the sticking and changes the effective collision rate (see Fig. 10 for the effect of the collision rate). The vibration frequency ν₀ has a negligible effect in the range of reasonable values. The most important microphysical parameter is the chemisorption energy T_chem. Figure 12 shows the strong influence of this parameter on the result. A lower chemisorption energy allows for smaller fluctuations to trigger desorption. Bigger grains are thus affected, resulting in a much lower formation rate (e.g., a reduction of more than 90% for T_chem = 3500 K).

5. Integration into the Meudon PDR code

We now investigate the effects of the changes in the formation rate found in the previous section on the overall structure and chemistry of an interstellar cloud with special attention to the observable line intensities.

Our code Fredholm, which computes the H₂ formation rate from the LH and ER mechanisms as described in Sects. 2 and 3, has been coupled to the Meudon PDR code described in Le Petit et al. (2006).

The Meudon PDR code models stationary interstellar clouds in a 1D geometry with a full treatment of the radiative transfer (including absorption and emission by both gas and dust), thermal balance, and chemical equilibrium done in a self-consistent way. In its standard settings, the formation of H₂ on dust grains is treated using rate equations and including both the LH and ER mechanism, as described in Le Bourlot et al. (2012). We have now enabled the possibility to compute the H₂ formation rate using our code Fredholm from the gas conditions and the radiation field sent by the PDR code at each position in the cloud, thus taking temperature fluctuations into account. We compare the results of those coupled runs with the simpler rate equation treatment. The dust population in the PDR models presented here is a single component following a MRN-like power-law distribution from 1 nm to 0.3 μm with an exponent of −3.5, using properties corresponding to a 70−30% mix of graphites and silicates above 3 nm and a progressive transition to PAH properties from 3 nm to 1 nm.

A direct comparison between the rate equation approximation and a full account of stochastic effects is thus possible. In addition to the radiation field illuminating the front side of the cloud, whose strength is indicated for each model, the back side of the cloud is always illuminated by a standard interstellar radiation field (G₀ = 1) in the following models.

5.1. Detailed analysis

5.1.1. Effects on PDR edges

We focus first on the effects on PDR edges, as the main effects of a change in formation rate come from the shift of the H / H₂ transition that is induced. We have seen in Sect. 4 that the full treatment of fluctuations reduces the efficiency of the ER mechanism while strongly increasing the efficiency of the LH mechanism (compared to a rate equation treatment). The resulting effect in full cloud models depends on the relative importance of the two mechanisms.

For low radiation field models, the full treatment increases the LH mechanism at the edge to a level comparable or even higher than the ER mechanism, as shown in Fig. 13. The LH formation rate is increased by more than one order of magnitude and becomes higher than the ER rate, while it was one order of magnitude below in the standard rate equation model. The net result is an increase of the total formation rate by a factor of 4. While we were expecting the ER mechanism to be slightly decreased by the full treatment, we can note that it is slightly increased here. This is due to the gas temperature sensitivity of the ER mechanism. In this low radiation field model, the heating of the gas due to H₂ formation is not negligible compared to the usually dominant photoelectric effect. A higher total H₂ formation rate increases the heating of the gas resulting in a higher temperature at the edge (see Fig. 14). The ER mechanism is thus made slightly more efficient.

The main effect is a shift of the H / H₂ transition from A_V = 4 × 10^-2 to A_V = 8 × 10^-3, as shown of Fig. 15. As a result, the column density of excited H₂ is increased, as shown in Table 4. As we are in low excitation conditions, the effect is mainly visible on the levels v = 0,J = 2 and v = 0,J = 3.

We also note a significant increase of the H₂ abundance before the transition, which results in increased CH⁺ formation at the edge. This effect and the higher temperature at the edge result in stronger excited lines of CH⁺ with differences of a factor of 2 on line (4−3). The excitation temperature computed from the observed intensities above (4−3) closely follow the change in gas temperature at the edge. The increase in intensity thus becomes larger as we go to larger levels, but those intensities remain most probably unobservable in this case.

The situation is different for high radiation field PDRs. The increase in the efficiency of the LH mechanism at the edge is not sufficient to make it comparable to the ER mechanism and the net result is a decrease of the total formation rate as the ER efficiency is reduced (see Fig. 16). This results in inverse effects with the H / H₂ transition slightly shifted away from the edge, but of smaller amplitude as the formation rate difference is much smaller.

We present a wider investigation of the effects on observable intensities in Sect. 5.2, covering a large domain of cloud conditions.

5.1.2. Effects on cloud cores

Unexpectedly, temperature fluctuations of small dust grains are found to have an effect on H₂ formation even in cloud cores where the UV field has been fully extinguished. As an example, Fig. 17 shows the atomic hydrogen abundance in a specific cloud (n_H = 10⁴ cm^-3, χ = 100, total A_V = 10) and compares the standard models to the coupled models with full treatment of the fluctuations for the different binding energy values of Table 1. The residual atomic H abundance in the core is found to be systematically higher when taking the fluctuations into account, which is a direct consequence of a lower formation efficiency (the efficiencies are shown for the same models in Fig. 18). We recall that the formation efficiency is defined as the ratio between the formation rate and half the collision rate between H atoms and grains (whose expression is ${\mathit{k}}_{\max }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{\mathit{v}}_{\mathrm{th}}{}^{\mathrm{\right(}}\mathrm{H}^{\mathrm{\left)}}\hspace{0.17em}\frac{{\mathrm{𝒮}}_{\mathrm{gr}}\hspace{0.17em}{\mathit{n}}_{\mathrm{H}}}{\mathrm{4}}\hspace{0.17em}\mathit{n}\mathrm{\left(}\mathrm{H}\mathrm{\right)}$), as this rate is the maximum formation rate that would occur if one molecule was formed every two collisions. We also notice that standard models tend to show a strong variability of the formation efficiency (and consequently of the H abundance) with a gap around A_V = 2 with very low efficiency, where neither the ER mechanism (the gas is not warm enough) nor the LH mechanism (the grains are too warm) are efficient. This variability disappears in the coupled models as fluctuations make the LH mechanism much less dependent on the grain equilibrium temperature (see Sect. 4.2), and the efficiency stays of the order of 10% (for amorphous carbon surfaces) across the cloud.

The formation efficiency at the center of the cloud is decreased by almost one order of magnitude when including the fluctuations. This unexpected effect is detailed in Fig. 19 for amorphous carbon surfaces: the formation parameter R_f for the LH mechanism is shown for the standard model (black dashed line) and the coupled model (black solid line), and two inset plots show the effect of fluctuations on a 1 nm grain at two positions (A_V = 0.5 and A_V = 5, marked, respectively, by the blue and green vertical arrows on the main plot). Each inset contains a plot similar to Fig. 5, displaying the temperature PDF (blue line), the equilibrium temperature (red line), the equilibrium efficiency curve (dashed green line), and the effective conditional efficiency curve when taking fluctuations into account (solid green line).

Close to the edge of the cloud (at A_V = 0.5), the situation is similar to the discussion of Fig. 5 (see Sect. 4.2). The main effect occurs because the temperature PDF has an important part inside the efficiency domain, while the equilibrium temperature is out of this domain. It results in increased efficiency. The dangers of using the equilibrium temperature are strikingly emphasized on this figure: when allowing fluctuations, a very small fraction of high temperature grains is sufficient to radiate away all the absorbed energy and most of the grains can be cold, while the grains must all be warm enough to balance the absorbed energy if we force them all to have a single temperature.

Deeper in the cloud (at A_V = 5), the situation is different: the equilibrium temperature now falls close to the actual peak of the PDF and inside the efficiency domain, so that the actual position and spread of the PDF does not change the result much. The dominant effect is now that the conditional efficiency with fluctuations is strongly reduced compared to the equilibrium efficiency. As discussed in Sect. 4.2, this reduction comes from the competition between the fluctuation timescale and the adsorption timescale. At this optical depth, the grains are mainly heated by IR photons emitted by the hot dust at the edge, and for 1 nm grains, the IR photons are sufficient to bring the grains temporarily out of the efficiency domain (above 18 K here) where they lose their surface population. On the other hand, adsorption of H atoms are extremely rare, as the gas is almost completely molecular. As a result, the grain rarely manages to have two atoms on its surface, and the formation efficiency is thus reduced.

The effect described here concern very small grains, and one might question the presence of such very small grains and PAHs deep in the cloud core (the grain population is independent of position in our model). However, secondary UV photons, which were not considered for grain heating in these models, could induce a similar effect for much bigger grains, and further reduce the formation efficiency at high optical depth inside the cloud. This mechanism tends to increase the residual atomic H fraction in the core and could be relevant in explaining the slightly higher than expected atomic H abundances deduced from observations of HI self-absorption in dark clouds (e.g., Li & Goldsmith 2003).

No detectable effects on intensities were found in these models but the inclusion of secondary UV photons may induce stronger effects.

5.2. Grid of models

A grid of constant density PDR models was run to check the effect on observable intensities over a wider parameter range. We explore gas densities from n_H = 10² cm^-3 to n_H = 3 × 10⁵ cm^-3 and scaling factors on the standard ISRF (Mathis et al. 1983) from 1 to 10⁴ (as explained previously, the scaling factor only affects the UV part of the radiation field).

The effects of the fluctuations are mostly found on tracers of the PDR edges. The formation efficiency at the edge of the cloud is shown as a color map in Fig. 20, which compares the standard models with the full models including fluctuations. In the upper part of the map, the ER formation is decreased, and the LH mechanism remains negligible, resulting in an overall (small) decrease in efficiency. In the lower part, the LH mechanism is not negligible anymore, and its strong increase due to fluctuations results in an overall increase in efficiency. The region around G₀ = 10 and n_H = 10³ − 10⁴cm^-3 is strikingly affected. In this region, the ER mechanism stops being efficient due to the low gas temperature at the edge. In standard models, the equilibrium temperature of the dust at the edge is still too high for the LH mechanism to work, giving a very low total formation efficiency. When taking into account fluctuations, the LH mechanism works efficiently, assuring a high formation efficiency. The effect is further amplified as gas heating by H₂ formation is dominant in this region. The efficient LH formation heats up the gas, thus keeping ER formation efficient as well.

Figure 21 shows the modification of the H₂ lower rotational intensities as the ratio between the result of the full computation, which includes dust temperature fluctuations, and the result of a standard PDR model. Higher rotational lines are affected very similarly to the S(3) line. Figure 22 shows the effect on the ro-vibrationnal line (1−0) S(1) in the same way. As expected from our discussion of Fig. 20, we see increased intensities for low radiation fields (roughly G₀< 200), especially in regions where the gas temperature is highly sensitive on H₂ formation heating and where LH formation can further trigger an amplification of ER formation. For higher radiation fields, the intensities are slightly decreased. The effect remains small, the increase for low radiation field models is at most a factor of 2.7 (on the S(1) line), and the decrease on the high radiation field models is at most 30%.

Another species affected by the more efficient formation of H₂ at the edge is CH⁺. It is formed by the reaction of C⁺ + H₂, which is allowed by the high gas temperature at the edge of the PDR and by the internal energy of the H₂ molecule, as described in Agúndez et al. (2010) (pumped by UV photons and by its formation process). A higher H₂ formation rate at the edge induces a higher H₂ fraction before the transition, which leads to a higher CH⁺ formation rate. Gas heating by H₂ formation can further enhance its formation. Figure 23 shows the effect on CH⁺ intensities. We again observe increased emission in the domain where the H₂ formation efficiency is enhanced, due to both higher abundances of CH⁺ and higher gas temperatures. The domain of strong enhancement is quite similar to the one for H₂ rotational lines, corresponding to the domain where gas heating is dominated by H₂ formation heating. This temperature effect is visible on all high excitation lines emitted before the H / H₂ transition (e.g., the most excited lines of HCO⁺, CS, HCN, ...).

6. Conclusion and perspectives

We presented a method to compute H₂ formation on dust grains taking fully the fluctuating temperature of the grain for both the ER and the LH mechanism into account.

As a side result of this work, we revisited the computation of dust temperature fluctuations. We confirm that the continuous cooling approximation (e.g., used in DustEM) is sufficient for an accurate computation of the mid- and near-infrared spectrum. However, the detailed treatment of surface processes requires a better estimation of the low temperature part of the dust temperature PDF, which is done in this paper by including the discrete nature of photon emission.

The coupled fluctuations problem was solved exactly for both the ER and the LH mechanisms. A numerically efficient way of computing the ER formation rate was found. The effect on the ER formation rate is limited to a factor of 2, and the effect on full cloud simulations is limited. The standard rate equation treatment of this ER mechanism is thus a good treatment for the purpose of complete PDR simulations.

With the exact solution for the LH mechanism being numerically intractable for big grains, we constructed an approximation based on simple physical arguments to compute the fluctuation effect on the LH mechanism and checked for small grains that the approximation yields accurate results. The effect is strong as fluctuating grains spend a significant part of their time in the low temperature range, in which formation is efficient, even when their average temperature is out of the rate equation efficiency domain. This results in a strong increase in LH formation efficiency in unshielded environments where the average dust temperature is high and only fluctuations can allow LH formation.

The LH mechanism is found to become of comparable importance with the ER mechanism at the edge of PDRs, except for strong radiation fields. The resulting increase of H₂ formation rate at the edge induces a shift in the H / H₂ transition toward the edge of the PDR. The resulting effect on observable lines emitted in this region remains limited and does not reach one order of magnitude in the strongest cases.

On the contrary, fluctuations of the smallest grains caused by IR photons reduce the formation efficiency by almost one order of magnitude in the core of the clouds, as collisions with H atoms in a completely molecular gas are much rarer than absorption of IR photons heating temporarily the grain above 18 K. No significant effects are found on intensities, but including secondary UV photons could induce stronger effects.

Overall, the variations induced in observational intensities are comparable to what can be expected from many other poorly known parameters. The interpretation of observations with standard PDR codes, thus, need not use the detailed treatment analyzed here and a simple rate coefficient approximation can be used most of the time. However, this cannot be blindly extended to any arbitrary process and must be checked in each specific case. In particular, the existence of temperature fluctuations of the smallest grains deep inside the cloud should be taken into account for other surface chemical processes. Possible consequences on the abundance of light organic molecules are currently being investigated.

---

Acknowledgments

This work was funded by grant ANR-09-BLAN-0231-01 from the French Agence Nationale de la Recherche as part of the SCHISM project. We thank Laurent Verstraete for fruitful discussions on dust properties, and Evelyne Roueff for many discussions. We thank our referee for his useful comments.

References

Appendix A: Effect of H₂ formation energy on the equilibrium rate equation calculation

We first consider the heating of the grain by the Eley-Rideal mechanism. Each formation reaction gives 1 eV to the grain.

The liberation of energy left in the grain by each formation reaction creates an additional coupling between the surface population and the grain temperature. The method described in Sect. 2.2.2, which was based on a one-way coupling between the two variables, is thus not possible here.

To evaluate this effect, we use the standard equilibrium temperature equation and chemical rate equation with a coupling term representing the heating of the grain by the formation reaction. We assume that each formation reaction gives an energy E_fH2 = 1 eV to the grain, which is instantly converted into thermal energy.

For ER formation on chemisorption sites, we thus have the system of equations: $\begin{array}{ccc}{\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}{\mathit{T}}_{\mathrm{eq}}\mathrm{\right)}& \mathrm{=}& {\mathit{N}}_{\mathrm{s}}\hspace{0.17em}\frac{\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{s}\mathrm{\left(}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}{\mathit{k}}_{\mathrm{des}}\mathrm{\left(}{\mathit{T}}_{\mathrm{eq}}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{coll}}}}\mathit{,}\\ & & {\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{U}\hspace{0.17em}{\mathit{P}}_{\mathrm{abs}}\mathrm{\left(}\mathit{U}\mathrm{\right)}\mathrm{+}{\mathit{E}}_{{\mathit{f}}_{{\mathrm{H}}_{\mathrm{2}}}}\hspace{0.17em}{\mathit{k}}_{\mathrm{coll}}\hspace{0.17em}\frac{{\mathit{n}}_{\mathrm{eq}}\mathrm{\left(}{\mathit{T}}_{\mathrm{eq}}\mathrm{\right)}}{\mathit{N}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{E}}_{\mathrm{grain}}}\mathrm{d}\mathit{U}\hspace{0.17em}{\mathit{P}}_{\mathrm{em}}\mathrm{\left(}\mathit{U,}{\mathit{T}}_{\mathrm{eq}}\mathrm{\right)}\mathit{,}\end{array}$which can be easily solved. However, we find that the effect of this average heating of the grain by formation reactions is completely negligible in the range of realistic values of parameters, resulting in a change of the grain temperature by less than 1%.

Independently of the detail of the formation mechanism, we can put an upper limit on the heating term caused by H₂ formation by assuming that all H atoms hitting the grain are turned into H₂. In this case, the formation rate is simply $\frac{\mathrm{1}}{\mathrm{2}}{\mathit{k}}_{\mathrm{coll}}$, and the heating power is $\frac{\mathrm{1}}{\mathrm{2}}{\mathit{k}}_{\mathrm{coll}}{\mathit{E}}_{{\mathrm{f}}_{{\mathrm{H}}_{\mathrm{2}}}}$. We can then compare this term to the absorbed power coming from photons.

The result is shown here in Fig. A.1 in an extreme case with a low radiation field (G₀ = 1) and a high collision rate (atomic gas with n = 10³ cm^-3 and T = 350 K). Despite assuming a total formation efficiency in a low radiation field, high collision rate environment, we can see that the heating term due to H₂ formation becomes barely significant only for the smallest grain sizes. We can thus safely neglect the grain heating coming from H₂ formation (through the ER mechanism, and all the more for the LH mechanism as the E_fH2 is probably much lower as physisorbed atoms are much more weakly bound) in most usual conditions. Moreover, under a radiation field with G₀ = 1, the gas would not be as warm as we assumed here and the ER efficiency would be very low.

Appendix B: Eigenvalues of the integral operator for the thermal energy PDF

The operator ℒ defined in equation 8 can be rewritten as $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{A}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}& & \end{array}$with $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{E,}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$where $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathit{,T}\mathrm{(}\mathit{E}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{·}& & \end{array}$The existence of a solution to our problem is equivalent to the existence of a positive eigenfunction associated with the eigenvalue 1.

This operator is continuous, compact, and strongly positive. The Krein-Rutman theorem (see Du 2006, Chap. 1) tells us that its spectral radius is a simple eigenvalue associated with a positive eigenfunction and that this eigenfunction is the only positive eigenfunction.

We first show that 1 is an eigenvalue of the adjoint operator ${\mathrm{ℒ}}^{\mathrm{\ast }}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{A}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,E}\mathrm{\right)}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}$ and that it is associated with a positive eigenfunction. We can easily find such a positive eigenfunction with the eigenvalue 1 for this adjoint operator: $\begin{array}{ccc}{\mathrm{ℒ}}^{\mathrm{\ast }}\mathrm{[}\mathrm{P}\mathrm{\left]}\mathrm{\right(}\mathrm{E}\mathrm{)}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{A}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathit{,E}\mathrm{\right)}\hspace{0.17em}\mathit{P}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\\ \hspace{0.17em}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{\mathit{E}}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\frac{{\mathit{R}}_{\mathrm{em}}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}^{\mathrm{\prime }}\mathit{,T}\mathrm{(}\mathit{E}\mathrm{\left)}\mathrm{\right)}}{\mathit{P}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}}\mathit{P}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\\ \hspace{0.17em}& \mathrm{+}& {\mathrm{\int }}_{\mathit{E}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}{\mathit{E}}^{\mathrm{\prime }}\frac{{\mathit{R}}_{\mathrm{abs}}\mathrm{(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{-}\mathit{E}\mathrm{)}}{\mathit{P}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}}\mathit{P}\mathrm{\left(}{\mathit{E}}^{\mathrm{\prime }}\mathrm{\right)}\\ \hspace{0.17em}& \mathrm{=}& \mathit{P}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathit{,}\end{array}$and P(E) is a positive function. Thus, 1 is the spectral radius of ℒ^∗ and hence of ℒ. This proves the existence and the unicity of the solution to our problem. Moreover, it proves that all the other eigenvalues λ are such that | λ | ≤ 1 and thus Re(λ) < 1. We can therefore iterate the application of the operator to an initial function and converge toward the solution, assuming that the initial function is not orthogonal to the solution. Starting with a positive initial function is sufficient.

Appendix C: Verification of the approximations

In this appendix, we evaluate the accuracy of our approximations by comparing their results to those of the exact methods presented in Sect. 2.2.

Appendix C.1: LH mechanism

For the LH mechanism, the exact method (cf. Sect. 2.2.3) is only tractable for a limited range of grain sizes. We compare the approximation (cf. Sect. 3.1) with it on this range and verify that the approximation joins the equilibrium result for bigger grains.

Figure C.1 shows the resulting curves for ⟨n|T⟩ for different grain sizes and compares the approximation to the exact result. The binding and barrier energies are taken for ice surfaces in Table 1. We see a very good match for small grains in the low temperature domain. Large discrepancies appear at high temperature but do not matter for the average formation rate on the grain as the contribution of this regime is negligible. For bigger grains, the approximation of the low temperature regime starts to be less accurate. However, we know that total formation rate is mostly influenced by the smallest grain sizes due to the −3.5 exponent of the size distribution.

The resulting formation rate per grain is shown in Fig. C.2 as a function of grain size for 5 models spanning the parameter domain. We again assumed here an ice surface. The match is very good on the domain on which the exact computation could be performed, and the approximation smoothly joins with the equilibrium result neglecting fluctuations for big grains. The relative errors on the final formation rate integrated over the size distribution are given for each model in the first part of Table C.1. The results are given for two size distributions. Both are MRN-like power-law distributions with exponent −3.5, but one starts at 5 Å and the other at 1 nm. Both go up to 0.3 μm. We find extremely accurate results for the approximation in most case with a relative error at most around 5%.

When using the energy values corresponding to an amorphous carbon surface (see Table 1), the same levels of accuracy are achieved. The formation rate as a function of grain size is shown in Fig. C.3. When integrated over the size distribution, the formation rate obtained from the approximation is again within 5% from exact result (see second part of Table C.1).

Appendix C.2: ER mechanism

We now perform the same verification for the ER mechanism. The approximation of Sect. 3.2 is compared to the exact method of Sect. 2.2.2. Figure C.4 shows the formation rate as a function of grain size. The qualitative effect is captured by the approximation, but large discrepancies appear locally for some sizes, due to a shifted transition from a low efficiency regime to a high efficiency one. However, the results when integrating over the dust size distribution remain reasonably accurate, as shown in Table C.2. We recall that this approximation is not used in this paper as the exact method is numerically efficient.

All Tables

All Figures