© ESO, 2015

1. Introduction

Dust grains play a vital role in the making and breaking of molecules in interstellar gas clouds because many chemical reactions proceed much faster on solid surfaces than reactions in the gas phase. As a result, dust can enhance molecular abundances in interstellar gas clouds by acting as a catalyst for the formation of complex molecules (e.g., recent studies by Meijerink et al. 2012; Le Bourlot et al. 2012; Gavilan et al. 2012). However, dust grains can also lock up gas-phase molecules through freeze-out in cold environments (e.g., Jones & Williams 1984; Lippok et al. 2013). This duality creates an intricate relationship between gas and dust particles. Their involvement in chemical reactions affects the thermodynamic properties of molecular clouds and from this their whole evolution (Goldsmith 2001; Cazaux et al. 2010). The kinetic temperature of a gas cloud is especially sensitive to changes in abundances of the dominant coolants, like CO (Hocuk et al. 2014). These findings highlight the importance of considering the formation of ices and depletion of gas-phase species in models and theories of cloud evolution and star formation.

Ices form on dust grains in the interstellar medium (ISM) during the evolution of interstellar clouds, the progenitors of star-forming regions. Diffuse clouds, where the dust-gas coupling and grain surface chemistry is still negligible, will evolve and undergo phase transitions to form molecular clouds, where freeze-out is effective, and eventually form dense clumps where dust-gas coupling becomes dominant. These are the critical phases that clouds undergo before a star finally forms and which, ultimately, determine the stellar masses at birth. During the evolutionary stages, gas-phase molecules are deposited on grain surfaces to form icy layers, thereby depleting the gaseous molecules. Eventually, the ices will grow thick mantels on dust grains and will remain there until the medium becomes warmer, allowing the species to evaporate back into the gas phase. Heating can be caused by radiation (UV photons, X-rays, or cosmic rays) penetrating the cloud and interacting with the gas, shock-waves injecting kinetic heating, or gravitational collapse where compressional heating ensues and radiation becomes trapped. Hot cores already revealed a glimpse of the rich organic chemistry locked up in the ice mantles (Cazaux et al. 2003), which suddenly became visible as a result of evaporation into the gas phase by protostellar heating.

Depletion of molecular species, such as CO, are also observed toward the formation of prestellar cores (Tafalla et al. 2006; Liu et al. 2013), with depletion factors of up to 80 in dense clumps (Fontani et al. 2012), indicating the presence of thick ice mantles. Observations of cold prestellar cores reported the lack of gas-phase H₂O, demonstrating a much more serious freeze-out than previously predicted (van Dishoeck et al. 2011). The presence of frozen water was corroborated by the detection of strong water emission from shocks in protostellar environments (van Dishoeck et al. 2011; Suutarinen et al. 2014). Recent observations of starburst galaxy M82 showed that CO₂ and H₂O ices are present for different physical conditions (Yamagishi et al. 2013). This is intriguing because it is demanding for models to predict high abundances of CO₂. The pathway of forming CO₂ is still very uncertain (Ioppolo et al. 2013), but the gas-phase route is widely accepted as inefficient. This leads to the notion that CO ice is not passive.

Observations of deeply embedded protostars show an anticorrelation between the abundance of CO (depleted) and methanol in the gas (Kristensen et al. 2010). This demonstrates that CO converts into methanol on the surface of dust grains and is subsequently released into the gas phase. The observed ortho-to-para ratio of water in the Orion Bar is also found to be inconsistent with only gas-phase water formation (Choi et al. 2013), indicating that nonthermal processes must be at work on dust grains. An important mechanism to supply the gas phase with ice species is the nonthermal process of photodesorption from either UV or cosmic-ray-induced UV photons. Photons with energies of mostly around 8 eV can directly desorb CO molecules into the gas phase (Fayolle et al. 2011). Photodesorption is considered as a possible important gas-phase supplier of CO (Muñoz Caro et al. 2010), methanol and formaldehyde (Yusef-Zadeh et al. 2013; Guzmán et al. 2013), water (Keto et al. 2014), CO₂ (Öberg et al. 2009; Bahr & Baragiola 2012), or O₂ (Zhen & Linnartz 2014). Observations by Yıldız et al. (2013) suggest, however, that photodesorption rates need to be increased by a factor of 2 to explain their abundances of O₂.

Recent experimental studies have unveiled another nonthermal mechanism, called chemical desorption, that directly converts species formed on dust surfaces into gas-phase species (Garrod et al. 2007; Dulieu et al. 2013). This mechanism occurs for exothermic reactions, where the products cannot thermalize with the dust surface. These findings show that the formation of species through surface reactions does not just lock up species from the gas phase, but will also directly enrich the gas-phase medium and is therefore an integral part of interstellar cloud evolution.

In this work we track the formation of ices during the evolution of an interstellar gas cloud, starting from a diffuse, fully atomic stage until the formation of dense clumps. We observe the formation of the first ice layers on the surfaces of dust grains, report their compositions, and determine the distribution of ices around a dense clump as well as their formation rates. For this purpose, we developed a chemical network using rate equations that incorporates grain surface reactions on two different substrates, bare grain (no ices) surface, and water ice surface. We also consider the processes of chemical and photodesorption (Cazaux et al. 2010; Noble et al. 2012a; Dulieu et al. 2013) and the most recent reactions through quantum tunneling (Oba et al. 2012; Minissale et al. 2013b). The paper is organized as follows: In Sect. 2 we present the code that is used in this work and describe our initial conditions. In Sect. 3 we explain the chemical processes on dust surfaces, describe each process separately, and present our equations. In Sect. 4 we outline the important thermal processes and report the heating and cooling terms used in our calculations. In Sect. 5 we present our results on the composition of the ice layers, show the dominant species formation rates, and give the distribution of ices around a dense clump. We also discuss the implications of our results. In Sect. 6 we conclude and discuss the caveats.

2. Numerical method

2.1. Numerical code

We performed the numerical simulations with the adaptive-mesh hydrodynamical code flash, version 4.0 (Fryxell et al. 2000; Dubey et al. 2009). Our work encompasses a broad range of physics, such as hydrodynamics, chemistry, thermodynamics (using time-dependent heating and cooling rates), turbulence, multispecies, gravity, and radiative transfer for UV. To solve the hydrodynamic equations, we applied the directionally split piecewise-parabolic method (PPM; Colella & Woodward 1984), which is well-suited to handling the type of calculations in this study. Our research captures the physics that act on small and large scales. These are the micro-sized scales required for the chemical reactions in the gas phase and on grain surfaces and the parsec-sized scales that are necessary to describe star formation processes.

To track multiple fluids, we employed the multispecies unit that is able to follow each species with its own properties. We applied the consistent multifluid advection scheme (Plewa & Müller 1999) to prevent overshoots in the mass fractions as a result of the PPM advection. The Poisson equations were solved with the Multigrid solver in which gravity is coupled to the Euler equations through the momentum and energy equations. The physics modules are well-tested and were either provided by flash or can be found in earlier works, for example, Hocuk & Spaans (2010, 2011) and Hocuk et al. (2014). A comprehensive chemistry and a thermodynamics module were created specifically for this work; they are explained in more detail in Sects. 3 and 4.

2.2. Initial conditions

We created a gravitationally bound diffuse gas cloud in which all hydrogen is in atomic form. Our model cloud starts with a uniform number density of n_H = 10 cm^-3 and an initial temperature of 100 K. n_H is defined as the total hydrogen nuclei number density, that is, n_H = ρ/m_H. The interstellar environment including the cloud surface has a temperature of 1000 K and a number density of n_H = 1 cm^-3. We allowed for a smooth, hyperbolic transition for the values from the cloud edge to the surrounding ISM. To follow its chemistry and temperature during its evolution, we placed our model cloud in a 3D cubic box of size 150 pc³. We applied periodic boundary conditions to our simulation domain. The spherical cloud has a radius of 42 pc and a total mass of 7.2 × 10⁴M_⊙. A graphical display of the initial conditions is given in Fig. 1.

The diffuse cloud was initiated with turbulent conditions that are representative of the ISM in the Milky Way, σ_turb = 1 km s^-1, with a power spectrum of P(k) ∝ k^-4 following the empirical laws for compressible fluids (Larson 1981; Myers & Gammie 1999; Heyer & Brunt 2004). This scaling is also known as Burgers turbulence (Burgers 1939; Bec & Khanin 2007). The turbulence is decaying and not driven. The cloud is also not supported by turbulence and will contract. The simulated cloud was placed in an environment with a background UV radiation flux of G₀ = 1 in terms of the Habing field (Habing 1968). This agrees with the ISM conditions of our Milky Way. The cloud center enjoys column densities of over 10²¹ cm^-3, which is equivalent to a visual extinction of A_V ~ 0.5. We do not consider magnetic fields in this work.

