© ESO, 2016

1. Introduction

Photodissociation regions (PDRs) consist of predominantly neutral gas and dust illuminated by far-ultraviolet (FUV) radiation (6 <hν< 13.6 eV). Dense PDRs are found in the vicinity of star-forming regions, since FUV photons usually arise from massive stars creating HII regions. Ultraviolet (UV) photons impinge on clouds of gas and dust and play an important role in the heating and chemistry of these irradiated regions; while the gas is heated to relatively high temperatures, the ionisation and, mainly, the photodissociation of different species are produced as the UV radiation penetrates into the region (Hollenbach et al. 1997). This, therefore, produces significant differences in the chemical composition of the cloud depending on the visual extinction.

To properly model the chemical composition of PDRs, it is key to take the role of dust in the chemistry of these regions into account. In environments powered by high radiation fields, the dust grains are mainly bare, since no ice mantles can form on their surfaces owing to radiation (Meijerink et al. 2012). Dust grains, when not covered by ice, provide an ideal place for chemical reactions to occur that directly enrich the gas phase. For higher extinctions, however, dust grains are mainly covered by ice mantles, which usually sublimate into the gas phase by star formation activities (Viti et al. 2004), enhancing the chemical composition of the gas phase as well.

The first numerical models of PDRs were developed by Hollenbach et al. (1971), Glassgold & Langer (1975), and Black & Dalgarno (1977), who assume a steady state to simulate the transitions from H to H₂ and from C⁺ to CO. Later, new models, such as those from van Dishoeck & Black (1988), Le Bourlot et al. (1993), Störzer et al. (1996), and Spaans (1996) were developed to focus on the chemical and thermal structure of clouds subject to an incident flux of FUV radiation. More recently, models considering time-dependent chemical networks (e.g. Bertoldi & Draine 1996), turbulence (e.g. Röllig et al. 2002; Bell et al. 2005) and 3D (Bisbas et al. 2012) have also been developed. All these models are mainly based on gas-phase chemistry. In fact, only a few PDR codes consider reactions on dust surfaces (apart from H₂ formation). These include the Meudon PDR code from Le Petit et al. (2006) and the PDR code from Hollenbach et al. (2009).

Many recent observations (e.g. Berné et al. 2007; Sellgren et al. 2010; Guzmán et al. 2013; Cuadrado et al. 2015) show the chemical richness of dense PDRs. This richness makes it evident that there is a need to consider a larger number of solid species in the current codes to understand the role of dust in the origin of the chemical complexity of PDRs, since surface chemical reactions can dramatically alter the gas-phase composition of these regions. It is, therefore, crucial to include a detailed treatment of the chemistry occurring on dust grains to really understand the link between dust chemistry and gas chemistry. This treatment should not only consider the most relevant physical processes for grains at large visual extinctions, but also study the presence or absence of ice mantles on their surfaces.

In this paper, we have significantly improved the chemical network of the Meijerink code (Meijerink & Spaans 2005), by considering 3050 new gas-phase chemical reactions and, mainly, by implementing dust grain chemistry to determine the different compositions in PDRs (gas/dust/ice). These chemical treatments include adsorption, thermal and chemical desorption, two-body reactions, photo processes, and cosmic-ray processes on dust grains, as well as the presence of some complex organic solid species, such as H₃CO and CH₃OH. In Sect. 2, we explain the new chemical treatment, especially for dust grains depending on the presence or absence of ice mantles, and we present the results in Sect. 3. In Sect. 4, we analyse the formation routes for ices of CO₂ and H₂O depending on the environmental conditions. We also study the impact of the presence of dust grains and of their type of substrate (bare or icy) on the chemical composition of regions powered by UV photons. In addition, we compare our results with observations from the Horsehead PDR and the Orion Bar PDR in Sect. 5. Finally, we summarise the main conclusions in Sect. 6.

2. The numerical code

2.1. Gas chemistry

We consider in our steady-state PDR code 7503 gas-phase chemical reactions from the Kinetic Database for Astrochemistry (KIDA; Wakelam et al. 2015)¹, including bimolecular reactions (A + B → C + D), charge-exchange reactions (A⁺ + B → A + B⁺), radiative associations (A + B → AB + photon), associative detachment (A⁻ + B → AB + e⁻), dissociative recombination (A⁺ + e⁻→ C + D), neutralisation reactions (A⁺ + B⁻→ A + B), ion-neutral reactions (A⁺ + B → C⁺ + D), ionisation or dissociation of neutral species by UV photons, and ionisation or dissociation of species by direct collision with cosmic-ray particles or by secondary UV photons following H₂ excitation. The initial gas-phase abundances (A_i) for the different elements that we consider (Jenkins 2004; Asplund et al. 2005, and Neufeld & Wolfire 2009) are listed in Table 1.

The thermal balance of PDRs is determined by different heating and cooling processes. As heating mechanisms, we consider the photoelectric effect on grains, carbon ionisation heating, H₂ photodissociation heating by UV photons, H₂ collisional de-excitation heating, gas-grain collisional heating, gas-grain viscous heating², and cosmic-ray heating. As cooling mechanisms, we consider fine-structure line cooling (being [CII] at 158 μm and [OI] at 63 μm and at 146 μm the most prominent cooling lines), metastable-line cooling (including lines of C, C⁺, Si, Si⁺, O, O⁺, S, S⁺, Fe, and Fe⁺), recombination cooling, and molecular cooling by H₂, CO, and H₂O (see Meijerink & Spaans 2005 for more details of each process).

2.2. Dust chemistry

While the gas-phase chemical network that we consider in our code is taken from KIDA 2014, we derived the surface chemical network from laboratory experiments (e.g. Dulieu et al. 2013; Minissale et al. 2015, 2016). The solid species that we consider in the code to model the dust chemistry are listed in Table 2. We included 117 chemical reactions occurring on grain surfaces covering different processes (adsorption, thermal and non-thermal desorption, two-body reactions, photo processes, and cosmic-ray processes). These processes are detailed in the following subsections.

2.2.1. Adsorption onto dust grains

