# Multilayer modeling of porous grain surface chemistry.

## I. The GRAINOBLE model

V. Taquet, C. Ceccarelli, C. Kahane

► To cite this version:

V. Taquet, C. Ceccarelli, C. Kahane. Multilayer modeling of porous grain surface chemistry. I. The GRAINOBLE model. Astronomy and Astrophysics - A&A, 2012, 538, 10.1051/0004-6361/201117802. insu-03612452

HAL Id: insu-03612452

https://insu.hal.science/insu-03612452

Submitted on 18 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multilayer modeling of porous grain surface chemistry

## I. The GRAINOBLE model

V. Taquet, C. Ceccarelli, and C. Kahane

UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 38041 Grenoble, France
e-mail: vianney.taquet@obs.ujf-grenoble.fr

Received 1 August 2011 / Accepted 10 November 2011

### ABSTRACT

**Context.** Mantles of iced water mixed with carbon monoxyde, formaldehyde, and methanol are formed during the so-called prestellar core phase. In addition, radicals are also thought to be formed on the grain surfaces, and to react to form complex organic molecules later on, during the so-called warm-up phase of the protostellar evolution.

**Aims.** We aim to study the formation of the grain mantles during the prestellar core phase and the abundance of formaldehyde, methanol, and radicals trapped in them.

**Methods.** We have developed a macroscopic statistic multilayer model that follows the formation of grain mantles with time and that includes two effects that may increase the number of radicals trapped in the mantles: i) during the mantle formation, only the surface layer is chemically active and not the entire bulk; and ii) the porous structure of grains allows the trapping reactive particles. The model considers a network of H, O, and CO forming neutral species such as water, CO, formaldehyde, and methanol, plus several radicals. We ran a large grid of models to study the impact of the mantle multilayer nature and grain porous structure. In addition, we explored how the uncertainty of other key parameters influences the mantle composition.

**Results.** Our model predicts relatively high abundances of radicals, especially of HCO and CH₃O (10⁻⁹–10⁻⁷). In addition, the multilayer approach enables us to follow the chemical differentiation within the grain mantle, showing that the mantles are far from being uniform. For example, methanol is mostly present in the outer layers of the mantles, whereas CO and other reactive species are trapped in the inner layers. The overall mantle composition depends on the density and age of the prestellar core as well as on some microscopic parameters, such as the diffusion energy and the hydrogenation reactions activation energy. Comparison with observations allows us to constrain the value of the last two parameters (0.5–0.65 and 1500 K, respectively) and provide some indications on the physical conditions during the formation of the ices.

**Key words.** astrochemistry – ISM: abundances – ISM: clouds – ISM: molecules – molecular processes – stars: formation

## 1. Introduction

Among the 150 molecules that have been detected in the interstellar medium (ISM) so far, a significant number are complex organic molecules (COMs), carbon bearing molecules with more than five atoms. Several COMs have been observed in large quantities for two decades in the warm and dense hot cores of massive protostars (Blake et al. 1987). They have received renewed interest in the last few years after the detection of abundant COMs in solar-type protostars, specifically in hot corinos (Cazaux et al. 2003; Bottinelli et al. 2007), and in the clouds of the Galactic Center (Requena-Torres et al. 2006).

Astrochemical models have shown that many COMs cannot be produced efficiently in the gas phase. Horn et al. (2004), for example, predicted methyl formate abundances to be less than 10⁻¹⁰ by considering formation pathways only in the gas phase, but this molecule has been detected with abundances up to 10⁻⁶ (see references in the caption of Fig. 1). On the other hand, several observational and experimental works have highlighted the catalytic behaviour of interstellar grains for chemical reactions. First, infrared (IR) observations of protostellar sources carried out with space or ground-based telescopes have shown that formaldehyde and methanol can be significant components of the grain mantles with fractional abundances up to 30% with respect to water in some cases (Gibb et al. 2004; Boogert et al. 2008). These mantles are believed to be formed during the cold and dense prestellar core phase mainly via hydrogenation and

oxygenation of CO and O. Second, several laboratory experiments simulating cold cloud conditions have confirmed the efficiency of the hydrogenation processes. Watanabe & Kouchi (2002) for instance have succeeded in producing solid formaldehyde and methanol at low temperature (~10 K) via CO hydrogenation reactions in a CO-H₂O ice mixture.

Modelling the chemistry on the grain surfaces is therefore crucial if we aim to understand the formation of COMs on the interstellar grains. Over the past 30 years, several numerical methods based on ab initio (Allen & Robinson 1977), Monte Carlo (Tielens & Hagen 1982; Charnley et al. 1992), or rate equations (Hasegawa et al. 1992; Caselli et al. 1998a) approaches have been developed. All these methods give a macroscopic description of the problem, meaning that they follow the overall behaviour of the mantle rather than the single particle. This method, based on the rate equations, though the least precise from a physical point of view, is the fastest and allows studies of the chemical evolution in objects with varying physical conditions as a function of the time. In this context, time-dependent models have been developed that predict the abundance of many COMs in hot cores/corinos by considering two phases: a dense and cold pre-collapse phase during which the grain mantles with simple hydrogenated species are formed, and a warm-up phase, caused by the collapse, during which heavier species, and particularly radicals, can react together on the grain surfaces before sublimating into the gas phase (Garrod & Herbst 2006; Garrod et al. 2008; Aikawa et al. 2008). A key point of these models is
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Gas phase methyl formate to methanol abundance ratio as function of the gas phase abundance of methanol. Observations of hot corinos with single-dish telescopes are represented by red crosses: NGC 1333-4A by Bottinelli et al. (2004), IRAS16293 by Cazaux et al. (2003), NGC 1333-2A and -4B by Bottinelli et al. (2007). Observations of hot cores with single-dish telescopes are represented by blue crosses: G34.3+0.15 by Mehringer & Snyder (1996), SgrB2(N) and SgrB2(M) by Nummelin et al. (2000), G327.3-0.6 by Gibb et al. (2000), OMC1 by Sutton et al. (1995), G34.3+0.2, SgrB2(N), DR23(OH), W51, Orion Hot-Core by Ikeda et al. (2001), AFGL2591, G24.78, G75.78, NGC 6334, NGC 7538, W3(H2O), W33A by Bisschop et al. (2007). Observations of hot cores with interferometers are represented by blue diamonds: G34.3+0.15 by MacDonald et al. (1996), G19.61-0.23 by Qin et al. (2010), Orion KL, G29.96 by Beuther et al. (2009), and G47.47+0.05 by Remijan et al. (2004). The red curves report the predictions of the methyl formate to methanol abundance ratio appropriate to the hot corinos case by Garrod et al. (2008, their Fig. 6) and Awad et al. (2010, their Fig. 5). The blue curves report the predictions appropriate to the hot cores case by Garrod et al. (2008, their Fig. 4) and Laas et al. (2011, their Fig. 5). The arrows represent the direction of the time in each model and the ticks refer to timescales.

therefore the presence of radicals (OH, HCO, CH₃O) trapped in the cold mantles. In these models, the radicals are assumed to be mainly produced by the UV photodissociation of the neutral species in the grain mantles. **Garrod & Herbst (2006)** tested this assumption by introducing reactions between OH, HCO, CH₃, CH₃O to produce methyl formate, dimethyl ether, and formic acid. **Garrod et al. (2008)** then expanded the chemical network of reactions between radicals to form other COMs such as ethanol, glycolaldehyde, or acetic acid. Finally, **Aikawa et al. (2008)** and **Awad et al. (2010)** have included the spatial distribution and evolution of a collapsing cloud with the aim to estimate the size of typical hot corinos such as the well-known source IRAS16293-2422.

Although the abundances of the simplest COMs, such as formaldehyde (H₂CO) or methanol (CH₃OH), are now quite well predicted, astrochemical models still fail to reproduce the abundance of more complex ones, such as methyl formate (HCOOCH₃), for example. This molecule is assumed to be mainly formed on grains via the reaction between HCO and CH₃O during the warm-up phase (Allen & Robinson 1977; **Garrod & Herbst 2006**), a reaction in competition with the hydrogenation reaction of CH₃O, which leads to methanol. Figure 1 shows the observed and the predicted methyl formate to methanol gas phase abundance ratios as function of the gas phase methanol abundance. The plotted observations refer to the abundance values estimated in the articles cited in the figure caption, where the sizes are either directly estimated by interferometric observations or indirectly by other considerations for single-dish observations. The plotted model predictions refer to gas phase abundances. The curves “Garrod-F” and “Garrod-S” have been obtained from Figs. 4 and 6 of **Garrod et al. (2008)** respectively,

the curve “Awad10” comes from Fig. 5 of Awad et al. (2010), and the curve “Laas11” from Fig. 5 of Laas et al. (2011). Each point of the curves represents a different time or temperature (depending on the model). In these models, the methanol abundance increases with time, except in Awad et al. (2010) who considered the destruction of the species in the gas phase. We remark that while the methanol abundances may suffer of the uncertainty on the source size and H₂ column density, the methanol to methyl formate abundance ratios are almost unaffected by these uncertainties. Therefore, unless the observed methanol abundance is always lower than 10⁻⁷, which is certainly not the case, no model can reproduce the totality of observations. Specifically, the observed ratios are roughly independent of the methanol abundance and are about 0.1 and 1 in hot cores and hot corinos respectively, regardless of the telescope used for the observations. State-of-art astrochemical models consequently underestimate the methyl formate to methanol abundance ratios and always predict a decreasing ratio with increasing methanol abundance (because methyl formate is in competition with methanol formation).

Motivated by this unsatisfactory situation, we developed a new model, called GRAINOBLE, whose ultimate goal is to simulate the synthesis of COMs on the grain surfaces. In this article, we focus on the first phase, namely the formation of the mantle during the dense and cold phase of the pre-collapse, the prestellar core phase. Our model differs from most previous published models because it treats the multilayer and porous structure of the iced mantles. In the mentioned previous models (Garrod & Herbst 2006; Garrod et al. 2008; Aikawa et al. 2008; Awad et al. 2010; Laas et al. 2011), COMs are synthesized during the warm-up phase thanks to the increased mobility of radicals, formed on the grain surfaces during the previous cold phase. The abundance of radicals in the mantles during the cold phase is, therefore, a key point. One possibility, explored by the above models, is that the radicals are formed by the UV photodissociation of frozen species. However, since the impact of the UV photons on the surface chemistry is poorly understood (see the detailed discussion in Sect. 2.1), we explore here the possibility that radicals are synthesized even without photolysis and are trapped in the grain mantles, thanks to the multilayer mantle treatment.

The article is structured as follows. In the next section, we review the models available in the literature, outlining their limits and pointing out possible improvement. The concluding paragraph of the section describes the improvements proposed by our new model GRAINOBLE. In Sect. 3 we describe the technical details of GRAINOBLE. Section 4 reports the results of the new computations, discussing a reference set of parameters and a large grid of runs where several input parameters are varied. In Sect. 5 we compare the obtained results with previous models with the twofold goal to validate our code and discuss the different results from the different adopted physical models and conditions. Section 6 compares the GRAINOBLE predictions with the observations of ices and discusses the constraints on the various parameters of the model. Finally, Sect. 7 concludes the article, emphasising the main results and perspectives.

## 2. Limits of the existing models and need for improvement

### 2.1. The physics of grain surface chemistry

#### a) Multilayer vers bulk chemistry

One of the main limits of many published models is that they do not distinguish the chemical processes that occur in the
mantle bulk and on the surface. In these models, species that are buried in the mantle can keep diffusing in the mantle and reacting with others even if they are covered by new layers. Furthermore, particles landing on the surface on top of the grain mantle can react with all particles, regardless of their depth within the mantle. Thus, heavy reactive particles such as radicals can continue to react with landing hydrogen atoms even though the radicals are buried deep in the mantle bulk. The “bulk chemistry method” may therefore (grossly) underestimate the abundance of radicals in grain mantles.

In contrast, the difference in the chemical behaviour of the mantle and of the surface may lead to the formation of non-homogeneous grain ices as suggested by the IR observations. Tielens et al. (1991), for instance, by observing two features in the CO bands, suggested a polar mixture enriched by water in the bottom of the mantle covered by a non-polar one mainly composed of CO located on the outermost layers. In addition, experiments have shown that UV and cosmic rays irradiation can alter the chemistry of astrophysically relevant ices composed of water, methanol, and other CO-bearing molecules even in dark cold cloud conditions (Gerakines et al. 1996; Bennett & Kaiser 2007). However, unlike highly energetic cosmic rays, which can cross through the entire grain mantle, UV photons only penetrate a limited number of monolayers, depending on the absorption cross-section of molecules that constitute the ice (see Gerakines et al. 2001, for an illustrative discussion). Modelling the photolytic processes in ices requires, therefore, a distinction between the surface and the bulk of the mantle.