We refined the computational grid that encloses the model cloud to a uniform resolution of 128³ cells. This yields a spatial resolution of 1.17 pc, which corresponds to a Jeans resolution of 44 grid cells per Jeans length for the initial state. The Jeans length is calculated as $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{J}}\mathrm{=}{\left(\frac{\mathit{\pi }{\mathrm{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{G}\mathit{\rho }}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{cm}\mathit{.}& & \end{array}$(1)Here c_s is the sound speed, which, for an ideal gas, can be formulated as ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma P}\mathit{/}\mathit{\rho }\mathrm{=}\mathit{\gamma }{\mathit{N}}_{\mathrm{A}}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathit{/}\mathit{\mu }$. In this, N_A and k_B are the Avogadro number and the Boltzmann constant, μ is the mean molecular mass, which is unity in the initial case, and the parameter γ depends on the equation of state (EOS). For a polytropic EOS, the pressure scales as P ∝ ρ^γ (Spaans & Silk 2000). The polytropic index γ will be affected by the thermal balance, that is, heating and cooling, of the cloud and therefore is not a fixed quantity. The sound speed of the diffuse cloud if it is isothermal (γ = 1) is c_s = 0.91 km s^-1.

We did not increase grid resolution during the simulation nor allowed for star formation to occur. The spatial resolution will therefore not be as high as contemporary advanced numerical studies of star formation, but our goal is not to resolve the fine details of cloud fragmentation, turbulence, or the direct spatial conditions prior to star formation. Our focus lies on describing the chemical composition (gas + dust) and thermal balance of an evolving gas cloud.

Because the formation of ices can take more than 10⁴ years (Cuppen & Herbst 2007) and is dependent on the changing conditions within the cloud, an equilibrium solution is not acceptable. Within our time-dependent solver, we also did not allow for variations in species densities of over 5% per iteration. To this end, we employed a very high time resolution with adaptive time-stepping that can extend to a time resolution of three days (0.01 yr) for a simulation that lasts over ten million years. In the most unfavorable case, it is necessary to iterate the chemistry routine more than 10 000 times within a hydrodynamical time step. The adaptive chemical time-stepping is handled within a subcycling loop of the hydrodynamical time step to avoid any unnecessary speed loss for the other routines. We let the cloud evolve until it reached 125% of its theoretical free-fall time, where ${\mathit{t}}_{\mathrm{ff}}\mathrm{=}\sqrt{\mathrm{3}\mathit{\pi }\mathit{/}\mathrm{32}\mathrm{G}\mathit{\rho }}$. Given the initial conditions of the model cloud, this adds up to a final simulation time of 2.0 × 10⁷ yr.

3. Analytical method: time-dependent chemistry

We performed time-dependent rate equations at each grid cell that include both gas and grain surface reactions to compute the chemical composition of the cloud. Our chemical model comprises 42 different species. Of these, 26 are in the gas phase and 16 are on dust grains. Our selection of species and their initial abundances (with respect to hydrogen atoms) is given in Table 1.

There are 257 relevant reactions. Gas-phase reactions with the corresponding equations and their coefficients were obtained from the Kinetic Database for Astrochemistry (KiDA; Wakelam et al. 2012). Surface reactions on dust grains were acquired from Cazaux et al. (2010) and comprise 89 reactions (see Appendix A).

3.1. Chemistry solver

Chemical reaction rates are solved using a fast and stable semi-implicit scheme, an improved scheme over the first-order backwards differencing (BDF) method developed by Anninos et al. (1997). The derivation of this scheme is given below.

The general expression for implicitly solving ordinary differential rate equations is defined as $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{C}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{-}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}\mathit{,}& & \end{array}$(2)where C_xi is the creation and D_xi is the destruction rate of species x_i at future, t + dt, time step. The first-order integration of this differential equation is also known as the backward Euler method. It differs from a forward Euler method in which C_xi and D_xi would be based on current, known values. Knowing that D_xi implicitly depends on the species x_i we can submit ${\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{\equiv }{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }}$, with ${\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }}$ the destruction rate coefficient (Schleicher et al. 2008). Following this, the semi-implicit scheme (Anninos et al. 1997) is obtained by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \frac{{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{new}}\mathrm{-}{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{old}}}{\mathrm{\Delta }\mathit{t}}\mathrm{=}{\mathit{C}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{-}{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{new}}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }}\mathit{,}\\ {\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{new}}& \mathrm{=}& \frac{{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{old}}\mathrm{+}{\mathit{C}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{\Delta }\mathit{t}}{\mathrm{1}\mathrm{+}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }}\mathrm{\Delta }\mathit{t}}\mathit{,}\end{array}$(3)where old refers to the values at the current time step and new points to the values at t + Δt, the future time step. Because the new rates, C_xi and ${\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }}$, are also not known at the current time step, which we ideally wish to find, they can be approximated using current, old values through iterating, applying a predictor-corrector method, or by using a mix between old and new values with the values that were previously evaluated. The latter does somewhat depend on the order in which the rates are calculated.

We devised a second-order variant of the semi-implicit method (Eq. (3)) to solve our equations. We applied the trapezoidal rule to integrate our differential equation to gain higher precision. In its derived form, the second-order semi-implicit scheme is presented as $\begin{array}{ccc}{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}^{\mathit{j}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }\mathit{j}}\mathrm{\Delta }\mathit{t}\right)\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{C}}_{{\mathit{x}}_{\mathit{i}}}^{\mathit{j}}\mathrm{+}{{\mathit{C}}_{{\mathit{x}}_{\mathit{i}}}^{\mathit{j}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}\mathrm{\Delta }\mathit{t}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{D}}_{{\mathit{x}}_{\mathit{i}}}^{\mathrm{\prime }\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{\Delta }\mathit{t}}\mathit{,}& & \end{array}$(4)where we have replaced the time step descriptor by j and j + 1 (previously, old and new). This scheme is easily computed while still being an improvement over the explicit schemes and the first-order semi-implicit scheme. The order of error does not drop with respect to a second-order fully implicit method, that is, the local error is O(h³), the global error O(h²). Semi-implicit schemes are symplectic integrators in nature and yield far better results than standard Euler methods. We note, however, that for all semi-implicit BDF methods, one must ensure mass conservation, and they are not time reversible.

A higher order BDF method will be more accurate for solving the chemistry, as proven by Bovino et al. (2013), but will also come at a higher computational price. Employing a high time resolution with this scheme is, in our experience, sufficient to have good accuracy while still performing the calculations at an acceptable speed.

3.2. Gas-phase chemistry

Gas-phase reactions were acquired from the Kinetic Database for Astrochemistry (Wakelam et al. 2012). We considered every possible reaction (within the scope of the database) that involves our selection of species as given in Table 1.

The 168 gas-phase reactions in our network include bimolecular reactions (e.g., A + B → C + D), charge-exchange reactions (e.g, A⁺ + B → A + B⁺), radiative associations (e.g., A + B → AB + photon), associative detachment (e.g., A⁻ + B → AB + e⁻), electronic recombination and attachment (e.g., AB⁺ + e⁻→ A + B), ionization or dissociation of neutral species by UV photons, ionization or dissociation of species by direct collision with cosmic-ray particles or by secondary UV photons following H₂ excitation.

3.3. Dust chemistry

The solid-phase reactions or the reaction rate coefficients were gathered from Cazaux et al. (2010). We have included several additional reactions in this work. The grain surface rate equations were simplified and set in easily accessible forms. To do this, we unified all the different reactions involving dust grains into five equation types. These are (A) adsorption of gas-phase species onto dust surfaces; (B) thermal desorption of ices; (C) two-body reactions on grain surfaces including chemical desorption; (D) cosmic-ray processes; and (E) photo-processes with UV photons that include photodissociation and photodesorption. The equations are explained in detail in the following five subsections.

3.3.1. Adsorption onto dust grains

Species in the gas phase can be accreted onto grain surfaces. This depends on the motion of the gas species relative to the dust particle. Since the motions are dominated by thermal velocity, the adsorption rate depends on the square root of the gas temperature as ${\mathit{v}}_{\mathrm{th}}\mathrm{=}\sqrt{\mathrm{8}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{g}}\mathit{/}\mathit{\pi m}}$. Once the gas species is in contact with the dust, there is a probability for it to stick on the surface of the grain. The sticking coefficient is calculated as $\begin{array}{ccc}\mathit{S}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\left(\mathrm{1}\mathrm{+}\mathrm{0.4}{\left(\frac{{\mathit{T}}_{\mathrm{g}}\mathrm{+}{\mathit{T}}_{\mathrm{d}}}{\mathrm{100}}\right)}^{\mathrm{0.5}}\mathrm{+}\mathrm{0.2}\frac{{\mathit{T}}_{\mathrm{g}}}{\mathrm{100}}\mathrm{+}\mathrm{0.08}{\left(\frac{{\mathit{T}}_{\mathrm{g}}}{\mathrm{100}}\right)}^{\mathrm{2}}\right)}^{-1}\mathit{,}& & \end{array}$(5)where T_g is the gas temperature and T_d is the dust temperature (Hollenbach & McKee 1979). We note that this coefficient is based on H atoms. Using this, the adsorption rate can be formulated as, $\begin{array}{ccc}{\mathit{k}}_{\mathrm{ads}}& \mathrm{=}& {\mathit{n}}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}{\mathit{v}}_{{\mathit{x}}_{\mathit{i}}}{\mathrm{s}}^{-1}\mathit{,}\\ {\mathit{R}}_{\mathrm{ads}}& \mathrm{=}& {\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{k}}_{\mathrm{ads}}\mathit{S}\mathrm{\left(}\mathit{T}\mathrm{\right)}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\end{array}$(6)where k_ads is the adsorption rate coefficient, R_ads is the adsorption rate, n_d is the number density of dust grains, σ_d is the cross section of the grain, v_xi is the thermal velocity of species x_i, that is, ${\mathit{v}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{=}\sqrt{\mathrm{8}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{g}}\mathit{/}\mathit{\pi }{\mathit{m}}_{{\mathit{x}}_{\mathit{i}}}}$, with m_xi the mass in grams, and n_xi is the number density of the engaging species. In this equation, n_dσ_d represents the total cross section of dust in cm^-1, which is obtained by integrating over the grain-size distribution. We adopted the grain-size distribution of Mathis et al. (1977), from here on MRN, with a value of α_MRN = ⟨ n_dσ_d/n_H ⟩ _MRN = 1 × 10^-21 cm². We chose this distribution rather than the one of Weingartner & Draine (2001), from here on WD, which has a total cross section larger by a factor three, that is, α_WD = 3 × 10^-21 cm², because MRN did not include poly-cyclic aromatic hydrocarbons (PAHs). The freeze-out of species on PAHs to form ices is not known.