Gas-phase species can be adsorbed on the dust grain surface. This adsorption is determined by the dynamics of the accretion processes acting on the dust in a determined region. In dense regions of the interstellar medium (ISM), accretion is favoured by the increased collision rates between gas-phase species and grains. Accretion efficiently depletes grains with a radii lower than 0.001 μm on a timescale of ≲10 Myr in solar-metallicity molecular clouds with densities n ~ 10³ cm^-3 (Hirashita 2000, 2012). The adsorption rate, R_ad (cm^-3 s^-1), of the species i is determined by the total cross section from dust³, σ_dn_d (cm^-1), the thermal velocity of the species i (v_thi), and the sticking coefficient S(T_g, T_d). The adsorption rate can be written as $\begin{array}{ccc}{\mathit{R}}_{\mathrm{ad}}\mathrm{=}{\mathit{n}}_{\mathit{i}}{\mathit{\sigma }}_{\mathrm{d}}{\mathit{n}}_{\mathrm{d}}{\mathrm{v}}_{{\mathrm{th}}_{\mathit{i}}}\mathit{S}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathit{i}}{\mathit{k}}_{\mathit{ads}}\mathit{S}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\mathit{,}& & \end{array}$(1)where k_ads is the adsorption rate coefficient, n_i is the number density of adsorbing species i, and the thermal velocity is written as $\begin{array}{ccc}{\mathrm{v}}_{{\mathrm{th}}_{\mathit{i}}}\mathrm{=}\sqrt{\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{g}}}{\mathit{\pi }{\mathit{m}}_{\mathit{i}}}}\mathit{,}& & \end{array}$(2)where m_i is the mass in grams of species i and k_B the Boltzmann constant. The sticking coefficient for all the species is given by $\begin{array}{ccc}\mathit{S}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\mathrm{=}{\left(\mathrm{1}\mathrm{+}\mathrm{0.4}\sqrt{\frac{{\mathit{T}}_{\mathrm{g}}\mathrm{+}{\mathit{T}}_{\mathrm{d}}}{\mathrm{100}}}\mathrm{+}\mathrm{0.2}\frac{{\mathit{T}}_{\mathrm{g}}}{\mathrm{100}}\mathrm{+}\mathrm{0.08}{\left(\frac{{\mathit{T}}_{\mathrm{g}}}{\mathrm{100}}\right)}^{\mathrm{2}}\right)}^{-1}& & \end{array}$(3)(Hollenbach & McKee 1979), where T_g and T_d are the gas and dust temperatures, respectively. In this paper, we use a mean grain cross section σ_MRN = ⟨ σ_dn_d/n_H ⟩ _MRN = 10^-21 cm² computed assuming a MRN grain size distribution (Mathis et al. 1977), with grain radius extending from ~50 Å to ~0.25 μm. In Table A.1, we list the adsorption reactions considered.

2.2.2. Thermal desorption