The desorption of grain mantles during the warm-up phase depends on their chemical composition. By carrying out temperature programmed desorption (TPD) experiments, Collings et al. (2004) demonstrated that CO ices desorb at ~20 K, whereas methanol and water ices desorb at temperatures higher than 100 K, for instance. Numerical values of desorption energies were deduced from these experiments and were incorporated in macroscopic rate equation models. Consequently, the model of Garrod et al. (2008) predicts a total desorption of the CO reservoir at 25 K, whereas more refractory species are desorbed at temperatures higher than 100 K. Actually, a significant amount of volatile ices can be trapped with highly polar molecules, such as water, in the mantle bulk. Thus, the desorption of such species is limited by the gradual evaporation of the mantle main matrix, which depends on other molecules.

The continuous time random walk (CTRW) Monte-Carlo method, introduced by Chang et al. (2005) and Cuppen & Herbst (2005), strongly differs from the macroscopic methods described above. This “microscopic” model considers the actual hopping of a discrete number of particles from one site to another but, unlike previous Monte Carlo methods, it takes into account the actual position of particles so that the spatial distribution is also included in the simulation. Cuppen & Herbst (2007) and Cuppen et al. (2009) have explained the growth of ice monolayers in dense clouds and showed a differentiation of the chemical composition in the mantle. More particularly, Cuppen et al. (2009) showed that under dense cloud conditions, CO can be trapped in the mantle by the accretion of new species above it before reacting with other particles. They also highlighted the survivability of very reactive particles in the mantle bulk by this process. Unfortunately, this method can only study simple chemical systems because of the CPU time required by the Monte-Carlo algorithm, and a direct link with the gas phase is not possible at present.

Hasegawa & Herbst (1993b) have attempted to treat separately the chemical processes within the bulk and on the surface

with a three-phase model. They considered the formation of an inert mantle caused by the accretion of particles onto it. They found that the differences between the three-phase and the two-phase models are small for stable species but become non-negligible for reactive species at long (≥10⁵ yr) timescales. Furthermore, the values of desorption and diffusion energies used in 1993 are lower than the values measured by the experiments carried out during the past decade. Because higher desorption and diffusion energies values tend to decrease the rate of the chemical processes on grain surfaces, the impact of the three-phase model on the frozen radical abundances would be significantly higher with the new values. However, at the time of the submission of this article, this three-phase model has not been pursued in other works.

### b) Porous grain surfaces

Astrochemical models used for predicting the COM formation have so far only considered perfectly spherical and smooth grains. However, observations (Mathis 1996), theoretical studies (Ossenkopf 1993; Ormel et al. 2009), and analyses of solar system bodies have provided evidence for the fluffy and porous structure of interstellar grains. For example, Mathis (1996) needed to use interstellar grains with a substantial fraction of vacuum to fit the observed extinction curves. Numerical simulations of dust coagulations carried out by Ossenkopf (1993) from a MRN-like grain size distribution have shown an increase of the size and the porosity of the dust in molecular cloud conditions. This fluffy structure of interstellar grains may alter the surface chemistry because grains can have micro-pores of small apertures compared to their volume in which species can be trapped. Thus, porous grains can increase the reactivity of chemical reactions by trapping volatile particles such as atomic hydrogen in their pores.

### c) Photolysis

As described in the introduction, the new generation of COMs formation models predicts that COMs are synthesized during the warm-up phase thanks to the increased mobility of radicals, formed on the grain surfaces during the previous cold phase (Garrod & Herbst 2006; Garrod et al. 2008; Aikawa et al. 2008; Awad et al. 2010; Laas et al. 2011). The abundance of radicals in the mantles during the cold phase is, therefore, a key point. In these works, UV photons from cosmic rays play a major role in synthesising radicals during this phase. However, as for any model, a series of inevitable assumptions are adopted to treat the process.

The first important assumption regards the photodissociation rates on the ices. They are assumed to be the same as those in the gas phase, computed by Sternberg et al. (1987); Gredel et al. (1989). However, there are various reasons to think that this may be a serious overestimate. First, grain surfaces could absorb part of the UV photon energy, ending up with a lower dissociation rate. Second, once the molecule is broken, if it is broken, being not in the gas phase but on a surface with almost no mobility, the photoproducts may recombine almost instantaneously. In addition, some products can even have enough energy to sublimate (see for example Andersson et al. 2006). Thus, the products of the photodissociation are likely different from those in the gas phase, and the overall formation rate of the assumed products (and radicals) may be severely lower than assumed.

A second point is related to the flux of CR induced UV photons in dense clouds. The exact value is fairly uncertain, because it largely depends on the primary cosmic ray energy spectrum, and on the grain extinction cross section, which are also
uncertain (see for example Padovani et al. 2009). For example, Shen et al. (2004) have shown that the typical uncertainty on the low-energy cosmic ray spectrum leads to a variation of the UV flux of more than one order of magnitude.

Another important assumption regards the penetration of the UV photons through the bulk of the mantle. Unfortunately, there exists almost no observational or theoretical ground for that, and it is not obvious how realistic this assumption is. Indeed, the experimental work by Gerakines et al. (2000) suggests that the UV photons can only penetrate a limited number of layers, depending on the optical properties of the ice.

Finally, even the branching ratios of species caused by the photodissociation in the gas phase are very poorly known (e.g., Laas et al. 2011). In this context, various authors omit the photolysis in their model (e.g., Chang et al. 2007; Cuppen et al. 2009).

## 2.2. The computational methods

The rate equations method allows us to study the evolution of the chemical composition of grain mantles via one temporal differential equation for each species. This approach allows us to use of complex chemical networks that include photolytic processes involving hundreds of species, and allows us to directly link the gas phase and the grain surfaces. However, rate equations are based on the evolution of densities of species in gas-phase and on grain surfaces, which are average values. Therefore this approach can be inaccurate when the number of particles on the grain becomes low because of the finite surface of grains. Caselli et al. (1998a) took into account the discrete nature of the individual composition of grain mantles by modifying the reaction rate coefficients, while Garrod et al. (2008) went further by modifying the functional form of the reaction rates. Biham et al. (2001) and Green et al. (2001) introduced a stochastic method based on the resolution of the master equations (ME), in which the system solves the time derivatives of the probability of each discrete mantle composition. In spite of its good accuracy for grain surfaces, this method can only be used for small chemical networks since the number of equations grows exponentially with the number of species. Lipshtat & Biham (2003) have thus introduced another method based on the time derivatives of the moments to study the formation of molecular hydrogen, while Barzel & Biham (2007) have extended this method for more complex chemical networks. The comparisons between the methods introduced here show that the rate equations method tends to overestimate the reaction rates and the abundance of the reactants, especially when they are in low quantities.

## 2.3. Need for improvement: the GRAINOBLE model

Although the past few years have seen a huge improvement in the modelling of grain surface chemistry, models are not yet able to accurately reproduce observations (e.g. Fig. 1). Therefore, something is still missing in those models. In the following, we review what are, in our opinion, the areas where improvements are needed and possible. Our new grain surface model, GRAINOBLE, has been developed to include the improvements described below.

### a) Multilayer structure of the ices

The first obvious limitation of the Herbst, Garrod, and collaborators class of models, the present state of the art, is that they do not take into account the multilayering formation of the ices, which leads to the differentiation observed in the

interstellar ices (Tielens et al. 1991; Pontoppidan et al. 2003) that is also predicted by the microscopic Monte-Carlo (MC) models of Cuppen and collaborators, and to the potential survival of reactive particles in the mantle bulk. On the other hand, the microscopic MC treatment is too cumbersome in computer time and difficult, if not impossible, to apply to model realistic cases. GRAINOBLE uses the rate equation method developed by Herbst and collaborators (which is at the base of the Garrod models) in a way that permits us to follow the multilayering structure of the ice. Therefore, it benefits from the advantages of the rate equation method (computer speed) and of the microscopic MC (multilayer) approach.

### b) Porosity of the grains

Laboratory experiments, numerical simulations of grain coagulation, mantle formation, and astronomical observations all demonstrate that interstellar grains are porous (Jones 2011). As previously discussed, the presence of pores creates traps for the molecules and atoms on the grain surfaces, leading to an enhancement of reactivity on the grains. GRAINOBLE includes pores, following the treatment of Perets & Biham (2006). The number and area of the pores are treated as parameters (Sect. 3).

### c) Multiparameter approach

Many parameters in the microphysics of the problem have highly uncertain values, like the various activation barriers of key chemical reactions, or the barriers against diffusion of adsorbed particules. Additionally, several macroscopic parameters are also either uncertain or vary depending on the astronomical object. In particular, it is worth mentioning that the densities and the temperatures of the prestellar cores depend on their age and their surroundings.

# 3. The GRAINOBLE model

## 3.1. General description

Following the work of Hasegawa et al. (1992), four main processes occuring on grains are taken into account.

1) Gas-phase particles can accrete onto the grains, which are considered spherical. For a given gas-phase species i, the accretion rate is a function of the thermal velocity of the gas-phase species $v(i)$, of the cross section of the grains $\sigma(a_d)$, of the density of the grains $n_d(a_d)$ where $a_d$ is the grain diameter, and of the sticking coefficient $S(i)$. We assumed the sticking coefficient estimated by Tielens (2005) for atomic hydrogen and a sticking coefficient equal to 1 for heavier particles, based on the TPD experimental work of Bisschop et al. (2006). For dark cloud conditions, particles constituting the ice are bound mainly via the physisorbed van de Waals interactions.

2) Once the particles are stuck onto the grains, they can diffuse along the surface via thermal hopping according to the Boltzmann law, which gives the hopping rate $R_{\text{hop}}$. Following experimental work of Katz et al. (1999), we neglected tunelling diffusion.

3) Physisorbed particles can react only via the Langmuir-Hinshelwood (Hinshelwood 1940) process, in which two species can react when they meet in the same site. The reaction rate $R_r$ is given by the product between the diffusion rate (i.e., the number of times per second that a species sweeps over a number of sites equal to the number of sites of the layer) and the probability of the chemical reaction, which is a function of its activation energy.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic view of a portion of the mantle layer as modelled by GRAINOBLE. All pores are assumed to be square, to have the same size, and no walls (see text). The edges sites, that connect the non-porous surface and the pores are depicted in blue. In this example, each pore is constituted by 16 sites, of which the 12 blue sites are edge sites, and the fraction of area occupied by the pores is 0.1.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic view of the multilayer approach adopted by the GRAINOBLE code. The inner layers of the mantle are chemically inert, whereas the layer at the surface is chemically reactive (see text).</figcaption></figure>

4) Surface species can desorb only via thermal processes, i) the “classic” thermal process caused by the thermal balance of the grain; ii) the cosmic ray induced desorption process in which cosmic rays heat the grains, as suggested by Leger et al. (1985), whose heating rate was computed by Hasegawa & Herbst (1993a). The sum of these two thermal desorption rates gives the total evaporation rate $R_{ev}$.

A full list of the symbols used in this work is provided in the Appendix.

## 3.2. Porosity

To model the impact of the porous structures of interstellar grains, we introduce two types of sites: the non-porous, and the porous sites, following Perets & Biham (2006). The non-porous sites form a smooth surface that is perfectly regular, where particles interact directly with the gas phase. The gas phase species can accrete onto the grain, and the adsorbates can desorb into the gas phase. The porous sites correspond to sites where no direct accretion and desorption are possible. They are filled up only by diffusion from non-porous sites at the edge of the pores. The porous sites cover a fraction $F_{por}$ of the grain surface. For the sake of simplicity, all the pores are assumed to be square and to have the same size, constituted by $N_{pore}$ sites, where $N_{pore} \ge 4$. The fraction of the grain surface occupied by the pores, $F_{por}$, is poorly known because, to date, no experimental measurements or numerical estimates have been attempted. Therefore we ran models with four values for $F_{por}$: 0, 0.3, 0.6, and 0.9, to cover a large range of possibilities,

The fraction of the edge porous sites, $F_{ed}$, on the grain depends then on $F_{por}$ and $N_{pore}$ as follows:

$$F_{ed} = 4 \times F_{por} \frac{\sqrt{N_{pore}} - 1}{N_{pore}}. \qquad (1)$$