3.3.2. Thermal desorption, evaporation

After species are bound on grain surfaces, they can evaporate back into the gas. The evaporation rate depends exponentially on the dust temperature and on the binding energies of the species with the substrate. The binding energy of each species differs according to the type of substrate. We consider two possible substrates in this study, bare surfaces (assuming carbon substrate) and water ice substrate, since ices are mostly made of water. Species adsorbed on water ice substrate have in most cases binding energies higher than on bare dust or other ices (e.g., see Cuppen & Herbst 2007). Species on top of CO, which can attain a significant coverage on dust, have binding energies that more closely resemble the binding energies of bare dust (e.g., Sandford & Allamandola 1988; Karssemeijer et al. 2014). Therefore, in this study we consider the binding energies on CO to be similar to those of bare dust.

We computed the fraction of the dust covered by (water) ice ℱ_ice and bare ℱ_bare to distinguish between the two. Together with the deposited amount of water ice, this fraction depends on the total number of possible attachable sites on grain surfaces per cubic cm of space, designated as n_dn_sites, which is defined as $\begin{array}{ccc}{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}\mathrm{=}{\mathit{n}}_{\mathrm{d}}\frac{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{d}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{pp}}^{\mathrm{2}}}\mathrm{=}{\mathit{n}}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\frac{\mathrm{4}}{{\mathit{a}}_{\mathrm{pp}}^{\mathrm{2}}}\mathrm{\simeq }\mathrm{4.44}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-6}{\mathit{n}}_{\mathrm{H}}\mathrm{c}{\mathrm{m}}^{-3}\hspace{0.17em}\mathrm{ml}{\mathrm{y}}^{-1}\mathit{,}& & \end{array}$(7)where the radius of dust is given by r_d and the typical separation between two adsorption sites on a grain surface is given by a_pp, which we assume to be 3 Å. A full monolayer (mly) is reached when all the possible sites on a grain surface are occupied by an atom or a molecule. To convert this from monolayers to number densities, one needs to multiply by n_dn_sites. For the dust-to-gas mass ratio ϵ_d, intrinsic to the dust number density n_d, we assumed the typical value of ϵ_d = 0.01.

We obtained the fractions, ℱ_ice and ℱ_bare, in the following manner: if the grain surface is covered by less than 1 mly of water ice, $\begin{array}{ccc}{\mathrm{ℱ}}_{\mathrm{ice}}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{\perp }{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}}\mathrm{·}& & \end{array}$(8)When the grain is covered by more than 1 layer of water ice, ℱ_ice = 1. The bare fraction of the dust is obtained by $\begin{array}{ccc}{\mathrm{ℱ}}_{\mathrm{bare}}\mathrm{=}\mathrm{1}\mathrm{-}{\mathrm{ℱ}}_{\mathrm{ice}}\mathit{.}& & \end{array}$(9)We note that at this stage we assumed that the bound species are homogeneously distributed and neglected the cases where stratified layers of species can form. This assumption can lead to an over- or underestimation of reaction rates because the species abundances on different substrates can differ. We expect, however, that since the situation ℱ_bare ~ ℱ_ice is not very common, the over- or underestimation will be marginal. With these definitions, we can formulate the evaporation rate as follows: $\begin{array}{ccc}{\mathit{k}}_{\mathrm{evap}}& \mathrm{=}& {\mathit{\nu }}_{\mathrm{0}}\left({\mathrm{ℱ}}_{\mathrm{bare}}\exp \left(\mathrm{-}\frac{{\mathit{E}}_{\mathrm{bare}\mathit{,i}}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}{\mathrm{ℱ}}_{\mathrm{ice}}\exp \left(\mathrm{-}\frac{{\mathit{E}}_{\mathrm{ice}\mathit{,i}}}{{\mathit{T}}_{\mathrm{d}}}\right)\right){\mathrm{s}}^{-1}\mathit{,}\\ {\mathit{R}}_{\mathrm{evap}}& \mathrm{=}& {\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{k}}_{\mathrm{evap}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\end{array}$(10)where k_evap is the evaporation rate coefficient, R_evap is the evaporation rate, ν₀ is the oscillation frequency, which is typically 10¹² s^-1 for physisorbed species, E_bare,i and E_ice,i are the binding energies of species x_i on bare grains and ices. The species specific binding energies can be found in Table 2.

3.3.3. Two-body reactions on dust grains

While species are attached to grain surfaces, they can move around by thermal diffusion and meet other species with which they can react to form new molecules. The mobility of the species depends on the oscillation frequency ν₀. We only consider physisorbed species in this work, except for neutral hydrogen, where we also take chemisorption into account. In addition to the mobility, the reaction rate depends on the specific binding energy with the substrate and on the dust temperature. When two species meet, they can immediately react if there is no (or little) reaction barrier. If the reaction barrier is high, however, it might be crossed by tunneling. The probability of overcoming the reaction barrier by tunneling is given by $\begin{array}{ccc}{\mathit{P}}_{\mathrm{reac}}\mathrm{=}\exp \left(\mathrm{-}\mathrm{(}\mathrm{2}\mathit{a}\mathit{/}\hslash \mathrm{)}\sqrt{\mathrm{2}\hspace{0.17em}{\mathit{m}}_{\mathrm{red}}\hspace{0.17em}{\mathit{k}}_{\mathrm{B}}\hspace{0.17em}{\mathit{E}}_{\mathrm{a}}}\hspace{0.17em}\right)\mathit{,}& & \end{array}$(11)where a = 1 Å, m_red is the reduced mass of the two engaging species, that is, m_red = (m_i × m_j)/(m_i + m_j), h is the Planck constant, and E_a is the energy of the barrier needed for the reaction to occur. The probability is P_reac = 1 if there is no barrier for the reaction to take place. We did not consider the reaction diffusion competition (Garrod & Pauly 2011). If this were included, the tunneling probabilities would become much higher.

For CO hydrogenation reactions that have a barrier, that is, ⊥H+⊥CO and ⊥H+⊥H₂CO, we deduced an “effective” tunneling barrier energy from the reaction rates obtained by Atkinson et al. (2004) and Fuchs et al. (2009) at low (~10 K) temperatures. The activation barriers we calculated are 600 K and 400 K, respectively. The height of these barriers are verified from recent experimental results (Minissale et al. priv. comm.). More precise measurements are being performed. The used activation barriers can be found in Appendix A.

When a reaction occurs, the product either remains on the surface or immediately desorbs into the gas phase because of the exothermicity of the reaction. The probabilities of desorption are given by δ_bare and δ_ice for our two substrates. The fraction that remains on the surface will be 1−δ. Non-exothermic reactions that do not desorb are by definition mutiplied by 1. With this information, we formulate the two-body reaction rate on grain surfaces as $\begin{array}{ccc}{\mathit{k}}_{\mathrm{2}\mathrm{body}}& \mathrm{=}& {\mathit{\nu }}_{\mathrm{0}}{\mathrm{ℱ}}_{\mathrm{bare}}\left(\exp \left(\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\frac{{\mathit{E}}_{\mathrm{bare}\mathit{,i}}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}\exp \left(\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\frac{{\mathit{E}}_{\mathrm{bare}\mathit{,j}}}{{\mathit{T}}_{\mathrm{d}}}\right)\right)\mathrm{\left(}\mathrm{1}\mathrm{-}{\mathit{\delta }}_{\mathrm{bare}}\mathrm{\right)}\\ & & \mathrm{+}{\mathit{\nu }}_{\mathrm{0}}{\mathrm{ℱ}}_{\mathrm{ice}}\left(\exp \left(\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\frac{{\mathit{E}}_{\mathrm{ice}\mathit{,i}}}{{\mathit{T}}_{\mathrm{d}}}\right)\mathrm{+}\exp \left(\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\frac{{\mathit{E}}_{\mathrm{ice}\mathit{,j}}}{{\mathit{T}}_{\mathrm{d}}}\right)\right)\mathrm{\left(}\mathrm{1}\mathrm{-}{\mathit{\delta }}_{\mathrm{ice}}\mathrm{\right)}\mathit{,}\\ {\mathit{R}}_{\mathrm{2}\mathrm{body}}& \mathrm{=}& \frac{{\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{n}}_{{\mathit{x}}_{\mathit{j}}}}{{\mathit{n}}_{\mathrm{d}}{\mathit{n}}_{\mathrm{sites}}}{\mathit{P}}_{\mathrm{reac}}{\mathit{k}}_{\mathrm{2}\mathrm{body}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\end{array}$(12)where k_2body is the two-body rate coefficient, R_2body is the two-body reaction rate. The exponent in this equation represents the diffusion of species on the surface, and we assume that diffusion occurs with a barrier of two-thirds of the binding energy (67% Dulieu et al. 2013; 40% Collings et al. 2003; 90% Barzel & Biham 2007). The desorption rate is obtained from the complement of 1δ of the rate coefficient.

The desorption probabilities δ_bare for the exothermic reactions ⊥ H + ⊥ O → OH, ⊥ H + ⊥ OH → H₂O, ⊥ O + ⊥ O → O₂ were adopted from Dulieu et al. (2013) and are δ_bare = 0.5,0.9,0.6. For the hydrogenation reactions ⊥ HCO + ⊥ H → H₂CO and ⊥ CH₃O + ⊥ H → CH₃OH, we assumed a desorption fraction of δ_bare = 5% and for ⊥ H₂CO + ⊥ H → CH₃Oδ_bare = 50%. In our next paper, we will report the different chemical desorption yields on different types of surfaces, estimated by laboratory experiments (Cazaux et al., in prep.). For the desorption probabilities on icy substrates δ_ice, we considered that chemical desorption is much weaker than on bare surfaces and assumed that δ_ice = δ_bare/ 5 (deduced from Dulieu et al. 2013). We did not take into account that chemical desorption has a dependence on surface coverage (Minissale & Dulieu 2014).

3.3.4. Cosmic-ray processes on grain surfaces

Cosmic-ray reaction rates on grain surfaces are assumed to be the same as the rates found in the gas phase. We included several of the cosmic-ray reactions for surface bound species. Cosmic-ray processes are usually inefficient destruction mechanisms, but can dominate the destruction rates deep inside the cloud. These reaction rates depend on the cosmic-ray ionization rate per H₂ molecule, which we adopted as ζ_H2 = 5 × 10^-17 s^-1 (Indriolo et al. 2007; Hocuk & Spaans 2011; Chaparro Molano & Kamp 2012). The cosmic-ray reaction rate is formulated as $\begin{array}{ccc}{\mathit{k}}_{\mathrm{CR}}& \mathrm{=}& {\mathit{z}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{\zeta }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathrm{s}}^{-1}\mathit{,}\\ {\mathit{R}}_{\mathrm{CR}}& \mathrm{=}& {\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{k}}_{\mathrm{CR}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\end{array}$(13)where k_CR is the cosmic-ray rate coefficient, R_CR is the cosmic-ray reaction rate, and z_xi the cosmic-ray ionization rate factor that is subject to the ionizing element (see KiDA database).

Cosmic-ray-induced UV (CRUV) photons were also considered within the same equation (Eq. (13)). In this case, z_xi is replaced by z_CRUV, the UV photon generation rate per cosmic-ray ionization. UV photons from cosmic rays do not suffer from radiation attenuation as normal UV photons do. Hence, the lack of dependence on optical depth in this formula.

3.3.5. Photo-processes on grain surfaces

When UV photons arrive on a dust particle, they can interact with the adsorbed species and either photodissociate or photoevaporate them. We used the same formula for both types of photo-processes. Photoreactions scale linearly with the local radiation flux (erg cm^-2 s^-1). The radiation field strength is necessarily a function of extinction, which is given by ξ_xiA_V, where ξ_xi is the extinction factor that is contingent on the relevant species. We obtain A_V by dividing the column density N_H over the scaling factor, that is, A_V = N_H/2.21 × 10²¹ mag (Güver & Özel 2009). The column densities were computed by integrating the density from the simulation boundaries to each point. Simply, this is N_H = Σ_in_H ds_i, with ds_i being the path length of the smallest resolution element in which the density remains constant. For this purpose, we constructed a ray-tracing algorithm using 14 equally weighted rays with long characteristics (traveling from inside to outside). The ray separation is set by our highest resolution. We assumed an isotropic UV radiation field with a flux of 1 G₀, where G₀ ≡ 1.6 × 10^-3 erg cm^-2 s^-1. The general photo-process rate equations are defined as $\begin{array}{ccc}{\mathit{k}}_{\mathrm{phot}}& \mathrm{=}& {\mathit{\alpha }}_{{\mathit{x}}_{\mathit{i}}}{\mathrm{e}}^{\mathrm{-}{\mathit{\xi }}_{{\mathit{x}}_{\mathit{i}}}{\mathit{A}}_{\mathit{V}}}{\mathrm{s}}^{-1}\mathit{,}\\ {\mathit{R}}_{\mathrm{phot}}& \mathrm{=}& {\mathit{n}}_{{\mathit{x}}_{\mathit{i}}}{\mathit{f}}_{\mathrm{ss}}{\mathit{k}}_{\mathrm{phot}}{\mathit{F}}_{\mathrm{UV}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}\end{array}$(14)where k_phot is the photo-process rate coefficient, R_phot is the photo-process reaction rate, α_xi is the unattenuated rate coefficient, f_ss is the self-shielding factor, and F_UV is the UV flux in units of 1.71 G₀, that is, G₀ = 1 gives F_UV = 0.58. The factor 1.71 arises from the conversion from the often used Draine field (Draine 1978) to the Habing field for the far-ultraviolet (FUV) intensity. We used the same α_xi, ξ_xi, and f_ss for the gas phase as we did for the surface reactions.

When UV photons arrive on an icy surface with multiple layers, we only allowed the first two layers to be penetrated by UV photons. As shown by Andersson et al. (2006), Arasa et al. (2010), and Muñoz Caro et al. (2010), only the uppermost few layers contribute to photodesorption. Photodissociation does seem to occur deeper into the ice, but trapping and recombination of species tend to dominate (Andersson & van Dishoeck 2008). This means that the highest number density that the photons can see is n_xi = min(n_xi,2n_dn_sites). This restriction is also enforced for reactions with CRUV photons, as given in Sect. 3.3.4.

Self-shielding, denoted as f_ss, was taken into account for H₂ and CO molecules. These molecules can shield the medium against photo-processes on grain surfaces as well as for species in the gas phase. The self-shielding factor for H₂ was obtained from Draine & Bertoldi (1996), Eq. (37), which is formulated as $\begin{array}{ccc}{\mathit{f}}_{\mathrm{ss}}\mathrm{=}\frac{\mathrm{0.965}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}\mathit{/}{\mathit{b}}_{\mathrm{5}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{0.035}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}{\mathrm{)}}^{\mathrm{0.5}}}\mathrm{\times }\exp \mathrm{[}\mathrm{-}\mathrm{8.5}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-4}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}{\mathrm{)}}^{\mathrm{0.5}}\mathrm{]}\mathit{,}& & \end{array}$(15)where x ≡ N_H2/ 5 × 10¹⁴ cm^-2, b₅ ≡ b/ 10⁵ cm s^-1, and b is the line broadening of H₂ lines, which we take as 3 km s^-1. This factor is only a function of H₂ column density N_H2. We used the column density algorithm to compute the H₂ column for this purpose. For CO molecules, self-shielding was achieved by incorporating the self-shielding tables from Visser et al. (2009) into our code. Given an H₂ column and a CO column, the factor f_ss was acquired from the tables. For all other species we took f_ss = 1.

