© ESO, 2012

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^-10 by considering formation pathways only in the gas phase, but this molecule has been detected with abundances up to 10^-6 (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 oxydation 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 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^-7, 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 the 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 modifiying 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 areconsidered 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 σ(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_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.

• 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 ≥ 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: ${\mathit{F}}_{\mathrm{ed}}\mathrm{=}\mathrm{4}\mathrm{\times }{\mathit{F}}_{\mathrm{por}}\frac{\sqrt{{\mathit{N}}_{\mathrm{pore}}}\mathrm{-}\mathrm{1}}{{\mathit{N}}_{\mathrm{pore}}}\mathrm{·}$(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 + Δt and a new number of particles on the grain is computed at t + Δ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 ${\mathit{d}}_{\mathrm{s}}^{\mathrm{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 intestellar 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₂, or methanol. Molecules such as NH₃ and CH₄ 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 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 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. 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: $\mathit{n}\mathrm{\left(}\mathit{H}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2.3}{\mathit{\zeta }}_{\mathrm{CR}}\mathit{n}\mathrm{\left(}{\mathit{H}}_{\mathrm{2}}\mathrm{\right)}}{\mathrm{2}\mathit{v}\mathrm{\left(}\mathit{H}\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{d}}{\mathit{X}}_{\mathrm{d}}\mathit{n}\mathrm{\left(}{\mathit{H}}_{\mathrm{2}}\mathrm{\right)}}$(2)where v(H) is thermal velocity of atomic hydrogen. The resulting H density is independent of the density and equal to 1.2 cm^-3 at 10 K, for a cosmic rays ionization rate ζ of 3 × 10^-17 (Caselli et al. 1998b), a grain abundance X_d of 1.33 × 10^-12 relative to H nuclei and for a grain radius r_d of 0.1 μm (${\mathit{\sigma }}_{\mathrm{d}}\mathrm{=}\mathit{\pi }{\mathit{r}}_{\mathrm{d}}^{\mathrm{2}}$), giving n_dσ_d/n_H = 4.2 × 10^-22 cm². 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 × 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 × 10^-5 or less. The initial abundance of O is therefore considered as a free parameter and can vary between 2 × 10^-4 and 2 × 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  ~10⁴ to  ~10⁶ 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 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^-2, depends on the structure of the surface and its composition. Most astrochemical models assume a value of s = 10¹⁵ cm^-2 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 (${\mathit{d}}_{\mathrm{s}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathit{s}}$), and is equal to 3.2 Å. However, Biham et al. (2001) found a value of s = 5 × 10¹³ cm^-2 (d_s = 14 Å) for an amorphous carbon surface and Perets & Biham (2006) assumed a value of s = 5 × 10¹⁵ cm^-3 (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 litterature (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: $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}& \mathrm{-}& \mathit{S}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{v}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{\sigma }\mathrm{\left(}{\mathit{a}}_{\mathrm{d}}\mathrm{\right)}{\mathit{n}}_{\mathrm{d}}\mathrm{\left(}{\mathit{a}}_{\mathrm{d}}\mathrm{\right)}\\ & \mathrm{+}& \\ \frac{\mathrm{d}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}& & \frac{\mathrm{1}}{\mathrm{4}}\mathit{S}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{d}}_{\mathrm{s}}^{\mathrm{2}}\mathit{v}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{-}{\mathit{R}}_{\mathrm{ev}}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\\ & \mathrm{-}& \frac{{\mathit{F}}_{\mathrm{ed}}}{{\mathit{F}}_{\mathrm{np}}}{\mathit{R}}_{\mathrm{hop}}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{+}\frac{{\mathit{F}}_{\mathrm{ed}}}{{\mathit{F}}_{\mathrm{por}}}{\mathit{R}}_{\mathrm{hop}}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{P}}_{\mathrm{por}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\\ & \mathrm{+}& \sum _{{\mathit{i}}_{\mathit{f}}}{\mathit{R}}_{\mathrm{r}\mathit{,}\mathrm{ext}}\mathrm{\left(}{\mathit{i}}_{\mathit{f}}\mathrm{\right)}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}{\mathit{i}}_{\mathit{r}\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}{\mathit{i}}_{\mathit{r}\mathrm{2}}\mathrm{\right)}\\ & \mathrm{-}& \\ \frac{\mathrm{d}{\mathit{P}}_{\mathrm{por}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}& \mathrm{+}& \frac{{\mathit{F}}_{\mathrm{ed}}}{{\mathit{F}}_{\mathrm{np}}}{\mathit{R}}_{\mathrm{hop}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{·}{\mathit{P}}_{\mathrm{np}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{-}\frac{{\mathit{F}}_{\mathrm{ed}}}{{\mathit{F}}_{\mathrm{por}}}{\mathit{R}}_{\mathrm{hop}}\mathrm{\left(}\mathit{i}\mathrm{\right)}{\mathit{P}}_{\mathrm{por}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\\ & \mathrm{+}& \sum _{{\mathit{i}}_{\mathit{f}}}{\mathit{R}}_{\mathrm{r}\mathit{,}\mathrm{in}}\mathrm{\left(}{\mathit{i}}_{\mathit{f}}\mathrm{\right)}{\mathit{P}}_{\mathrm{por}}\mathrm{\left(}{\mathit{i}}_{\mathit{r}\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathrm{por}}\mathrm{\left(}{\mathit{i}}_{\mathit{r}\mathrm{2}}\mathrm{\right)}\\ & \mathrm{-}& \end{array}$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 17 496 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 × 10⁵ 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₃ (ozone), OH and H₂O₂ (hydrogen peroxyde). Indeed, the high O to H initial gas phase abundance ratio (~6) and the relatively high temperature (15 K) allow a significant amout 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 × 10⁴ 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₃O, reach abundances between  ~5  ×  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 × 10⁵ 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⁵ yr because they have a fast formation rate (each layer is formed in less than  ~2 × 10³ yr). The formation time of layers increases sharply at  ~10⁵ 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. Indeed, “intermediate” molecules such as O₃, 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₂O, remain relatively constant throughout the mantle up to the last outermost layers, where 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 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 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)_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₂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)_ini increases between 2 × 10^-5 and 2 × 10^-4. At 10⁵ 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₂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  ≥ 10⁵ yr. Conversely, higher values of E_a give uniform distributions of H₂CO and CH₃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₃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₂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 ${\mathit{d}}_{\mathrm{s}}^{\mathrm{2}}$, whereas the latter is proportional to ${\mathit{d}}_{\mathrm{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_part and the number of grain sites N_s (proportional to ${\mathit{d}}_{\mathrm{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 ${\mathit{d}}_{\mathrm{s}}^{\mathrm{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 mantlecomposition are the approach to model the chemical behaviour(multilayer versus bulk), the density, the diffusion energy and theactivation energy of CO and H₂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_b(i) = 14E_i − CO), 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^-3 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_d = 2  ×  10¹²n_h ; b) as a) but n_d = 1  ×  10¹²n_h; 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_por = 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 agglomerate together to form islands where they can be fixed onto the grains more efficiently. Our multilayer method is a macroscpic approach where segration effects are not taken into account. However, the pores can play the same role as the islands, because particles that enter into the pores cannot desorb 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 necessarly 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 observed 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 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⁵ cm^-3) 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⁴ cm^-3 gives similar constraints: E_a = 1450 K implies E_d/E_b = 0.5 and E_a = 400 K, E_d/E_b ~ 0.65. The case 10⁶ cm^-3 does not provide more constraints, it is indeed excluded by the observations (see below).

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⁴, 10⁵ and 10⁶ cm^-3) 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⁴ and 10⁵ cm^-3. 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⁵ cm^-3, the duration of the pre-collapse phase must have lasted at least 10⁵ yr. Lower densities would allow shorter pre-collapse duration times. But densities higher than 10⁵ 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 × 10⁴ yr.

One important result is that the CH₃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₃OH/CO abundance ratio. Note that the CH₃OH/H₂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 ~ 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₂, O₃, N₂, or CO₂ (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 thespecies within the mantle, with the innermost layers beingrich in CO, the intermediate layers in formaldehyde, and theoutermost layers in methanol, reflecting the different for-mation time of each species in the mantle. This differentia-tion will likely lead to a differentiation in the deuteration offormaldehyde and methanol on the ices. Because the deuter-ation increases with the CO depletion in the gas phase (Robertset al. 2003; Ceccarelli& Dominik 2005), methanol will likely be moredeuterated than formaldehyde, as indeed observed in Pariseet al. (2006). A forthcoming paperwill 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₃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₃OH/CO abundance ratio) suggest that intermediate – and high – mass protostars evolved from condensations less dense than about 10⁵ 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.

Acknowledgments

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).

References

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.

Appendix B: List of symbols used in this work

All Tables

All Figures