A posteriori, the actual value of $N_{pore}$ does not influence the results significantly because what matters is mostly the fraction of area occupied by the pores and not the size or the number of pores. In the following, we assumed $N_{pore} = 9$.

The approach that we used for modelling the grain porosity is sketched in Fig. 2.

## 3.3. Multilayer approach

In order to distinguish the chemistry on the surface and in the bulk of the mantle, we used a multilayer approach. In practice, the code follows the chemical processes on a single layer, the one at the surface of the mantle. In this layer, particles can accrete, diffuse, react with each other and desorb, depending on whether they are on non-porous or porous sites. The mantle layer growth is treated at the same time as the chemistry. Specifically, at each time $t$, the chemistry is followed by solving the set of differential equations described in Sect. 3.9 between $t$ and $t + \Delta t$ and a new number of particles on the grain is computed at $t + \Delta t$. The layer is considered chemically inert as soon as the number of particles that are on the layer is equal to the number $N_s(a_d)$ of (porous plus not-porous) sites of the layer, where $N_s(a_d)$ is given by the ratio between the grain surface area and the site area $d_s^2$ (see Sect. 3.7). Because the porous sites are populated only via diffusion (Sect. 3.2), they fill up more slowly than the non-porous sites. It is therefore possible that the layer is considered inert when some pores are still empty. In this case, the number of non-porous sites is artificially increased. Although this is somewhat arbitrary, a posteriori this assumption has almost no impact on the results, as described in Sect. 4.3. When the layer becomes inert, the code memorizes its composition and a new and reactive layer is started. Note that no chemical exchange is allowed between layers. Finally, the code also takes into account the growth of the grain size and of the number of sites of the new layer by assuming that the thickness of the layer is equal to the site size $d_s$. All layers have, therefore, the same thickness, equal to the site size. Figure 3 gives a schematic view of this approach.

## 3.4. Chemical network

We adopted the network based on the work by Allen & Robinson (1977) and Tielens & Hagen (1982), updated for the new theoretical and experimental results and adapted to model the formation of the main constituents of interstellar ices. As previously mentioned in the introduction, observations of ices have shown that grain mantles are mainly composed of O- and CO- bearing molecules such as water, CO, $CO_2$, or methanol. Molecules such as $NH_3$ and $CH_4$ produced by the hydrogenation of N and C can also be found, but only in smaller quantities (Gibb et al. 2004). In this work, we therefore neglected the formation of ammonia and
<table><caption>**Table 1.** Species and binding energies, assuming that the ice matrix is mainly water ice.</caption><thead><tr><th>Species</th><th>E<sub>b</sub>(water ice) (K)</th></tr></thead><tbody><tr><td>H</td><td>450<sup>a</sup></td></tr><tr><td>O</td><td>800<sup>b</sup></td></tr><tr><td>CO</td><td>1150<sup>c</sup></td></tr><tr><td>OH</td><td>2820<sup>d</sup></td></tr><tr><td>H<sub>2</sub>O</td><td>5640<sup>e</sup></td></tr><tr><td>O<sub>2</sub></td><td>1000<sup>f</sup></td></tr><tr><td>O<sub>3</sub></td><td>1800<sup>f</sup></td></tr><tr><td>HO<sub>2</sub></td><td>5570<sup>f</sup></td></tr><tr><td>H<sub>2</sub>O<sub>2</sub></td><td>6300<sup>f</sup></td></tr><tr><td>HCO</td><td>1600<sup>f</sup></td></tr><tr><td>H<sub>2</sub>CO</td><td>2050<sup>f</sup></td></tr><tr><td>CH<sub>3</sub>O</td><td>5080<sup>f</sup></td></tr><tr><td>CH<sub>3</sub>OH</td><td>5530<sup>f</sup></td></tr></tbody></table>

Notes. (a) Hollenbach & Salpeter (1970); (b) Tielens & Hagen (1982); (c) Collings et al. (2004); (d) Avgul & Kiselev (1970); (e) Speedy et al. (1996); (f) Garrod & Herbst (2006).

methane. In addition, for the reasons explained in Sect. 2.1, we also do not consider photolytic processes. Our model, therefore, considers the formation of stable species (e.g., formaldehyde and methanol) and radicals (e.g., HCO, and CH₃O) via hydrogenation of frozen CO. Note that we did not consider CH₂OH because, in absence of photolytic processes, it is formed via the hydrogenation of H₂CO. Theoretical studies have shown that the reaction H₂CO + H → CH₂OH has an activation barrier (10.5 kcal/mol) more than twice the activation barrier leading to CH₃O (4–6 kcal/mol) (Sosa & Schlegel 1986; Walch 1993; Woon 2002). We therefore neglected CH₂OH in the present study. For the H₂O formation, we adopted the new reaction scheme based on the experimental works by Dulieu et al. (2010), Ioppolo et al. (2010) and Cuppen et al. (2010), who estimated the energy barriers of several reactions. The recombination reaction of hydrogen atoms, forming molecular hydrogen, is also taken into account. Since H₂ has a low binding energy with respect to water ice (Tielens & Allamandola 1987; Perets et al. 2005), and since its reaction of formation releases chemical energy (Hollenbach & Salpeter 1970), we assumed that the H₂ molecules are desorbed in the gas phase as soon as they are formed. Finally, IR observations have shown that CO₂ may be a major component of ices with abundances reaching 30% with respect to water in some cases (Dartois et al. 1999; Pontoppidan et al. 2008). However, the formation pathways of carbon dioxide, either by a direct surface route or by UV photolysis, in such cold conditions are still not well understood. Accordingly, we decided to neglect the formation of this molecule.

Table 1 lists the chemical species considered in our network with their binding energies relative to an amorphous water ice substrate. Indeed, interstellar grains located in prestellar cores are already covered by a polar ice mainly composed of water, formed during an earlier phase when atomic abundances are still high, before the depletion of CO. The network used for the grain surface chemistry is reported in Table 2.

### 3.5. Gas phase initial abundances

The observed and predicted gas-phase abundance of atomic carbon is much lower than CO or atomic oxygen abundances.

<table><caption>**Table 2.** Surface reactions and activation energies.</caption><thead><tr><th>Number</th><th>Reaction</th><th>E<sub>a</sub> (K)</th></tr></thead><tbody><tr><td>1</td><td>O + O → O<sub>2</sub></td><td>0</td></tr><tr><td>2</td><td>O + O<sub>2</sub> → O<sub>3</sub></td><td>0</td></tr><tr><td>3</td><td>H + H → H<sub>2</sub>.gas</td><td>0</td></tr><tr><td>4</td><td>H + O → OH</td><td>0</td></tr><tr><td>5</td><td>H + OH → H<sub>2</sub>O</td><td>0</td></tr><tr><td>6</td><td>H + O<sub>2</sub> → HO<sub>2</sub></td><td>0</td></tr><tr><td>7</td><td>H + HO<sub>2</sub> → H<sub>2</sub>O<sub>2</sub></td><td>0 (36%)<sup>a</sup></td></tr><tr><td>8</td><td>H + HO<sub>2</sub> → OH + OH</td><td>0 (57%)<sup>a</sup></td></tr><tr><td>9</td><td>H + H<sub>2</sub>O<sub>2</sub> → H<sub>2</sub>O + OH</td><td>1400<sup>b</sup></td></tr><tr><td>10</td><td>OH + OH → H<sub>2</sub>O<sub>2</sub></td><td>0</td></tr><tr><td>11</td><td>H + O<sub>3</sub> → O<sub>2</sub> + OH</td><td>0</td></tr><tr><td>12</td><td>H + CO → HCO</td><td>400<sup>c</sup> &lt; E<sub>a</sub> &lt; 2500<sup>d</sup></td></tr><tr><td>13</td><td>H + HCO → H<sub>2</sub>CO</td><td>0</td></tr><tr><td>14</td><td>H + H<sub>2</sub>CO → CH<sub>3</sub>O</td><td>400<sup>c</sup> &lt; E<sub>a</sub> &lt; 2500<sup>d</sup></td></tr><tr><td>15</td><td>H + CH<sub>3</sub>O → CH<sub>3</sub>OH</td><td>0</td></tr></tbody></table>

Notes. (a) Branching ratios measured by Cuppen et al. (2010); (b) Klemm et al. (1975); (c) Fuchs et al. (2009); (d) Woon (2002).

Therefore, we considered only the accretion of gaseous H, O, and CO onto the grain mantles.

In dark cloud conditions, the abundance of atomic hydrogen in the gas phase results from a balance between the cosmic rays ionization of H₂ on the one hand and the accretion on dust grains (which react to form H₂ or heavier iced species) on the other hand. At steady state, the density of H is therefore given by the ratio between these two processes (Tielens 2005), as follows:

$$n(H) = \frac{2.3\zeta_{CR}n(H_2)}{2\nu(H)\sigma_d X_d n(H_2)} \quad (2)$$

where $v(H)$ is thermal velocity of atomic hydrogen. The resulting H density is independent of the density and equal to $1.2 \text{ cm}^{-3}$ at 10 K, for a cosmic rays ionization rate $\zeta$ of $3 \times 10^{-17}$ (Caselli et al. 1998b), a grain abundance $X_d$ of $1.33 \times 10^{-12}$ relative to H nuclei and for a grain radius $r_d$ of $0.1 \text{ µm}$ ($\sigma_d = \pi r_d^2$), giving $n_d \sigma_d / n_H = 4.2 \times 10^{-22} \text{ cm}^2$. Observations and astrochemical models show that CO is a very stable species and the most abundant molecule after H₂ in dark clouds. Following the estimates of Frerking et al. (1982), we assumed an initial abundance of CO of $4.75 \times 10^{-5}$ with respect to H nuclei. This abundance decreases with time because of the freeze-out onto the grain mantles. Finally, the abundance of atomic oxygen in the gas phase remains poorly constrained. Observations of the [OI] 63 µm line obtained by Caux et al. (1999); Vastel et al. (2000); Lis et al. (2001) with the Infrared Space Observatory (ISO) have measured abundances higher than $10^{-4}$, suggesting that all the gaseous oxygen is in atomic form. But astrochemical models predict an atomic oxygen abundances of $6 \times 10^{-5}$ or less. The initial abundance of O is therefore considered as a free parameter and can vary between $2 \times 10^{-4}$ and $2 \times 10^{-5}$.

### 3.6. Physical conditions

As supported by observational and theoretical arguments, the physical properties of prestellar cores vary substantially with the age and the mass of the object (see di Francesco et al. 2007; Bergin & Tafalla 2007). The densities $n_H$ evolve typically from $\sim 10^4$ to $\sim 10^6 \text{ cm}^{-3}$, while the gas kinetic temperature (equal to the grain temperature in most cases) varies between 7 and 20 K,
depending on the type of the object. Indeed, massive stars are suspected to be formed from warmer prestellar cores than low mass ones. Consequently, we included several values of the density and the temperature, listed in Table 3, in our parameter grid.

Additionally, as shown by numerical simulations (Ossenkopf 1993; Ormel et al. 2009) and observations (Stepnik et al. 2003; Steinacker et al. 2010), the coagulation of grains occuring in molecular clouds and in prestellar cores tend to increase their size. That is why the grain size distribution, which is assumed to follow a broad nearly power-law in diffuse interstellar clouds (Mathis et al. 1977), tends to evolve with time, the bigger grains dominating most of the dust mass. We assumed all interstellar grains to have the same diameter $a_d$, and we considered three values: 0.1, 0.2, and 0.3 µm respectively. The grain abundance $X_d$ changes in accordance with the grain size, to keep a dust-to-gas mass ratio of 1%.

### 3.7. Other key parameters of GRAINOBLE

In addition to the above described parameters, astrochemical models of grain chemistry are based on several other parameters related to the particle-surface, surface particle-particle interaction and surface properties. We list below three additional key parameters of the GRAINOBLE model whose values are poorly constrained because of the difficulties involved in measuring or computing them. We considered them as free parameters in a range defined by the different values found in the literature.

#### a) Diffusion energy to binding energy ratio $E_d/E_b$