Photodesorption is only implemented here for ⊥CO molecules. Since CO molecules in the gas phase can strongly affect the thermal balance of collapsing molecular clouds, as we have shown in an earlier work (Hocuk et al. 2014), detailed processes were taken into account to correctly determine gas-phase abundances of CO. This encompasses both the normal and the CRUV photons. Photodesorption is not expected to be the dominant destruction mechanism of ⊥CO, but is a route to desorb surface-bound CO molecules directly into the gas phase. The two implemented reactions are $\begin{array}{ccc}\mathrm{\perp }\mathrm{CO}& \mathrm{+}& \mathrm{UV}\hspace{0.17em}\mathrm{photon}\mathrm{\to }\mathrm{CO}\mathit{,}\\ \mathrm{\perp }\mathrm{CO}& \mathrm{+}& \mathrm{CRUV}\hspace{0.17em}\mathrm{photon}\mathrm{\to }\mathrm{CO}\mathit{.}\end{array}$For these two photo-chemical reactions α_xi of Eq. (14) represents the number of CO photodesorptions per second per unit radiation flux (Hollenbach et al. 2009; Chaparro Molano & Kamp 2012). This variable depends on the photodesorption yield, which was experimentally obtained and adopted by us as Y_CO = 1.0 × 10^-2 (Fayolle et al. 2011). With respect to previous values (Öberg et al. 2009), these yields are higher by about a factor 4. In addition, only the first two monolayers were assumed to be penetrable by UV photons for these rates.

4. Analytical method: thermodynamics

To address the thermodynamics of the gas cloud, we calculated time-dependent heating and cooling rates that complemented the chemistry calculations. In this way, we obtained the gas and dust temperatures by solving the thermal balance.

4.1. Heating and cooling

We included the most prominent heating processes that are relevant to our work. The non-equilibrium heating processes include

• photoelectric heating;

• H₂ photodissociation heating;

• H₂ collisional de-excitation heating;

• cosmic-ray heating;

• gas-grain collisional heating (when T_gas<T_dust).

Compressional heating and shock heating are by default taken into account through the hydrodynamics, the EOS, and the shock-detection routines of our code, flash.

The cooling of the gas is ensured by several different types of non-equilibrium processes. These processes are

• electron recombination with PAHs cooling;

• electron impact with H (i.e., Ly-α) cooling;

• metastable transition of [OI]-630 nm cooling;

• fine-structure line cooling of [OI]-63 μm and [CII]-158 μm;

• molecular ro-vibrational cooling by H₂, CO, OH, and H₂O;

• and gas-grain collisional cooling (when T_gas>T_dust).

Cooling by adiabatic expansion was, again, handled by the standard EOS routines of flash. The heating and cooling functions and their rates are described in Hocuk et al. (2014).

4.2. Dust temperature