Once species are depleted on dust grains, they can evaporate back into the gas. This process depends on the dust temperature and on whether the surfaces of the grains are covered by ice (icy grains) or not (bare grains), since the binding energies in both cases are different. We calculate the fraction of the dust that is bare, f_bare, and icy, f_ice, taking into account the density of sites, n_dn_sites (cm^-3), where species are locked on the dust. The density of sites can be written as $\begin{array}{ccc}{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}\mathrm{=}{\mathit{n}}_{\mathrm{d}}\frac{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{d}}^{\mathrm{2}}}{\mathrm{(}{\mathit{a}}_{\mathit{pp}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{=}{\mathit{n}}_{\mathrm{d}}\frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{d}}}{\mathrm{(}{\mathit{a}}_{\mathit{pp}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\simeq }\mathrm{4.44}\mathrm{\times }{\mathrm{10}}^{-6}{\mathit{n}}_{\mathrm{H}}\mathit{,}& & \end{array}$(4)where r is the radius of dust and a_pp is the distance between two sites that we assume to be 3 Å. This is the typical size between sites that should be considered to obtain a typical site density of ~10¹⁵ sites/cm². A full monolayer is formed when all the possible sites on a grain surface are occupied by an atom or molecule. It means that for solid species, as abundances are higher than 4.44 × 10^-6, more than one monolayer is reached. When the grain surface is covered by less than one monolayer, we calculate f_ice as $\begin{array}{ccc}{\mathit{f}}_{\mathrm{ice}}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{J}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{\right)}}}{{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}}\mathit{,}& & \end{array}$(5)where n _J(H2O) is the number density of solid H ₂ O. If the grain is covered by more than one layer of water ice (i.e. when n_J(H2O)>n_dn_sites), f_ice = 1. The expression for f_bare can be obtained by $\begin{array}{ccc}{\mathit{f}}_{\mathrm{bare}}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{ice}}\mathit{.}& & \end{array}$(6)As previously mentioned, the fraction of bare or icy dust has important consequences for the binding energies of the species. We, therefore, consider two types of binding energies: on bare dust (E_b) and on icy surfaces (E_i). The desorption rate, R_des (cm^-3 s^-1), can be written as $\begin{array}{ccc}{\mathit{R}}_{\mathrm{des}}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}{\mathit{n}}_{\mathit{i}}\left[{\mathit{f}}_{\mathrm{bare}}\exp \left(\frac{\mathrm{-}{\mathit{E}}_{\mathrm{b}}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}{\mathit{f}}_{\mathrm{ice}}\exp \left(\frac{\mathrm{-}{\mathit{E}}_{\mathrm{i}}}{{\mathit{T}}_{\mathrm{d}}}\right)\right]\mathrm{·}& & \end{array}$(7)ν₀ is the oscillation frequency that is determined by $\begin{array}{ccc}{\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\sqrt{\frac{\mathrm{2}{\mathit{N}}_{\mathrm{s}}\mathit{E}}{{\mathit{\pi }}^{\mathrm{2}}\mathit{m}}}\mathit{,}& & \end{array}$(8)where N_s is the surface number density of sites on the grain, m is the mass of the species, and E is the energy of the site where the species is bound. ν₀ is typically ~10¹² s^-1 for physisorbed species. In Tables A.2 and A.3, we list the desorption reactions considered and the binding energies for each species, respectively.

2.2.3. Two-body reactions on grain surfaces

Once gas species are adsorbed onto dust grains, they can also move around the grain surfaces. There are two different types of interaction between the species and the surface: physisorption and chemisorption. The physisorption is a weak interaction due to Van der Waals forces between the adsorbed atom and the surface (dipole-dipole interaction). The typical depth of wells associated with physisorption are of the order of 0.01–0.2 eV (Vidali et al. 1991). The chemisorption is adsorption in which the forces involved are similar to valence forces, and the interaction potential depends not only on the distance from the surface, but also on the position on the surface. The typical binding energies for chemisorption are of the order of ~1 eV (Barlow & Silk 1976; Zangwill 1988).

The surface of a dust grain is mainly irregular, with the presence of peaks and valleys. The valleys represent the physisorbed (highest) and chemisorbed (deepest) wells and they are separated by saddle points (diffusion barriers). The mobility of adsorbed species is associated with transfer across these barriers, which can occur through thermal diffusion. Dulieu et al. (2013) and Collings et al. (2003) find experimentally that diffusion occurs with a barrier of 67% and 40%, resepctively, of the binding energy. Theoretical results show percentages of 90% (Barzel & Biham 2007) and 30% (Karssemeijer & Cuppen 2014a). In this paper, we assume that diffusion occurs with a barrier of 2/3 of the binding energy. The adsorbed species i and j have different probabilities of mobility depending on whether the grain is bare (P_bare) or icy (P_ice) as follows: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{bare}}\mathrm{=}{\mathit{f}}_{\mathrm{bare}}\left[\exp \left(\frac{\mathrm{-}\mathrm{2}{\mathit{E}}_{\mathrm{b}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathrm{3}{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}\exp \left(\frac{\mathrm{-}\mathrm{2}{\mathit{E}}_{\mathrm{b}}\mathrm{\left(}\mathit{j}\mathrm{\right)}}{\mathrm{3}{\mathit{T}}_{\mathrm{d}}}\right)\right]\mathit{,}& & \end{array}$(9)and $\begin{array}{ccc}{\mathit{P}}_{\mathrm{ice}}\mathrm{=}{\mathit{f}}_{\mathrm{ice}}\left[\exp \left(\frac{\mathrm{-}\mathrm{2}{\mathit{E}}_{\mathrm{i}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathrm{3}{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}\exp \left(\frac{\mathrm{-}\mathrm{2}{\mathit{E}}_{\mathrm{i}}\mathrm{\left(}\mathit{j}\mathrm{\right)}}{\mathrm{3}{\mathit{T}}_{\mathrm{d}}}\right)\right]\mathrm{·}& & \end{array}$(10)The two-body reaction rate, R_2body (cm^-3 s^-1), is therefore given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{2}\mathrm{body}}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}\left(\frac{{\mathit{n}}_{\mathit{i}}{\mathit{n}}_{\mathit{j}}}{{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}}\right){\mathit{k}}_{\mathrm{2}\mathrm{body}}{\mathit{P}}_{\mathrm{react}}\mathit{,}& & \end{array}$(11)where ν₀ is the oscillation factor, k_2body is the two-body rate coefficient $\begin{array}{ccc}{\mathit{k}}_{\mathrm{2}\mathrm{body}}\mathrm{=}{\mathit{P}}_{\mathrm{bare}}{\mathit{\delta }}_{\mathrm{bare}}\mathrm{+}{\mathit{P}}_{\mathrm{ice}}{\mathit{\delta }}_{\mathrm{ice}}\mathit{,}& & \end{array}$(12)and δ_bare and δ_ice are the theoretical probabilities of desorption upon reaction for bare and icy⁴ substrates, respectively (Minissale et al. 2016). To determine the probability, P_react, of overcoming a reaction barrier with energy ϵ (K), we need to consider thermal diffusion, P_react(therm), and tunneling, P_react(tunn) as follows: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{react}}\mathrm{\left(}\mathit{therm}\mathrm{\right)}\mathrm{=}\exp \left(\frac{\mathrm{-}\mathit{ϵ}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathit{,}& & \end{array}$(13)and $\begin{array}{ccc}{\mathit{P}}_{\mathrm{react}}\mathrm{\left(}\mathit{tunn}\mathrm{\right)}\mathrm{=}\exp \left(\mathrm{-}\mathit{a}\sqrt{\frac{\mathrm{2}{\mathit{m}}_{\mathrm{red}}{\mathit{k}}_{\mathrm{B}}\mathit{ϵ}}{{ħ}^{\mathrm{2}}}}\right)\mathit{,}& & \end{array}$(14)where a is the width of the barrier of 1 Å (for a square barrier; Tielens & Hagen 1982; Hasegawa et al. 1992), ħ is the Planck constant divided by 2π, and m_red is the reduced mass of the reaction between two species i and j, m_red = (m_i × m_j)/(m_i + m_j). The probability P_react can be written as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{react}}\mathrm{=}\frac{{\mathit{P}}_{\max }}{{\mathit{P}}_{\max }\mathrm{+}{\mathit{P}}_{\mathrm{bare}}\mathrm{+}{\mathit{P}}_{\mathrm{ice}}}\mathit{,}& & \end{array}$(15)with P_max the maximum between P_react(tunn) and P_react(therm). Table A.4 lists the considered reactions occurring on grain surfaces and the parameters δ_bare and δ_ice. We also consider reactions on dust grains between physisorbed and chemisorbed hydrogen to form molecular hydrogen. See Appendix B for more details.

2.2.4. Photo processes on dust grains

In the outskirts of molecular clouds and in the vicinity of high-mass stars, FUV photons can dominate the chemistry. These photons usually impinge on dust grains located in the surroundings of the forming stars and can lead to the photodesorption and/or photodissociation of the species adsorbed onto dust grains. Since recent results (Muñoz-Caro et al. 2010; Fayolle et al. 2011; Chen et al. 2014) indicate that photons adsorbed deeper than a few monolayers have no effect on the photodesorption because they cannot transfer their energy to the uppermost monolayers, in our code we consider that incident photons can only interact the first two layers of ice and produce photodesorption.

Since photoreactions scale linearly with the local radiation flux (erg cm^-2 s^-1) and the radiation field strength is a function of extinction (ξ_iA_V), where ξ_i is the extinction factor for the relevant species, the photo-process reaction rate, R_photo (cm^-3 s^-1), can be written as $\begin{array}{ccc}{\mathit{k}}_{\mathrm{photo}}& \mathrm{=}& {\mathit{\alpha }}_{\mathit{i}}{\mathrm{e}}^{\mathrm{-}{\mathit{\xi }}_{\mathit{i}}{\mathit{A}}_{\mathit{V}}}\mathit{,}\\ {\mathit{R}}_{\mathrm{photo}}& \mathrm{=}& {\mathit{n}}_{\mathit{i}}{\mathit{f}}_{\mathrm{ss}}{\mathit{k}}_{\mathrm{photo}}{\mathit{F}}_{\mathrm{UV}}\mathit{.}\end{array}$The parameter n_i is the number density of the photodissociated species, k_photo is the photo-process rate coefficient, α_i is the unattenuated rate coefficient, f_ss is the self-shielding factor, and F_UV is the UV flux in units of 1.71G₀ (G₀ = 1 gives F_UV = 0.58). The factor 1.71 arises from the conversion of the often used Draine field (Draine 1978) to the Habing field for the FUV intensity. We consider self-shielding for H₂ and CO molecules (van Dishoeck & Black 1988). For the rest of the species, we assumed f_ss = 1. Table A.5 lists the reactions occurred on grain surfaces due to photon impacts and the parameters α_i and ξ_i. In particular, we consider in our code direct photodesorption of CO, H₂O, and H₂CO. We do not include direct photodesorption of CH₃OH, given that recent laboratory results (Bertin et al. 2016) conclude that this mechanism is not very efficient to release methanol into the gas phase.

2.2.5. Cosmic-ray processes on dust grains

Cosmic rays impacting on dust grains can provide non-thermal energy to desorb molecules frozen on grain surfaces. In addition, given that cosmic rays have larger penetrating power than UV photons and X-rays, they can have a greater impact on the chemistry in well-shielded regions. In particular, cosmic rays can heat dust grains partially or completely leading to explosive desorption (d’Hendecourt et al. 1982; Léger et al. 1985; Ivlev et al. 2015). The grain temperature increase due to cosmic rays depends on their flux, the projected area of the grain, and the energy lost by a cosmic ray as it passes through the grain. Since some electrons excited by cosmic rays have large energies and escape from the grains leading to a reduction of the effective heating, cosmic rays need to pass through a sufficiently long path in the grains to deposit enough energy to produce desorption. Desorption is not, however, the only consequence from the interaction between cosmic rays and dust grains. In fact, recent experiments with interstellar ices (Dartois et al. 2015) show that cosmic irradiation can also alter the ice mantle state on dust grains.

The cosmic-ray reaction rates on grain surfaces, R_CR (cm^-3 s^-1), are assumed to be the same as the rates for gas-phase reactions. They are determined by $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{CR}}& \mathrm{=}& {\mathit{z}}_{\mathit{i}}{\mathit{\zeta }}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{,}\\ {\mathit{R}}_{\mathrm{CR}}& \mathrm{=}& {\mathit{n}}_{\mathit{i}}{\mathit{\kappa }}_{\mathrm{CR}}\mathit{,}\end{array}$where κ_CR (s^-1) is the cosmic-ray rate coefficient, n_i is the number density of the photodissociated species, z_i is the cosmic-ray ionisation rate factor, which depends on the ionising element (see KIDA database and Table A.6), and ζ_H2 is the cosmic-ray ionisation rate per H₂ molecule (ζ_H2 = 5 × 10^-17 s^-1; Indriolo et al. 2007; Hocuk & Spaans 2011). The cosmic-ray reactions on dust grains considered in this code are listed in Table A.6.

2.3. Dust temperature

The dust temperature is a key parameter in the thermal balance calculation, since it influences the gas temperature, through heating and cooling rates along with chemical reaction rates. In addition, the dust temperature value is also crucial for the formation of ice mantles on grain surfaces. There are several expressions for dust temperature in the literature, such as those derived by Werner & Salpeter (1969), Hollenbach et al. (1991), Zucconi et al. (2001), and Garrod & Pauli (2011). See also Hocuk et al. (in prep.), who show a fit of dust temperature observations as a function of visual extinction, using the different analytical expressions of T_d previously mentioned.

In our code, we considered the most recent dust temperature expression derived by Garrod & Pauli (2011), but with an adaptation (R. T. Garrod 2015, priv. comm.), that includes a dependence on radiation field, since the original expression only depends on the visual extinction (A_V). The final expression is $\begin{array}{ccc}{\mathit{T}}_{\mathrm{d}}& \mathrm{=}& \mathrm{18.67}\mathrm{-}\mathrm{1.637}\mathrm{(}\begin{array}{c}\\ \end{array}\mathrm{)}\\ & & \mathrm{+}\mathrm{0.07518}{\mathrm{(}\begin{array}{c}\\ \end{array}\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{0.001492}{\mathrm{(}\begin{array}{c}\\ \end{array}\mathrm{)}}^{\mathrm{3}}\mathit{.}\end{array}$(20)Dust grains in strong radiation field environments present not only high temperatures, but also grain temperature fluctuations, as derived by Cuppen et al. (2006) and Iqbal et al. (2014) using Monte Carlo simulations, and by Bron et al. (2014) using an analytical approach. In this paper, however, we have not included a formalism taking temperature fluctuations of a grain size distribution into account, since we have large grains (≳50 Å because of the considered MRN distribution) and large fluctuations mainly occur for smaller grains (Draine & Li 2001)⁵.

3. Results

In this section, we discuss the results for three models in which we varied the radiation field⁶ (G₀) and the density (n_H) for a semi-infinite slab geometry and irradiation from one side without geometrical dilution. The adopted model parameters are listed in Table 3 (see also Fig. 1, which indicates the regions of the G₀ − n_H space occupied by the three models and various astrophysical objects). In Model 1, we study a typical PDR (e.g. the Orion Bar) characterised by high density and strong radiation field conditions. In Model 2 (typical conditions of an extreme starburst) and Model 3, we study the consequences of varying the density and intensity of the radiation field, respectively.

3.1. Heating and cooling

In Sect. 2.1, we listed the different heating and cooling mechanisms considered in our code. The results for the most relevant heating rates for Models 1-3 are shown in Fig. 2 as a function of the visual extinction. For both radiation fields and densities considered (see Table 3), the dominant heating source for visual extinctions A_V< 5 mag is photoelectric emission from grains. At A_V< 0.5 mag, if the intensity of the radiation field is high (Model 1 and 2), we also obtain a high contribution from the H₂ photodissociation heating. By contrast, if G₀ is low (Model 3) the second highest contribution to the heating arises from carbon ionisation. Deeper in the cloud, A_V> 2 mag, the heating by photoelectric emission, although still dominant, becomes progressively less efficient. For a model with high density and low G₀ (Model 3), cosmic rays dominate the heating of the region together with photoelectric emission processes at A_V ~ 5 mag. We also find that viscous heating never contributes significantly to the heating and that gas-grain collisions act as a cooling mechanism.

In Fig. 3, we show the cooling rates for Models 1–3 as a function of the visual extinction. For a PDR with intermediate density (n_H = 10⁴ cm^-3, Model 1), the cooling is dominated by [OI] 63 μm at A_V ≲ 5 mag. For a higher density PDR (Models 2 and 3), the cooling is dominated by [OI] 63 μm and gas-grain collisions up to A_V ~ 2 mag and ~0.5 mag, respectively. For higher extinctions, however, [OI] 63 μm becomes inefficient and gas-grain collisions represent the main coolant. Other processes, such as recombination of electrons with grains, represent minor coolants, especially at A_V> 1 mag.

3.2. Chemical structure

As stated previously, the physical and chemical processes in PDRs are dominated by interactions with photons, whose timescales are much lower than those for dynamical processes occurring in the opaque interiors of the clouds (Hollenbach et al. 2009). At very large visual extinctions, steady-state chemical codes, such as our PDR code, do not apply because of the need of considering the time-dependent chemical network, since certain chemical timescales are comparable to cloud lifetimes (~10⁷ yrs). Therefore, to analyse the chemical structure of different types of PDRs, we consider results for low visual extinctions (A_V ≲ 1 mag), where physical and chemical processes are purely dominated by interactions with photons, and for extinctions of translucent clouds (A_V ~ 1–5 mag).

3.2.1. Gas-phase species

In Fig. 4, we show the fractional abundances of several gas-phase species as a function of the visual extinction for Models 1, 2, and 3. We also plot the dust temperature (T_d) and the gas temperature (T_g) for each model. We obtain a H→H₂ transition with its location varying significantly depending on the radiation field and density. The atomic hydrogen is converted into H₂ at deeper locations in the cloud for high radiation fields (Models 1 and 2), since the photodissociation rates are larger. With a low G₀ and high density model (Model 3), however, the transition occurs closer to the cloud surface, since the chemical rates depend on n². We also observe that the H-H₂ transition becomes sharper as the intensity of the radiation field increases and the density decreases, in agreement with analytical results obtained by Sternberg et al. (2014).

Unlike results obtained by Meijerink & Spaans (2005) and with other PDR codes (which only consider gas-phase chemistry), where C⁺ presents high abundances at low visual extinctions, C at intermediate extinctions, and CO at high extinctions, here the transition C⁺ → C → CO is no longer well-defined. While the abundance of C⁺ is high at low extinctions for all models and it decreases as the cloud becomes denser, at intermediate extinctions (A_V ~ 2 mag) for Models 1 and 2, and at A_V ~ 0.5 mag for Model 3, most of the atomic carbon is rapidly converted into CO leading to low C abundances. Similar results for atomic carbon are found by Hollenbach et al. (2009) in their steady-state PDR code modelling the formation of CO and H₂O ices, which suggests that surface chemistry is accelerating the formation of CO. In Models 2 and 3, we also observe how the O abundance dramatically decreases by several orders of magnitude at A_V ~ 3 mag. These decreases are due to the oxygen depletion, which is locked in water ice. This fast increase of the solid water abundance can be seen in Fig. 5 (Sect. 3.2.2).

3.2.2. Dust-phase species

In Fig. 5, we show the fractional abundances of several solid species as a function of the visual extinction for Models 1, 2, and 3. The dash-dotted black line represents the number of possible attachable sites on grain surfaces per cm³ of space (see expression 4), i.e. the limit to reach a full monolayer of ice. We observe that the solid CO₂ and H₂O abundances in all of the models become larger than this limit, leading to the formation of several ice layers of these two species (see Fig. 6, which describes the exact number of ice layers of CO₂ and H₂O that are formed).

The radiation field strength is one of the most important factors governing ice formation, since it determines the dust temperature. Stronger radiation fields lead to higher dust temperature and prevent ice formation until large depths. We observe this effect comparing results from Models 2 and 3. In Model 3, with a low G₀ (10²), the dust temperature at A_V ~ 5 mag is <15 K, while in Model 2, T_d ~ 25 K owing to the higher intensity of the radiation field in this case (G₀ = 10⁴). This leads to the CO₂ ice formation at lower A_V (~0.7 mag) in Model 3 compared to Model 2 (A_V ~ 2.5 mag). Similar results are found for H₂O ices; the formation of a full water ice monolayer takes places at ~2 mag lower in Model 3 (low G₀) than in Model 2 (high G₀).

The density also plays an important role in the ice formation processes, since it affects the number of ice layers that are formed for water (see Fig. 6, Models 1 and 2). In particular, in the PDR with the lowest density (Model 1), ~6 ice monolayers of H₂O are formed at A_V ≤ 5 mag against the ~15 ice monolayers formed when the density is increased by two orders of magnitude. Comparing Models 2 and 3 in Fig. 6, we also conclude that the variation of the radiation field intensity barely affects the number of formed monolayers of CO₂ and H₂O ice. In a low radiation and high density PDR (Model 3, bottom pannel of Fig. 6), we find the formation of the first monolayer of CO₂ and H₂O ice at very low extinctions (A_V ≲ 1.5 mag). This is due to the shielding effect produced by the high density of the region, which prevents ice destruction by the impact of UV photons. In particular, in Model 3 (G₀ = 100 and n = 10⁶ cm^-3) we find the H₂O ice threshold extinction at A_V ~ 1.5 mag. Recent observational results of the IC 5146 dark coud (n ~ 10⁵ cm^-3) and of low-mass young stellar objects show the H₂O ice threshold extinction at A_V ~ 3 mag, which is equivalent to that found for the Taurus dark cloud (Chiar et al. 2011; Noble et al. 2013). Taking into account our results showing that a decrease of density in two orders of magnitude leads to an increase of the threshold extinction to form ices of ~1.5 mag, our model results are in agreement with observations. For CO, the ice threshold is significantly higher (A_V ~ 5–11 mag), according to observations of Taurus and ρ Ophiuchi (Whittet et al. 1989; Shuping et al. 2000; Velusamy et al. 2005).

In the previous Section, it was mentioned that the abundance of solid water increases as the gas-phase oxygen is depleted. In Fig. 5 (Models 2 and 3), we also observe that at a given high visual extinction (different for each model), the water ice growth slows down, while the gas-phase oxygen abundance keeps decreasing. This is because the ice abundance saturation after most of the O nuclei are locked in water ice. From results in Fig. 5, we also deduce that a low radiation field promotes the formation of solid methanol, since the surface reaction between H and H₃CO becomes more efficient as the visual extinction increases. Regarding to solid H₂O₂, we obtain the highest fractional abundances to be ~10^-14 for a low G₀ PDR and ~10^-10 for a high G₀ PDR in agreement with results from Ioppolo et al. (2008). After reaching its maximum abundances, we observe that solid H₂O₂ decreases sharply, mainly owing to its destruction by reacting with solid H to form solid water.

4. Discussion

4.1. Ice species formation rates

4.1.1. H₂O ice

In Fig. 7, we show rates for the main surface reactions forming solid water. We obtain that H₂O₂ on dust grains is an important intermediate in the formation of solid water for a high G₀ PDR (Models 1 and 2) at A_V ≲ 3 mag, in agreement with Du et al. (2012). However, the main surface reaction leading to water ice monolayers at larger extinctions is the reaction between solid H and OH. By contrast, for a low G₀ PDR (Model 3), the rates of this reaction decreases as the visual extinction increases, until the reaction between solid OH and H₂ becomes equally important in the water ice formation. Other surface reactions present rates that are too low (mainly due to the low abundances of some of the solid reactant) to significantly contribute to the formation of water ice at A_V ≲ 5 mag. This is, for example, the case of the surface reaction OH+CH₃OH→H₃CO+H₂O. We also find a small contribution to solid water formation from water depletion in all the models.

4.1.2. CO₂ ice

The rates of the main chemical reactions forming CO₂ on dust grains are shown in Fig. 8. We obtain that the surface reaction between solid O and solid CO dominates the formation of solid CO₂ before reaching its maximum number of monolayers (at A_V = 4.5, 3.5, and 1.5 mag for Models 1, 2, and 3, respectively). For larger extinctions, the CO₂ formation becomes dominated by the reaction between solid OH and solid CO. We also find that CO₂ depletion barely contributes to the formation of CO₂ on dust grains at A_V ≲ 5 mag.

4.2. Desorption probabilities

Table A.4 lists the desorption coefficients considering bare grains (δ_bare) and icy grains (δ_ice) for each two-body surface reaction. As previously mentioned in Sect. 2.2.3, the δ_bare coefficients were obtained theoretically (see expression (2) in Minissale et al. 2016), while the δ_ice coefficients were obtained considering the 10% of δ_bare. In order to analyse the effect of varying this desorption coefficient on the gas-phase abundances of different species, we also considered the case where δ_ice = 0 to compare it with the case where δ_ice = 0.1δ_bare. The results are shown in Fig. 9. We observe that the difference between not considering desorption when the two-body reaction takes place on icy grains and to consider a small percentage with respect to the desorption when the grains are bare becomes significant at A_V ≳ 4.5 mag. In particular, this variation in δ_ice implies small differences in the gas-phase abundances of several molecules, such as O₂ and H₂O, but differences of up to three orders of magnitude for other species, such as CH₃OH. We also find that these differences can be larger as the visual extinction increases, demonstrating the importance of considering desorption processes not only with bare substrate, but also with icy grains.

4.3. Bare versus icy grains

Although bare grains represent an ideal place for surface chemical reactions to occur, the presence of ice mantles can significantly enrich the gas phase when these ices are desorbed either by thermal or non-thermal processes. In Model 1 (G₀ = 10⁴ and n_H = 10⁴ cm^-3), we considered the case in which all grains are bare (f_ice = 0) and the case in which all grains are icy (f_bare = 0) in order to study in more detail the role of ice mantles on dust grains and to quantify the impact on the gas-phase and dust-phase abundances. A comparison of both results is shown in Fig. 10. For gas-phase water, unlike other species such as H, O₂, and CO, the presence of ice mantles leads to significant differences in its abundance at low extinctions (A_V< 1 mag). In particular, the water abundance increases up to about three orders of magnitude when grains are icy. This could be due to efficient thermal desorption at these low extinctions, since the water ice reservoir on grains is also significantly larger with the presence of ice mantles for A_V< 3 mag. For other species, such as CO, H₂CO, and CH₃OH, the differences between both cases are significant only for A_V> 3 mag, with a decrease in their gas-phase abundances of up to about two orders of magnitude when grains are icy at 3 <A_V< 6 mag.

Regarding the abundances of solid species (right panel of Fig. 10), we obtain large differences for several species at any extinction when considering that grains are bare or icy because of the different binding energies for each case. In addition, we also observe that the presence of ice mantles not only determines the abundance of each species, but also the visual extinction at which a full monolayer of ice is formed. In the case of CO₂ and H₂O, the formation of the first full monolayer occurs at ~1 mag lower when surface chemistry takes place on icy grains rather than on bare grains. With respect to the abundances, the largest differences are found for solid H₂O₂, with bare grains promoting its formation. CO, CO₂, and H₂O also present significant differences between both cases, although these differences become smaller as visual extinction increases, especially for CO₂ and water. When all grains are icy, the abundances of solid methanol also increase by about six orders of magnitude at A_V< 6 mag.

4.4. Effect of dust on the chemical composition of PDRs

We ran Model 1 (characteristic of starbursts) again considering that H₂ formation is the only reaction taking place on dust grains to analyse how the implementation of dust chemistry affects the abundances of gas-phase species. The results for different gas-phase species are shown in Fig. 11. In the top panel, we present results for the gas-phase abundances of H, O, CO, CO₂, and H₂O. For low visual extinctions, we barely find differences in the abundances of most of these species, independently of whether we consider dust chemistry or not. For A_V ≳ 4 mag, however, we clearly observe a decrease (up to three orders of magnitude) in all their abundances as a result of the formation of ices, when adsorption, desorption, and two-body processes on grain surfaces, along with the incidence of UV photons and cosmic rays on dust grains, are considered in the chemical network.

The bottom panel of Fig. 11 shows the abundances as a function of A_V for H₂CO and CH₃OH in the gas phase. We again observe large differences between both cases at A_V ≳ 4 mag. In particular, at A_V ~ 5 mag, we observe a change in the trend of H₂CO and CH₃OH gas-phase abundances; while their abundances start to decrease in the case without dust chemistry, they start to increase by several orders of magnitude when dust chemistry is implemented in the code, since they are mainly formed on grain surfaces (Chutjian et al. 2009). A similar trend is found for HCN, revealing the importance of considering grain chemical processes to explain the enhanced abundance of these molecules in the gas phase (Lintott & Viti 2006; Akimkin et al. 2013). In Fig. 12, we show the gas-phase HCN/HCO⁺ ratio obtained when surface chemistry is considered in the models with the highest radiation field intensities (Models 1 and 2). Although we find an increasing HCN/HCO⁺ ratio for A_V ≳ 0.1 mag in both cases, this increase is particularly pronounced at A_V ≳ 3.5 mag because of the presence of surface chemistry on grains.

4.5. Comparison with the original Meijerink PDR code

We also compared some of our results with those obtained by Meijerink & Spaans (2005) using the original version⁷ of the Meijerink code, in which H₂ formation was the only chemistry considered on dust grains. In Fig. 13 (left panel), we observe that in the original Meijerink code the decrease of the H abundance to form H₂ is slightly sharper than in the current version of the code. At A_V< 4 mag, while atomic oxygen does not present significant differences between both models, the update of the gas chemical network is mainly responsible for the increase of the abundances of CO, H₂O, and O₂ by several orders of magnitude. For higher extinctions, we also obtain large differences between both versions of the Meijerink code (higher abundances for H, O, and CO, and lower for O₂ and H₂O in the original code), but in this case, these differences are mainly due to the implementation of dust chemistry. For other molecules, such as HCO⁺, H₂CO, and CH₃OH (right panel), we observe the same effect at A_V ≲ 4 mag as that observed in the molecules of the top panel. At A_V> 4 mag, the differences observed between both versions of the code for HCO⁺ and CH₃OH are due to both the implementation of more than 3000 new gas-phase reactions and to surface chemistry, leading to a decrease and increase, respectively, of their gas-phase abundances of several orders of magnitude.

With respect to the thermal balance, we considered a more recent analytical expression for the dust temperature (obtained from Garrod & Pauli 2011, Sect. 2.3) than that considered in the original version of the Meijerink code (from Hollenbach et al. 1991). As we see in Fig. 14 (right panel), both expressions lead to very distinct values of T_d, especially for a high G₀ model (Model 1) at A_V< 0.5 mag with a difference of more than 30 K. For 0.5 <A_V< 5 mag, the differences for the dust temperature between both codes become smaller for the high G₀ PDR and larger for the low G₀ and high density PDR. Deep in the cloud, different dust temperature expressions significantly affect the gas temperature, determining how its profile decreases as the visual extinction increases. In the left panel of Fig. 14, we obtain differences of a few hundreds degrees between the two cases at A_V ~ 1 mag in Model 1, and of a few tens of degrees at A_V ≳ 2 mag in Model 3. These differences are mainly due to the influence of T_d on T_g through the heating and cooling rates. Variations in the reaction rates (for example, produced when considering a different chemical network with updated parameters for each reaction) can originate large differences in the number density of species with a direct impact on several heating and cooling processes (e.g. molecular cooling) and, therefore, in the final gas temperature of each model as well.

5. Comparison to observations

5.1. The Horsehead case

The Horsehead nebula PDR (d ~ 400 pc) is located at the western edge of the molecular cloud L1630 illuminated by the O9.5V star σ Ori (Habart et al. 2005). The far-UV intensity of the incident radiation field illuminating the Horsehead nebula is χ ~ 60 in Draine⁸ units (G₀ ~ 100; Habart et al. 2005; Goicoechea et al. 2007), which is moderate compared to those of classical PDRs illuminated by O stars (generally χ ~ 10⁴–10⁵; e.g. Tielens et al. 1993). In the last decade, several molecular observations were carried out in this region to study, for example, sulphur chemistry (Goicoechea et al. 2006), the presence of ions, such as C⁺ and CF⁺ (Guzmán et al. 2012), and the presence of organic molecules, such as CCH, C₄H, H₂CO, and CH₃OH (Goicoechea et al. 2007; Guzmán et al. 2011, 2013). In particular, recent results from a theoretical analysis of the chemistry forming H₂CO and CH₃OH indicate that the observations of these molecules in the PDR position of the Horsehead cannot be reproduced when only gas-phase chemistry is considered in the model (Guzmán et al. 2013).

We considered models with a radiation field intensity of G₀ = 100 (Models 4 and 5 with densities of 10⁴ and 10⁵ cm^-3, respectively) and G₀ = 65 (Model 6) to reproduce the conditions of the Horsehead PDR to analyse the origin of H₂CO and CH₃OH forming on dust grain surfaces. The results for the abundances of H₂CO and CH₃OH as a function of the visual extinction are shown in Fig. 15. We also included the value range for the observed abundances for both molecules (magenta lines) in the plots. These observations are taken from Guzmán et al. (2013). The left panel of Fig. 15 shows that the observations of H₂CO are reproduced by the three models at A_V ~ 2.5–3 mag and the right panel of Fig. 15 shows the model results for CH₃OH. In this case, we observe that only models with high density (n = 10⁵ cm^-3, Models 5 and 6) reproduce the observations at A_V ~ 6 mag. We also obtain that the difference in the radiation field intensity mainly determines the visual extinction at which the observations are reproduced (with a difference of ~0.5 mag).

Figure 16 shows the rates of the surface reactions leading to H₂CO (left panel) and CH₃OH (right panel) gas that we considered in our code. The results are for Models 5 and 6; those reproducing the abundances of both molecules, H₂CO and CH₃OH. In particular, for the case of H₂CO, we considered thermal desorption, direct cosmic-ray and photodesorption, and formation of H₂CO gas through chemical desorption coming from the reaction between solid hydrogen and solid HCO. For the CH₃OH formation, we also considered thermal desorption, cosmic-ray desorption, and the surface reaction between solid hydrogen and solid H₃CO. We did not include, however, direct photodesorption, since laboratory experiments from Bertin et al. (2016) show its minor role in the formation of gas methanol, as mentioned in Sect. 2.2.4. In both models we obtain that H₂CO gas is mainly formed through the surface reaction between solid hydrogen and solid HCO, whose rates are more than two orders of magnitude higher than those for direct photodesorption at A_V ≲ 4.5 mag. We find that direct photodesorption has a high contribution in the formation of H₂CO gas only at A_V> 4.5 mag. For the same models (5 and 6), we also obtain that methanol is mainly formed at A_V ~ 6 mag through chemical desorption upon the surface reaction between solid H and solid H₃CO. These findings disagree with recent results from Guzmán et al. (2013), who concluded that direct photodesorption is the main process of releasing methanol into the gas phase where thermal desorption of ice mantles is ineffective. This discrepancy is due to the overestimation of their photodesorption rates (Bertin et al. 2016).

5.2. The Orion Bar case

The Orion Bar (d ~ 414 pc) is located between the Trapezium cluster and the Orion molecular cloud (Menten et al. Menten2007). It is the prototypical warm PDR with a FUV radiation field at the ionisation front of G₀ ~ (1−4) × 10⁴ (Marconi et al. 1998) and densities of n ~ 10⁴–10⁶ cm^-3 (Leurini et al. 2010).

We also analyse the origin of H₂CO and CH₃OH gas in the interclump of the Orion Bar. For that, we considered observations from Leurini et al. (2010) and Models 1 and 2, as well as a new model with G₀ = 10⁴ and n ~ 10⁵ cm^-3. Figure 17 shows the model results for the H₂CO and CH₃OH gas abundances and the range of the observational abundances (magenta lines). We note that observations of H₂CO and CH₃OH are only reproduced at A_V> 7.5 mag and >10 mag, respectively. Since our code is a steady-state chemical code, we would need the implementation of time dependence to study the chemistry occurring in this visual extinction range properly.

6. Summary and conclusions

We have updated the Meijerink PDR code by implementing 3050 new gas-phase reactions and surface chemistry. In particular, we include 117 surface reactions including adsorption, thermal and non-thermal desorption, two-body reactions, photo processes and cosmic-ray processes on dust grains. In total, we used 7621 chemical reactions in this new version of the Meijerink code. In our dust treatment, we also distinguish between icy and bare grains. We analysed in detail how the variation of density and the intensity of the radiation field influence the thermal balance and gas/surface chemistry. In addition, we also studied how chemistry is affected by the presence/absence of ice mantles on dust grains. Our code also gives an indication of the origin of the formation of ices, such as CO₂ and water, that are present in different environments and, by comparing with observations, it shows the most efficient surface processes forming H₂CO and CH₃OH gas. The main results can be summarised as follows:

• Ice formation: the formation of a large number of H₂O ice monolayers is favoured by high densities. Varying the radiation field intensity mainly determines the threshold extinction to form ices of H₂O and, especially, of CO₂ (ice formation at lower A_V for low G₀). Our models also indicate that a low G₀ promotes the formation of CH₃OH and CO ices. In particular:

  • (a)

    H₂O ice formation: this formation is mainly dominated by thesurface reaction between solid H and solid OH atA_V ≲ 5 mag. H₂O₂ on dust grains is also an important intermediate in the formation of solid water for a high G₀ PDR.

  • (b)

    CO₂ ice formation: this formation is mainly dominated by the surface reaction between solid O and solid CO. Once the maximum number of monolayers is reached, the main reaction producing CO₂ on dust grains is the reaction between solid OH and solid CO.

• Bare versus icy grains: the presence or absence of ice mantles on dust grains determines the (gas and solid) abundances of many species leading to differences of several orders of magnitude for many of them, especially at A_V ≳ 3 mag. Some of the most affected species are water, methanol, CO, and H₂O₂ with differences of at least two orders of magnitude. The type of substrate (bare or icy) also varies the extinction at which ices of CO₂ and H₂O are formed with differences of up to 1 mag.

• Desorption probabilities: the difference between not considering chemical desorption with icy grains and to consider a small (10%) probability of desorption with respect to that for bare grains leads to variations of up to two-three orders of magnitude in the gas-phase abundances of molecules, such as methanol.

• PDR observations: the main mechanism releasing H₂CO and CH₃OH into the gas phase of the Horsehead PDR is chemical desorption upon the surface reaction between solid H with solid HCO and H₃CO, respectively, rather than photodesorption.

  These results demonstrate the important role of dust grains when studying chemistry in molecular clouds exposed to different UV radiation fields. They also show the need to implement time dependence in order to properly analyse the chemical composition at large visual extinctions.

---

Acknowledgments

This work is supported by the Netherlands Organization for Scientific Research (NWO). S.C. is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.017) and by the European Research Council (ERC; project PALs 320620). P.C. and M.S. acknowledge the financial support of the European Research Council (ERC; project PALs 320620). We thank S. Hocuk and R. Garrod for stimulating discussions on dust temperature. We also thank the referee F. Le Petit for his helpful comments.

References

Appendix A: Additional tables

We list in Tables A.1–A.6 all the surface reactions considered in the code and explained in Sect. 2.2. In addition, we show the binding energies (Table A.3) for each species depending on the type of the grain substrate (bare or icy). In Sect. B, we show the rates of the chemisorption reactions forming H₂.

Appendix B: Chemisorption reactions

We also consider in the code the chemical reaction JH→JH_c and the following reactions forming H₂ from chemisorbed hydrogen (H_c): $\begin{array}{ccc}\mathrm{H}\mathrm{+}\mathrm{J}\mathrm{\left(}{\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}& & \end{array}$(B.1)and $\begin{array}{ccc}\mathrm{J}\mathrm{\left(}\mathrm{H}\mathrm{\right)}\mathrm{+}\mathrm{J}\mathrm{\left(}{\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}\mathit{,}& & \end{array}$(B.2)where J(X) means solid X. The transmission coefficient through a barrier width a, with energy ϵ(K), is: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{chem}}\mathrm{=}\exp \left(\mathrm{-}\mathit{a}\sqrt{\frac{\mathrm{2}{\mathit{m}}_{\mathrm{red}}{\mathit{k}}_{\mathrm{B}}\mathit{ϵ}}{{ħ}^{\mathrm{2}}}}\right)\mathrm{+}\exp \left(\frac{\mathrm{-}\mathit{ϵ}}{{\mathit{T}}_{\mathrm{g}}}\right)\mathit{,}& & \end{array}$(B.3)where ħ is the Planck constant divided by 2π, m_red is the reduced mass of the reactants, and k_B is the Boltzmann constant. The chemisorption reaction rate, R_chem (cm^-3s^-1), for (B.1) is given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{chem}}\mathrm{=}{\mathit{k}}_{\mathrm{ads}}{\mathit{P}}_{\mathrm{chem}}{\mathit{f}}_{\mathrm{bare}}\mathit{S}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\sqrt{\frac{{\mathit{T}}_{\mathrm{g}}}{\mathrm{100}}}\mathit{,}& & \end{array}$(B.4)where k_ads is the adsorption rate coefficient (see Sect. 2.2.1) and S(T_g, T_d) the sticking coefficient.

When the formation of H₂ is through a chemical reaction between solid physisorbed and chemisorbed hydrogen (reaction B.2), it is necessary to consider the binding energies of H and H_c for bare grains (E_bare). In this case, the chemisorption reaction rate, R_chem (cm^-3 s^-1) is given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{chem}}\mathrm{=}{\mathit{\nu }}_{\mathrm{0}}\left(\frac{{\mathit{n}}_{\mathit{i}}{\mathit{n}}_{\mathit{j}}}{{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}}\right){\mathit{f}}_{\mathrm{bare}}{\mathit{k}}_{\mathrm{chem}}\mathit{.}& & \end{array}$(B.5)ν₀ is the oscillation factor and k_chem determines the probability of mobility by tunneling (m_tunn) and thermal diffusion (m_diff) $\begin{array}{ccc}{\mathit{k}}_{\mathrm{chem}}\mathrm{=}{\mathit{m}}_{\mathrm{tunn}}\mathrm{+}{\mathit{m}}_{\mathrm{tunn}}\mathrm{=}\mathrm{\left[}\mathrm{8}\sqrt{\mathit{\pi }{\mathit{T}}_{\mathrm{d}}}\mathrm{\right(}\mathit{A}\mathrm{\times }\mathit{B}\mathrm{\left)}\mathrm{\right]}\mathrm{+}\mathrm{\left[}\mathrm{4}\mathit{C}\mathrm{\right]}\mathit{,}& & \end{array}$(B.6)with $\begin{array}{ccc}\mathit{A}\mathrm{=}\frac{\sqrt{{\mathit{E}}_{\mathrm{bare}}\mathrm{\left(}{\mathrm{J}\mathrm{(}\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{)}\mathrm{-}{\mathit{E}}_{\mathrm{bare}}\mathrm{(}\mathrm{J}\mathrm{\left(}\mathrm{H}\mathrm{\right)}\mathrm{)}}}{{\mathit{E}}_{\mathrm{bare}}\mathrm{\left(}{\mathrm{J}\mathrm{(}\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{)}\mathrm{+}\mathit{ϵ}}\mathit{,}& & \\ \mathit{B}\mathrm{=}\exp \mathrm{(}\mathrm{-}\mathrm{0.406}\mathrm{\times }\mathit{d}\mathrm{\times }{\sqrt{{\mathit{E}}_{\mathrm{bare}}\mathrm{\left(}{\mathrm{J}\mathrm{(}\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{)}\mathrm{+}\mathit{ϵ}}}^{\mathrm{)}}\mathit{,}& & \end{array}$where d is the distance between physisorbed and chemisorbed sites (we assume 3 Å), and $\begin{array}{ccc}\mathit{C}\mathrm{=}\sqrt{\frac{{\mathit{E}}_{\mathrm{bare}}\mathrm{\left(}\mathrm{J}\mathrm{\left(}\mathrm{H}\mathrm{\right)}\mathrm{\right)}\mathrm{+}\mathit{ϵ}}{{\mathit{E}}_{\mathrm{bare}}\mathrm{\left(}{\mathrm{J}\mathrm{(}\mathrm{H}}_{\mathrm{c}}\mathrm{\right)}\mathrm{)}\mathrm{+}\mathit{ϵ}}}\mathrm{\times }\exp \left(\frac{\mathrm{-}\mathrm{\left[}{\mathit{E}}_{\mathrm{bare}}\mathrm{\right(}\mathrm{J}\mathrm{\left(}\mathrm{H}\mathrm{\right)}\mathrm{)}\mathrm{+}\mathit{ϵ}\mathrm{]}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathit{.}& & \end{array}$(B.9)See Cazaux & Tielens (2004) for more details about these rates.

All Tables

All Figures