Theoretical calculations of barriers against diffusion for physical adsorption $E_d$ on perfectly smooth surfaces have been carried out by Jaycock & Parfitt (1986) and have shown that $E_d$ is typically about 30% of the binding energy $E_b$. However, “real” surfaces, which show defects, irregularities, or steps, tend to increase this ratio (Ehrlich & Hudda 1966). Accordingly, several experimental studies have been carried out to constrain the diffusion barrier on astrophysically relevant surfaces. By fitting experimental formation of molecular hydrogen on refractory bare surfaces and on amorphous water ices with a rate equation model, Katz et al. (1999) and Perets et al. (2005) experimentally found high $E_d/E_b$ ratios (0.77 and ~0.85 respectively). However, Collings et al. (2003), Ulbricht et al. (2002), and Matar et al. (2008) estimated lower energy ratios (of about 0.5) for CO and D on water ice, and for atomic oxygen on carbon nanotubes. Clearly, the energy $E_d/E_b$ ratio is highly uncertain and strongly depends on the composition and the structure of the substrate. Three values of the $E_d/E_b$ energy ratio have been considered here: 0.5, 0.65, and 0.8.

#### b) Activation energy $E_a$

The activation energy values of the CO and H₂CO hydrogenation reactions are also very uncertain. Theoretical studies of these reactions on an icy mantles via calculations of quantum chemistry carried out by Woon (2002) and Goumans (2011), showed that the activation energies are about 2000–2500 K, and depend slightly on the presence of water molecules around the reactants. In contrast to these theoretical values, experiments show that these reactions take place at low temperature (Watanabe & Kouchi 2002). Fuchs et al. (2009), by modelling the experiments of a hydrogen deposition on a CO ice with a CTRW method, found values of about 400–500 K for both reactions. We therefore consider these barriers as a free parameter. Three values of activation energies, assumed to be equal

<table><caption>**Table 3.** Key free parameters of GRAINOBLE and value range.</caption><thead><tr><th>Parameter</th><th>Values</th></tr></thead><tbody><tr><td>Density n<sub>H</sub></td><td>10<sup>4</sup> – <b>10<sup>5</sup></b> – 10<sup>6</sup> cm<sup>-3</sup></td></tr><tr><td>Temperature T<sub>g</sub> = T<sub>d</sub></td><td>10 – <b>15</b> – 20 K</td></tr><tr><td>Initial oxygen abundance X(O)<sub>ini</sub></td><td>2 × 10<sup>-5</sup> – <b>6</b> × 10<sup>-5</sup> – 2 × 10<sup>-4</sup></td></tr><tr><td>Grain size a<sub>d</sub></td><td>0.1 – <b>0.2</b> – 0.3 µm</td></tr><tr><td>Energy ratio E<sub>d</sub>/E<sub>b</sub></td><td>0.5 – <b>0.65</b> – 0.8</td></tr><tr><td>Porosity factor F<sub>por</sub></td><td><b>0</b> – 0.3 – 0.6 – 0.9</td></tr><tr><td>Site size d<sub>s</sub></td><td>1.4 – <b>4.2</b> – 7 Å</td></tr><tr><td>Activation energy E<sub>a</sub></td><td><b>400</b> – 1450 – 2500 K</td></tr></tbody></table>

**Notes.** Bold values mark the values adopted in the reference model (see text, Sect. 4).

for the H+CO and H+H₂CO reactions, have been considered: $E_a$ = 400, 1450, 2500 K.

#### c) Site size $d_s$

The surface density of sites $s$, i.e. the number of sites per surface unit given in sites x cm⁻², depends on the structure of the surface and its composition. Most astrochemical models assume a value of $s = 10^{15}$ cm⁻² which has been measured experimentally by Jenniskens et al. (1995) for a high-density amorphous water ice. By assuming that the distance between two sites $d_s$ is constant along the surface, $d_s$ can be easily deduced from the surface density $s$ ($d_s = 1/\sqrt{s}$), and is equal to 3.2 Å. However, Biham et al. (2001) found a value of $s = 5 \times 10^{13}$ cm⁻² ($d_s = 14$ Å) for an amorphous carbon surface and Perets & Biham (2006) assumed a value of $s = 5 \times 10^{15}$ cm⁻³ ($d_s = 1.4$ Å) in their model of porous grains. We therefore treat $d_s$ as a free parameter, with values equal to 1.4, 4.2 and 7 Å.

### 3.8. Multi-parameter approach

To summarize, the GRAINOBLE model depends on eight key parameters presented in the previous paragraphs and listed in Table 3 along with the range and values considered in this study for each of them. The values reported in bold are used in the “reference” model (Sect. 4). They correspond to the average physical conditions of prestellar cores (density, temperature, and grain size), to the average values found in the literature (site size, diffusion energy to binding energy ratio, initial abundance of atomic oxygen), to the values of parameters used by most of previous astrochemical models (porosity factor) or to the most recent values measured by experiments (activation energy).

In total, we ran a grid of 17 496 models varying these eight free parameters.

### 3.9. Computational aspects

The code solves three sets of differential equations: the first one describes the evolution of the densities of species (other than H, whose density is given by Eq. (2)) in the gas phase, the other two describe the chemical composition of the mantle outermost layer on the non-porous surface and within the pores. The equations giving the evolution of the chemical composition of the mantle layer do not depend on the grain size and are expressed in monolayers/sec (MLs/s) where a monolayer is given by the number of particles of the considered species divided by the number of sites
of a layer. The equations are the following:

$$ \frac{dn_g(i)}{dt} = -S(i)v(i)n_g(i)\sigma(a_d)n_d(a_d) + R_{ev}(i)P_{np}(i)N_s(a_d)n_d(a_d) \quad (3) $$

$$ \frac{dP_{np}(i)}{dt} = \frac{1}{4}S(i)d_s^2v(i)n_g(i) - R_{ev}(i)P_{np}(i) - \frac{F_{ed}}{F_{np}}R_{hop}(i)P_{np}(i) + \frac{F_{ed}}{F_{por}}R_{hop}(i)P_{por}(i) + \sum_{i_f} R_{r,ext}(i_f)P_{np}(i_{r1})P_{np}(i_{r2}) - \sum_{i_d} R_{r,ext}(i_d)P_{np}(i_{r1})P_{np}(i_{r2}) \quad (4) $$

$$ \frac{dP_{por}(i)}{dt} = + \frac{F_{ed}}{F_{np}}R_{hop}(i) \cdot P_{np}(i) - \frac{F_{ed}}{F_{por}}R_{hop}(i)P_{por}(i) + \sum_{i_f} R_{r,in}(i_f)P_{por}(i_{r1})P_{por}(i_{r2}) - \sum_{i_d} R_{r,in}(i_d)P_{por}(i_{r1})P_{por}(i_{r2}) \quad (5) $$

where $P_{np}(i)$ and $P_{por}(i)$ are the surface populations of the species $i$ (in MLs) on the non-porous surface and within the pores respectively, and $n_g(i)$ its gas-phase density. The first and second terms in Eq. (3) describe the accretion rate and the evaporation rate. The first and second terms of Eq. (4) describe the accretion and the evaporation rates. The third and fourth terms describe the exchange rate between the non-porous surface and the pores via diffusion. The fifth and sixth terms describe the production and the destruction rates of $i$ via reactions on the non-porous surface. Similarly, the first and second terms of Eq. (5) describe the exchange rate between the non-porous surface and the pores via diffusion. The third and fourth terms describe the production and the destruction rates of $i$ via reactions in the pores. Note that the evolution of composition within all the pores can be treated with a single equation because all pores are assumed to have the same size. The equations giving the evolution of each pore can be deduced from Eq. (5) by multiplying it with a simple factor of proportionality.

### 3.10. What the model does not consider

We aim to study the impact of the multilayer behaviour and the porous structure of the grains on the mantle formation and its chemical composition. To this end, we considered thermal processes only. Therefore we did not consider the desorption processes caused by exothermic surface reactions (Garrod et al. 2007), and the photolytic desorption (Öberg et al. 2009b,a; Arasa et al. 2010) in spite of their relative possible importance. Neither did we account for the photodissociation processes on the ice in this version of GRAINOBLE owing to the reasons given in Sect. 2.1.

## 4. Results

We ran a grid of 17496 models, varying the eight free parameters listed in Table 3. In this section, we describe the results of our model. We start by describing the results relative to the reference model, and then we describe how the results change when each of the parameters changes. The analysis is based on the distribution of the obtained mantle abundances when all parameters are varied, except for the parameter under consideration.

This allows us to understand if and by how much the considered parameter influences the results. The plot of the multilayer and bulk distributions is shown in Sect. 4.2, while the distributions of other parameters (using the multilayer approach only) are shown in the Appendix. Finally, Table 4 summarises the main effects on the resulting mantle abundances caused by each parameter.

### 4.1. The reference model

This section describes the results obtained using the reference set of parameters (the boldface values in Table 3). Below we refer to it as the “reference model”. Figure 4 presents the evolution with time of the abundance of each species in the gas phase and in the grain mantles.

The evolution of the chemical composition of the mantle depends on the initial gas phase abundances. In the reference model, O and CO abundances are approximately five times higher than the H abundance. Consequently, H atoms cannot hydrogenate all O and CO molecules that accrete onto the mantles. Besides, because the CO hydrogenation reaction has a barrier of 400 K, H atoms react preferentially with O and its products via barrierless reactions, ending up in water molecules. In $2 \times 10^5$ years, 80% of the CO reservoir is trapped within the mantle bulk, whereas ~60% of O atoms are contained in water. The remaining O atoms are shared between $O_3$ (ozone), OH and $H_2O_2$ (hydrogen peroxyde). Indeed, the high O to H initial gas phase abundance ratio (~6) and the relatively high temperature (15 K) allow a significant amount of O atoms to diffuse on the surface and to react with each other before meeting hydrogen atoms.

However, gas phase CO and O abundance ratios relative to H gradually decrease with time (because they freeze-out onto the mantles), increasing the possibility of hydrogenation reactions. After ~5 $\times 10^4$ years, the formation rates of formaldehyde and methanol start to increase while the formation of ozone stops, corresponding to the time when H, O, and CO have similar abundances. The percentage of ozone, formaldehyde, and methanol in the ices therefore strongly depends on the age of the core. Young cores would show high abundances of ozone and of hydrogen peroxyde with respect to water, whereas older ones would show high abundances of formaldehyde and methanol.

Another important result of the model is that a significant amount of radicals is trapped in the bulk because of the multilayer treatment. Radicals assumed to be the precursors of COMs such as OH, HCO or $CH_3O$, reach abundances between ~5 $\times 10^{-6}$ and $10^{-9}$. OH is more abundant because of its barrierless formation reaction, and its formation mostly occurs at shorter times.

The evolution of the chemical composition can also be studied “spatially”. Indeed, our multilayer approach allows us to study the composition of each monolayer within the grain mantle. Figure 5 shows the evolution of the formation time of each layer and the evolution of the mantle thickness with time. In the reference model, a grain mantle of 77 layers is formed in $2-3 \times 10^5$ yr. Assuming that the layer thickness is equal to the site size (here 4.2 Å), it gives a total mantle thickness of ~0.04 µm, namely 32% of the grain radius. The first 60 layers are created in less than $10^5$ yr because they have a fast formation rate (each layer is formed in less than ~2 $\times 10^3$ yr). The formation time of layers increases sharply at ~$10^5$ yr, because of the drop of the gas phase abundances.