The dust temperature is a crucial parameter that not only influences the gas temperature through the heating and the cooling rates, but also affects the chemical reaction rates. This eventually drives the formation and build-up of ices. We follow Hollenbach et al. (1991), also mentioned in Latif et al. (2012), but with some adaptations of our own. The initial estimate of the dust temperature, T_d,i, which incorporates the attenuated incident UV radiation field, the cosmic microwave background temperature T_CMB = 2.725 K, and the infrared dust emission, is defined as $\begin{array}{ccc}{\mathit{T}}_{\mathrm{d}\mathit{,i}}^{\mathrm{5}}& \mathrm{=}& \mathrm{8.9}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-11}{\mathit{\nu }}_{\mathrm{0}}{\mathit{G}}_{\mathrm{0}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathrm{1.8}{\mathit{A}}_{\mathit{V}}\mathrm{)}\mathrm{+}{\mathit{T}}_{\mathrm{CMB}}^{\mathrm{5}}\mathrm{+}\mathrm{3.4}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-2}\\ & & \hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{[}\mathrm{0.42}\mathrm{-}\ln \mathrm{(}\mathrm{3.5}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-2}{\mathit{\tau }}_{\mathrm{100}}{\mathit{T}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{]}}{\mathit{\tau }}_{\mathrm{100}}{\mathit{T}}_{\mathrm{0}}^{\mathrm{6}}\mathit{,}\end{array}$(16)where the adopted value ν₀ = 3 × 10¹⁵ s^-1 represents the most efficient absorbing frequency over the visual and UV wavelengths, τ₁₀₀ is the emission optical depth at 100 μm, and T₀ is the equilibrium dust temperature at the cloud surface due to unattenuated incident FUV field alone (Hollenbach et al. 1991). T₀ in this case equates to ${\mathit{T}}_{\mathrm{0}}\mathrm{=}\mathrm{12.17}{\mathit{G}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathrm{5}}$ K. If it is assumed that the incident FUV flux equals the outgoing flux of dust radiation from T₀, then ${\mathit{\tau }}_{\mathrm{100}}\mathrm{=}\mathrm{2.7}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{3}}{\mathit{G}}_{\mathrm{0}}{\mathit{T}}_{\mathrm{0}}^{-5}$ (Hollenbach et al. 1991). Knowing T₀ fixes τ₁₀₀ to a value of 0.001. This makes it independent of optical depth. We display this by plotting the dust temperature as a function of visual extinction in Fig. 2, which also shows other dust temperature calculations.

To accommodate for the depth dependence of the infrared emission, we calculated τ₁₀₀ in a different fashion. Since A_V = 2.5log ₁₀(e) τ_V ≃ 1.086 τ_V, where τ_V is the opacity at optical (predominantly 550 nm) wavelengths, and because $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{\lambda }}\mathrm{=}{\mathit{\tau }}_{\mathrm{100}}{\left(\frac{\mathrm{100}\mathit{\mu }\mathrm{m}}{\mathit{\lambda }}\right)}^{{\mathit{\beta }}_{\mathrm{sed}}}\mathit{,}& & \end{array}$(17)with λ being the wavelength and β_sed the sub-mm slope of the spectral energy distribution known as the spectral emissivity index, we can rewrite τ₁₀₀ as a function of A_V in the form $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{100}}\mathrm{=}{\mathit{\tau }}_{\mathrm{550}\mathrm{nm}}{\left(\frac{\mathrm{550}\mathrm{nm}}{\mathrm{100}\mathit{\mu }\mathrm{m}}\right)}^{\mathrm{1.3}}\mathrm{=}\mathrm{9.4}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-4}{\mathit{A}}_{\mathit{V}}\mathit{.}& & \end{array}$(18)For β_sed we take the value of 1.3, which gives the same τ₁₀₀ at A_V ≃ 1 as Hollenbach et al. (1991) advocated. A spectral emissivity index between β_sed = 1−2 is typically found for the Milky Way (Miettinen et al. 2012; Arab et al. 2012). However, we now have a higher value of τ₁₀₀ at higher A_V to account for the depth dependence.

We also considered the heating of dust grains by the gas through the gas-grain collisional heat exchange. Since dust grains have a larger heat capacity, the heating of dust grains will be considerably weaker than the cooling of the gas, but this process might still slightly increase the dust temperature. Assuming that there is a temperature equilibrium in which the equilibrium timescale is much shorter than a free-fall time, the energy balance to reach a stable dust temperature is given by $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{gg}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{rad}}\mathit{,}& & \end{array}$(19)where Λ_gg is the gas-grain collisional heating and Γ_rad is the radiative losses due to blackbody radiation. Note that these are the losses from raising the dust temperature to a higher value by gas-grain collisional heat exchange than the equilibrium temperature given in Eq. (16). The losses can be described as $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{rad}}\mathrm{=}\mathrm{4}{\mathit{\sigma }}_{\mathrm{SB}}\mathrm{(}{\mathit{T}}_{\mathrm{d}}^{\mathrm{4}}\mathrm{-}{{\mathit{T}}_{\mathrm{d}\mathit{,i}}^{\mathrm{4}}}^{\mathrm{)}}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\kappa }}_{\mathrm{P}}{\mathrm{erg}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}\hspace{0.17em}\mathrm{s}}^{-1}\mathit{,}& & \end{array}$(20)with σ_SB = 5.67 × 10^-5 erg cm^-2 s^-1 K^-4 the Stefan-Boltzmann constant, ρ_d the mass density of dust, and κ_P the Planck mean opacity. Here, we also make use of our initial estimate of the dust temperature T_d,i and expect T_d,i to be at a stable equilibrium dust temperature when there is no additional heating source. We adopt the relation $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{P}}\mathrm{=}\mathrm{4}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-4}{\mathit{T}}_{\mathrm{d}}^{\mathrm{2}}{\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\mathrm{g}}^{-1}& & \end{array}$(21)as presented by Omukai (2000) for a typical molecular cloud composition (Pollack et al. 1994) and for T_d ≲ 50 K, but we did not consider Rosseland mean opacity in our solution (Stamatellos et al. 2007; Dopcke et al. 2011). Combining Eqs. (19) and (20), we obtain $\begin{array}{ccc}{\mathit{T}}_{\mathrm{d}}^{\mathrm{4}}\mathrm{=}\frac{{\mathrm{\Lambda }}_{\mathrm{gg}}}{\mathrm{4}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\kappa }}_{\mathrm{P}}}\mathrm{+}{\mathit{T}}_{\mathrm{d}\mathit{,i}}^{\mathrm{4}}\mathit{.}& & \end{array}$(22)Replacing κ_P with Eq. (21) and by using ρ_d = ϵ_dρ for the dust mass density, with ϵ_d = 0.01, we can formulate $\begin{array}{ccc}{\mathit{T}}_{\mathrm{d}}^{\mathrm{6}}\mathrm{=}\frac{{\mathrm{\Lambda }}_{\mathrm{gg}}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}\mathrm{\ast }}\mathrm{\right)}}{\mathrm{9.072}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-10}\mathit{\rho }}\mathrm{+}{\mathit{T}}_{\mathrm{d}\mathit{,i}}^{\mathrm{4}}{\mathit{T}}_{\mathrm{d}\mathrm{\ast }}^{\mathrm{2}}\mathit{,}& & \end{array}$(23)where T_{d ∗} is a prediction for the real dust temperature. This so-called sextic equation is a transcendental equation and unsolvable analytically for T_{d ∗} = T_d, but we can solve it numerically by way of iteration while taking the first-order estimate of the dust temperature as T_{d ∗} = T_d,i.

In the end, we have a slightly higher, more accurate dust temperature than our initial estimate that increases with density as a result of better gas-grain coupling, and with optical depth for infrared radiation. To highlight the difference in dust temperature by our adjustments, the three different dust temperatures are displayed together in Fig. 2.

5. Results

When we started our simulation, the turbulence created a filamentary structure of the gas inside the cloud. It slowly evolved into a more clumpy structure with dense clumps that were gravitationally unstable. Ten million years after the initial phases, we chose the densest part of our cloud, a collapsing clump of radius r_clump = 1 pc to follow its evolution. The density of our clump grew from n_H ~ 10² to 10^4.5 cm^-3 during this time. The clump was located close to the cloud center. In Fig. 3 we display this clump after nearly 13 Myr of simulation.

To investigate how the chemistry adapts, we recorded the species abundances, their formation rates, and the thermodynamic quantities in time and in space. The results we report as a function of time are based on the mean value of the cells located inside the chosen clump. We also present plots of the clump conditions pertaining to a region of space at a fixed time. In this case, we show the growth of ice layers as a function of visual extinction and present maps in Cartesian coordinates. To also expose the outer regions, we then examine a greater radial distance, with r_map = 5 pc.

5.1. Evolution of the clump

The density of the clump after 10 million years of cloud evolution is n_H = 200 cm^-3, corresponding to an A_V of near unity. This increases to a 3 × 10⁴ cm^-3 after another 10 million years, see Fig. 4. Both the gas temperature and the dust temperature of the clump within the same time interval (10–20 Myr) initially decrease as the clump becomes denser, and they rise as a result of compressional heating after rapid collapse sets in. At densities of n_H = 2 × 10⁴ cm^-3 and above, gas-dust coupling allows the two physical states to enjoy the same temperature, with T ≳ 12 K. This occurs around an A_V of 20. The optical depth is high because the clump is deeply embedded inside the cloud, close to its center. The collapse of our clump is not delayed, nor does it fragment during the course of the simulation.

5.2. Time evolution of species abundances

We plot in Fig. 5 the time evolution of the species abundances in the gas phase and on grain surfaces. These are the results of the abundances within our 1 pc clump.

Gas-phase species abundances In the upper panel of Fig. 5, where we plot the gas-phase species, the amount of ionized carbon decreases with passing time and as the cloud becomes denser, while in the meantime the neutral carbon abundance increases. Carbon is readily converted into CO when the cloud enters a translucent stage at the clump density of n_H = 10³ cm^-3 and at A_V = 3. CO peaks at t = 1.6 × 10⁷ yr, with a peak abundance of n_CO/n_H = 1.28 × 10^-4, which amounts to 99% of all the carbon. This occurs at a density of n_H = 3 × 10³ cm^-3 while the visual extinction has reached A_V = 5. The only significant way to deplete CO from the gas phase beyond this point is through CO freeze-out on grain surfaces. We can see this happening by the decrease in CO abundance after t = 1.65 × 10⁷ yr. A full layer of water ice has covered the dust surface at this point. See Fig. 5 bottom panel where the black (water) line crosses the dotted line. The CO abundance drops to ~2 × 10^-5 near the end of our simulation. We perceive that it continues to drop to ~2 × 10^-6 at t = 2.2 × 10⁷ yr, which extends beyond the plot range. Oxygen is also increasingly depleted from the gas phase when the dust is enveloped by a mantle of water ice. Neutral oxygen has a steeper decline than CO, with an abundance of n_O/n_H = 1.5 × 10^-6 at t = 2 × 10⁷ yr. This means that oxygen is also being locked up in something other than CO, which we know to be mostly water ice.