The differentiation of the composition between the inner and the outer layers of the mantle can also be seen in Fig. 6.
<table><caption>**Table 4.** Summary of the effect of each free parameter on the surface chemistry and on the distribution of results (shown in the Appendix).</caption><thead><tr><th>Parameter</th><th>Effect on the surface chemistry</th><th>Distributions of abundances for stable species</th><th>Distributions of abundances for radicals</th></tr></thead><tbody><tr><td>Multilayer versus bulk approach</td><td>The multilayer approach<br/>– decreases the formation of stable species<br/>– traps reactive species in the mantle</td><td>– X(CO) &gt; 10<sup>-5</sup> for 88% (ML), 50% (bulk) cases<br/>– X(H<sub>2</sub>O) &gt; 10<sup>-5</sup> for 60% (ML), 93% (bulk) cases<br/>– X(CH<sub>3</sub>OH) &gt; 10<sup>-5</sup> for 15% (ML), 55% (bulk) cases</td><td>– X(OH) &gt; 10<sup>-9</sup> for 75% (ML), 1% (bulk) cases<br/>– X(HCO) &gt; 10<sup>-9</sup> for 25% (ML), 2% (bulk) cases</td></tr><tr><td>Porosity</td><td>The increase of porosity<br/>– traps volatile and diffusing particles<br/>– increases the formation rate of stable species</td><td>– X(H<sub>2</sub>CO) &gt; 10<sup>-7</sup> for 52% (smooth), 61% (porous) cases<br/>– X(CH<sub>3</sub>OH) &gt; 10<sup>-7</sup> for 51% (smooth), 61% (porous) cases</td><td>– X(OH) &gt; 10<sup>-6</sup> for 55% (smooth), 20% (porous) cases</td></tr><tr><td>Density</td><td>The increase of density<br/>– decreases the formation of stable species formed by hydrogenation reactions</td><td>– X(H<sub>2</sub>O) &gt; 10<sup>-5</sup> for 2% (dense), 85% (sparse) cases with X(H<sub>2</sub>O) &gt; 2 &times; 10<sup>-4</sup><br/>– X(CH<sub>3</sub>OH) &gt; 10<sup>-5</sup> for 0% (dense), 35% (sparse) cases with X(CH<sub>3</sub>OH) &gt; 9 &times; 10<sup>-5</sup></td><td>– X(OH) &gt; 10<sup>-6</sup> for 0% (sparse), 75% (dense) cases</td></tr><tr><td>Temperature</td><td>The increase of temperature<br/>– strongly increases the desorption of H and O<br/>– decreases the rates of hydrogenation reactions – decreases the abundance of stable species formed by these reactions</td><td>– X(CH<sub>3</sub>OH) &gt; 10<sup>-7</sup> for 65% (cold), 45% (warm) cases<br/>– X(H<sub>2</sub>CO) &gt; 10<sup>-7</sup> for 60% (cold), 45% (warm) cases</td><td>– X(OH) &gt; 10<sup>-6</sup> for 55% (cold), 48% (warm) cases</td></tr><tr><td>Grain size a<sub>d</sub></td><td>The increase of the grain size<br/>– increases the mantle thickness<br/>– does not influence the integrated mantle composition</td><td>Same evolution of distributions for all stable species</td><td>Same evolution of distributions for all radicals</td></tr><tr><td>Initial abundance of atomic oxygen</td><td>The increase of X<sub>ini</sub>(O):<br/>– increases the formation of water (and other molecules formed from reactions involving O)<br/>– slightly decreases the formation of molecules formed from CO</td><td>– X(H<sub>2</sub>O) &gt; 2 &times; 10<sup>-5</sup> for 0% (low-O), 65% (high-O) cases with X(H<sub>2</sub>O) &gt; 2 &times; 10<sup>-4</sup><br/>– No influence is seen for CO, H<sub>2</sub>CO and CH<sub>3</sub>OH</td><td>– X(OH) &gt; 4 &times; 10<sup>-6</sup> for 0% (low-O), 50% (high-O) cases with X(OH) &gt; 4 &times; 10<sup>-5</sup></td></tr><tr><td>Diffusion energy</td><td>The increase of the diffusion energy<br/>– decreases the diffusion rate of mobile species<br/>– decreases the formation rate of stable species<br/>– slightly increases the survival of radicals</td><td>– X(H<sub>2</sub>O) &gt; 5 &times; 10<sup>-6</sup> for 87% (fast), 74% (slow) cases<br/>– X(CH<sub>3</sub>OH) &gt; 10<sup>-7</sup> for 80% (fast), 30% (slow) cases<br/>– X(H<sub>2</sub>CO) &gt; 10<sup>-7</sup> for 75% (fast), 40% (slow) cases</td><td>Same evolution of distributions for all radicals</td></tr><tr><td>Activation energies</td><td>The increase of activation energies<br/>– strongly decreases the formation rate of H<sub>2</sub>CO and CH<sub>3</sub>OH<br/>– strongly decreases the survival of radicals HCO and CH<sub>3</sub>O<br/>– slightly increases the reaction rates of reactions involving O<br/>– slightly increases the formation of water</td><td>– X(CO) &gt; 10<sup>-5</sup> for 35% (low-Ea), 98% (high-Ea) cases<br/>– X(H<sub>2</sub>CO) &gt; 10<sup>-7</sup> for 83% (low-Ea), 25% (high-Ea) cases<br/>– X(CH<sub>3</sub>OH) &gt; 10<sup>-6</sup> for 80% (low-Ea), 5% (high-Ea) cases</td><td>– X(HCO) &gt; 10<sup>-10</sup> for 90% (low-Ea), 0% (high-Ea) cases<br/>– X(CH<sub>3</sub>O) &gt; 10<sup>-10</sup> for 70% (low-Ea), 0% (high-Ea) cases</td></tr><tr><td>Site size</td><td>The increase of the site sizes<br/>– decreases the formation time of layers<br/>– increases the mantle thickness<br/>– does not influence the integrated mantle composition</td><td>Same evolution of distributions for all stable species</td><td>Same evolution of distributions for all radicals</td></tr></tbody></table>

Indeed, “intermediate” molecules such as O<sub>3</sub>, which is created from atomic oxygen, are rapidly formed and are therefore abundant in the inner layers of the mantle (the first 50 layers), whereas formaldehyde and methanol are mostly formed

later, their fractional composition showing significant abundances only in the 20 outermost layers. The abundances of the two main species, CO and H<sub>2</sub>O, remain relatively constant throughout the mantle up to the last outermost layers, where
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Abundance of species in interstellar grain mantles (solid) and in gas phase (in dashed) as a function of time for the reference model (boldface values in Table 3). **Left panel:** the most abundant stable species. **Middle panel:** reactive species. **Right panel:** less abundant stable species.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Formation time of each monolayer (blue line) and mantle formation time versus its thickness expressed in monolayers (red line) for the reference model (boldface values in Table 3).</figcaption></figure>

they decrease in favour of H₂CO and CH₃OH. The composition of radicals within the mantle follows the same evolution as their precursor molecules: the HCO abundance stays relatively constant until the drop in the last layers like CO, while the OH abundance decreases following the O abundance behaviour. Finally, the CH₃O abundance increases approaching the surface, as H₂CO.

## 4.2. Multilayer versus bulk approach

The introduction of the multilayer approach modifies the chemical behaviour of grain mantles and thus their overall chemical composition. Figure 7 shows the evolution of the mantle composition with time for the reference model, adopting a multilayer (described in Sect. 3.3) or a bulk (old method described in 2.1) approach, respectively. The figure shows that at short timescales (≤5 × 10⁴ yr), both approaches give similar abundances for the main species because the H abundance is initially low compared to the CO and O abundances (Sect. 4.1). In contrast, O and its associated “intermediate” molecules (OH, O₂, and HO₂) are more efficiently burned in O₃, H₂O₂, and H₂O in the bulk compared to the multilayer approach.

At later times (≥5 × 10⁴ yr), reactive species, which are trapped by the multilayer approach, continue to react in the bulk method. Consequently, in the bulk approach, radicals (OH, HCO, HO₂, CH₃O) and “intermediate” stable species (O, O₃, H₂O₂, O₂) are totally burned in less than ~10⁶ yr. In contrast, the abundances of these species do not evolve after ~3 × 10⁵ yr with the multilayer approach. Finally, the predicted methanol and formaldehyde abundances substantially diverge in the two approaches at times longer than ~10⁵ yr, namely the presumed ages of the prestellar cores. The bulk method predicts a final (at 10⁶ yr) methanol abundance higher by factor 6 compared to the multilayer approach, and a formaldehyde abundance lower by a factor higher than 10⁴.

To understand how robust these results are, we considered the whole grid of models described in Sect. 3.8 and built the distribution of the predicted abundances. Figure 8 shows this distribution for key stable species (CO, H₂O, H₂CO, and CH₃OH) and the radicals (OH, HCO, and CH₃O) for the bulk and multilayer approaches for two times (10⁵ and 10⁶ yrs). The plots of Fig. 8 tell us that at short times (10⁵ yr) the predicted species abundance distributions are quite similar in the bulk and multilayer approaches, with the exceptions of H₂CO and CH₃O, for which the distribution are slightly different. In contrast, at long times (10⁶ yr), the distribution of all radicals are substantially different in the two approaches, demonstrating that the age is a key parameter. The plots also provide another important information. For example, considering the multilayer approach results, CO and OH have quite peaked distributions and their predicted abundances depend slightly if at all, on the assumed values of the parameters: they are robust predictions. On the contrary, H₂CO, CH₃OH, HCO and CH₃O have broad abundance distributions implying that their predicted abundances are very sensitive to the values of the other model parameters.

## 4.3. Porous versus non-porous grains

The influence of the grain porosity (Sect. 3.2) is shown in Fig. 7 for the reference set of parameters. The presence of the porosity in the grain has the effect of slightly speeding up and increasing (by less than a factor 4) the formation of formaldehyde, methanol, and CH₃O, as well as to increase the destruction of
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Fractional composition of each mantle monolayer for the reference model (boldface values in Table 3).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Abundance of species in interstellar grain mantles as a function of time for the reference model (solid lines), the bulk approach model (Sect. 4.2; dashed lines), and the porous grain model (Sect. 4.3; dotted lines). Left panel: the most abundant stable species. Middle panel: reactive species. Right panel: less abundant stable species.</figcaption></figure>

HO₂, O₂ and O. The reason is that porous grains trap H atoms more efficiently, increasing their number on the grains, but in the end, the H atom numbers is still limited by their low gas phase abundance (Fig. 7). Furthermore, the porosity only slightly modifies the distribution of the results shown in Fig. A.1. Therefore, the porosity only enhances the formation of hydrogenated molecules by a few factors at most, regardless of the values of other parameters.

## 4.4. Influence of physical conditions

### Density $n_H$

The density plays an important role for two reasons. First, increasing the density reduces the time needed to form each monolayer (because the accretion rate is proportional to the square of the density). Second, the ratio between the number of H atoms and the number of heavy particles (CO and O) that land on the grain is inversely proportional to the density (see Sect. 3.5). As a consequence, the higher the density, the lower the abundance of stable species (H₂O, H₂CO, CH₃OH) created by the hydrogenation reactions, especially at long times

(≥10⁵ yr) (Fig. A.2). In contrast and for the same reasons, mantle CO and OH are more abundant for higher density.

### Temperature $T_g = T_d$

The grain temperature has a moderate influence on the predicted mantle composition, although the diffusion (and thus the reaction) and the desorption rates depend exponentially on it. However, the desorption rate increases much more quickly than the diffusion rate with temperature (by factors of 10²–10⁶ depending on the $E_d/E_b$ ratio). Consequently, with increasing grain temperature, H atoms have a higher probability to desorb into the gas phase before encountering another particle. The final mantle abundance of molecules created by hydrogenation reactions (H₂CO and CH₃OH) is therefore lower at any time at higher temperatures, as shown by Fig. A.3. The mantle abundance of radicals, on the other hand, is not affected by the grain temperature.

### Grain size $a_d$

While the mantle thickness depends roughly quadratically on the grain size, the final mantle composition depends on it
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Distribution of the predicted mantle abundance $X_{\text{mantle}}$ of the key species at $10^5$ (left panels) and $10^6$ yr (right panels). The thin blue and thick red lines refer to the bulk and multilayer approach. The distribution has been built by considering all the ~18 000 runs of the grid (Sect. 3.8).</figcaption></figure>

very little. The reason is that the particle accretion rate (inversely proportional to the grain size square) and the formation time of each layer compensate each other. The result is that the mantle composition does not depend on the assumed grain size, whether 0.1 or 0.3 µm.

### 4.5. Influence of other key parameters of GRAINOBLE

#### Initial abundance of atomic oxygen $X(O)_{\text{ini}}$

Not surprisingly, the initial gas phase abundance of atomic oxygen is an important parameter for the final mantle water abundance. The higher the O abundance, the higher the iced H$_2$O. Less evidently, oxygenation is in competition with hydrogenation. The initial O abundance also affects the final mantle abundance of formaldehyde and, even more, methanol, because the elemental C/O ratio decreases between 0.7 and 0.2 when $X(O)_{\text{ini}}$ increases between $2 \times 10^{-5}$ and $2 \times 10^{-4}$. At $10^5$ yr, the predicted methanol abundance is a factor ten higher when the O abundance is a factor ten lower.

#### Diffusion energy to binding energy ratio $E_d/E_b$

The diffusion energy is an important parameter in the predicted mantle abundance of formaldehyde and methanol.

Because the diffusion rates depend exponentially on the diffusion energy, the mantle abundance of stable species, like formaldehyde and methanol, decreases with increasing $E_d$. Conversely, the diffusion energy does not affect the mantle abundance of radicals, because their formation and destruction rates compensate each other.

#### Activation energy $E_a$

The activation energy of CO and H$_2$CO hydrogenation reactions strongly influences the formation rates of formaldehyde and methanol because their reaction rates depend exponentially on $E_a$. Most of the runs with a low value of $E_a$ (400 K) predict mantle abundances of formaldehyde and methanol higher than $10^{-6}$ at $\ge 10^5$ yr. Conversely, higher values of $E_a$ give uniform distributions of H$_2$CO and CH$_3$OH abundances, between $10^{-18}$ and $10^{-6}$. Only ~5% of the runs using a activation energy equal to 2500 K predict formaldehyde and methanol abundances higher than $10^{-6}$. The final abundances of formaldehyde and methanol, consequently strongly depend on the other model parameters for high $E_a$ values. A significant amount of radicals, like HCO and CH$_3$O, can survive in the mantle only if the activation energies are low. Runs with activation energies of 400 K predict radical mantle abundances between $10^{-10}$ and $10^{-8}$ for both radicals, while runs with $E_a$ = 2500 K predict abundances less than ~$10^{-11}$. Not surprisingly, the activation energy of CO and H$_2$CO hydrogenation reactions are, therefore, critical parameters for the survival of radicals in the mantle.

#### Site size $d_s$

The reaction rates and the formation time of each monolayer are both functions of the site size $d_s$: the former is proportional to $d_s^2$, whereas the latter is proportional to $d_s^{-2}$. Therefore, the two processes, which are in competition, cancel each other, and consequently, the total mantle abundance of stable species and radicals do not sensitively depend on the size of sites.

However, the thickness of the grain mantle strongly depends on the site size. Indeed, the maximum number of monolayers is given by the ratio between the number of particles on the grain $N_{\text{part}}$ and the number of grain sites $N_s$ (proportional to $d_s^{-2}$). Furthermore, if we assume that the thickness of a layer is equal to the distance between two sites, the thickness of the mantle is proportional to $d_s^3$. Thus, the bigger are the sites, the faster the grains grow. Consequently, the depletion is more efficiently on grains with big sites, increasing the abundance of methanol and decreasing the abundance of CO for example.

### 4.6. Concluding remarks

Based on the previous paragraphs, we can identify three classes of parameters, depending on their influence on the results:

- The parameters that have a *strong* influence on the mantle composition are the approach to model the chemical behaviour (multilayer versus bulk), the density, the diffusion energy and the activation energy of CO and H$_2$CO hydrogenation reactions. Modifying the values of these parameters within their chosen range drastically changes the distribution of the predicted mantle abundances, both in the shape and in the values.

- The grain porosity and the temperature as well as the initial abundance of atomic oxygen have only a *moderate* impact on the predicted mantle composition.

- The grain and the site sizes have a *negligible* influence on the total mantle composition. However, these parameters strongly modify the ice thickness, which could have an
impact on photolytic processes and for the desorption time of
the mantle during the warm-up phase in protostar envelopes.

**5. Comparisons with previous microscopic models**

The Monte-Carlo continuous time random walk (CTRW) method was introduced by Chang et al. (2005) and Cuppen & Herbst (2005) for the hydrogen recombination system, and was then extended to more complex networks in Chang et al. (2007), Cuppen & Herbst (2007), Cuppen et al. (2009). The natural output of this model is the multilayer structure of the mantle. This section compares our results with the model of Cuppen et al. (2009) with the twofold aim of validating the GRAINOBLE code and of highlighting the differences between the two models.

**5.1. Validation of the GRAINOBLE code**

Cuppen et al. (2009) considered the accretion of H, CO, and H₂ on smooth surfaces as described by Cuppen & Herbst (2005), leading to the formation of formaldehyde, methanol, and their associated radicals HCO and CH₃O via the Langmuir-Hinshelwood and the Eley-Rideal mechanisms. A major difference between GRAINOBLE and the Cuppen et al. model is that they consider an α-CO ice, whereas we assumed that the bulk of the mantle is formed by iced water. To reproduce the interactions of molecules within an α-CO ice, we modified the set of binding energies in our model accordingly. Second, the CTRW model takes into account the individual interaction energies between each particle, while our model only considers the total binding energy between the adsorbate and the substrate. To simulate their approach, we assumed for each species *i* a binding energy equal to fourteen times the interaction between the species *i* and a CO molecule (*E*<sub>b</sub>(*i*) = 14*E*<sub>i-CO</sub>), corresponding to the interaction between an adsorbate and a porous multi-layered CO ice. The binding energy of H is equal to 450 K, while the binding energy of CO is 880 K. The diffusion to desorption energy ratio is assumed to be 0.8, as in Cuppen et al. (2009) for a layer that is initially smooth. The activation energies for CO and H₂CO hydrogenation reactions are those measured by Fuchs et al. (2009), and slightly depend on the grain temperature.

To validate our code with respect to the model of Cuppen
et al. (2009), we run our model with a density of 10⁵ cm⁻³ and
considered the results at 2 × 10⁵ yr, for the following four cases
(as in Part 3.3 of Cuppen et al. 2009): a) grain temperature 12 K,
n(H)/n(CO) = 1 and n<sub>d</sub> = 2 × 10<sup>12</sup>n<sub>h</sub> ; b) as a) but n<sub>d</sub> = 1 ×
10<sup>12</sup>n<sub>h</sub> ; c) as a) but temperature 15 K; d) as b) but n(H)/n(CO) =
0.5. Figure 9 shows the chemical composition of each monolayer
for the four cases, for porous grains (F<sub>por</sub> = 0.7: this is the value
that agrees best for our and Cuppen’s model). The comparisons
between Fig. 4 of Cuppen et al. (2009) and our Fig. 9 shows that
the two codes give very similar results, which validates our code.

**5.2. Differences between the two approaches**

We highlight a few critical differences between the model of
Cuppen and ours.

Unlike our multilayer approach, the CTRW method allows
us to model the segregation effects. Heavier particles can ag-
glomerate together to form islands where they can be fixed
onto the grains more efficiently. Our multilayer method is a
macroscopic approach where segregation effects are not taken into
account. However, the pores can play the same role as the is-
lands, because particles that enter into the pores cannot desorb

<figure><img src="image_10.png" /></figure>