Within our densest clump, we never reach a situation with a C/O abundance ratio of unity. A C/O ratio above unity might lead to interesting carbon chemistry in the gas phase and is attained by models using varying local interstellar background radiation fields, G₀ = 1−1000 (Hollenbach et al. 2009). We remind that we use G₀ = 1. The ratio we obtain always remains below unity and consistently below 0.01 above t = 1.55 × 10⁷ yr, when more than half of the dust surface is covered.

The gas-phase abundances of water (black), methanol (light blue), and formaldehyde (dark blue) rise above a value of 10^-7 after 15 Myr (Fig. 5), but remain relatively constant with abundances that, in the same order, linger around 3.2, 4.7, 7.5 × 10^-7 at the end of our simulation. The H₂CO/CH₃OH ratio at that moment is 1.6, which is a consequence of the high chemical desorption of formaldehyde. The abundances of species inside the clump after 20 Myr of cloud evolution are listed in Table 3.

We reach a high amount of methanol ice at the end of our simulation. This is partially because we did not incorporate species larger than methanol in our network. In addition to this, we have not considered diffusion through tunneling for oxygen atoms (Minissale et al. 2013a), which should lead to a lower methanol and a higher CO₂ abundance.

Ice species abundances The lower panel of Fig. 5 displays the abundance of frozen species in terms of ice layers covering the surface of the dust, which are given in units of monolayers along the Y-axis. Dust grains within the clump grow thicker ice mantles as a function of time, which starts to level off around t = 1.8 × 10⁷ yr. The maximum number of ice layers is mainly limited by the total surface area and the availability of the oxygen atoms, n_O/n_H = 2.9 × 10^-4 (Table 1). We reach a total of 59 ice layers at the end of our simulation. This means that eventually 97% of the oxygen freezes out on dust and that only 3% resides in the gas phase locked in CO. The majority of ices are in the form of ⊥H₂O, ⊥CO, ⊥H₂CO (formaldehyde), or ⊥CH₃OH (methanol).

When the cloud is still in a diffuse stage early on in cloud evolution, ⊥CO is the predominant surface-bound species, with frozen water only second to ⊥CO. This means that the initial ice mantle is well mixed. This contradicts the idea that the first ice layer should consist of water ice alone. The ⊥CO/⊥H₂O ratio is about unity when half the dust surface is covered. This ratio decreases quickly over time due to the repeated hydrogenation of ⊥CO to eventually form formaldehyde and methanol. The amount of ⊥CO steeply decreases after the first ice layer has formed, while water ice grows more rapidly with increasing density. CO and water each dominate the surface coverage at different epochs, but water ice always dominates the ice composition when more than one layer of ice covers the dust. The dominance of solid CO on surfaces lasts for about 3 Myr in our simulation, which is only a fraction of a cloud dynamical time, whereafter water ice becomes the dominant ice species. We expect that frozen CO can only be observed mixed with ices.

As stated earlier, the first mly of ice is formed at t = 1.65 × 10⁷ yr after we start our simulation. More importantly than time, the first ice layer is formed when the cloud is dense and cold enough, at which point it enters a molecular stage. Before that, we find that solid species start to grow rapidly, that is, beyond 0.1 mly, above a density of n_H = 360 cm^-3, with n_H2/n_HI = 2, and below a gas temperature of T_g = 16 K. The cloud environment where these frozen species form has an A_V ≥ 1.8 with a dust temperature of 7.5 K. Water ice becomes the main constituent of the icy surface when the gas density rises above n_H = 2.0 × 10³ and the gas temperature falls below 11.5 K. The visual extinction of our clump now reaches upto 5.2. This also marks the turnover point where gas-phase CO starts to become depleted and the gas and dust temperatures start rising. More than two-thirds of the dust surface is covered by ice at this time. This in turn shifts the adsorption energies of the species from bare surface to icy surface binding energies, see Table 2. In most cases, the binding energies increase when the substrate becomes water ice. These changes affect the thermal evaporation and the reaction rates on dust surfaces, as can be seen by Eqs. (10) and (12). Where the former one causes a build-up of ices by inhibiting evaporation, the latter slows down the formation of species by reducing the mobility. Moreover, the chemical desorption probabilities also drop by a factor of 5 when the dust is covered by ice, adding to the ice build-up. The combined effects work in favor of constructing more ices. This is shown by the sharp rise in the curves in Fig. 5 around the point where one mly of ice has formed.

5.3. Time evolution of species formation rates

By following the formation rates of the species in our model, we can understand the preferred formation pathways for each of our species. Here we report the formation rates of the species H₂O, CO, H₂CO, and CH₃OH during collapse. In the following images, Figs. 6 (gas species) and 7 (ice species), we plot the main formation rates for each of these species.

Gas-phase formation rates In the top left panel of Fig. 6, we see that gas-phase water has two main formation channels. One through the reaction of H₃O⁺ with e⁻ in the gas phase, the other through chemical desorption by the reaction of ⊥ OH with ⊥ H on dust surfaces from which 90% of the reactions on bare grain surfaces desorb into the gas phase. This drops to 18% when the surface is covered by ice. This means that if the dust formation channel is not considered, the gas-phase water formation rate will typically be underestimated by a factor of ~2 with the conditions used in our simulations. Other channels, such as OH with H, H₂ with O⁻, and the second important desorption channel ⊥ H₂O₂ with ⊥ H are less significant routes to form gas-phase water. The latter channel may become much more relevant if reaction-diffusion competition is considered for tunneling because of the high (E_a = 1000 K) barrier.

In the top right panel of Fig. 6, we address the formation rates of CO in the gas phase. CO mainly forms through the dissociative recombination reaction HCO⁺ + e⁻. This recombination rate is quite high as long as there are enough electrons around. The CO formation rate is also strongly dependent on the HCO supply in the gas phase, in which, as we show below, surface chemistry plays a major role. Below n_H = 10³ cm^-3, HCO⁺ is primarily formed by the reaction C⁺ with H₂O, but at higher densities HCO⁺ is mainly produced through the ionization of chemically desorbed HCO. The second important formation route to form CO is through C + OH, but only early on in cloud evolution, that is, when the cloud is still in a diffuse stage. The CO formation rate by C + OH sharply decreases after it peaks around t = 1.5 × 10⁷ yr. This is because atomic carbon is becoming depleted as more and more carbon is converted into CO and HCO. Once atomic carbon is depleted, the remainder of O and OH can be channeled into other reactions. We can see, for example, an upturn in the water formation rate through the reaction of OH + H at the same moment in time. H + HCO becomes the second-most important channel at later stages (t> 1.55 × 10⁷ yr) as we enter the translucent cloud stage (n_H> 10³ cm^-3). Even though these are all gas-phase reactions, they are heavily affected by the production of HCO and OH on grain surfaces that are subsequently being released into the gas phase. Without the catalyzed species formation on dust, gas-phase CO formation rates will be constrained. We note that since hydrocarbon chemistry is beyond the scope of this work, an important CO formation channel in diffuse regions and at edges of molecular clouds, i.e., CH₂ + O, is omitted (Tielens & Hollenbach 1985; Keto & Caselli 2008). Long carbon chain species, such as C₉H⁺, which do not take part here, may also influence CO formation rates (Ruffle et al. 2002).

After the formation of the first ice layer at t = 1.65 × 10⁷ yr, H and HCO are more strongly bound to the surface, see Table 2 for adsorption energies. This reduces the mobility of atoms on the dust surface as well as the desorption probabilities because they both depend on the binding energy, which explains the acute momentary decline in the second-most important rate at this time, R_{H + HCO}. We also see that photodesorption is not an important CO producer during the whole evolution of the cloud, from diffuse conditions to the first core formation. Even without our two-ice-layer penetration restriction for the UV photons, the photodesorption rates would still be an order of magnitude lower than the main CO formation rate. Since CO photodesorption occurs most efficiently by photons with energies of around 8–9 eV (Muñoz Caro et al. 2010; Fayolle et al. 2011), while CO is photodissociated in the gas phase by UV photons with energies of >11 eV, we might be underestimating the photodesorption rate when taking into account CO self-shielding for photodesorption. A test run without any self-shielding showed that the photodesorption rate is then much higher, and almost equal to the rate R_{C + OH}, but is still lower by more than one order of magnitude than the main CO formation rate R_{HCO⁺ + e⁻}. This did not affect our results. Lastly, we can see that gaseous CO is not enhanced by thermal evaporation as the dust temperature is too cold T_dust = 7−8 K for ⊥CO to be released. In fact, this will eventually lead to the depletion of CO from the gas phase as accretion starts to become more important with increasing density.

In the bottom left panel of Fig. 6, we display the formation rates of formaldehyde in the gas phase. The most dominant pathway to form formaldehyde in diffuse clouds (t< 1.5 × 10⁷ yr, n_H< 10³ cm^-3) is by direct chemical desorption of ⊥H₂CO following the reaction ⊥HCO + ⊥H. Only 1% (ice) to 5% (bare) of this reaction results in the formation of H₂CO (gas), but it still emerges as the dominant rate. Guzmán et al. (2013) also concluded that the gas-phase formaldehyde formation is due to nonthermal desorption, but attributed its production to photodesorption. Van der Wiel et al. (2009) affirmed the need of a continuous supply of formaldehyde, for instance, from grain surface chemistry, to explain the abundances in their observations of the Orion Bar. Above t = 1.55 × 10⁷ yr, chemical desorption following the reaction ⊥H₂CO + ⊥H, together with the photodissociation of methanol in the gas by CRUV photons, overtakes the former (R_{⊥ HCO + ⊥ H}) desorption rate. In 10% (ice) to 50% (bare) of the cases, the reaction ⊥H₂CO with ⊥H results in the desorption of the products. This reaction leads to the formation of a CH₃O radical that may be released into the gas only to quickly form H₂CO + H. The in-between steps are omitted here, since they are relatively fast. Thermal desorption of ⊥H₂CO is not a significant gas-phase supplier for formaldehyde either.

In the bottom right panel of Fig. 6, we show that methanol has only one essential formation mechanism. This is by chemical desorption from the grain surface reaction of ⊥CH₃O + ⊥H. The rate remains below 10^-14 cm^-3 s^-1, which means that it will not result in much methanol into the gas phase in a cloud lifetime of ~10¹⁴ s. Most of the methanol will remain frozen onto grain surfaces. The gas can only be enriched by methanol through evaporation of the frozen species at higher temperatures inside the cloud.

Ice species formation rates In the top left panel of Fig. 7, we present the formation rates of water on surfaces. Water ice primarily forms by the reaction ⊥ OH + ⊥ H. This exothermic reaction occurs without a barrier and has a high chance to chemically desorb. The percentage that remains on the dust surface is 10%, which increases to 82% for an icy substrate. This effect is directly visible in the figure from the jump in reaction rates at t = 1.65 × 10⁷ yr when the dust is covered by an ice mantle. The reaction ⊥ OH + ⊥ H₂ also increases greatly upon first ice layer formation. With 100% of the products remaining on dust, this reaction has a barrier of E_a = 2100 K that needs to be overcome. Despite this, it dominates the water-ice formation rate at later times when the cloud density is higher, that is, n_H ≥ 10⁴ cm^-3 at the dense molecular stage. This is mainly due to the increasing n_H2/n_HI ratio during cloud evolution. Accretion of gaseous water becomes important at densities of n_H ≃ 4 × 10³ cm^-3, with gas and dust temperatures of around 10 K. By forming H₂O on dust, releasing it into the gas and reaccreting it back onto the dust will make the ice more porous. After fully covering the dust surface with water ice, the accretion rate decreases like the gas-phase water formation did as a result of the decline in the desorption probability of the main formation rate.

In the top right panel of Fig. 7, we give the formation rates of ⊥CO. CO ice grows on grain surfaces predominantly by accretion. The high CO abundance in the gas phase and the low gas and dust temperatures that allow for a high sticking coefficient result in a high accretion rate. The CO ice formation mechanism becomes somewhat circular above t = 15 Myr, because CO on grain surfaces are hydrogenated to form HCO, thereupon to be chemically desorbed into the gas phase. The desorbed HCO molecule is quickly ionized through photoionization or ion exchange reactions to form HCO⁺, making it the primary route to form HCO⁺. This ion, as we know, dissociates into CO in the gas phase only to be reaccreted on dust grains where the whole cycle restarts. A demonstration of the cycle is given below. $\begin{array}{ccc}\mathrm{C}{\mathrm{O}}_{\mathrm{gas}}\stackrel{\mathrm{a}}{\mathrm{\to }}\mathrm{C}{\mathrm{O}}_{\mathrm{ice}}\stackrel{\mathrm{b}}{\mathrm{\to }}\mathrm{HC}{\mathrm{O}}_{\mathrm{ice}}\stackrel{\mathrm{c}}{\mathrm{\to }}\mathrm{HC}{\mathrm{O}}_{\mathrm{gas}}\stackrel{\mathrm{d}}{\mathrm{\to }}\mathrm{HC}{\mathrm{O}}_{\mathrm{gas}}^{\mathrm{+}}\stackrel{\mathrm{e}}{\mathrm{\to }}\mathrm{C}{\mathrm{O}}_{\mathrm{gas}}& & \end{array}$a = accretion b = hydrogenation c = chemical desorption d = ionization e = dissociative recombination.

The chicken-and-egg problem is circumvented since CO initially forms in the gas phase through a channel independent of itself, namely H₂O + C⁺ → HCO⁺ + H at n_H< 10³ cm^-3. This route is still associated with surface chemistry to a certain degree because gas-phase water formation is enhanced by the reactions on dust grains. A relatively high CO abundance is sustained in the gas phase by the continued supply through the cycle at low, T_d< 10 K, temperatures, which would otherwise result in the rapid freeze-out of CO. Photodissociation by UV photons at low extinction A_V< 3, and the dissociative reaction ⊥ HCO + ⊥ H at high extinction A_V> 5 are other, but minor producers of ⊥ CO. Since the CO freeze-out occurs throughout cloud evolution, it is expected that CO ice will be present and well mixed within every ice layer covering the grain surfaces. The CO in the upper layers will subsequently become more and more hydrogenated as the medium density rises. This will decrease the amount of CO ice, but the first ice layer(s) should still have pristine CO ice mixed with water ice.

In the bottom left panel of Fig. 7, we report the rates needed to form ⊥H₂CO. Formaldehyde ice mainly forms after two successive hydrogenations of ⊥ CO. The reaction ⊥ HCO + ⊥ H is, therefore, the main producer of formaldehyde ice throughout cloud evolution. This reaction has a 95% (bare) to 99% (ice) probability to let the product remain on the surface of the dust grains. Reaccretion of formaldehyde by the chemically desorbed ⊥ H₂CO from the this reaction, as this is one of the main gas-phase suppliers, is low because of the low desorption percentages. De-hydrogenation (removing an H atom) of ⊥ CH₃O also produces by about two orders of magnitude less formaldehyde ice than the main rate. This is because the reaction ⊥ CH₃O + ⊥ H → ⊥ H₂CO + H₂ has a reaction barrier, albeit low, of E_a = 150 K that needs to be overcome.

The bottom right panel of Fig. 7 shows the formation rates of CH₃OH on surfaces. Methanol ice essentially forms by four repeated hydrogenations of CO. In a similar fashion as formaldehyde, methanol ice forms through the exothermic reaction ⊥ CH₃O + ⊥ H with a 95% (bare) to 99% (ice) probability to make ⊥ CH₃OH. The formation rate is initially low because first ⊥ CH₃O has to be created in ample amounts. The rate increases after t = 1.33 × 10⁷ yr when ⊥ CH₃O is steadily formed, and this rate is similar to that of formaldehyde. Methanol will still be more abundant than formaldehyde since it is a more stable, strongly bound molecule. To break up methanol with a neutral hydrogen atom, a barrier of E_a = 3000 K needs to be overcome, which allows it to survive and build up over time. Methanol is also being reaccreted by the molecules that were initially created on dust and were immediately desorbed into the gas.

A caveat in our chemical network is that we did not include species larger than methanol, nor did we treat all the associated ions of methanol in the gas phase. Even though the gas-phase chemistry is not expected to be as important as grain surface chemistry, this will still create a sink out of methanol and therefore make us overestimate it, especially at later times.

Destruction rates The most important destruction rates for the species discussed in this section are displayed in Fig. 8.

For the gas phase, only the accretion rate of CO is presented in this figure to highlight the competition between the destruction rate of ⊥CO. Other accretion rates are shown in Fig. 7. For ice species, the dominant destruction rates are given with the addition of CO photodesorption and methanol dehydrogenation.

The CO that is accreted is quickly hydrogenated to form HCO. At t = 16.5 Myr the hydrogenation rate diminishes somewhat, which will result in the increase of ⊥CO. The accretion rate of CO also slows down shortly thereafter because the chemical desorption of HCO that influences the CO abundance in the gas phase is also reduced. In the end, a new balance is reached. Like CO, formaldehyde ice is mainly destroyed by hydrogenation. The rate is initially lower than that of CO due to the lower abundance of ⊥ H₂CO. That the two hydrogenation rates are equivalent at later times suggests that a balance is reached between formation and destruction. Methanol ice is mostly destroyed by radiation. Only at very late times, dehydrogenation of methanol becomes the strongest destruction mechanism. This is due to the high activation barrier (E_a = 3000 K) in the dehydrogenation of methanol.

The main point is that all rates exhibit a twist in their curves at the time when ice completely covers the dust surface (t = 16.5 Myr). Following the transition in surface binding energies of species at this time, the reaction rates on grain surfaces change, that is, the mobility is reduced when the binding energies increase, while the chemical desorption probabilities decrease when the surface is covered by water ice.

5.4. Ice distribution around the clump

We also examined the abundance profiles of the ices within a radius of 5 pc from the clump center. Inside this region, we followed the growth of ices on dust surfaces as a function of optical depth and present column density maps. We inspected the optical depth behavior for two different epochs, at a time of t = 16.6 Myr and t = 18.2 Myr after we started our simulation. In the period between these two time intervals, the transition from CO ice to methanol ice occurs inside the clump. The density of the clump center for the two snapshots is n_H = 4 × 10³ cm^-3 and 10⁴ cm^-3, respectively, while the gas temperature in both time frames lingers around 11 K (see Fig. 4). The gas temperature in the border region (A_V ≲ 1) is higher T_g ≥ 20 K, while the dust temperature is low T_d ~ 7 K. Our resolution limits us in mapping the inner pc of our clump, and we did not have many data points at the center. This makes our curves appear rather smooth. A least-squares fifth-order polynomial was fitted to our data points and with this, the ice distributions for the two epochs are displayed in Fig. 9.