**Fig. 9.** Chemical composition of the mantle layers as a function of the monolayers at a time of 2 × 10⁵ years, for a density *n*<sub>H</sub> = 10⁵ cm⁻³, as in Cuppen et al. (2009) (see text): CO (red), formaldehyde (blue) and methanol (green). The four panels refer to the four cases of Fig. 4 of Cuppen et al. (2009) (see text): **a)** grain temperature 12 K, *n*(H)/*n*(CO) = 1 and *n*<sub>d</sub> = 2 × 10<sup>12</sup>*n*<sub>h</sub>; **b)** as **a**) but *n*<sub>d</sub> = 1 × 10<sup>12</sup>*n*<sub>h</sub>; **c)** as **a)** but temperature 15 K; **d)** as **b)** but *n*(H)/*n*(CO) = 0.5.

either. This is supported by the fact that our porous grains case
agrees best with Cuppen’s model.

An important difference between ours and Cuppen et al.'s model is the inclusion of the water ice formation. Because Cuppen et al. did not consider the formation of water, their ices are formed by pure α-CO ice. Consequently, the interaction energies are lower than those used in our model. Therefore in Cuppen's model ices can form only at temperatures lower than 16–18 K, whereas in our model iced CO can still survive until 21–23 K. In our model, then, the formation of formaldehyde and methanol can occur in a wider range of temperatures than in Cuppen et al. (2009). Consequently, the abundance of formaldehyde and methanol formed on grain mantles decreases much more efficiently with the increase of the temperature in Cuppen et al.'s model than in ours. For example, using the same physical parameters, the CTRW model predicts a decrease of the mantle thickness of 85% (~40 to 6 monolayers), whereas our model does not predict any thickness decreases between 15 and 16.5 K.

**6. Comparisons with observations**

**6.1. *The observations***

In order to compare our model predictions with observations, we used the data pertaining to the ices only, because the abundance of the sublimated mantle species does not necessarily reflect the chemical composition of their precursor ices. First, it can be altered by reactions in the gas phase. Second, the different sublimation temperatures may introduce errors when comparing the abundance ratios of different species, because they may refer to different regions.

Unfortunately, solid CO, H₂CO, and CH₃OH have been ob-
served simultaneously only towards a very small sample of high-
and intermediate-mass protostars (Gibb et al. 2004; Pontoppidan
et al. 2004). To increase the statistics and to include low-mass
<table><caption>**Table 5.** Mantle abundances of CO and CH₃OH with respect to water along with the relative CH₃OH/CO abundance ratio, as observed towards low-(LM), intermediate-(IM), and high-mass (HM) protostars.</caption><thead><tr><th>Source</th><th>X(CO)<br>% wrt H₂O</th><th>X(CH₃OH)<br>% wrt H₂O</th><th>X(CH₃OH)<br>% wrt to CO</th></tr></thead><tbody><tr><td><strong>Low-Mass</strong></td><td></td><td></td><td></td></tr><tr><td>L1448 IRS 1</td><td>45.5<sup>a-e</sup></td><td>&lt;14.9<sup>d</sup></td><td>&lt;32.7</td></tr><tr><td>L1455 SMM 1</td><td>5.6<sup>a-e</sup></td><td>&lt;13.5<sup>d</sup></td><td>&lt;241.1</td></tr><tr><td>RNO 15</td><td>13.6<sup>a-e</sup></td><td>&lt;5.0<sup>d</sup></td><td>&lt;36.8</td></tr><tr><td>IRAS 03254</td><td>13.6<sup>a-e</sup></td><td>&lt;4.6<sup>d</sup></td><td>&lt;33.8</td></tr><tr><td>IRAS 03271</td><td>8.2<sup>a-e</sup></td><td>&lt;5.6<sup>d</sup></td><td>&lt;68.3</td></tr><tr><td>B1-a</td><td>13<sup>a-e</sup></td><td>&lt;1.9<sup>d</sup></td><td>&lt;14.6</td></tr><tr><td>B1-c</td><td>28.6<sup>a-e</sup></td><td>&lt;7.1<sup>d</sup></td><td>&lt;24.8</td></tr><tr><td>L1489</td><td>15.2<sup>a-e</sup></td><td>4.9<sup>d</sup></td><td>32.3</td></tr><tr><td>DG Tau B</td><td>20.5<sup>a-e</sup></td><td>&lt;5.7<sup>d</sup></td><td>&lt;27.8</td></tr><tr><td>IRAS 12553</td><td>12.6<sup>a-e</sup></td><td>&lt;3.0<sup>d</sup></td><td>&lt;23.8</td></tr><tr><td>IRAS 13546</td><td>22.7<sup>a-e</sup></td><td>&lt;3.9<sup>d</sup></td><td>&lt;17.2</td></tr><tr><td>IRAS 15398</td><td>40.3<sup>a-e</sup></td><td>10.3<sup>d</sup></td><td>25.6</td></tr><tr><td>CRBR 2422.8-3423</td><td>11.2<sup>a-e</sup></td><td>&lt;9.3<sup>d</sup></td><td>&lt;83.0</td></tr><tr><td>RNO 91</td><td>19.0<sup>a-e</sup></td><td>&lt;5.6<sup>d</sup></td><td>&lt;29.5</td></tr><tr><td>IRAS 23238</td><td>4.0<sup>a-e</sup></td><td>&lt;3.6<sup>d</sup></td><td>90.0</td></tr><tr><td><strong>Intermediate-Mass</strong></td><td></td><td></td><td></td></tr><tr><td>AFGL989</td><td>19.0<sup>b</sup></td><td>1.7<sup>b</sup></td><td>8.9</td></tr><tr><td>SMM4</td><td>30.0<sup>c</sup></td><td>28.0<sup>c</sup></td><td>93.3</td></tr><tr><td><strong>High-Mass</strong></td><td></td><td></td><td></td></tr><tr><td>W33A</td><td>7.4<sup>a-e</sup></td><td>14.7<sup>d</sup></td><td>198.6</td></tr><tr><td>GL 2136</td><td>10.2<sup>a-e</sup></td><td>8.5<sup>d</sup></td><td>83.3</td></tr><tr><td>S140</td><td>16.6<sup>a-e</sup></td><td>&lt;3.0<sup>d</sup></td><td>&lt;18.1</td></tr><tr><td>NGC 7538 IRS 9</td><td>16.5<sup>a-e</sup></td><td>7.5<sup>d</sup></td><td>45.5</td></tr><tr><td>AFGL2136</td><td>5.0<sup>b</sup></td><td>5.0<sup>b</sup></td><td>100.0</td></tr><tr><td>AFGL7009</td><td>16.0<sup>b</sup></td><td>33.0<sup>b</sup></td><td>206.3</td></tr></tbody></table>

**References.** (a) Pontoppidan et al. (2003); (b) Gibb et al. (2004); (c) Pontoppidan et al. (2004); (d) Boogert et al. (2008); (e) Pontoppidan et al. (2008).

protostars, too, we therefore restricted the comparison of our model to the observed solid CO and CH₃OH abundances only. Fortunately, a posteriori, when considering the model predictions, we found that using the H₂CO would not provide (substantially) more constraints. Table 5 reports the compilation of observations that we used for the comparison.

## 6.2. Constraining the activation and diffusion energies

As a first step, we compared the model predictions with the observations with the goal to constrain the two microphysics parameters, which have a high impact on the model predictions (Sect. 4.6), namely $E_d/E_b$ and $E_a$. Figure 10 shows the CH₃OH/CO abundance ratio as a function of time obtained assuming reference parameters (and more particularly a density $10^5$ cm⁻³) and for different values of the activation energy $E_a$ and diffusion energy $E_d$. The comparison of the model predictions with the values observed towards intermediate- and high-mass protostars clearly excludes an activation energy higher than 1450 K and a diffusion energy higher than 0.65 times the binding energy $E_b$. More stringently, if $E_a = 1450$ K, $E_d/E_b$ has to be equal to 0.5. The analysis of a similar plot for a density of $10^4$ cm⁻³ gives similar constraints: $E_a = 1450$ K implies $E_d/E_b = 0.5$ and $E_a = 400$ K, $E_d/E_b \sim 0.65$. The case $10^6$ cm⁻³ does not provide more constraints, it is indeed excluded by the observations (see below).

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Mantle CH₃OH/CO abundance ratio versus time. The blue, green, and red lines refer to $E_d/E_b = 0.65$, 0.5, and 0.8, respectively. Solid, dotted, and dashed lines refer to $E_a = 400$ K, 1450 K, and 2500 K, respectively. In these computations the density is $10^5$ cm⁻³. The box with hatching shows the interval of CH₃OH/CO abundance ratio observed towards intermediate- and high-mass protostars (Table 5). Because low-mass protostars provide only less stringent upper limits, they are not reported in the plot.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Mantle CH₃OH/CO abundance ratio versus time. The red lines refer to Model 1 ($E_a = 400$ K and $E_d/E_b = 0.65$), the blue lines refer to Model 2 ($E_a = 1450$ K and $E_d/E_b = 0.5$) (see text). Solid, dotted and dashed lines refer to densities $10^4$, $10^5$ and $10^6$ cm⁻³ respectively. The grey box with hatching shows the interval of CH₃OH/CO abundance ratio observed towards high- and intermediate-mass protostars, while the green dashed box shows the values (upper limits) observed in low-mass protostars (Table 5).</figcaption></figure>

## 6.3. Constraining the pre-collapse phase duration and density

Once we limited the range of possible values of $E_a$ and $E_b$ consistent with the observations, we can attempt to constrain the duration of the pre-collapse phase and the density. To this end, Fig. 11 shows the curves of the CH₃OH/CO abundance ratio versus time for three different densities ($10^4$, $10^5$ and $10^6$ cm⁻³) and computed for the two following sets of $E_a$ and $E_d$ values: 1) $E_a = 400$ K and $E_d/E_b = 0.65$ (and the reference parameters); 2) $E_a = 1450$ K and $E_d/E_b = 0.5$ (and the reference parameters). The comparison between the predicted and observed CH₃OH/CO abundance ratio in intermediate- and high-mass protostars suggests that the bulk of methanol has been formed when the pre-collapse condensation had densities between $10^4$ and $10^5$ cm⁻³. Higher densities would produce an insufficiently CH₃OH/CO abundance ratio compared to the values observed in intermediate- and high-mass protostars. If the densities are
around $10^5$ cm$^{-3}$, the duration of the pre-collapse phase must have lasted at least $10^5$ yr. Lower densities would allow shorter pre-collapse duration times. But densities higher than $10^5$ cm$^{-3}$ are not ruled out for the condensations that produce the low-mass protostars. In that case, the duration of the pre-collapse phase lasted more than $5 \times 10^4$ yr.

One important result is that the CH$_3$OH/CO abundance ratio does not provide an estimate of the duration of the pre-collapse phase, but only a lower limit to it. It neither provides an estimate of the density at the time of the collapse but only a range of possible densities for any given observed CH$_3$OH/CO abundance ratio. Note that the CH$_3$OH/H$_2$CO abundance ratio does not allow us to arrive at better constraints, because the theoretical curves are very similar to those of Fig. 11.

## 6.4. Chemical differentiation within grain mantles

As briefly mentioned in Sect. 2, IR observations suggest a chemical differentiation of the grain mantles with a polar matrix formed already at $A_V \sim 3$ mag (Whittet et al. 1988), and a non-polar CO matrix formed at higher densities (Tielens et al. 1991). In between, CO can be mixed either with water, methanol, or both at the same time (Bisschop et al. 2007), as well as mixed with O$_2$, O$_3$, N$_2$, or CO$_2$ (Pontoppidan 2006). Actually, this observational fact is one of the motivations of our work. As discussed in Sect. 4 and, for example, shown in Fig. 6, our code predicts inhomogeneous grain mantle composition. For the reference model, for example, CO is mixed with water in the inner layers while a mixture between CO and methanol is predicted in the outer layers. Given the lack of a dataset to compare our results with, it is difficult to carry out quantitative comparisons between observations and predictions. However, it is worth noticing that the exact grain mantle composition and stratification depends on the overall evolution and structure of the pre-collapse condensation, because the chemical differentiation depends on the density, the temperature and time.

## 7. Conclusions and perspectives

We presented a new model, GRAINOBLE, which computes the chemical composition of interstellar grain mantles during the cold and dense phase of the pre-collapse. The model presents two main differences from most other published codes: 1) it assumes that only the outermost mantle layer is chemically active, while the mantle bulk is not; 2) it considers porous grains.

We run GRAINOBLE for a large set of parameters, which are either unconstrained (density, temperature and grain size of the prestellar condensation) or poorly known (gas atomic oxygen abundance, diffusion energy, hydrogenation reactions activation energy, grain site size). We obtained a grid of about 18 000 models and built distribution plots of the predicted mantle abundances. This allows us to study the influence of each of these parameters on the predicted abundances and, consequently, the robustness of the predictions.

Finally, we compared the GRAINOBLE predictions with those obtained by the microscopic CTRW model of Cuppen et al. (2009), obtaining a fair agreement between the two models. The comparison clearly validates our treatment and our code (Sect. 5). The most important results of the GRAINOBLE model are the following.

- The multilayer treatment shows a differentiation of the species within the mantle, with the innermost layers being rich in CO, the intermediate layers in formaldehyde, and the

outermost layers in methanol, reflecting the different formation time of each species in the mantle. This differentiation will likely lead to a differentiation in the deuteration of formaldehyde and methanol on the ices. Because the deuteration increases with the CO depletion in the gas phase (Roberts et al. 2003; Ceccarelli & Dominik 2005), methanol will likely be more deuterated than formaldehyde, as indeed observed in Parise et al. (2006). A forthcoming paper will focus on the detailed modeling of the deuteration process.

- The multilayer treatment predicts a relatively high abundance of radicals trapped in the mantle. For example, HCO and CH$_3$O show abundances of $10^{-9}-10^{-7}$ with respect to H nuclei. As suggested by Garrod & Herbst (2006), these radicals can react to form COMs when the grain temperature increases during the formation of the protostar. At this stage, we cannot say whether the predicted radical abundances are sufficient to explain the present observations (e.g. Fig. 1), and, unfortunately, comparisons with previous models are not possible because of the lack of specific information. We will explore this aspect, the formation of COMs, in a forthcoming paper.

- The presence of porosity in the grains only moderately influences the mantle chemical composition, causing an enhanced abundance of formaldehyde and methanol by less than a factor four.

- The chemical composition of grain mantles strongly depends on the physical conditions of the prestellar condensation, particularly on its density as well as its age. Therefore, it will be important to model the evolution of the condensation for realistic predictions of the mantle composition. Conversely, the observed mantle composition can provide valuable information on the past history of protostars, and will be the focus of a forthcoming paper. Comparisons of the present predictions with the ices observations (specifically the CH$_3$OH/CO abundance ratio) suggest that intermediate- and high-mass protostars evolved from condensations less dense than about $10^5$ cm$^{-3}$, while no stringent constraints can be put on the prestellar condensations of low-mass protostars.

- The predicted mantle composition critically depends on the values of the diffusion energy and activation energy of the hydrogenation reactions. Comparing our model predictions with observations of ices allows to constrain the values of these two important parameters: the diffusion to binding energy ratio has to be around 0.5-0.65, and the activation energy has to be less than 1450 K.

- Other parameters of the GRAINOBLE model, like the gas atomic oxygen abundance and the site size, do not substantially influence the predicted mantle abundances.

In summary, GRAINOBLE is a versatile and fast code to model the grain surface chemistry. Its main improvement with respect to previous similar models is the treatment of the layer-by-layer chemistry, which provides a more realistic stratified mantle composition. It predicts that the mantles are indeed stratified and that radicals are trapped inside the mantles. Because it is not expensive from the point of view of computing time, modeling of complex and more realistic cases, like the evolution of a prestellar core, are now feasible.

*Acknowledgements.* The authors would like to thank O. Biham, P. Caselli, T. Hasegawa, and V. Wakelam for useful exchanges on grain chemistry modelling, P. Peters and L. Wiesenfeld for discussions about the physical and chemical processes on grain surfaces, and X. Tielens for helpful discussions on the abundances in the gas phase. This work has been supported by l'Agence Nationale pour la Recherche (ANR), France (project FORCOMS, contracts ANR-08-BLAN-022).
**Appendix A: Distributions of abundances**

As introduced in Sect. 3.8, we have built a model grid in which eight input parameters vary. To study the results of this grid, we computed the distributions of selected species abundances on grain mantles $X_m(i)$, namely the number of results giving $X_m(i)$ in each abundance interval. The following figures show the distributions of abundances for each input parameter and for two extremal values, using the multilayer approach only.

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Distribution of abundances on grain mantles at $t = 10^5$ and $10^6$ years for two porosity values: $F_{por} = 0$ (thick red) and 0.9 (narrow blue).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. A.2.** Distribution of abundances on grain mantles at $t = 10^5$ and $10^6$ years for two densities: $n_H = 10^4$ (thick red) and $10^6$ cm$^{-3}$ (narrow blue).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.3.** Distribution of abundances on grain mantles at $t = 10^5$ and $10^6$ years for two temperatures: $T = 10$ (thick red) and 20 K (narrow blue).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. A.4.** Distribution of abundances on grain mantles at $t = 10^5$ and $10^6$ years for two grain sizes: $a_d = 0.1$ (thick red) and 0.3 µm (narrow blue).</figcaption></figure>
<figure><img src="image_17.png" /></figure>

**Fig. A.5.** Distribution of abundances on grain mantles at *t* = 10⁵ and 10⁶ years for two initial oxygen abundances: X ini (O) = 2 × 10⁻⁵ (thick red) and 2 × 10⁻⁴ (narrow blue).

<figure><img src="image_18.png" /></figure>

**Fig. A.7.** Distribution of abundances on grain mantles at *t* = 10⁵ and 10⁶ years for two site sizes: *d* s = 1.4 (thick red) and 7 Å (narrow blue).

<figure><img src="image_19.png" /></figure>

**Fig. A.6.** Distribution of abundances on grain mantles at *t* = 10⁵ and 10⁶ years for two energy ratios: *E* d/*E* b = 0.5 (thick red) and 0.8 (narrow blue).

<figure><img src="image_20.png" /></figure>

**Fig. A.8.** Distribution of abundances on grain mantles at *t* = 10⁵ and 10⁶ years for two CO and H₂CO hydrogenation reactions: *E* a = 400 (thick red) and 2500 K (narrow blue).
## Appendix B: List of symbols used in this work

<table><caption>Table B.1. List of symbols and their explanation.</caption><thead><tr><th>Symbol</th><th>Parameter</th></tr></thead><tbody><tr><td>a<sub>d</sub></td><td>diameter of interstellar grains (in µm)</td></tr><tr><td>d<sub>s</sub></td><td>size of each site (in Å)</td></tr><tr><td>E<sub>d</sub></td><td>barrier energy against diffusion of a considered species (in K)</td></tr><tr><td>E<sub>b</sub></td><td>barrier energy against desorption of a considered species (in K)</td></tr><tr><td>E<sub>a</sub></td><td>activation energy of the CO and H<sub>2</sub>CO hydrogenation reactions (in K)</td></tr><tr><td>F<sub>por</sub></td><td>fraction of the grain surface occupied by the pores</td></tr><tr><td>F<sub>np</sub></td><td>fraction of the grain surface occupied by the non-porous surface</td></tr><tr><td>F<sub>ed</sub></td><td>fraction of the grain surface occupied by the edge porous sites</td></tr><tr><td>r<sub>d</sub></td><td>radius of interstellar grains (in µm)</td></tr><tr><td>n<sub>d</sub></td><td>density of interstellar grains (in cm<sup>-3</sup>)</td></tr><tr><td>N<sub>pore</sub></td><td>number of sites occupied by each pore</td></tr><tr><td>N<sub>s</sub></td><td>number of sites of the considered layer</td></tr><tr><td>n<sub>H</sub></td><td>total density of H nuclei (in cm<sup>-3</sup>)</td></tr><tr><td>n<sub>g</sub></td><td>gas phase density of a considered species (in cm<sup>-3</sup>)</td></tr><tr><td>P<sub>np</sub></td><td>surface population of a considered species in the non-porous surface (in monolayers)</td></tr><tr><td>P<sub>por</sub></td><td>surface population of a considered species in the pores (in monolayers)</td></tr><tr><td>R<sub>acc</sub></td><td>accretion rate of a considered species (in s<sup>-1</sup>)</td></tr><tr><td>R<sub>hop</sub></td><td>hopping rate from site to site of a considered species (in s<sup>-1</sup>)</td></tr><tr><td>R<sub>diff</sub></td><td>diffusion rate of a considered species (in s<sup>-1</sup>)</td></tr><tr><td>R<sub>r np</sub></td><td>reaction rate of a considered reaction in the non-porous surface (in s<sup>-1</sup>)</td></tr><tr><td>R<sub>r por</sub></td><td>reaction rate of a considered reaction in the pores (in s<sup>-1</sup>)</td></tr><tr><td>R<sub>ev</sub></td><td>evaporation rate of a considered species (in s<sup>-1</sup>)</td></tr><tr><td>S</td><td>sticking coefficient of a considered species</td></tr><tr><td>σ<sub>d</sub></td><td>cross section of interstellar grains (in cm<sup>2</sup>)</td></tr><tr><td>T<sub>g</sub></td><td>gas temperature (in K)</td></tr><tr><td>T<sub>d</sub></td><td>dust temperature (in K)</td></tr><tr><td>v</td><td>thermal velocity of a considered species in the gas phase (in cm s<sup>-1</sup>)</td></tr><tr><td>X<sub>d</sub></td><td>abundance of interstellar grains relative to H nuclei</td></tr><tr><td>X<sub>g</sub></td><td>abundance relative to H nuclei of a considered species in the gas phase</td></tr><tr><td>X<sub>m</sub></td><td>abundance relative to H nuclei of a considered species on grain mantles</td></tr><tr><td>ζ</td><td>cosmic rays ionization rate (in s<sup>-1</sup>)</td></tr></tbody></table>

## References

Aikawa, Y., Wakelam, V., Garrod, R. T., & Herbst, E. 2008, ApJ, 674, 984

Allen, M., & Robinson, G. W. 1977, ApJ, 212, 396

Andersson, S., Al-Halabi, A., Kroes, G.-J., & van Dishoeck, E. F. 2006, J. Chem. Phys., 124, 064715

Arasa, C., Andersson, S., Cuppen, H. M., van Dishoeck, E. F., & Kroes, G.-J. 2010, J. Chem. Phys., 132, 184510

Avgul, N. N., & Kiselev, A. V. 1970, Chem. Phys. Carbon, 6, 1

Awad, Z., Viti, S., Collings, M. P., & Williams, D. A. 2010, MNRAS, 407, 2511

Barzel, B., & Biham, O. 2007, J. Chem. Phys., 127, 144703

Bennett, C. J., & Kaiser, R. I. 2007, ApJ, 661, 899

Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339

Beuther, H., Zhang, Q., Bergin, E. A., & K., S. T. 2009, ApJ, 137, 406

Biham, O., Furman, I., Pirronello, V., & Vidali, G. 2001, ApJ, 553, 595

Bisschop, S. E., Fraser, H. J., Öberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297

Bisschop, S. E., Jrgensen, J. K., van Dishoeck, E. F., & de Wachter E. B. M. 2007, A&A, 465, 913

Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621

Boogert, A. C. A., Pontoppidan, K. M., Knez, C., et al. 2008, ApJ, 678, 985

Bottinelli, S., Ceccarelli, C., Lefloch, B., et al. 2004, ApJ, 615, 354

Bottinelli, S., Ceccarelli, C., Williams, J. P., & Lefloch, B. 2007, A&A, 463, 601

Caselli, P., Hasegawa, T. I., & Herbst, E. 1998a, ApJ, 495, 309

Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998b, ApJ, 499, 234

Caux, E., Ceccarelli, C., Castets, A., et al. 1999, A&A, 347, L1

Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 593, L51

Ceccarelli, C., & Dominik, C. 2005, A&A, 440, 583

Chang, Q., Cuppen, H. M., & Herbst, E. 2005, A&A, 434, 599

Chang, Q., Cuppen, H. M., & Herbst, E. 2007, A&A, 469, 973

Charnley, S. B., Tielens, A. G. G. M., & Millar, T. J. 1992, ApJ, 399, L71

Collings, M. P., Dever, J. W., Fraser, H. J., & McCoustra, M. R. S. 2003, Ap&SS, 285, 633

Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133

Cuppen, H. M., & Herbst, E. 2005, MNRAS, 361, 565

Cuppen, H. M., & Herbst, E. 2007, ApJ, 668, 294

Cuppen, H. M., Ioppolo, S., Romanzin, C., & Linnartz, H. 2010, Phys. Chem. Chem. Phys. (Incorporating Faraday Transactions), 12, 12077

Cuppen, H. M., van Dishoeck, E. F., Herbst, E., & Tielens, A. G. G. M. 2009, A&A, 508, 275

Dartois, E., Schutte, W., Geballe, T. R., et al. 1999, A&A, 342, L32

di Francesco, J., Evans, II, N. J., Caselli, P., et al. 2007, Protostars and Planets V, 17

Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, A&A, 512, A30

Ehrlich, G., & Hudda, F. G. 1966, J. Chem. Phys., 44, 1039

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629

Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927

Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103

Garrod, R. T., Weaver, S. L. W., & Herbst, E. 2008, ApJ, 682, 283

Gerakines, P. A., Schutte, W. A., & Ehrenfreund, P. 1996, A&A, 312, 289

Gerakines, P. A., Moore, M. H., & Hudson, R. L. 2000, A&A, 357, 793

Gerakines, P. A., Moore, M. H., & Hudson, R. L. 2001, J. Geophys. Res., 106, 33381

Gibb, E., Nummelin, A., Irvine, W. M., Whittet, D. C. B., & Bergman, P. 2000, ApJ, 545, 309

Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJS, 151, 35

Goumans, T. P. M. 2011, MNRAS, 262

Gredel, R., Lepp, S., Dalgarno, A., & Herbst, E. 1989, ApJ, 347, 289

Green, N. J. B., Toniazzo, T., Pilling, M. J., et al. 2001, A&A, 375, 1111

Hasegawa, T. I., & Herbst, E. 1993a, MNRAS, 261, 83

Hasegawa, T. I., & Herbst, E. 1993b, MNRAS, 263, 589

Hasegawa, T. I., Herbst, E., & Leung, M. C. 1992, Ap&SS, 82, 167

Hinshelwood, C. N. 1940, The Kinetics Of Chemical Change (The Clarendon Press)

Hollenbach, D. & Salpeter, E. E. 1970, ApJ, 163, 155

Horn, A., Møllendal, H., Sekiguchi, O., et al. 2004, ApJ, 611, 605

Ikeda, M., Ohishi, M., Nummelin, A., et al. 2001, ApJ, 560, 792

Ioppolo, S., Cuppen, H. M., Romanzin, C., van Dishoeck, E. F., & Linnartz, H. 2010, Phys. Chem. Chem. Phys. (Incorporating Faraday Transactions), 12, 12065
Jaycock, M. J., & Parfitt, G. D. 1986, Chemistry of Interfaces (New York: Wiley and Sons)

Jenniskens, P., Blake, D. F., Wilson, M. A., & Pohorille, A. 1995, ApJ, 455, 389

Jones, A. P. 2011, A&A, 528, A98

Katz, N., Furman, I., Biham, O., Pirronello, V., & Vidali, G. 1999, ApJ, 522, 505

Klemm, R. B., Payne, W. A., & Stief, L. J. 1975, in Chemical kinetic data for the upper and lower atmosphere, ed. S. W. Benson (New York: Wiley)

Laas, J. C., Garrod, R. T., Herbst, E., & Widicus Weaver, S. L. 2011, ApJ, 728, 71

Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147

Lipshtat, A., & Biham, O. 2003, A&A, 400, 585

Lis, D. C., Keene, J., Phillips, T. G., et al. 2001, ApJ, 561, 823

MacDonald, G. H., Gibb, A. G., Habing, R. J., & Millar, T. J. 1996, A&AS, 119, 333

Matar, E., Congiu, E., Dulieu, F., Momeni, A., & Lemaire, J. L. 2008, A&A, 492, L17

Mathis, J. S. 1996, ApJ, 472, 643

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Mehringer, D. M., & Snyder, L. E. 1996, ApJ, 471, 897

Nummelin, A., Bergman, P., Hjalmarson, A., et al. 2000, ApJS, 128, 213

Öberg, K. I., Linnartz, H., Visser, R., & van Dishoeck, E. F. 2009a, ApJ, 693, 1209

Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009b, A&A, 496, 281

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

Ossenkopf, V. 1993, A&A, 280, 617

Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619

Parise, B., Ceccarelli, C., Tielens, A. G. G. M., et al. 2006, A&A, 453, 949

Perets, H. B., & Biham, O. 2006, MNRAS, 365, 801

Perets, H. B., Biham, O., Manicó, G., et al. 2005, ApJ, 627, 850

Pontoppidan, K. M. 2006, A&A, 453, L47

Pontoppidan, K. M., Fraser, H. J., Dartois, E., et al. 2003, A&A, 408, 981

Pontoppidan, K. M., van Dishoeck, E. F., & Dartois, E. 2004, A&A, 426, 925

Pontoppidan, K. M., Boogert, A. C. A., Fraser, H. J., et al. 2008, ApJ, 678, 1005

Qin, S., Wu, Y., Huang, M., et al. 2010, ApJ, 711, 399

Remijan, A., Shiao, Y., Friedel, D. N., Meier, D. S., & Snyder, L. E. 2004, ApJ, 617, 384

Requena-Torres, M. A., Martín-Pintado, J., Rodríguez-Franco, A., et al. 2006, A&A, 455, 971

Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41

Shen, C. J., Greenberg, J. M., Schutte, W. A., & van Dishoeck, E. F. 2004, A&A, 415, 203

Sosa, A., & Schlegel, B. H. 1986, Int. J. Quant. Chem., 29, 1001

Speedy, R. J., Debenedetti, P. G., Scott Smith, R., Huang, C., & Kay, B. D. 1996, J. Chem. Phys., 105, 240

Steinacker, J., Pagani, L., Bacmann, A., & Guieu, S. 2010, A&A, 511, A9

Stepnik, B., Abergel, A., Bernard, J., et al. 2003, A&A, 398, 551

Sternberg, A., Dalgarno, A., & Lepp, S. 1987, ApJ, 320, 676

Sutton, E. C., Peng, R., Danchi, W. C., et al. 1995, Ap&SS, 97, 455

Tielens, A. G. G. M. 2005, Phys. Chem. Interstellar Medium (Cambridge University Press)

Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245

Tielens, A. G. G. M., & Allamandola, L. J. 1987, in Interstellar Processes, ed. D. J. Hollenbach, & H. A. Thronson, 397

Tielens, A. G. G. M., Tokunaga, A. T., Geballe, T. R., & Baas, F. 1991, ApJ, 381, 181

Ulbricht, H., Moos, G., & Hertel, T. 2002, Phys. Rev. B, 66

Vastel, C., Caux, E., Ceccarelli, C., et al. 2000, A&A, 357, 994

Walch, S. P. 1993, J. Chem. Phys., 98, 3076

Watanabe, N., & Kouchi, A. 2002, ApJ, 567, 651

Whittet, D. C. B., Bode, M. F., Longmore, A. J., et al. 1988, MNRAS, 233, 321

Woon, D. E. 2002, ApJ, 569, 541