A thicker layer of ice covers the surface of the dust in the inner parts of the clump where A_V and density are higher. Solid CO is more extended than solid water at a surface coverage below one mly. ⊥CO, however, also strongly contributes at the center of the clump. This establishes that water and CO ices are well mixed everywhere in the clump during the first ice layer formation. Clearly water ice becomes the main ice constituent of the icy mantle when one mly of water ice has formed. The transition at which solid water becomes more important than solid CO occurs at A_V = 4.8 during the translucent cloud stage (see left panel of Fig. 9), while 1.6 Myr later, at a dense molecular stage, this transition takes place at A_V = 2.3 (see right panel of Fig. 9). Deeper inside the core, methanol and formaldehyde ice is rapidly being formed at the loss of ⊥CO. The methanol ice surface coverage surpasses the CO ice coverage around an A_V of 7 at t = 1.82 × 10⁷ yr, but this transition also shifts to lower A_V at later times. For example, the transition occurs at A_V = 4.5 at t = 2 × 10⁷ yr.

From the difference between two time frames, we can understand that the composition of the ice mantles do not strongly depend on the optical depth, since we have different compositions at different times for a given value of A_V. We can infer from this that photo-processes with a strong dependence on optical depth, such as photodissociation and photodesorption, are not significant factors to the ice composition within the conditions used in this work.

We also present the distribution of ices in spatial dimensions. In Fig. 10 we show maps of the species as they would appear on the sky.

Here we plot the line-of-sight column densities of the species in and around the dense clump. These maps were created at the time of t = 1.82 × 10⁷ yr after simulation start, which is the same time as in the right panel of Fig. 9. The column densities range from 10¹⁰ cm^-2 to 10²⁰ cm^-2 for all species. Water has the strongest peak at the center, followed by methanol. ⊥CO has a weaker central feature, but is more extended. Most of the cells are colored green with columns of around N_CO = 10¹⁶ cm^-2. This tells us that CO ices will be present for a wide range of environmental conditions. Methanol, on the other hand, is very centralized with a high peak at the center, but with very little methanol ice extending outward. At this stage, water ice and methanol ice emission lines should be the easiest to detect; they have column densities of 10¹⁹ cm^-2 and above inside the clump.

6. Summary of conclusions and discussion

We performed for the first time hydrodynamical simulations of a collapsing gas cloud with detailed gas and grain surface chemistry in which the interplay between gas and dust is interlinked with the thermodynamics of the cloud. We presented our results on the impact of dust chemistry on the formation of ices and the enrichment of gas-phase species during the evolution of a gas cloud. We also revealed the dominant formation routes for several species at different stages of cloud evolution.

To obtain our results, we followed a gas cloud from a fully diffuse atomic stage (n_H = 10 cm^-3), which contracted and underwent phase transitions to finally form molecular clumps (n_H> 10⁴ cm^-3). During the simulation, thick, n_mly = 59, ice layers formed within the densest clump. We find that the first ice layer covering the surface of the dust has a strong presence of adsorbed CO. This mainly comes from accretion of CO from the gas phase, since gaseous CO is in ample supply. However, to form CO in the gas phase, the gas cloud needs to be supplied with H₂O, HCO, and OH in the first place (CH₂ not included in this work). During the translucent and molecular cloud stages, these species are mostly formed on dust grains to be subsequently released into the gas phase through the process of chemical desorption. Without the enhanced HCO formation in translucent clouds through the exothermic reaction ⊥H + ⊥CO at n_H ≥ 10³ cm^-3, CO in the gas phase would form at a much lower rate. We conclude from this that gas-phase CO formation requires grain surface reactions to be effective when starting from fully atomic conditions. We also conclude that in this chemical desorption is essential in supplying HCO to the gas, but note that the desorption rate is sensitive to the activation barrier of the reaction.

Our results show us that the first ice layer is formed during the translucent cloud stage at a density of n_H = 4 × 10³ cm^-3. CO ice is well mixed at this stage with water ice. After one mly of ice has formed, freeze-out occurs more rapidly as a result of the change in binding energy of the species with respect to an icy surface. This influences the mobility of species on surfaces as well as the chemical desorption probabilities, which eventually causes the formation of more water ice. We also see that CO ice gradually decreases with increasing density and opacity because it increasingly is more hydrogenated through the successive reactions with ⊥H to form formaldehyde and methanol.

From the distribution of ices in the region surrounding the clump we see that most of the ices are at the clump center. CO ice is more extended toward the outer regions than water ice, and it is consumed in the clump center to form formaldehyde and methanol. Methanol favors the high-density regions, where it becomes the second-most abundant ice after water ice.

We conclude that grain surface chemistry strongly affects the abundances of species in the gas phase by creating strong, sometimes dominant, pathways for forming key molecules, such as water, formaldehyde, and methanol. Gas-phase species are also depleted by freeze-out onto dust surfaces at different cloud evolutionary stages. Water and especially CO freeze-out is seen to occur, which is confirmed by rates. Dust chemistry also allowed us to follow the formation and build-up of ices during the evolution of a gas cloud. Water ice becomes the dominant constituent of the ice mantle after 15.5 Myr of cloud evolution when the density has risen above 2 × 10³ cm^-3, and methanol ice has a strong ~43% presence inside dense clumps. These chemical influences affect the thermodynamic properties of the progenitors of star-forming regions. Moreover, we note that these changes occur and have an impact at early, translucent cloud stages that evolve into molecular clouds. These are critical moments for cloud fragmentation.

In summary, our results are that

• in the first ice layer(s), CO ice is well mixed with H₂O ice;

• freeze-out of species greatly increases after one mly of ice covers the dust surface;

• chemical desorption from grain surfaces can be important in supplying the gas phase with CO and formaldehyde;

• methanol and formaldehyde (gas) is seen in higher density regions (n_H> 4 × 10³ cm^-3), while formaldehyde ice and CO ice are more likely to be found in the surrounding area; and finally;

• surface chemistry alters the species abundances in translucent clouds, which will affect the whole cloud evolution.

6.1. Discussion

Our chemical network consists of 42 species and 257 reactions. Although this is a respectable amount given the complexity of this work, the network is still restricted. In this, methanol is the largest molecule that we allow to form in our network. This will inevitably create a sink out of this molecule and makes us overestimate the methanol abundance since there are no routes to form higher order molecules or ions.

Chemical desorption following the hydrogenation of CO and formaldehyde is sensitive to the activation barriers we used. The activation barriers for the reactions CO + H and H₂CO + H are sometimes found to be high, 2000–2500 K (e.g., Garrod et al. 2008; Peters et al. 2013). Such high activation barriers would cause much less CO hydrogenation and chemical desorption to occur. Monte Carlo simulations often use these high activation barriers, but also naturally consider the reaction-diffusion competition in their kinematic models. The reaction-diffusion competition makes the reactions much more likely to occur. Because we did not incorporate the reaction-diffusion competition, we used an “effective” activation barrier to account for this process. We adopted “effective” barriers of 600 K for CO + H and 400 K for H₂CO + H (Awad et al. 2005; Fuchs et al. 2009) instead of the high values considered for quantum tunneling when the reaction-diffusion competition is taken into account (5400–9600 K, Awad et al. 2005). However, recent experiments by Minissale et al. (priv. comm.) tend to find barriers on the order of 1000 K. The barriers we used in that case would overestimate the hydrogenation of CO. The exact content of the ices at long time scales will be re-addressed in a future work, once the barriers for the hydrogenation of CO and H2CO are more clearly defined for the solid state.

Higher resolution simulations will allow us to probe the inner pc of a clump, which was not the prime target of this work. We resolved the Jeans length by 44 cells for the initial state, but our Jeans resolution drops to 1.2 cells at the final stages of our simulation (taking into account the increase in mean mololecular mass). This is lower by a factor three than the Truelove criterion (Truelove et al. 1997), but artificial fragmentation at this stage is not a concern, because our simulation concludes shortly after we drop below 4 cells. The Jeans resolution at the final stages of our simulation should not affect the results on which we focused, that is, the chemical reactions.

Higher resolution will also enable obtaining more clumps with a wider range of clump densities and opacities. A statistical study can be performed with such a set of clumps to determine wether our results hold true for the general case. We note that although the chemistry calculations are a local phenomenon, they reflect the cloud history and thus are affected by the evolution of the cloud, which makes the initial conditions (initial abundances, density profile, temperatures, background radiation field) and the large-scale physics (gravity, turbulence, feedback, magnetic fields) important for the results. However, because of its nature, grain surface chemistry acts as a balancing factor to the thermodynamics of the gas by the processes of freeze-out and by its catalytic nature, which causes it to form complex molecules. This, in the end, self-regulate the evolution of the cloud to some extent. To which degree this is true can be realized by a parameter study of simulations with a set of different initial and physical conditions.

A problem arising from our non-conserving numerical solver is that the ion-electron ratio is not conserved. We find that his occurs when the differential equation becomes stiff at higher, n_H ≥ 4 × 10³ cm^-3, densities. The electron abundance diverges after 17 Myrs and is too high. This problem will have negligible effects on the thermal balance, while for the chemical evolution of the species other than electrons, we believe that the remaining simulation time is too short to have a strong impact.

Acknowledgments

This work is supported by the Netherlands Organization for Scientific Research (NWO). The authors thank P. Caselli, M. Spaans, M. Minissale, and F. Dulieu for all the useful discussions on various aspects of this work. The software used in this work was developed in part by the DOE NNSA ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago. The simulations have been run on the dedicated special purpose machines at the Kapteyn Astronomical Institute and at the Donald Smits Center for Information Technology using the Millipede Cluster, University of Groningen. Some of the kinetic data implemented in this work have been downloaded from the online database KiDA (Wakelam et al. 2012, http://kida.obs.u-bordeaux1.fr).

References

Appendix A: Reaction tables

Here we list in tables the gas-phase and the surface reactions. Since adsorption reactions (Sect. 3.3.1) and thermal desorption reactions (Sect. 3.3.2) are relatively straightforward, they are omitted from this list. We also supply the adopted activation barriers, E_a, in Kelvin within the relevant tables.

All Tables

All Figures