Grain Surface Models and Data for Astrochemistry

Abstract

The cross-disciplinary field of astrochemistry exists to understand the formation, destruction, and survival of molecules in astrophysical environments. Molecules in space are synthesized via a large variety of gas-phase reactions, and reactions on dust-grain surfaces, where the surface acts as a catalyst. A broad consensus has been reached in the astrochemistry community on how to suitably treat gas-phase processes in models, and also on how to present the necessary reaction data in databases; however, no such consensus has yet been reached for grain-surface processes. A team of ${\sim}25$ experts covering observational, laboratory and theoretical (astro)chemistry met in summer of 2014 at the Lorentz Center in Leiden with the aim to provide solutions for this problem and to review the current state-of-the-art of grain surface models, both in terms of technical implementation into models as well as the most up-to-date information available from experiments and chemical computations. This review builds on the results of this workshop and gives an outlook for future directions.

1 Introduction

The very presence of anything but atoms and obscuring minuscule dust grains in the interstellar medium (ISM) was inconceivable by astronomers merely a hundred years ago. Even the brightest minds of the time, such as Sir Arthur Eddington, were doubtful about the existence of molecules in the vast interstellar void. In his Bakerian lecture he pointed out that “…it is difficult to admit the existence of molecules in interstellar space because when once a molecule becomes dissociated there seems no chance of the atoms joining up again” (Eddington 1926).

However, around one decade later, absorption electronic transitions of the first interstellar molecular species, CN, CH, and $\mbox{CH}^{+}$, were identified (Swings and Rosenfeld 1937; McKellar 1940; Douglas and Herzberg 1941). The rapid development of radio and infrared detectors following World War II has since allowed the discovery of ${\sim} 190$ molecules in the ISM, as of March 2016 (see http://www.astro.uni-koeln.de/cdms/molecules). These interstellar species have a multitude of orbital electronic configurations and include stable molecules, radicals, open-shell molecules, cations, and anions.

Many interstellar molecules are recognizable from terrestrial and atmospheric chemistry. Among those are relatively stable species, e.g., water (H₂O), molecular hydrogen and nitrogen (H₂ and N₂), and carbon dioxide (CO₂), all of which consist of just a few atoms. More complex, hydrogen-rich saturated organic molecules are also present in space, e.g., formaldehyde (H₂CO), glycolaldehyde (HCOCH₂OH), methanol (CH₃OH), formic acid (HCOOH), and dimethyl ether (CH₃OCH₃) (Ehrenfreund and Charnley 2000; Herbst and van Dishoeck 2009). These species “bridge the gap” between the simple species listed previously and those considered of prebiotic and biological importance, e.g., amino acids. Other interstellar molecules are more exotic and unique to space. These include highly-unsaturated carbon chains and cages, e.g., HC₁₁N (Bell et al. 1997), and the fullerenes, C₆₀, ${\mbox{C}_{60}}^{+}$, and C₇₀ (Cami et al. 2010; Berné et al. 2013; Campbell et al. 2015), the latter of which are also the largest molecular species discovered to date in the ISM. Even larger macromolecules, polyaromatic hydrocarbons (PAHs), consisting of between tens and hundreds of carbon atoms, are identifiable in space as a distinct class of species through their characteristic infrared bands (see the review by Tielens 2008). In summary, it is now known that the interstellar matter out of which stars and planets form has a substantial molecular component, which plays a pivotal role in the thermal balance of the ISM and its evolution (Tielens 2010).

The first theoretical models that successfully explained the presence and abundances of early observed molecular species were developed by Bates and Spitzer (1951), Watson and Salpeter (1972b), Herbst and Klemperer (1973), and Watson (1974), and thereafter significantly extended. The common perception in modern astrophysics is that many interstellar molecules, including complex unsaturated molecules, can be readily formed through purely gas-phase kinetics. Ion-molecule reactions and dissociative recombination reactions are of particular importance. Such processes typically do not have activation barriers and thus the rate coefficients are large, and possibly even enhanced, at low temperatures (∼10–20 K, Adams et al. 1985; Herbst and Leung 1986). However, gas-phase chemistry alone cannot efficiently synthesize saturated organic species. The available reaction pathways typically require high temperatures and/or three-body reactions—conditions that are not usually met in the ISM.

Another efficient route towards increasing molecular complexity in the ISM is the chemical kinetics that occurs on dust-grain surfaces. Intriguingly, the most abundant molecule in space, molecular hydrogen, is formed almost exclusively via surface chemistry in the local universe (Gould and Salpeter 1963; Hollenbach and Salpeter 1971; Watson and Salpeter 1972a). The dust-grain surface has several roles. Firstly, the surface serves as a local “meeting point” for molecules or atoms that become bound to the dust grain via electrostatic or van der Waals forces, so-called physisorption, or by forming chemical bonds with its surface, so-called chemisorption. Secondly, the dust grain lattice can accommodate a portion of the excess energy usually generated during surface-mediated association reactions, stabilizing the product, and thus allowing large polyatomic species to be efficiently synthesized. Thirdly, in dense ($\gtrsim 10^{4}~\mbox{cm}^{-3}$) and cold ($\ll 100~\mbox{K}$) interstellar environments, thick ice mantles can grow on dust grains (${\sim} 100$ monolayers). Molecular species trapped in the ice mantle under these conditions are protected from further processing by the FUV interstellar radiation field (ISRF), although they are exposed to the significantly lower strength ambient radiation field generated internally by the interaction of cosmic rays with molecular hydrogen. Eventually, over long timescales ($\gtrsim 10^{4}~\mbox{yr}$) these pristine ices can be processed by the heating and enhanced irradiation associated with the star-formation process, potentially forming even more complex volatile and refractory organic compounds, including amino acids (Kvenvolden et al. 1970; Ehrenfreund and Charnley 2000; Elsila et al. 2009). These molecules may then be delivered to young protoplanets and planets via accretion early in the evolution of the planetary system, or at a later time via bombardment by pristine icy bodies (Anders 1989; Chyba et al. 1990; Cooper et al. 2001).

Most modern astrochemical models of ISM chemistry simulate dust-grain surface chemical kinetics processes with various degrees of complexity (Tielens and Hagen 1982; Hasegawa and Herbst 1993a; Garrod 2013). Models using large reaction networks typically adopt the rate-equation approach as is done for the gas-phase chemistry where the time evolution of surface species’ abundances is described by a set of coupled ordinary differential equations, and the abundances considered “averaged” over the entire dust grain population, i.e., the mean-field approximation. One of the major challenges in these models, is the accurate treatment of the stochasticity of diffusive surface processes. This becomes critical when abundances of reactants on the dust-grain surface becomes very low, i.e., $\ll 1$ reactant per dust grain, and fluctuates with time, thus rendering the rate-equation approach unfeasible (Gillespie 1976; Green et al. 2001; Charnley 2001). This is the so-called “accretion-limited” case. A number of approximate or precise micro- and macroscopic Monte Carlo techniques have been proposed to overcome this issue (e.g., Biham et al. 2001; Charnley 2001; Lipshtat and Biham 2004; Stantcheva and Herbst 2004; Chang et al. 2005; Garrod 2008; Vasyunin and Herbst 2013). Another challenge is to account for the multilayered nature of dust-grain ice mantles, and to take all relevant processes into account in the modeling, e.g., inter-lattice diffusion, mobility/immobility of reactants, desorption, porosity trapping (see Cuppen and Herbst 2007; Chang et al. 2007; Kalvāns and Shmeld 2010; Wolff et al. 2011; Taquet et al. 2012; Vasyunin and Herbst 2013). A further obstacle in both approaches (rate equation and stochastic) is the lack of appropriate laboratory data on binding energies and desorption efficiencies of molecular ices of astrophysical interest, as well as the energy barriers and branching ratios for surface reactions. Large reaction networks can treat up to a few hundred different surface species; however, only a handful of reaction systems and molecules have been theoretically or experimentally studied. Moreover, the underlying molecular mechanism is not always fully understood, which makes it hard to scale up to astrophysically relevant timescales.

A wealth of evidence suggests that dust-grain surface processes are important over a wide range of interstellar conditions and star-formation environments, while models and observations are rapidly advancing to trace this chemical evolution through to at least the protoplanetary disk phase (Henning and Semenov 2013; Dutrey et al. 2014; Walsh et al. 2014). The new Atacama Large Millimeter Array (ALMA), with its orders-of-magnitude increase in sensitivity and resolving power, is expected to give us an unprecedented view of potentially pre-biotic and biologically-relevant molecules in various astrophysical environments over the coming years. The analysis of these new data will require much more elaborate, and more diverse, gas-grain astrochemical models than have been developed so far. Unfortunately, a major stumbling-block in our understanding of pre-biotic chemistry in the ISM is the lack of a standardized and comprehensive approach to simulate grain-surface chemistry. In the case of gas-phase chemistry, several publicly available databases with reactions and the corresponding rate data exist, of which the UMIST Database for Astrochemistry (UDfA) and the KInetic Database for Astrochemistry (KIDA) are the most widely used (McElroy et al. 2013; Wakelam et al. 2012). Consensus on how these data should be used, including how the rate coefficient is calculated and over which temperature ranges it is viable, has been reached, and the quality of the data in these databases is regularly reviewed. However, for grain-surface chemistry, this is not yet the case. Modelers often compile their own grain-surface reaction networks, and most are not publicly available, primarily due to the lack of an agreed and standardized approach.

Fortunately, many of the assumptions within the models can now be tested using surface science techniques with interstellar ice analogs. Over the past few years, substantial progress has been made on the understanding of various grain-surface reaction systems, including which processes are dominant and under what conditions, as well as the underlying mechanisms. In the summer of 2014, astronomers, experimentalists and theoretical chemists came together during a Lorentz Center Workshop (“Grain-Surface Networks and Data for Astrochemistry”) to identify the needs of modelers for their models, the appropriate formalisms to use, and to identify how recent experimental techniques and results can help to test and improve the models. In this paper we summarize the key findings of this workshop and relay our recommendations for the treatment of grain-surface chemistry to the astrochemical community. First we describe the outline of a typical gas-grain model (Sect. 2) and in Sects. 3 to 8 we discuss, in turn, the various processes that need to be considered in astrochemical models of surface chemistry: accretion, desorption, surface reactions, diffusion (thermal diffusion in the surface and bulk and quantum tunneling), and photoprocesses. We also address the more technical aspects of writing and executing code such as numerical precision in chemical models (Sect. 9) and finally end with a discussion on a test case of CO hydrogenation to form complex molecules (Sect. 10) and a future outlook (Sect. 11).

2 Outline of a Generic Gas-Grain Code

Gas-grain astrochemical models typically use the rate-equation approach to describe both the gas-phase and grain-surface chemistry using chemical kinetics. This generates a set of stiff ordinary differential equations that can be numerically solved using a multi-step integrator, e.g., via Runge-Kutta or Adams algorithms. In chemistry, rate equations are often applied to describe macroscopic experimental effects and account for many-body effects with a mean-field approach. As we mentioned above, this may not be the case on a dust-grain surface under particular conditions; hence, using the rate-equation approach to describe interstellar surface chemistry can lead to large errors when compared with results using more realistic stochastic techniques. The main reason why modelers persist with such a method is the convenience, stability, and the rather fast numerical performance of the pure chemical kinetics codes, even for reaction networks which consist of thousands of reactions involving hundreds of molecules (e.g., Dalgarno and Black 1976; Leung et al. 1984; D’Hendecourt et al. 1985; Brown and Charnley 1990; Hasegawa et al. 1992; Bergin et al. 1995; Millar et al. 1997; Aikawa et al. 1996; Willacy et al. 1998; Semenov et al. 2010; Wakelam et al. 2010; Agúndez and Wakelam 2013; Albertsson et al. 2013; McElroy et al. 2013; Grassi et al. 2014). As an indication, rate equations require CPU time of ${\sim}1\mbox{--}60$ seconds. Typically, the addition of the modified rate approach to the rate equation model, which use the same numerical scheme, slows it by a factor of several due to the performance penalty for accounting for probabilities of reactants to be on the grain surface (Garrod et al. 2009). A multiphase model (Furuya et al. 2016) without bulk ice chemistry, but with swapping only, takes typically ∼30–60 minutes per trajectory, using a full gas-ice network with deuterium chemistry. It is hence computationally feasible to use such a model to simulate a collapsing core model; tracing 35000 parcels from the prestellar core phase to the circumstellar disk phase results in ∼35000 CPU hours or ${\sim}2$ weeks on a ∼100-core machine. When bulk ice chemistry is included, the CPU time increases by a factor of 10–100. Adding in bulk chemistry increases this to months and hence a multi-phase model with bulk ice chemistry coupled with 2-D/3-D physical models remains a computationally challenging problem. On the other hand, macroscopic Monte Carlo models like presented in Vasyunin et al. (2009) require much more CPU time, from hours to days, for a simulation of a TMC-1 type cloud. What is more important, Monte Carlo models usually have a rather limited range of physical conditions that can be considered due to their slow performance. Microscopic Monte Carlo models (Lamberts et al. 2014b; Cuppen et al. 2009; Chang and Herbst 2014, 2016) are restricted to an even smaller chemical network and require days to weeks. The method of choice is hence highly dependent on the available computer power, the problem that one would like to address which dictates the level of detail in the grain description required, and the heterogeneity and complexity of the astrophysical object that one is interested in. In recent times, efforts have been made to simulate both laboratory and astrophysical conditions with the same model, thus using the laboratory simulation as a benchmark (e.g., Lamberts et al. 2013). Especially for these cases, microscopic Monte Carlo methods are worth the extra computational effort since they allow to include more surface complexity that might be crucial to gain insight in the physical and chemical processes occurring in the experiments. At the same time, laboratory environments typically deal with well-constrained physical conditions and a limited chemical network.

Here we present a recipe for the construction of a chemical kinetics model based on the rate-equation approach (see also Semenov et al. 2010). The chemical system consists of two major phases: the gas-phase and the dust-grain surface ice mantle. If all the necessary kinetic data are provided (e.g., from a database) and the initial abundances are assumed (e.g., from observations) or generated (e.g., using a similar model), a chemical kinetics code numerically solves the equations of first- and second-order kinetics and returns time-dependent molecular concentrations. Under typical ISM conditions, i.e., low densities, three-body reactive collisions are usually irrelevant and hence ignored. Here, we focus solely on the grain-surface chemistry aspect of the code. The treatment of gas-phase chemistry has been described in a number of papers, including McElroy et al. (2013) and Wakelam et al. (2010).

As schematically shown in Fig. 1, species on a grain surface generally experience four types of processes: (i) accretion (or adsorption) onto the surface, (ii) desorption from the surface, (iii) diffusion across the surface or on/within the ice mantle, and (iv) reaction. When grain-surface ice mantles are still exposed to far-UV radiation, species contained within can also be photodissociated. This leads to the following expression for the change in surface abundance:

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d} t} n_{\mathrm{s}}(\mathrm{A},t) =&\sum _{i, j} f_{\mathrm{react}}(i + j\longrightarrow \mathrm{A})- \sum_{i} f_{\mathrm{react}}(i + \mathrm{A} \longrightarrow) \\ &{} + \sum_{i} f_{\mathrm{UV}\ \mathrm{diss},i} - f_{\mathrm{UV}\ \mathrm{diss},\mathrm{A}} \\ &{}+ f_{\mathrm{acc},\mathrm{A}}- [f_{\mathrm{evap},\mathrm{A}} + f_{\mathrm{nonth},\mathrm{A}}]. \end{aligned}$$

(1)

The first four terms in this expression account for the gain and loss of species A due to grain-surface reactions or photodissociation reactions, respectively. The fifth term expresses the accretion of species A from the gas phase onto the grain, and the final term denotes the desorption of species A from the grain back into the gas phase. This latter process can occur via thermal desorption or by non-thermal processes, whereby desorption is trigged by the input of external energy in the form of far-UV or X-ray photons or high energy particles or by energy released during in-situ exothermic reactions. In the subsequent Sections we will discuss in detail the functional forms usually adopted for each of these chemical processes.

Fig. 1

Overview of the most important grain surface processes that will be covered in this review

As previously mentioned, for low surface abundances, the mean-field assumption inherent in the rate-equation approach breaks down, and several stochastic methods have been developed to overcome this issue. Although the description of the chemistry is intrinsically more accurate, stochastic models are computationally much more demanding than rate equations, and for the purpose of this review we will limit ourselves to rate-equation models. Modifications to the rate-equation approach can be made to better treat the surface chemistry in the accretion-limited case. Caselli et al. (1998) were the first to propose such an adjustment. They applied a semi-empirical approach to scale down the reaction rates for those cases where the surface migration of atomic hydrogen is significantly faster than its accretion rate onto grains (Caselli et al. 1998). This method gave good agreement with stochastic methods for a number of cases; however, it was not clear how applicable the method was outside of the tested regime. More recently, a new modified-rate approach was suggested by Garrod (2008) which improves upon the original.

3 Accretion

The accretion term $f_{\mathrm{acc},\mathrm{A}}$ in Eq. (1) accounts for the adsorption of gas-phase species onto the dust grains. It is determined by the collisional frequency of a gas-phase species with a grain, times a sticking efficiency, $S_{\mathrm{A}}$:

$$ f_{\mathrm{acc},\mathrm{A}}= S_{\mathrm{A}} v_{\mathrm{A}} n_{\mathrm{grain}} \pi r_{\mathrm{grain}}^{2} n_{\mathrm{g}}(\mathrm{A}), $$

(2)

where $n_{g}(\mathrm{A})$ is the number density of species A in the gas phase, $r_{\mathrm{grain}}$ is the average radius of a dust grain (${\sim} 0.1~\upmu \mbox{m}$ for ISM-like grains) with number density, $n_{\mathrm{grain}}$, and $v_{\mathrm{A}}$ is the average gas-phase thermal velocity,

$$ v_{\mathrm{A}} = \sqrt{\frac{8 k T_{\mathrm{gas}}}{\pi m_{\mathrm{A}}}}. $$

(3)

This in turn depends on the gas temperature, $T_{\mathrm{gas}}$, the mass of the species, $m_{\mathrm{A}}$, and Boltzmann’s constant $k$. The sticking efficiency or probability, $S_{\mathrm{A}}$, of species A to the surface is determined by how well it can dissipate its kinetic energy. This depends on the dust-grain and gas temperature, on the relative masses of the substrate molecules and the incoming species, and on the presence of a barrier for sticking, typically restricted to chemisorption. For most species at low gas and grain temperatures, this results in a sticking fraction near unity, with the exception of hydrogen. Computationally, sticking fractions have been determined by Molecular Dynamics (Buch and Czerminski 1991; Al-Halabi et al. 2002, 2003, 2004; Batista et al. 2005; Veeraghattam et al. 2014), perturbation and effective Hamiltonian theories, close coupling wavepacket, and reduced density matrix approaches (Lepetit et al. 2011), or by the much-more-approximate soft-cube method (Logan and Keck 1968; Burke and Hollenbach 1983). These studies typically focus on the accretion of a single atom or molecule on an otherwise bare surface, whereas the sticking coefficient could be coverage dependent, especially for chemisorption where for high coverage there are simply fewer sites available for sticking. Experimentally it was found that the sticking coefficient of physisorbed H₂ increases linearly with the number of deuterium molecules already adsorbed on the surface (Amiaud et al. 2007; Chaabouni et al. 2012).

4 Desorption

The desorption term in Eq. (1) represents the desorption of the species from the grain surface back into the gas phase. Various desorption processes are possible and usually a particular distinction is made between thermal desorption and non-thermal desorption. For the latter process, a multitude of different mechanisms is possible, such as photodesorption (Westley et al. 1995b; Öberg et al. 2007, 2009a,b), sputtering by cosmic rays or grain heating by cosmic rays (Hasegawa and Herbst 1993a; Herbst and Cuppen 2006), and reactive or chemical desorption, where the excess heat generated upon reaction allows desorption of the products (Garrod et al. 2006, 2007; Dulieu et al. 2013). Here, thermal and reactive desorption are briefly discussed. Photodesorption is discussed in a separate section on photoprocesses since photodesorption and photodissociation are parallel processes which require a different treatment.

4.1 Thermal Desorption

The residence time of a species on the dust-grain surface is predominantly determined by its desorption rate. The thermal desorption rate, in turn, depends on the binding energy of the species to the surface, $E_{\mathrm{bind},\mathrm{A}}$,

$$ k_{\mathrm{evap},\mathrm{A}} = \nu \exp \biggl(-\frac{E_{\mathrm{bind},\mathrm{A}}}{kT} \biggr), $$

(4)

where $\nu$ is a characteristic attempt frequency. Here $\nu$ and $E_{\mathrm{bind},\mathrm{A}}$ are important input parameters for astrochemical models. Usually the following equation for the characteristic frequency (Tielens and Allamandola 1987) is assumed,

$$ \nu = \sqrt{\frac{2 N_{s} E_{\mathrm{bind},\mathrm{A}}}{\pi^{2} m_{A}}}, $$

(5)

where $N_{s}$ is the surface density of binding sites and $m_{A}$ is the mass of species $\mbox{A}$. Tielens and Allamandola (1987) derived this expression assuming that the vibrational frequency perpendicular to the surface equals the vibrational frequency parallel to the surface and that the binding can be described by a harmonic potential, which might not be an accurate assumption for a physisorbed species. They also derived an expression including rotational degrees of freedom and one for the frequency of a free particle $\nu = \frac{kT}{h}$. Together this leads to an estimation of $\nu = 10^{12}~\mbox{s}^{-1}$, in accordance with all three approaches.

The binding energy, $E_{\mathrm{bind},\mathrm{A}}$, and the thermal desorption rate, $k_{\mathrm{evap},\mathrm{A}}$, can be experimentally obtained using Temperature Programmed Desorption (TPD). These experiments are usually performed under ultra-high-vacuum conditions (base pressure better than ${\sim}10^{-9}~\mbox{mbar}$) coupled with a quadrupole mass spectrometer. The temperature of the substrate can be carefully controlled using a cryostat. A TPD experiment consists of two phases: (i) the substrate is brought to a constant low temperature and a known quantity of one or more species is deposited, and (ii) the temperature is linearly increased and the desorption monitored using the mass spectrometer.

Different types of analysis methods can be applied to obtain kinetic parameters, such as desorption energy, desorption order, and the so-called “prefactor”, which is analogous to (although not entirely equivalent to) the characteristic frequency, $\nu$. Sometimes the latter two parameters are assumed and only the first is obtained from the analysis, other groups use for instance “leading edge fitting” to obtain all three simultaneously. Whichever method is applied, the three parameters are not completely independent and therefore a desorption energy derived from experiment should be used in combination with its corresponding prefactor. In most gas-grain codes, the computationally convenient description of the prefactor by Eq. (5) is adopted, although sometimes not physically accurate. It is recommended that experimentalists always quote the desorption energy derived in combination with the prefactor, and a fixed integer value for the desorption order, so that the binding energies are used in an appropriate manner in astrochemical models.

The desorption order is an important consideration worthy of further discussion. Zeroth-order desorption, i.e., a constant desorption rate, generally occurs when multiple layers of the same species are deposited. The number of surface species available for desorption (limited to the top few monolayers) remains the same; hence, the desorption rate is independent of the number of total species on the surface. In the sub-monolayer regime, first-order desorption is observed. Second-order desorption, i.e., a quadratic dependence of the desorption rate on the number of surface species, is also seen. This can occur in two cases: (i) when the surface exhibits a distribution of binding sites, and (ii) through chemical desorption of species that are formed via a second-order surface reaction.

In many astrochemical models, the first-order thermal desorption rate is assumed,

$$ f_{\mathrm{thermal},\mathrm{A}} = k_{\mathrm{evap},\mathrm{A}} n_{\mathrm{s}}(\mathrm{A}) $$

(6)

where $n_{\mathrm{s}}(\mathrm{A})$ is the number density of species A adsorbed on the grain surface. As mentioned above, this only strictly occurs in the sub-monolayer regime. In two-phase gas-grain astrochemical models there is no positional information on the various species, so it is not known which species occupy the top layers of the ice mantle. However, it is possible to apply a fix to the thermal desorption rate to account for the fractional composition of the ice mantle, as well as treating thermal desorption as a zeroth-order process in the multilayer regime. This involves counting the number of monolayers present within the ice mantle,

$$ N_{\mathrm{mono}} = \frac{\sum_{i} n_{\mathrm{s,i}}}{N_{\mathrm{s}}\sigma_{g} n_{\mathrm{grain}}}, $$

(7)

where the numerator is the total number density of surface species per unit volume, and the denominator is the total number of surface sites available per unit volume. For $N_{\mathrm{mono}} > N_{\mathrm{act}}$, where, $N_{\mathrm{act}}$ is the assumed number of “active” monolayers (typically ${\sim} 2\mbox{--}4$), the thermal desorption rate is given by,

$$ f_{\mathrm{thermal},\mathrm{A}} = k_{\mathrm{evap},\mathrm{A}} N_{\mathrm{act}} \chi_{\mathrm{A}} N_{\mathrm{s}} \sigma_{g} n_{\mathrm{grain}}, $$

(8)

where $\chi_{\mathrm{A}} = n_{\mathrm{s}}(\mathrm{A})/\sum_{i}n_{\mathrm{s,i}}$ is the fractional abundance of species A in the ice mantle. For $N_{\mathrm{mono}} \le N_{\mathrm{act}}$, the rate switches to the first-order desorption rate.

Table 1 lists the binding energies of a wide collection of stable species that have been determined using the TPD technique and are therefore relatively well constrained. The binding energies have been mostly determined for the desorption of pure ices from different substrates. The differences between the different substrates are rather small and become negligible in the multilayer regime (Green et al. 2009). The uncertainties on the binding energies quoted in Table 1 can have different origins: experimental errors, errors in the fit, or—especially for the amorphous silicate surfaces—they can represent a range of binding energies which is an intrinsic property of the substrate.

Table 1 List of experimentally determined binding energies

Desorption rates depend exponentially on binding energies and uncertainties in these binding energies can have a large effect, even at dark cloud conditions where the temperature is well below the desorption temperature of the vast majority of the surface species (Penteado et al. 2016). Since in most systems the diffusion barrier is calculated by taking a fixed ratio with the binding energy, changing binding energies not only affects the temperature at which species desorb, i.e., the temperature at which species cannot participate in the grain surface chemistry, but also the onset temperature at which species start to diffuse. A sensitivity analysis of grain surface chemistry under dark cloud conditions to binding energies of ice species showed that the model results appear particularly sensitive to the binding energy of H₂ (Penteado et al. 2016). The dust temperatures in molecular cloud cores are relatively well constrained to precision of ∼1 K by a number of Herschel studies (e.g., Stutz et al. 2010; Launhardt et al. 2013; Lippok et al. 2016).

The experiments show that the molecules indeed desorb with a (close to) zeroth-order rate in the multilayer regime whereas they desorb with a (close to) first-order rate in the monolayer regime, as explained above (Fraser et al. 2001). It is difficult to experimentally obtain similar results for radical species due to their high reactivity (and correspondingly short lifetime). Binding energies for radicals can only be obtained in an indirect manner, usually involving the simulation of experimental data, and an exploration of the possible parameter space, using stochastic chemical models. However, there are recent experimental results reporting the experimental determination of the binding energy of atomic oxygen on a range of surfaces (Dulieu et al. 2013; He et al. 2014, 2015), showing that for some species, direct measurements are possible.

Most TPD experiments are performed using pure ices to allow an unambiguous interpretation of the results and to minimize the chance of contamination. Some studies on mixed and layered ice have been done to better mimic the composition of interstellar ice mantles. The introduction of more species in the ice immediately increases the complexity of the desorption process. The binding energy of individual species will vary depending on its surrounding material, and the dominant ice-mantle species can prevent other species from desorbing. Collings et al. (2004) showed, for instance, that a fraction of molecules like CO and CO₂ can become trapped in an ice mantle which consists predominantly of water ice. The trapped CO and CO₂ are then released at the higher temperatures expected for water desorption. However, laboratory timescales are significantly shorter than those in the ISM; hence, trapped species may have sufficient time to escape the ice mantle since they will also have had sufficient time to segregate (Öberg et al. 2009d). This process depends on a large number of parameters including surface temperature, ice composition and mixing ratio. Two-phase astrochemical models can include the effects of trapping in a somewhat empirical manner by allowing a fraction of volatile species such as CO to have a binding energy similar to the species within which they are trapped (e.g., CO₂ or H₂O, see, e.g., Viti et al. 2004). Three-phase and multilayer models can simulate the effects of trapping by allowing diffusion of surface species into the bulk ice mantle (and vice versa). We discuss the treatment of bulk diffusion in Sect. 7.

4.2 Reactive/Chemical Desorption

Chemical desorption is desorption of reaction products from the grain surface by excess reaction energy. This type of desorption is also referred to as reactive desorption. Garrod et al. (2006) were the first to suggest this mechanism to explain, e.g., the gas-phase detection of methanol in cold dark clouds. They based their initial model on the Rice-Ramsperger-Kessel (RRK) theory, which relates the excess energy and the binding energy of species to a desorption probability. They modified this theory by adding an unconstrained $a$ parameter which they chose to be 0.1. In a follow-up study, Garrod et al. (2007) showed that chemical desorption may play an important role in explaining the observed abundances of different gas-phase chemical species, particularly in dark molecular clouds. Later Cazaux et al. (2010) came to similar conclusions when they included this mechanism in their model for water formation on grains.

The first constraints on the probability of this mechanism were obtained using Molecular Dynamics simulations (Andersson et al. 2006; Arasa et al. 2010, 2011). In these studies, the fate of photoproducts of water ice photodissociation—OH and H—were monitored in time. In some cases, the photoproducts were found to recombine to form water which subsequently escaped from the ice mantle: this can loosely be described as reactive desorption driven by photodissociation. However, as will be discussed later, this can also be thought of as “photodesorption” (see Sect. 8). In these simulations, the desorption probability was highly dependent on the location of the dissociated molecule in the ice mantle. Recombination events in the fourth layer of the ice or further below almost exclusively resulted in trapping of the reformed water molecule. These results are limited to a water-rich environment and they may not be applicable to the formation of the first monolayers of the water ice mantle.

What remains to be quantified is the efficiency of reactive desorption which is not driven by photoprocessing. The first experimental study by Dulieu et al. (2013) measured the chemical desorption of reaction products through sequential O₂ hydrogenation experiments on an amorphous silicate or a graphite surface, where the amount of deposited O₂ remained in the (sub)monolayer regime. They find substantial desorption of the formed H₂O molecules, which is caused, at least in part, by the lack of binding opportunities with surrounding molecules. Moving to the multilayer regime, they find that the desorption probability for the $\mbox{O}+\mbox{O}$ recombination reaction drops to negligible values (Minissale and Dulieu 2014). Expanding their studies to other reaction systems (e.g., $\mbox{CO}+\mbox{H}$, $\mbox{H}_{2}\mbox{CO}+\mbox{H}$), they determined relatively low reactive desorption rates in the (close-to) sub-monolayer regime (≲10 %, Minissale et al. 2016).

Despite the lack of conclusive experimental evidence for chemical desorption driven purely by exothermicity of reactions (and not by photoprocessing), especially in the multi-layer regime, astrochemical models typically still account for such a process by implementing the Garrod et al. (2007) prescription

$$ P = \frac{a\varXi}{1+a\varXi} $$

(9)

with

$$ \varXi = \biggl[1 - \frac{E_{\mathrm{bind}}}{E_{\mathrm{exo}}} \biggr]^{s-1} $$

(10)

where $E_{\mathrm{exo}}$ is the exothermicity of the reaction, and $s$ is the number of vibrational modes in the molecule/surface-bond system. This number is $s = 2$ for diatomic species; for all others, $s = 3 N-5$, where $N$ is the number of atoms in the molecule, which is assumed to be non-linear and forms an extra “bond” to the surface. The efficiency parameter $a$ is not well constrained and is generally used as universal input parameter with a value between 0.01 and 0.1. Figure 2 shows the sensitivity of gas phase and ice abundances to this parameter in a laminar solar nebula model by changing $a$ from 0.05 to 0.01. The figure shows that changes can locally be several orders of magnitude, but integrated over the height of the disk the changes are relatively small for several species. Ices in disks concentrate around the midplane and are nearly absent in upper layers due to thermal evaporation and photodesorption. The changes in column densities are hence mainly determined by the changes in abundance around the midplane. Interestingly, the abundances of CH₃OH and CH₄ increase both in the gas phase and in the ice by lowering the reactive desorption. This is presumably since these species are formed in several steps and a lower reactive desorption efficiency keeps the intermediate species on the grain, enabling the full reaction route to proceed.

Fig. 2

The change in molecular abundance of a selection of species in a laminar solar nebula at 1 Myr, when using a chemical desorption efficiency of 0.01 ($N_{1}$) instead of 0.05 ($N_{2}$). The log of relative ratios are given both as function of location (height $z$ and radius $r$) and integrated over $z$ as a function of $r$. Relative gas phase abundances are in the left panels, the corresponding ice abundances in the panels on the right-hand side

5 Reactions

Generally, surface reactions are thought to occur via one of three mechanisms: the diffusive Langmuir-Hinshelwood mechanism, where both species move over the surface and react upon meeting, the stick-and-hit Eley-Rideal mechanism where one (stationary) reactant is hit by another species from the gas phase, and the hot-atom mechanism (which is a combination of both) where non-thermalized species travel some distance over the surface before finding a fellow reactant. Photodissociation is usually treated separately and will be discussed in Sect. 8. Under astrophysical conditions where ices are abundant, the gas and dust grains typically have similar temperatures, and the chemical timescales tend to be significantly longer than the thermalization timescale; hence, the hot-atom mechanism is often considered not important.

The analytical expression to describe the Langmuir-Hinshelwood mechanism on a surface is

$$ f_{\mathrm{react}, \mathrm{LH}}(i+j\longrightarrow \mathrm{A}) = P_{\mathrm{react}, \mathrm{LH},i,j} (k_{\mathrm{scan}, i}+ k_{\mathrm{scan}, j} ) n_{\mathrm{s}}(i,t) n_{\mathrm{s}}(j,t) $$

(11)

where $n_{\mathrm{s}}(i,t)$ is the number of species $i$ present on the surface at time $t$ and $k_{\mathrm{scan}, i}$ the rate by which species $i$ scans the grain surface. The scanning rate is given by

$$ k_{\mathrm{scan}} = \frac{k_{\mathrm{hop}}}{N_{\mathrm{sites}}} $$

(12)

where $N_{\mathrm{sites}}$ is the number of binding sites per grain. For amorphous solid water with a density of $0.94~\mbox{g}\,\mbox{cm}^{-3}$, the site density is $1\times 10^{15}$ molecules $\mbox{cm}^{-2}$. Simulations of CO₂, which is a bulkier molecule, on top of water ice showed a site density of $6\times 10^{14}$ molecules $\mbox{cm}^{-2}$ (Karssemeijer et al. 2014a). Both lead to approximately $\approx 10^{6}$ binding sites per grain for a standard grain of 0.1 μm.

The scanning rate determines the meeting frequency of the two particles $i$ and $j$ due to the mobility of one, or both, reactant(s). The Langmuir-Hinshelwood mechanism is dependent upon the abundances of both reactant and hence is a second-order process. The hopping rate, $k_{\mathrm{hop}}$, will be discuss in Sect. 6.

The rate coefficient, $P_{\mathrm{react}, \mathrm{LH},i,j}$, accounts for the probability that a possible reaction barrier will be crossed during the encounter. This probability is assumed to be 1 for a reaction with zero activation energy, and 0.5 if the two reactants are the same species. For reactions with an activation barrier, $E_{a}$, that occur on a dust-grain surface with a temperature, $T$, this probability is

$$ P_{\mathrm{react}, \mathrm{LH},i,j} = \exp \biggl(-\frac{E_{a}}{kT} \biggr) $$

(13)

when the barrier is crossed through thermal activation as schematically depicted in Fig. 3 for the reaction H + H₂CO. In this case there is a clear transition state that determines the rate limiting energy barrier. Tunneling through the reaction barrier is also possible, greatly increasing the probability of reaction (see, e.g., Hasegawa et al. 1992). This occurs through delocalization of the transition state. As can be seen for the $\mbox{H--H}_{2}\mbox{CO}$ complex, light species are much more delocalized and quantum-mechanical tunneling is hence of main importance for reactions where light species, e.g., H, D, are involved in the bond breaking or forming. Tunneling is discussed in more detail in Sect. 5.2.

Fig. 3

Schematic representation of crossing a reactive barrier either through tunneling or through thermal activation. The ${\mbox{H} + \mbox{H}_{2}\mbox{CO} \longrightarrow \mbox{H}_{3}\mbox{CO}}$ reaction is used as an example

Although conceptually simple, in reality the situation is more complex. First, a surface reaction may have several exit channels leading to a number of various products, similar to reactions in the gas phase. For most examples each of these channels will have its own transition state and corresponding (temperature-dependent) rate. The reaction constant does not need to include a special scaling to account for this effect, but the branching ratios $\alpha$ are a natural outcome of the model

$$ \alpha_{l} = \frac{k_{l}}{\sum_{j=1,m} k_{j}} $$

(14)

for the all possible $m$ reaction channels. In constructing a reaction network, one should be very careful when adding new product channels especially when the reaction rates come from very different sources (surface vs. gas phase experiments, computations) since some product channels might be heavily suppressed. For some reactions, only a destruction rate is known and branching ratios are determined separately. Individual product rate should in this case be adjusted accordingly.

Second, for a diffusive surface reaction to happen, the two molecules must remain adsorbed in close vicinity until they react, otherwise they can migrate away from each other or even desorb as schematically depicted in Fig. 4. Therefore, a surface reaction process is in competition with diffusion and desorption (Herbst and Millar 2008). Consequently, the reaction constant for product channel $k$ takes the following expression (see Equation 6 in Garrod and Pauly 2011),

$$ P^{\mathrm{full},k}_{\mathrm{react}, \mathrm{LH},i,j} = P_{\mathrm{react}, \mathrm{LH},i,j}^{k} \frac{\nu \sum_{l=1,m}P_{\mathrm{react}, \mathrm{LH},i,j}^{l}}{\nu \sum_{l=1,m} P_{\mathrm{react}, \mathrm{LH},i,j}^{l} + k_{\mathrm{hop},i} + k_{\mathrm{hop},j} + k_{\mathrm{evap},i} + k_{\mathrm{evap},j}}, $$

(15)

where the $k_{\mathrm{hop}}$ and $k_{\mathrm{evap}}$ are the thermal hopping (scanning) and evaporation rates for the reactants $i$ and $j$, consequently. In the majority of astrophysically relevant situations the evaporation terms are small in magnitude compared with the hopping terms and can be safely neglected.

Fig. 4

Surface reactions can be attempted as long as the reactions are each others vicinity. Hence, reaction competes with diffusion and desorption. Faster diffusion will lead to a higher meeting rate but also in a shorter reaction time

The Eley-Rideal mechanism is considered to be important only for high surface coverages or low surface mobility (Ruffle and Herbst 2001). An example where this can become important is during catastrophic freeze-out of CO in prestellar cores. During this phase, a layer of reactive CO ice forms on the grains (Pontoppidan 2006). Under these circumstances, Eley-Rideal could be an important mechanism in the formation of methanol. It can be included in models by using the following expressions

$$\begin{aligned} f_{\mathrm{react}, \mathrm{ER}}(i+j\longrightarrow \mathrm{A}) = P_{\mathrm{react}, \mathrm{ER},i,j} f_{\mathrm{acc},i} n_{\mathrm{s}}(j,t) + P_{\mathrm{react}, \mathrm{ER},i,j} f_{\mathrm{acc},j} n_{\mathrm{s}}(i,t). \end{aligned}$$

(16)

The reaction constant is different for the Eley-Rideal mechanism. Here the two reactants have only one attempt to cross the reaction barrier $E_{a}$, so diffusion and desorption competition is of no importance. The corresponding rate coefficient is much simpler than in the LH-case (Eq. (15)):

$$ P_{\mathrm{react}, \mathrm{ER},i,j} = \alpha_{l}\exp \biggl(- \frac{E_{a}}{kT} \biggr). $$

(17)

5.1 Surface Experiments

Surface reactions can be monitored in the laboratory using an ultra-high vacuum setup (better than ${\sim} 10^{-9}~\mbox{mbar}$ and H₂ as main gas residue in the chamber) and experiments are generally performed in two ways. Either the reactants are deposited in sequence, referred to as pre-deposition experiments, or in tandem in a so-called co-deposition experiment. The first gives a better control over the total dose and the predeposited amount can either be in the monolayer regime on top of an astrophysically relevant surface (e.g., a silicate or carbonaceous substrate), or a thicker ice if the reactant is a stable species, e.g., CO. For pre-deposition ice experiments, the final yield of the newly formed species for a selected radical fluence and ice temperature is largely limited by the penetration depth of the reactants in the ice. For the case of hydrogenation of CO ice, the maximal penetration depth of H atoms is four monolayers; therefore, only the upper layers of the ice are hydrogenated. Co-deposition experiments generally give a higher signal as they do not suffer from such penetration effects. Moreover, they are particularly useful in experiments involving radical species other than H atoms. Radicals are generally formed in a microwave discharge or RF (radio frequency) source in which the stable precursor is injected and subsequently dissociated. The dissociation products are then piped to the substrate, usually through a cooling pipe that thermalizes or cools the species to room temperature or even below. Generally, the dissociation is not 100 % efficient and recombinations of the radical species can occur in the cooling pipe. Depending on the initial stable precursor gas that is used, the desired radical might not be the only radical formed during the discharge which can lead to a beam of mixed composition. For hydrogen gas, this is H and H₂ where the latter sticks to the surface only in the monolayer regime. For oxygen gas, this is O and O₂ and in a predeposition experiment where first an ice is deposited which is then exposed to O atoms, O₂ has been seen to form a thick layer on top of the predeposited species depending on the surface temperature, which can inhibit any reaction.

Surface reactions are generally monitored by means of Fourier transform infrared spectroscopy and gas-phase molecules, generated by performing a post-experiment TPD, by mass spectrometry using a quadrupole mass spectrometer. With an in-situ technique like infrared spectroscopy, the amount of species formed can be obtained as a function of time if the band strength of the species is known. If the dose is also known, a formation rate in the experiment can then be determined. However, this depends on a combination of processes and their competition, i.e., diffusion, desorption, and reaction. Moreover, often several reactive species are present and a multitude of different reactions are possible at the same time. Thus, the reaction pathways and associated branching ratios can become nontrivial to disentangle as the systems studied increase in complexity. In addition, the characteristics of the substrate—its composition and surface structure—can affect reaction pathways and rates in the sub-monolayer regime. However, these data (rates, pathways, and branching ratios) are critical for advancing surface-chemical networks for use in astrochemistry. Simulating laboratory conditions can aid in extracting or constraining reaction data.

Although quantitative information on the reaction barrier or rate is not easily obtainable from surface experiments, experiments are extremely useful in detecting whether specific reactions can proceed under circumstances close to those present in the ISM and this is the aim of most experimental studies. It is also not trivial to calculate rates of surface reactions by quantum chemical methods. Surface reactions require many atoms in the calculation, which unequivocally increases the computational time. One of the approximations that can be made is the use of a model surface, e.g., a coronene molecule to model a carbonaceous grain surface (Adriaens et al. 2010), or a water cluster to mimic a thick ice layer (Rimola et al. 2014). Another promising method that can be used to include the effect of the surface on reaction rates is quantum mechanics/molecular mechanics (QM/MM, Goumans et al. 2009). Although such methods are robust for studying the effect of the surface on the reaction barrier height and shape, they cannot account for the adsorption of energy through phonon excitation. Therefore, they can result in different dynamical behavior of the products and thus, different rates.

Much more data is available for gas-phase reactions and a reasonable assumption for the construction of a surface reaction network is to lend from gas-phase data. As far as we are aware, there have been no reports of reactions that are efficient on the surface and not in the gas phase (or vice versa), except for association reactions of the form, ${\mathrm{A} + \mbox{B} \longrightarrow \mbox{AB}}$. This type of reaction is very inefficient in the gas phase without the presence of a third body, since it has to proceed through radiative association. For surface reactions, the grain acts as a third body. For reactions on clusters, the efficiency of association reactions likely lies somewhere in between. Examples of reactions that have a high barrier in the gas phase but nevertheless proceed in the solid phase at low temperatures are ${\mbox{CO} + \mbox{H}}$, ${\mbox{H}_{2}\mbox{CO} + \mbox{H}}$, ${\mbox{H}_{2} + \mbox{OH}}$ and ${\mbox{H} + \mbox{H}_{2}\mbox{O}_{2}}$ (see Tables 2 and 3). This is thanks to the possibility of many crossing attempts, energy dissipation for reactions of the type ${\mathrm{A} + \mbox{B} \longrightarrow \mbox{AB}}$, and tunneling; although the latter is relevant in the gas phase as well.

Table 2 Astrochemically relevant surface reactions for which tunneling has been confirmed experimentally

Table 3 Astrochemically relevant surface reactions for which tunneling has been studied theoretically

5.2 Tunneling

Experiments have shown that a number of surface reactions proceed through quantum-mechanical tunneling, a small summary of which is shown in Table 2. An extensive explanation can be found in Hama and Watanabe (2013). Briefly, both hydrogen-addition and hydrogen-abstraction reactions are important processes that can occur through tunneling.

When tunneling is involved in “crossing” the reaction barrier, the description of the rate coefficient changes. To accurately describe tunneling, a full quantum-mechanical calculation is preferred and there are numerous methods available to calculate tunneling rates (summarized recently by Kästner 2014). However, such methods typically lie beyond the scope of most astrochemical models. Instead, the most common way to account for tunneling in models is by using the Wentzel-Kramers-Brillouin approximation and the (crude) assumption of a rectangular barrier, which leads to a rate constant,

$$ k_{\mathrm{react},i,j} = \nu \exp \biggl( -\frac{2a}{h} \sqrt{2 \mu E_{\mathrm{a}}} \biggr), $$

(18)

where $\nu$ is an attempt frequency (as described previously), $a$ is the width of the barrier, $h$ is Planck’s constant, $\mu$ is the effective mass of the system and $E_{\mathrm{a}}$ is the reaction barrier. Temperature no longer plays a role, in contrast with the rate coefficient for thermally-activated reactions, and the reaction probability increases with decreasing reduced mass and barrier width.

For a particular reaction system, below the so-called cross-over temperature, the contribution of quantum-mechanical tunneling to the reaction rate dominates over the thermal contribution. Tunneling can only occur through the exothermic part of the barrier. For an endothermic reaction, the rate coefficient described by Eq. (18) should be scaled by a Boltzmann factor accounting for the difference in energy between the reactant and product states, following arguments of detailed balance

$$ k_{\mathrm{react},j,i} = \nu \exp \biggl( -\frac{2a}{h} \sqrt{2 \mu E_{\mathrm{a}}} \biggr) \exp \biggl( -\frac{\Delta E}{kT} \biggr) $$

(19)

where $E_{\mathrm{a}}$ and $\Delta E$ are defined in Fig. 5. This decreases the reaction rate considerably, and can potentially significantly alter the outcome of large astrochemical models (e.g., Lamberts et al. 2014a).

Fig. 5

Energy profile of a reaction the forward reaction is exothermic; the backward reaction is endothermic and its rate should be described by Eq. (19)

First of all, the width of the barrier, $a$, in the simple approximation mentioned above, is often not known and therefore commonly approximated as 1–2 Å (Garrod and Herbst 2006; Furuya et al. 2013; Walsh et al. 2014), initiated by older literature such as Tielens and Hagen (1982) and Hasegawa et al. (1992). However, the shape of the barrier can also have a significant effect on the tunneling transmission coefficient. Eq. (18) assumes a rectangular barrier, i.e., with equal initial and final energies, which is not the case for most chemical reaction systems. A computationally cheap method to improve on this rectangular shape is using the Eckart model for the shape of the potential (Eckart 1930). However, this requires knowledge of the forward and backward reaction barriers, as well as the imaginary frequency of the transition state. Therefore, (quantum) calculations need to be available, which is typically the case only for the gas phase. Recently, Taquet et al. (2013) demonstrated the significance of using the Eckart model versus a rectangular barrier on the calculation of the rates for a set of surface reactions, using gas-phase theoretical data as input.

Secondly, let us consider the mass dependence, which can be intuitively understood by considering the uncertainty principle $\Delta x \Delta p > \hbar/2$. A large mass will result in a large $\Delta p$, causing the accuracy of the position, $\Delta x$, to be better defined. Hence, a particle is more localized and is less likely to “leak” some of its probability density into or through a barrier. Eq. (18) requires the effective mass, $\mu$, along the reaction coordinate, which for a simple bond forming reaction can be approximated by,

$$ {\mu} = \biggl( \frac{1}{m_{\mathrm{1}}} + \frac{1}{m_{\mathrm{2}}} \biggr)^{-1} = \frac{m_{\mathrm{1}} m_{\mathrm{2}}}{m_{\mathrm{1}} + m_{\mathrm{2}}}, $$

(20)

where $m_{1}$ and $m_{2}$ are the masses of the individual atoms forming the bond. However, in astrochemical models, generally the reduced mass of the system is used,

$$ {\mu} = \biggl( \frac{1}{m_{\mathrm{A}}} + \frac{1}{m_{\mathrm{B}}} \biggr)^{-1} = \frac{m_{\mathrm{A}} m_{\mathrm{B}}}{m_{\mathrm{B}} + m_{\mathrm{A}}}, $$

(21)

where are A and B are the reacting species, regardless of the reaction mechanism. For abstraction reactions ${\mathrm{A} + \mbox{XB} \longrightarrow \mbox{AX} + \mbox{B}}$, the situation becomes more complex and depends on the incoming angle. The expression that can be applied to calculate the reduced mass for linear systems A–X–B is

$$ \mu = \frac{m_{\mathrm{A}}m_{\mathrm{X}}(1+c)^{2} + m_{\mathrm{X}}(m_{\mathrm{A}} + c^{2} m_{\mathrm{B}})}{ (m_{\mathrm{A}} + m_{\mathrm{X}} + m_{\mathrm{B}})(1+c^{2})} . $$

(22)

Here, $c$ is the ratio between the infinitesimal change in bond length of $R_{\mathrm{AX}}$ and $R_{\mathrm{XB}}$. For reactions where species A and B have equal mass, or both have masses much larger than that of X, $\mu$ is between $\frac{1}{2}m_{\mathrm{X}}$ and $m_{\mathrm{X}}$. In general, however, this depends on the forces acting on the system. More extensive theoretical literature can be found in Bell (1980).

As an example, the reduced mass for abstraction reactions by OH, e.g.,

$$ {\mbox{A--H} + \mbox{OH} \longrightarrow \mathrm{A}\cdot + \mbox{H}_{2} \mbox{O}}, $$

(23)

should reflect that the tunneling species is the hydrogen atom and the effective mass should be $\mu \approx m_{\mathrm{H}}$. However, using Eq. (21) will yield a much higher mass and this has a large influence on the reaction rate. Consider, for example, a reaction with a barrier $E_{\mathrm{a}}$ of 2000 K and a width of 1 Å. Using Eq. (18) with two different values for $\mu$ ($\mu_{1}=10$ and $\mu_{2}=1~\mbox{a.u.}$) yields $k_{1}/k_{2} = 2 \times 10^{-3}$.

Typically, gas grain codes use information on barrier height and the masses of the products to calculate the rates. We recommend to include the tunneling reaction rate as obtained from high-level gas-phase calculations as an additional input parameter in an astrochemical model. When such rates are unavailable, the original approximate Eq. (18) can be reinstalled.

The importance of tunneling can be experimentally confirmed in two ways: either by studying the temperature dependence—to determine whether this confirms the predicted behavior of either Eq. (15), (17) or (18)—or by making use of the mass dependence of tunneling. In this latter case, an experiment is performed (at least) twice, using different isotopes of the same species between the two (or more) runs. Since isotopes are chemically equivalent, but have a different mass, this is a convenient way to detect the difference in tunneling efficiency. If the reaction products scale with the effective mass of the reaction, this is usually interpreted as a proof of tunneling and is called the Kinetic Isotope Effect (KIE). For instance, a reaction with a hydrogen atom, ${\mbox{H} + \mbox{X} \longrightarrow \mbox{HX}}$, can be compared to the same reaction with a deuterium atom, ${\mbox{D} + \mbox{X} \longrightarrow \mbox{DX}}$. In both cases, the experiments are performed at temperatures as low as possible, which is typically setup-dependent.

In order to study the influence of tunneling on a reaction with a relatively high precision from a computational perspective, it is necessary to use methods that encompass more than only an Eckart calculation of the barrier, or a Wigner correction (Wigner 1938) to the rate. For instance, a variety of Transition State Theories (TST) have been applied, such as Variational TST (or Canonical Variational Theory, CVT) combined with a Multidimensional Tunneling correction (VTST/MT), Semi-Classical TST (SCTST), and Harmonic Quantum TST (HQTST). Other methods, such as Quantum-Reaction Path Hamiltonian method (Q-RPH), Ring Polymer Molecular Dynamics (RPMD), and Free Energy Instanton Theory (FEIT) have also been employed. Several recent papers comment on the differences and accuracy of these various methods with respect to each other (Nyman 2014; Zhang et al. 2014; Hele et al. 2015). In Table 3, a number of gas-phase studies are listed that focus on reactions that may be important in surface astrochemistry. However, this is by no means an exhaustive list and focuses only on recent papers. Note that we have excluded the reaction ${\mbox{H}_{2} + \mbox{O} \longrightarrow \mbox{OH} + \mbox{H}}$ from this table, based on its endothermicity (Lamberts et al. 2014a) and on the high barrier and, consequently, low reaction rate (Nguyen and Stanton 2014).

5.3 Temperature Dependence of Surface Reactions

Although many reaction rates are temperature independent, because they are either barrierless or occur through quantum tunneling, reactions still have a temperature window within which they are most effective. This is a consequence of the temperature dependence of the diffusion and desorption of reactants. Diffusion is required for reactants to meet in a Langmuir-Hinshelwood reaction and determines the lower bound of the temperature window; desorption, on the other hand, sets the upper bound. As a result, at low temperatures, surface chemistry will mostly be dominated by hydrogenation reactions,

$$ {\mbox{X}\cdot + \mbox{H} \longrightarrow \mbox{X--H}}, $$

(24)

effectively converting all radical species to more stable species such as H₂O, CH₄, NH₃, and CH₃OH. In parallel, hydrogen abstraction can recreate surface radicals,

$$ {\mbox{X--H} + \mbox{Y}\cdot \longrightarrow \mbox{X}\cdot + \mbox{Y--H}}, $$

(25)

but most likely with a much lower efficiency. Whereas hydrogenation reactions have been extensively studied, for hydrogen abstraction reactions, the reaction rates and the importance of tunneling remain unknown.

At higher temperatures, the residence time of hydrogen atoms on the surface becomes too short for hydrogenation reactions to dominate, and other surface reaction types increase in importance. For example, high-mass star-forming cores, that are in the so-called “hot core” phase, are found to be rich in gas-phase complex organic molecules, including alcohols, aldehydes, carboxylic acids, esters, and ethers (Bisschop et al. 2007b). While several of these classes of molecules may have viable gas-phase formation mechanisms, the larger species are thought to be formed primarily (or exclusively) on the surfaces of dust grains, or within/upon dust-grain ice mantles (Garrod and Herbst 2006; Garrod et al. 2008). The main reaction mechanism for the formation of complex organics on the grains is the creation of functional-group radicals on/within the ice mantle, such as CH₃ and CH₃O, both of which may be formed via the photodissociation of methanol (CH₃OH). As the temperature in the core increases to above $20~\mbox{K}$ or so, these radicals become mobile, diffuse across or through the ice mantle, thereby allowing radical-radical association reactions to become competitive with hydrogenation of radicals by abundant H atoms. As temperatures rise yet further ($\gtrsim 100~\mbox{K}$), the grain-surface-formed molecules desorb into the gas phase, where they are detected typically via (sub)mm rotational spectroscopy.

Ubiquitous molecules, including methyl formate (HCOOCH₃) and dimethyl ether (CH₃OCH₃), appear to have efficient dust-grain surface formation routes. Astronomical observations and astrochemical modeling (e.g., Belloche et al. 2014) indicate that surface/ice-mantle formation mechanisms may also be efficient for very large molecules like propyl cyanide (C₃H₇CN); however, in such cases, the earlier formation of smaller homologues (e.g. ethyl cyanide, C₂H₅CN) is usually required. The removal of a hydrogen atom from these molecules produces the necessary large radicals, to which other functional groups may be added, further increasing chemical complexity.

The most abundant molecule in the ices is water itself, which may be photodissociated to form the highly reactive OH radical. At relatively low temperatures, most OH produced recombines with H, or with mobile, abundant radicals like CH₃ and HCO. However, models suggest that above around 60 K (Garrod 2013), OH becomes sufficiently mobile to find and react with large stable molecules, to abstract an H atom and produce a large radical. At this stage, H-abstraction by OH (along with NH₂, which is also formed by H-abstraction from ammonia) becomes the dominant formation mechanism for the molecular radicals that are the precursors to even larger species on the grains.

While relatively few of the H-abstraction reactions of OH invoked in this dust-grain chemistry have been directly measured, even in the gas phase, those for which rates have been determined display very small activation energy barriers to H abstraction by OH, typically less than $1000~\mbox{K}$. Crucially, barriers of this size are comparable to the expected diffusion barrier for surface OH, meaning that for these reactions, the surface diffusion of OH is the rate-limiting step. Thus, as soon as temperatures are sufficiently high for OH diffusion, OH becomes the key instigator of radical production.

Unfortunately, few of these barriers to H-abstraction by OH from large, saturated molecules have been determined, and the calculation of rates in models is further complicated by the fact that the abstracted hydrogen atom may be able to tunnel through the barrier, meaning that information about the barrier shape is required. Alternatively, if tunneling is efficient, then in many cases the key quantity needed in astrochemical gas-grain calculations will be the diffusion barrier for OH. However, an accurate determination of each one of these quantities will be necessary for a full understanding of how complex organic molecules form in different temperature regimes in interstellar regions.

The rather elegant “warm-up” scenario for the origin of so-called “hot core” molecules has been muddied by recent observational and laboratory results. High sensitivity observations have shown that complex organic molecules are also present in the gas phase in cold environments (∼ 10 K), albeit at a low level with respect to molecular hydrogen (Öberg et al. 2010, 2011; Cernicharo et al. 2012; Bacmann et al. 2012; Vastel et al. 2014). Bisschop et al. (2007a) investigated the laboratory hydrogenation of solid acetaldehyde, CH₃CHO, under interstellar relevant conditions, using both RAIRS and TPD to analyze the results. The experiments showed that the hydrogenation of CH₃CHO leads to the formation of 20 % of ethanol, C₂H₅OH, showing for the first time that surface hydrogenation of unsaturated complex species can be responsible for the abundances of more complex saturated species detected in dense interstellar clouds. In addition, it has now been demonstrated experimentally that molecules such as glycolaldehyde (HOCH₂CHO) and ethylene glycol (HOCH₂CH₂OH) can form at 10 K without the need for irradiation (Fedoseev et al. 2015). Here, the reaction scheme also relies mainly on hydrogenation reactions which are known to be efficient at low temperatures; however, the scheme also involves dimerization of the HCO radical which allow the formation of the C–C backbone (see also Woods et al. 2013). These results suggest that complex molecules are already present in the ice mantles prior to warm up in the environs of young stars.

A few alternatives to the UV photo-induced surface chemistry hypotheses involve gas-phase reactions (Rawlings et al. 2013; Balucani et al. 2015). The so called “Rapid-Radical Association” mechanism proposes that COMs are formed by three-body gas-phase reactions between radicals in warm high density gas (Rawlings et al. 2013). This environment exists, for a very short period of time, following the sudden and total sublimation of dust-grain ice mantles driven by the catastrophic recombination of trapped hydrogen atoms, and other radicals, in the ice. More recently, Balucani et al. (2015) proposed a new gas-phase model to form some of the COMs, such as dimethyl ether and methyl formate, starting from methanol in the gas phase. In the proposed scheme, dimethyl ether is the precursor of methyl formate via an efficient reaction overlooked by previous models. Very recently, electronic structure and kinetic calculations showed that the gas-phase reaction ${\mbox{NH}_{2} + \mbox{H}_{2}\mbox{CO}}$ leading to formamide is barrierless. Hence, for some species, there is no need to invoke grain-surface chemistry provided that the necessary precursors are available in the gas phase (Barone et al. 2015). The main complicating factor in deciphering the chemistry of COMs is that their abundances and abundance ratios usually vary significantly from source to source. Moreover, not all the detected COMs and their abundances can be explained by gas-phase reactions alone. Therefore, it is likely that surface reactions on ice grains still play an important role in the dense regions of the ISM.

Figure 6 taken from Linnartz et al. (2015) summarizes the extensive laboratory work that investigated the surface formation of several simple and more complex species through atom-addition and radical-radical recombination reactions under dense cloud conditions. The figure shows how the surface reaction routes to H₂O, CO₂, HCOOH, CH₃OH, (CH₂OH)₂, HNCO, NH₃, and NH₂OH, are interlinked and are initiated by the H/O/N-atom addition to simpler precursor species, like CO molecules. The efficiency of these reaction routes, and their contribution to the molecular abundances observed in space, depends on a number of physico-chemical parameters. For instance, in atom addition/abstraction reactions, the height of an activation barrier determines whether or not a reaction can proceed at 10 K. Radical-radical recombination reactions are barrierless and therefore their inefficient thermal diffusion at 10 K is the limiting factor. Another important parameter is the molecular environment. In the case of exothermic reactions, polar (water-rich) ices can promote surface chemistry through the dissipation of extra energy in their H-bond network.

Fig. 6

A summary of non-energetic surface chemistry through atom-addition and radical-radical recombination reactions based on experimental evidence. The arrows indicate possible pathways, but other (energetic) processes are at play as well. The figure clearly shows the complexity of non-energetic ice chemistry and the possibility for this type of chemistry to create complex molecules without additional energy input. Figure reproduced from Linnartz et al. (2015)

6 Diffusion

As discussed in the previous section, diffusion rates are key to determining the rate for the Langmuir-Hinshelwood reaction mechanism, because they regulate the meeting frequency between reactants. Models require as input diffusion barriers and binding energies for each of the surface species included in the reaction network to determine the hopping rate

$$ k_{\mathrm{hop}, \mathrm{th}, i} = \nu \exp \biggl(- \frac{E_{\mathrm{diff}}}{kT} \biggr). $$

(26)

The scanning rate $k_{\mathrm{scan}}$ is determined by either thermal hopping or quantum mechanical tunneling, depending on the mass of the diffusing species.

In early models, H atoms were assumed to diffuse via quantum tunneling. After experimental studies showed that the diffusion of H atoms is rather slow, thermal hopping became favored in models (Pirronello et al. 1997b,a, 1999; Katz et al. 1999). More recently, the discussion on the nature of the diffusion mechanism for atoms was reopened, with experiments suggesting H atoms can quantum tunnel under particular conditions (Watanabe et al. 2010; Kuwahata et al. 2015). It has also been postulated that atoms as heavy as oxygen can diffuse via quantum tunneling (Congiu et al. 2014). Heavier species are most likely to diffuse through thermal hopping.

Obtaining diffusion barriers experimentally remains challenging. Diffusion rates cannot be measured directly and have to be inferred from experiments using a model. Often diffusion barriers are measured through reaction, where a known dose of both mobile reactants and relatively immobile (stationary) reactants are deposited. If the reaction between species is diffusion limited, i.e., there exists no reaction barrier and in the regime of low surface coverage, the diffusion rate can be inferred from the disappearance of the reactants as a function of temperature, dose, and/or time. This method is limited to reactive species and, because one works within the submonolayer regime, sensitivity of detection is an issue. Usually products are measured using mass spectrometry during TPD. Another method is to deposit a layer of ASW (amorphous solid water) on top of an ice consisting of the species of interest. If the temperature is raised above the desorption temperature of this species, the rate limiting step for desorption is the diffusion of the species through the ASW layer. In the case of porous ASW, pore-wall diffusion is probed which mimics surface diffusion. For compact ASW, bulk diffusion is more likely studied. Here, the disappearance of the diffusing species is monitored, usually by IR spectroscopic techniques. This method is limited to stable species with a desorption temperature below that of ASW. Surface diffusion rates for CO molecules have been studied this way (Mispelaer et al. 2013; Karssemeijer et al. 2014b; Lauck et al. 2015).

For both methods of study, the model that is used to extract the diffusion rates is crucial. In the case of reaction, one has to ensure that no other reactions play a role. Radical beams are often not 100 % pure and other species will also be present on the surface. Moreover, species might not be instantaneously thermalized when deposited on the substrate, so that they are able to move some distance superthermally. One has to be sure that the second reactant is indeed stationary. Furthermore, the effect of the warm-up phase during TPD needs to be considered which can thermally enhance diffusion (and subsequent reaction). In the case of the two-layer experiments, care should be taken when considering to which type of diffusion the system is limited: bulk diffusion, pore wall diffusion, or surface diffusion upon a layer of the same species adsorbed on the ASW (see Fig. 7).

Fig. 7

Schematic of the mixing process. (a) Layered system at $t = 0$. (b) Occurrence of mixing, with the inset showing the CO molecules diffusing along the micropore surfaces into the strong-binding nanopore sites. (c) After some period of time, the layer becomes fully mixed. Figure taken from Lauck et al. (2015)

The species that has been most studied to determine its surface diffusion is atomic hydrogen. Because it is a radical, it is studied through reaction, either with itself to form H₂, with D to form HD, or with O₂ to from H₂O₂, and ultimately H₂O. However, the results between different experimental groups are not in agreement (see Hama and Watanabe 2013, for a review). The latest results on an ASW surface show that very shallow-potential binding sites ($E_{\mathrm{diff}} \leq 18 \pm 2~\mbox{meV}$) are dominant, while there are also middle- ($E_{\mathrm{diff}} = 22 \pm 1~\mbox{meV}$) and deep- ($E_{\mathrm{diff}} = 30~\mbox{meV}$) potential binding sites (Watanabe et al. 2010; Hama et al. 2012). These particular results were obtained with a new technique for the detection of product species. The remaining species are photodesorbed by a laser as a function of delay time between deposition and laser desorption and the detection occurs through REMPI (resonance-enhanced multi-photon ionization). This eliminates the effect of warm-up which can be an issue for TPD experiments. These experiments support the following picture (Watanabe et al. 2010; Hama et al. 2012; Kuwahata et al. 2015): H atoms land on the ASW surface, diffuse rapidly using the abundant shallow and middle sites, and are finally trapped in the deepest sites. Once a significant number of H atoms become trapped, the subsequent H atoms recombine with the trapped H atoms. Therefore, the effective diffusion rate becomes dependent upon coverage and is also affected by the number of deep binding sites present on the ice surface.

This picture is also observed at an atomistic level through adaptive kinetic Monte Carlo simulations of CO diffusion on an ASW surface (Karssemeijer et al. 2014b). Here, for a single CO molecule on an ASW surface, the diffusion is limited by diffusion out of the deep binding sites ($E_{\mathrm{diff}}^{\mathrm{CO}} = 84\mbox{--}114~\mbox{meV}$), whereas when the deep sites are filled with CO molecules, the diffusion rate is increased ($E_{\mathrm{diff}}^{\mathrm{CO}} = 48\mbox{--}79~\mbox{meV}$). In the simulations, a large spread in diffusion barriers between different surfaces is found, in addition to the large spread observed for each specific surface. However, given a specific surface and surface coverage, the temperature-dependent diffusion constant follows an Arrhenius-like behavior, i.e., it can be described by a single diffusion barrier, which is consistent the diffusion barrier for moving out of the deepest available binding site. This is good news for astrochemical gas-grain models because, for simplicity, they typically do not account for local chemical or structural differences on the grain. Instead, these models aim to provide a macroscopic view. For this purpose, a single diffusion barrier per species likely suffices. A collection of different barriers and binding energies obtained in this way is presented in Table 4.

Table 4 Desorption and diffusion energy barriers for CO and CO₂ on the proton-disordered and the ordered Fletcher phase of ice Ih

Table 4 is far from complete, and to obtain the remaining diffusion barriers required for the astrochemical models, the diffusion energy of a species is assumed to be a universal, fixed fraction, $f$, of the desorption energy: $E_{\mathrm{diff},i} = f E_{\mathrm{bind},i}$. There is no fundamental physical argument for such a universal ratio to exist and it is used solely due to the lack of data. This ratio is most likely dependent upon the species, the substrate material and structure, and on the surface coverage (which determines the likely binding “partner”). There already exist problems with the concept of using a single diffusion and a single desorption barrier even for a single species because they both vary strongly from site to site, especially for amorphous ices. As shown by Vasyunin and Herbst (2013) the value of $f$ influences the outcome of the models. Due to the lack of diffusion information, the existence and possible value of fraction $f$ is very poorly constrained and values between 0.3 and 0.8 are used by the modeling community and is often treated as a free parameter within these typical bounds (Hasegawa et al. 1992; Ruffle and Herbst 2000; Cuppen et al. 2009). An accurate value of $f$ is also important because it affects the conditions under which the accretion limit is reached and thereby whether or not modelers should turn to stochastic models. The data in Table 4 suggests that, at least for stable species, there is a more or less constant ratio and that this ratio, $f$, is more likely to lie around $0.3\mbox{--}0.4$. It is not possible, currently, to give a definitive recommended value. We encourage the experimental and theoretical community to continue to work towards filling this large gap in the necessary input data for astrochemical models, especially for radical species, which, to date, have remained largely unstudied.

7 Bulk Processes

Although the chemistry on interstellar grains starts out with processes on bare carbonaceous or silicate grains, when the ice thickness increases to more than a few monolayers, it becomes important to differentiate between surface and bulk processes. The earliest models of grain-surface chemistry primarily considered the grain surface as a substrate for the formation of molecular hydrogen or other simple species, which could then rapidly desorb back into the gas phase upon recombination (e.g., Watson and Salpeter 1972a; Allen and Robinson 1977; Tielens and Hagen 1982). Hence, originally, ice chemistry was modeled using rate equations that consider only the averaged abundance of a species throughout the entire mantle (Hasegawa et al. 1992). This is a rather crude approximation because bulk species cannot diffuse as easily as those on the surface. “Bulk ice” implies that the species involved in “bulk” processes are fully surrounded by neighboring molecules and are therefore rather tightly bound, leading to the assumption that diffusion within the bulk ice at low temperature is inefficient and therefore, chemistry is also inhibited. The ice mantle can thus be conceptually divided into an ice surface and bulk ice. To this end, three-phase models have been introduced, where the distinction is made between gas-phase species, reactive surface species within the top (few) monolayer(s), and fully inert bulk species in the core of the ice mantle, with terms that allow surface material to be incorporated into the bulk (or vice versa), as the ice thickness increases or decreases (Hasegawa and Herbst 1993b; Garrod and Pauly 2011). Each phase is treated mathematically as an independent entity. More recent models, informed by the discovery of substantial ice mantles on interstellar dust grains via infrared absorption observations, have involved a more concerted effort to treat the build-up of surface ices (e.g., Hasegawa et al. 1992). It is in the extension of these models to the formation of multiple layers of ice that a significant conceptual error should be identified. The absolute reaction rates for surface processes are commonly formulated thus:

$$ f_{\mathrm{react}}(A+B) = n_{\mathrm{s}}(A) n_{\mathrm{s}}(B) \bigl[k_{\mathrm{hop}}(A) + k_{\mathrm{hop}}(B) \bigr] / N_{\mathrm{sites}} = n_{\mathrm{s}}(A) n_{\mathrm{s}}(B) k_{\mathrm{scan}}(AB) $$

(27)

where $k_{\mathrm{hop}}$ is the site-to-site hopping rate for either species A or B, and $N_{\mathrm{sites}}$ is the number of surface sites present on the grain surface. $n_{\mathrm{s}}(A)$ and $n_{\mathrm{s}}(B)$ are used here as a shorthand to represent the mean populations of each species, $\langle n_{\mathrm{s}}(A)\rangle$ and $\langle n_{\mathrm{s}}(B)\rangle$. This formulation is often further simplified, absorbing all coefficients into a single rate, $k_{\mathrm{scan}}(AB)$. The efficiency factor used for reactions involving activation energy barriers is omitted here for simplicity.

Equation (27) is valid in the case where the sum of all surface reactants, which we label $n_{\mathrm{all}}$, is smaller than $N_{\mathrm{sites}}$, i.e. that there is less than one layer of particles on the grain surface. However, if $n_{\mathrm{all}} > N_{\mathrm{sites}}$, the equation is necessarily invalid. This may be seen more clearly if Eq. (27) is re-arranged:

$$ f_{\mathrm{react}}(A+B) = n_{\mathrm{s}}(A) k_{\mathrm{hop}}(A) n_{\mathrm{s}}(B) / N_{\mathrm{sites}} + n_{\mathrm{s}}(B) k_{\mathrm{hop}}(B) n_{\mathrm{s}}(A) / N_{\mathrm{sites}} $$

(28)

It may be seen that the reaction rate is composed of two analogous parts. The first represents the rate at which all surface species of type A may hop into an adjacent site, multiplied by the probability that the neighboring site is occupied by a species of type B. The second term may be described similarly. This latter probability is the simple ratio $n_{\mathrm{s}}(B)/N_{\mathrm{sites}}$. Clearly, if $n_{\mathrm{s}}(B)$ exceeds $N_{\mathrm{sites}}$, such a probability is meaningless. More rigorously, this is also true where $n_{\mathrm{all}} > N_{\mathrm{sites}}$.

A majority of two-phase models (i.e. gas-phase and grain surface) retain the usage of $N_{\mathrm{sites}}$ in the reaction rates. In the case of an ice mantle with a thickness of, say, 100 monolayers, this could lead to a rate that is inaccurate by two to four orders of magnitude.

The simplest fix for this problem is to exchange $N_{\mathrm{sites}}$ for $n_{\mathrm{all}}$, which may be calculated easily within the model code, for cases where $n_{\mathrm{all}} > N_{\mathrm{sites}}$. This approach removes the immediate inconsistency in the probabilities; however, it retains another assumption implicit in the equations above, namely that not only are the species in all layers of the ice chemically active, but that all such reactants may diffuse within the bulk ice at the surface diffusion rate. Neither of these assumptions should necessarily be used, as explained elsewhere in this paper.

To remove the latter problem while retaining the structure of a two-phase model, one may adjust Eq. (28) further, by retaining $N_{\mathrm{sites}}$, but substituting $n_{\mathrm{s}}(A)$ for $n_{\max}(A)$, where:

$$\begin{aligned} n_{\max}(A) &= n_{\mathrm{s}}(A); \quad\mbox{for } n_{\mathrm{all}} \le N_{\mathrm{sites}} \end{aligned}$$

(29)

$$\begin{aligned} n_{\max}(A) &= n_{\mathrm{s}}(A) N_{\mathrm{sites}}/n_{\mathrm{all}}; \quad\mbox{for } n_{\mathrm{all}} > N_{\mathrm{sites}} \end{aligned}$$

(30)

And similarly for $n_{\max}(B)$, giving

$$ f_{\mathrm{react}}(A+B) = n_{\max}(A) n_{\max}(B) \bigl[k_{\mathrm{hop}}(A) + k_{\mathrm{hop}}(B)\bigr] / N_{\mathrm{sites}} $$

(31)

This solution removes the error in the probabilities, while limiting the reactive portion of the ice mantle to a single layer, on the assumption that a layer is composed of a total of $N_{\mathrm{sites}}$ particles. Such an approach may be considered a quasi-three-phase model, although the chemically-active portion of the mantle, which is assumed to be the surface layer itself, nevertheless represents the composition of the entire mantle. Only a fully realized three-phase model, consisting of entirely separate phases for the gas, ice mantle, and ice-surface layer can resolve this point. However, the added complexity of the physical treatment, combined with technical challenges (see below), makes three-phase modeling of astrophysically complex objects such as protoplanetary disks including bulk chemistry a more distant prospect, as outlined in Sect. 2.

Another approximation that can be made is the multilayer approach by Taquet et al. (2012, 2013), where the abundances of each single layer are stored in the simulation and these compositions are used to determine the binding energy of the species in the overlying layer. Again only the top monolayer is considered chemically active. Furthermore, recently another multilayer model (using a Monte Carlo approach) has been reported where the chemically active surface is extended to the uppermost four monolayers (Vasyunin and Herbst 2013). Both models account for a correct description of thermal desorption, but still leaves the bulk reactively inert. In microscopic Kinetic Monte Carlo routines, on the other hand, bulk reactivity can be included, for instance through allowing two neighboring molecules to swap position (Öberg et al. 2009d) or incorporating interstitial positions between the molecules in the lattice (Lamberts et al. 2013; Chang and Herbst 2014).

A number of experimental and theoretical studies indicate that bulk processes can be important, both at low and high ice temperatures. Here we make the distinction between bulk diffusion and pore-wall diffusion, since the latter is essentially surface diffusion. Experimentally it is not always straightforward to distinguish the two effects. Bulk diffusion appears to play a role in the segregation of mixed H₂O:CO₂ and H₂O:CO ices which leads to the occurrence of “islands” of both mixture components at temperatures of 30 K for CO and 60 K for CO₂ (Öberg et al. 2009d). Exposure of an O₂ ice to hydrogen atoms yields large quantities of H₂O₂ and H₂O, corresponding to tens of monolayers depending on the temperature (Ioppolo et al. 2010). For hydrogenation of CO, a maximum exposure of four monolayers was found experimentally (Fuchs et al. 2009). Furthermore, Molecular Dynamics simulations performed by Andersson et al. (2006) and Arasa et al. (2010) of photodissociation and photodesorption of water ice at low temperature (10–90 K) have shown that the desorption of species (H, OH, H₂O) depends strongly on the distance of the excited molecule to the top monolayer of the ice. The first 4–5 monolayers are, however, taking part in the desorption process, which indicates again that the assumption of a single reactive top monolayer cannot be valid. The desorption depth dependence has been observed for a range of excitation energies, ice temperatures and isotope compositions (Andersson et al. 2006; Andersson and van Dishoeck 2008; Arasa et al. 2010; Koning et al. 2013; Arasa et al. 2015). The studies further show that radical species which are created through photodissociation remain in the bulk of the ice, can use their excitation energy to move some short distance before they thermalize. This will most likely generally hold for radicals formed through an energetic process like photodissociation, exothermic reactions and through cosmic rays.

An analog for hydrogen diffusion at low temperature can also be brought about by reactions that have a net neutral effect, of the type ${\mbox{H--A} + \mathrm{A} \longrightarrow \mathrm{A} + \mbox{H--A}}$ as mentioned during the Lorentz Center Workshop by Fedoseev. This could be of interest, for instance, in the case of water ices where photodissociation can create OH radicals that may react with neighboring H₂O molecules, ${\mbox{H}_{2}\mbox{O} + \mbox{OH} \longrightarrow \mbox{OH} + \mbox{H}_{2}\mbox{O}}$ (Nguyen and Stanton 2013; Lamberts et al. 2016).

Recently, the discovery of a new class of thin films, spontelectrics, may also have interesting implications for diffusion and ordering in ISM ices (Field et al. 2013). Upon deposition onto a substrate at low temperatures, dipolar molecules were seen to spontaneously align, creating electric fields $> 10^{8}~\mbox{eV}\,\mbox{m}^{-2}$. The impact of such spontaneous ordering within ISM-like ices is yet to be quantified. An alternative mechanism that could work at high temperatures, is to allow the reactants to reach each other before the ice is close to evaporating, leaving the mantle molecules more freedom to diffuse.

Diffusion through the walls of macropores is less relevant for interstellar ices with respect to laboratory analogs. Interstellar water ice is mostly formed through reactions instead of deposition resulting in compact ice (Oba et al. 2009) by the exothermicity of surface reactions (Accolla et al. 2011). Moreover, laboratory experiments of deposited ices indicate that the macropores present in these ices collapse at least to a certain extent (Bossa et al. 2012; Mispelaer et al. 2013; Isokoski et al. 2014), making pores less relevant. Moreover, there is a class of reactions that are not (merely) diffusion limited, since one of the reactants is one of the main mantle species, but (mainly) thermally activated at higher temperatures. This is possible in particular for reactions between two neighboring species, for instance, in a water-rich environment. Examples of such reactions are ${\mbox{H}_{2}\mbox{O} + \mbox{D}_{2}\mbox{O} \longrightarrow 2 \mbox{HDO}}$ (Lamberts et al. 2015), and ${\mbox{H}_{2}\mbox{O} + \mbox{HNCO} \longrightarrow \mbox{H}_{3}\mbox{O}^{+} + \mbox{OCN}^{-} \longrightarrow \mbox{HOCN} + \mbox{H}_{2}\mbox{O}}$ (Theule et al. 2011). Other thermally activated reactions in an ice are possible, leading to the formation of salts (e.g., Bossa et al. 2009; Noble et al. 2014; Bergner et al. 2016). An example is the following network of equilibrium reactions

$${\text{NH}}_3+{\text{CO}}_2\rightleftharpoons {\text{NH}}_3+{\text{COO}}^{-}$$

(32)

$${\text{NH}}_3+{\text{NH}}_3+{\text{COO}}^{-}\rightleftharpoons {\text{NH}}_4+{\text{NH}}_2{\text{COO}}^{-}$$

(33)

$${\text{NH}}_4+{\text{NH}}_2{\text{COO}}^{-}\rightleftharpoons {\text{NH}}_3+{\text{NH}}_2\text{COOH}$$

(34)

leading to the formation of carbamic acid in a protic environment, i.e. a $\mbox{NH}_{3-}$ or H₂O-dominated ice, where the reaction barrier is sufficiently low for proton transfer to occur thermally without any external source of processing. However, experimentally only mixtures with a NH₃:H₂O ratio greater than 1 were considered. In a more dilute environment, which is more representative for interstellar ices these processes become bulk diffusion limited and they require a suitable bulk diffusion mechanism in astrochemical gas-grain models.

Observations of interstellar ices show absorption bands attributed to various cations and anions (OCN⁻, HCOO⁻, $\mbox{NH}_{4}^{+}$, Boogert et al. 2015), and hence these processes should occur. Moreover, it is known that interstellar (and circumstellar) grains are likely charged (e.g., Draine 2003). Grains can become the main carrier of negative charge, for instance in dark and dense midplanes, where the gas is weakly ionized. As grain charge affects the chemistry, it is important to consider. Ionic chemistry in most grain surface models is currently limited to cation-grain recombination, where an ion recombines with a free electron and dissociates rather than sticking to the grain, applying gas-phase products and branching ratios analogous to those from dissociative recombination

$$ {\mbox{XY}^{+} + \mbox{GRAIN}^{-} \longrightarrow \mbox{X} + \mbox{Y} + \mbox{GRAIN}}. $$

(35)

The rate of this process is determined by the cation-grain collision rate and is analogous to the accretion rate (Eq. (2)) but can include an enhancement (e.g., Draine and Sutin 1987).

8 Photoprocesses

It has long been known experimentally that FUV irradiation of interstellar ice analogues induces complex chemistry (e.g., Hagen et al. 1979; D’Hendecourt et al. 1986; Allamandola et al. 1988; Grim et al. 1989; Gerakines et al. 1996; Öberg et al. 2009c). The interaction of FUV photons with ices has both primary and secondary chemical effects which influence both the ice and gas composition. Which of these effects dominates depends upon the FUV flux and shape of the FUV spectrum, the photoexcitation/photodissociation cross sections of the molecular components which make up the ice mixture, and the temperature and structure of the ice. One of these effects is termed “photodesorption”, whereby absorption of FUV (or X-ray) photons by an ice leads to desorption (or sublimation) of ice species into the gas phase. Photodesorption was originally thought to occur as a single-step process: the FUV photon excites the molecule into a state which is antibonding (or repulsive) with the surface and the species is ejected. As will be discussed later, it is now known that it can also occur as a multi-step process where the excited species decays back to the ground state and imparts kinetic energy to a surrounding molecule which is then “kicked out”. Another effect is photodissociation, whereby the absorption of a FUV photon by a molecule induces dissociation, e.g., ${\mbox{H}_{2}\mbox{O} + \mbox{h}\nu \longrightarrow \mbox{OH} + \mbox{H}}$. It is this latter process which induces chemical changes in FUV-irradiated ices. As might be expected, photodesorption and photodissociation can occur in parallel. Photodissociation has primarily been induced to explain the origin of complex molecules in FUV-irradiated ices with an initial interstellar-like composition (e.g., CO, CO₂, H₂O, CH₄, N₂, NH₃, and CH₃OH).

Photodesorption is now considered an important process for releasing molecules into the gas phase in cold (≲ 100 K) astrophysical regions where they would otherwise be frozen out onto icy grain mantles (and thus wholly depleted from the gas, e.g., Hartquist and Williams 1990). The origin of the FUV radiation depends upon the nature of the region. In cold and dark molecular clouds, FUV radiation is generated via the excitation of H₂ molecules by cosmic rays (Prasad and Tarafdar 1983), whereas in weak photodissociation regions (PDRs), the origin is the external, often enhanced, interstellar radiation field (ISRF, Hollenbach and Tielens 1997). In the envelopes surrounding forming stars, the innermost regions are sufficiently hot that thermal desorption is thought to dominate; however, molecules can be released back into the gas in the extremities of the envelope by photodesorption by the ISRF, similar to PDRs (e.g., van Dishoeck and Blake 1998). Recent high-sensitivity observations have suggested that the origin of gas-phase molecules at offset positions from the location of the central forming star may be due to photodesorption by stellar FUV photons along the outflow cavity walls (e.g., Öberg et al. 2010). In protoplanetary disks around young stars, there are multiple origins, including FUV photons from the central star and the external interstellar radiation field, and those generated internally via the excitation of H₂ by cosmic rays or X-rays (e.g., Willacy and Langer 2000). Which source dominates the FUV radiation field depends on the radial distance away from the central star and the vertical distance away from the shielded disk midplane.

Given the potential importance of photodesorption when interpreting observations, this process has routinely been included in astrochemical models for some decades. However, it is only relatively recently that experiments have successfully quantified this process, and also that the underlying physical mechanisms at work have begun to be understood. Photoprocesses have only been included in models relatively recently once it was demonstrated that molecules, such as methanol (CH₃OH), were unable to form efficiently via gas-phase chemistry alone.

In this Section, we give a brief overview of the experimental studies specific to ice photodesorption and photodissociation (Sect. 8.1), and those theoretical calculations which have offered insight into the mechanisms at play (Sect. 8.2). We also provide a recommendation for the implementation of experimental and theoretical data into the computation of photodesorption and photodissociation rates in astrochemical models (Sect. 8.3).

8.1 Experimental Studies

The measured photodesorption yields for pure ices are listed in Table 5. Early experiments on the photodesorption of molecules of astrophysical interest resulted in yields spanning from ${\sim}10^{-6}~\mbox{photon}^{-1}$ (Greenberg 1973), to ${\sim}10^{-1}~\mbox{photon}^{-1}$ (Nishi et al. 1984; Hartquist and Williams 1990). The large range of values is likely due to the wide array of setups adopted and the challenges of direct detection of the photodesorbed molecules. The former experiments were conducted in the 100–275 nm range at 77 K using a Hg-Xe high pressure arc lamp, whereas the latter were conducted using an ArF 193 nm laser at relatively high fluence between 90 and 135 K. It was not until the experiments by Westley et al. (1995a) and Westley et al. (1995b) that a robust quantification of water ice photodesorption at low temperatures ($35\mbox{--}100~\mbox{K}$) was reported. They determined a yield of $(3\mbox{--}8\times10^{-3})~\mbox{photon}^{-1}$ at Lyman-$\alpha$ (121.6 nm). Photodesorption was revisited much later by (Öberg et al. 2007, 2009a,b) who investigated the photodesorption yields of pure CO, N₂, CO₂, and H₂O ices under ultra-high vacuum conditions (base pressure better than $10^{-9}~\mbox{mbar}$) at low temperatures using a broadband hydrogen microwave discharge lamp ($120\mbox{--}170~\mbox{nm}$). CO, CO₂, and H₂O were found to desorb from thick ices ($\gg10$ monolayers) at low temperatures ($< 20~\mbox{K}$) with similar yields: ($2.7\pm1.7\times10^{-3}$), ($2.3\pm1.4\times10^{-3}$), and $(1.3\pm0.8\times10^{-3})~\mbox{photon}^{-1}$. Pure N₂ ice, on the other hand, had a significantly lower yield $(<2\times10^{-4})~\mbox{photon}^{-1}$. Öberg et al. (2009a) also provided evidence that photodesorption occurs primarily from the top 1–2 monolayers of the ice mantle. For a CO ice covered in a thin layer of N₂, the signal for CO photodesorption is decreased whilst that for N₂ is enhanced, and above that expected for pure N₂ ice. This was the first evidence of co-desorption which we now term “kick out”. Öberg et al. (2009a,b) were also able to determine that photodesorption is a parallel process to photodissociation. For the pure H₂O ice experiments, they observe a product yield for H₂O:OH of 0.7:1.0 at 20 K to 1.7:1.0 at 100 K. The temperature dependence arises from the increased mobility of the photodissociation products at a higher temperature allowing more efficient reformation of the water ice. For pure CO₂ ice, they find CO₂ desorbs as CO with a yield between 20–50 % that of CO₂. Again, at higher temperatures, the ratio of CO₂:CO increases due to increased mobility of the fragments. Hence, for both H₂O and CO₂, the primary process is photodissociation, followed by (i) desorption of the energetic fragments (OH or CO), (ii) recombination of the fragments which have excess energy to desorb, or (iii) “kick out’’’ of a neighboring molecule by an energetic fragment. Pure O₂ ice was also recently studied by Zhen and Linnartz (2014) under similar experimental conditions, and was found to have a desorption yield $(6.0\pm2.0\times10^{-4})~\mbox{photon}^{-1}$. Ozone (O₃) was also detected (arising from photodissociation of O₂ and reaction between O₂ and O) and similar to CO₂, the yield had a weak temperature dependence where the yield increased at higher temperatures.

Table 5 Measured photodesorption yields for pure ice species and products (where determined)

Later experimental efforts moved towards wavelength-dependent studies. These measurements are more challenging and are conducted at synchrotron facilities where the FUV flux is lower, requiring longer irradiation times to match the fluence typical of interstellar conditions. On the other hand, microwave discharge lamps have a significantly higher flux; however, independent measurements of CO photodesorption have shown that the derived yield is dependent on having a well-calibrated lamp spectrum. Muñoz Caro et al. (2010) determined a desorption yield for pure CO ice of $(3.5\pm0.5\times10^{-2})~\mbox{photon}^{-1}$ which is one order of magnitude higher than that from Öberg et al. (2009a). Recent investigations showed that the FUV spectrum, in particular the ratio of Lyman-$\alpha$ (121.6 nm) to H₂ molecular emission bands (110–180 nm), is sensitive to the gas pressure, the type of lamp, the gas composition (pure H₂ versus a mixture of H₂ and He) and the window properties (Chen et al. 2014; Ligterink et al. 2015). Chen et al. (2014) determined CO desorption yields as high as $(2.12\pm0.03\times10^{-1})~\mbox{photon}^{-1}$ for a lamp setup with pure H₂ and a window giving a ratio of 20 % flux of Lyman-$\alpha$ relative to the molecular emission bands. This showed that desorption of CO at Lyman-$\alpha$ is less efficient than at wavelengths covering the H₂ emission bands and turned out to be the first step in quantifying the wavelength dependence of photodesorption. The same group have recently reported VUV (120–160 nm) absorption cross-sections for a plethora pure ices (Cruz-Diaz et al. 2014a,b). All data are publicly available for download.¹

Fayolle et al. (2011) and Bertin et al. (2012) presented the first wavelength-dependent measurement of CO photodesorption (between 90 and 160 nm) conducted at the SOLEIL² synchrotron facility. The results showed a strong wavelength dependence corresponding to the vibronic transition spectrum of CO ice. CO was thus confirmed to desorb via the DIET mechanism: Desorption Induced by Electronic Transitions. The results also confirmed the low desorption yield at Lyman-$\alpha$ wavelengths. Peak yields approaching ${\sim}5\times10^{-2}~\mbox{photon}^{-1}$ were found to correspond with the $A^{1}\varPi - X^{1}\varSigma^{+}$ electronic transitions of CO (${\sim}130\mbox{--}160~\mbox{nm}$). Followup studies on pure N₂, O₂, and CO₂ ices showed interesting results (Fayolle et al. 2013; Fillion et al. 2014). The pure N₂ photodesorption yield was found to peak at significantly shorter wavelengths ($\lesssim100~\mbox{nm}$) with a wavelength dependence corresponding to the $b^{1}\varPi_{u}-X^{1}\varSigma^{+}_{g}$ electronic transition. A peak yield of $3.77\pm1.13\times10^{-2}~\mbox{photon}^{-1}$ was found at 100.8 nm. Photodesorption yields at longer wavelengths were consistent with previous measurements using broad-band hydrogen lamps (Öberg et al. 2009a). Pure O₂ ice photodesorption was also studied and, in contrast to CO and N₂, was found to have a non-negligible contribution due to photodissociation (Fayolle et al. 2013), as also seen later in the broad-band photodesorption experiments (Zhen and Linnartz 2014). The peak photodesorption rates (${\sim}0.5\times10^{-2}~\mbox{photon}^{-1}$) overlapped with the Schumann-Runge band for O₂ ($\tilde{X}^{3}\varSigma_{g}^{-}-\tilde{B}\varSigma_{u}^{-}$) which leads to dissociation into $\mbox{O}({}^{3}\mbox{P})$ and $\mbox{O}({}^{1}\mbox{D})$ below ${\sim}176~\mbox{nm}$ (Parker 2000; Mason et al. 2006). Fillion et al. (2014) reported results from pure CO₂ ice. CO₂ showed a more complicated behavior than the more volatile species discussed above. Both CO and CO₂ molecules were found to photodesorb at wavelengths below 118 nm only, with yields ranging from ${\sim}1\mbox{--}3\times10^{-3}~\mbox{photon}^{-1}$. The steep increase in photodesorption signal at 118 nm is associated with the first dipole-allowed transition ($(3p\sigma_{u})^{1}\varSigma_{u}^{+}- ^{1}\varSigma_{g}^{+}$) and thereafter follows the photoabsorption continuum. That both CO and CO₂ were seen with similar photodesorption yields suggests that the mechanism likely has multiple paths: (i) photoexcitation of CO₂ which “kick out” surface CO₂ (or CO) molecules (DIET mechanism), (ii) photodissociation of CO₂ generating energetic fragments ($\mbox{CO}^{*} + \mbox{O}$) which either “kick out” surface CO₂ (or CO) molecules (DIET mechanism), or desorb directly, or recombine to reform CO₂ which then desorbs. In all cases, the photodesorption signature of CO₂ and CO were similar and as described above. However, for pre-irradiated CO₂ ice, the photodesorption signatures of CO₂ and CO were found to match that of pure CO ice (Fillion et al. 2014). Here, the DIET mechanism is driven by photoexcitation of CO (as opposed to photoexcitation or photodissociation of CO₂).

Wavelength-dependent results from binary N₂:CO ice mixtures confirmed the desorption of surface species via sub-surface excitation (Bertin et al. 2013). The resulting photodesorption signature of surface molecules was found to correspond with the dominant component of the bulk ice composition, e.g., a thin surface layer of N₂ on top of CO ice photodesorbs with a signature of CO photodesorption whereas for a 1:1 mixture, both molecules photodesorb with the signature of both species. Also, Bertin et al. (2013) found that a thin layer (${\sim}3$ monolayers) of N₂/CO on top of CO/N₂ bulk ice was sufficient to completely quench the CO/N₂ photodesorption signal.

The quantification and understanding of the process of photodesorption has made leaps and bounds in the past few years and for species other than CO and N₂, it has been shown that irradiation also induces photodissociation, that is, photodesorption and photodissociation are parallel processes. For pure (and relatively) simple ices, the resultant products and reaction pathways are straightforward to ascertain, as described above. Many experiments have shown that FUV irradiation of mixed ices induce chemistry; however, these experiments often couple irradiation with thermal processing and it can be difficult to disentangle reaction pathways specific to either FUV or heat processing (e.g., Hagen et al. 1979; D’Hendecourt et al. 1986; Allamandola et al. 1988).

Several experiments have conducted in-situ monitoring of ice composition via infrared spectroscopy during irradiation at low temperatures. Gerakines et al. (1996) irradiated numerous pure ices composed of known interstellar molecules with a broad-band hydrogen lamp and monitored changes in the ice composition with a Transmission Fourier-transform infrared spectrometer. Gerakines et al. (1996) reported a multitude of identified products, even for relatively simple ices (e.g., H₂O, NH₃, and CH₄). They also studied more complex ices: formaldehyde (H₂CO) and methanol (CH₃OH). In the former experiment, products such as polyoxymethylene (POM, $(\mbox{--H}_{2}\mbox{CO--})_{n}$), CO, HCO, CO₂, CH₄, CH₃OH, and methyl formate (HCOOCH₃), were identified after 1 hour of irradiation. For pure CH₃OH irradiation, similar products were identified (excluding POM). Gerakines et al. (1996) were also the first to extract quantitative information from the experiments and provided photodissociation cross sections (averaged over the lamp spectrum) for those species for which photodissociation was deemed efficient, and proposed potential reaction pathways.

The irradiation of pure and mixed methanol ices was revisited by Öberg et al. (2009c). They also studied the influence of temperature and FUV fluence and monitored the composition of the ice in-situ via RAIRS (Reflection-Absorption Infrared Spectroscopy). In all experiments conducted, complex molecules were observed to form and included CH₃CHO (acetaldehyde), CH₃CH₂OH (ethanol), CH₃OCH₃ (dimethyl ether), HCOOCH₃ (methyl formate), HOCH₂CHO (glycolaldehyde), and (CH₂OH)₂ (ethylene glycol). The exact products and relative ratios were found to have a weak dependence on temperature during irradiation and a stronger dependence on ice composition (pure CH₃OH versus CH₃OH mixed with CO or CH₄). For example, CO:CH₃OH ices were more enriched with $\mbox{--HCO}$ molecules while CH₄:CH₃OH ice were more enriched in $\mbox{--CH}_{3}$ molecules. Öberg et al. (2009c) also observed a large increase in the abundances of several complex molecules (around an order of magnitude) upon warming the irradiated sample from 20 to 70 K, confirming that radical mobility is temperature dependent and is important for building additional complexity. An analysis of the products allowed them to identify potential reaction pathways and formation cross sections. They also estimated the branching ratios for CH₃OH ice photodissociation at 20 K as $5:1:<1$ for ${\mbox{CH}_{2}\mbox{OH}} + \mbox{H}$, ${\mbox{CH}_{3}\mbox{O}}+\mbox{H}$, and $\mbox{CH}_{3} + \mbox{OH}$. The authors plan an extraction of fundamental data from the experimental results, for example, photodissociation rates, products, and branching ratios, diffusion rates, and binding energies of radical species. At the time of writing (March 2016), this challenging analysis was still ongoing.

As a final note in this Section, Öberg et al. (2009c) also report a photodesorption yield for CH₃OH of $(2.1\pm1.0\times10^{-3})~\mbox{photon}^{-1}$ which was determined by modelling the loss of methanol ice during irradiation as a combination of photodesorption (a zeroth-order process) and photodissociation (a first-order process). Very recently, several groups have revisited this system in light of recent revelations regarding the importance of robust lamp characterization for accurate determinations of photodesorption yields (Chen et al. 2014). The general consensus from several independent groups, is that pure methanol ice does not photodesorb intact. Photofragments only were observed, with a deduced upper limit to intact methanol photodesorption ${\sim} 10^{-6} - 10^{-5}$ molecules per incident photon, i.e., at least two orders of magnitude less than that found for water ice (Cruz-Diaz et al. 2016; Bertin et al. 2016). This low photodesorption yield was found for both pure and mixed (with CO) methanol ice, showing also that the DIET mechanism appears not to work for large molecules such as methanol. Note that this low photodesorption yield for methanol now calls into question the origin of gas-phase methanol in both cold and irradiated environments (e.g., dark clouds, PDRs, and protoplanetary disks), providing an intriguing challenge for molecular astrophysicists and astrochemical modelers.

8.2 Theoretical Calculations

Photodesorption has been studied theoretically via classical molecular dynamics (MD) simulations. Most work to date has been conducted on water ice (Andersson et al. 2005, 2006; Andersson and van Dishoeck 2008; Arasa et al. 2010), including its deuterated analogues, HDO and D₂O (Koning et al. 2013). Andersson et al. (2005, 2006) studied the photodesorption mechanisms of H₂O in both crystalline and amorphous water ice at 10 K following photoexcitation of H₂O into the first excited state with photons with wavelengths between $130\mbox{--}150~\mbox{nm}$. The MD calculations predicted a H₂O photodesorption efficiency between 0.5–1 % per absorbed photon in top few monolayers of an amorphous water ice surface (Andersson et al. 2005, 2006), which was in reasonable agreement with experimental data available at that time (e.g., Westley et al. 1995a). Andersson and van Dishoeck (2008) went a step further than earlier work and determined the relative contributions of all possible outcomes as function of depth into the ice. Upon photodissociation into $\mbox{H}+\mbox{OH}$, the following outcomes are possible: (i) H desorbs and OH is trapped, (ii) H is trapped and OH desorbs, (iii) H and OH are trapped, (iv) H and OH recombine and desorb as H₂O, and (v) H and OH recombine and are trapped as H₂O. There is also a parallel process whereby the excited H product had sufficient energy to “kick out” a neighboring water molecule. Outcomes (i)–(iv) and the “kick out” processes all have the potential to affect the chemical balance of both the ice and gas. Andersson and van Dishoeck (2008) determined that outcome (i) has the highest probability in the surface monolayer (${\sim} 90~\%$), and decreases in efficiency rapidly with depth into the ice reaching ${\sim}10~\%$ efficiency by the sixth monolayer. On the other hand, trapping of either the fragments or the reformed molecule increased in probability with depth, becoming the most probable outcome by the forth monolayer. The MD calculations also provided insight into the behavior of the photodissociation fragments post dissociation. For outcomes in which no desorption occurred, H and OH had sufficient excess energy to travel up to 70 and 60 Å, respectively, on or within the uppermost monolayers of an amorphous ice surface (Andersson and van Dishoeck 2008). Thus photodissociation of ices also creates mobile radicals at 10 K which can take part in further surface reactions.

Arasa et al. (2010) extended the calculations to higher temperatures (up to 90 K). Similar results were found: photodesorption has the highest efficiency in the top few monolayers, and trapping increases in probability with depth into the ice. The efficiency of OH and H₂O photodesorption was found to have weak dependence on temperature with an increase in probability of around $30~\%$. The probability of the parallel “kick out” mechanism was found to have a more complex behavior with temperature which the authors attribute to the amorphous nature of the ice. Koning et al. (2013) studied isotopic effects, and looked at the astrophysically relevant cases of HDO and D₂O photodissociation in an amorphous H₂O water ice. D atoms were found to photodesorb with a slightly lower efficiency than H atoms. On the other hand, the “kick out” probability of D atoms was found to be higher than that for H atoms due to more efficient energy transfer by the higher mass atom to the surrounding water molecules. Arasa et al. (2015) revisited these results to determine whether photodissociation of HDO (into either $\mbox{OH} + \mbox{D or OD} + \mbox{H}$) could affect the deuterium fractionation of water in the gas phase. The conclusion was that the effects are likely minimal; however, the lower probability for D atom photodesorption may lead to fractionation effects within the ice. This latter prediction remains to be tested in astrochemical models. Probabilities for each outcome as a function of monolayer and temperature have now been made available for use in astrochemical models (Arasa et al. 2015). Also available are values averaged over the top few monolayers for use in those models which are unable to treat the ice chemistry in a stochastic manner.

Very recent MD calculations have also attempted to understand the underlying mechanism for CO photodesorption via the DIET process (van Hemert et al. 2015). van Hemert et al. (2015) excited selected CO molecules with $8.7\mbox{--}9.5~\mbox{eV}$ photons contained within both crystalline and amorphous CO clusters (containing up to 1200 CO molecules). They determined two potential mechanisms. In the first, the excited CO molecule decays to the ground state in an unfavorable configuration and subsequently ejected, and the second corresponds to the aforementioned “kick out” mechanism. Photodesorption yields of $(4.0\times10^{-3})~\mbox{photon}^{-1}$ agreed well with the lower end of the range of values currently constrained by experiment.

8.3 Implementation in Astrochemical Models

To date, the implementation of experimentally-determined photodesorption yields ($Y_{\mathrm{X}}$, typically expressed in units of molecules $\mbox{photon}^{-1}$) in astrochemical models has been fairly pragmatic. The photodesorption rate of a molecule is dependent upon the flux of FUV photons striking the grain surface, $F_{\mathrm{UV}}~\mbox{photons}\,\mbox{cm}^{-2}\,\mbox{s}^{-1}$, the cross-sectional area of grain surface available per unit volume, $\sigma_{g} n_{g}~\mbox{cm}^{2}\,\mbox{cm}^{-3}$, and the surface coverage fraction of the specific molecule in the top 1–2 monolayers of the ice mantle, $\theta_{\mathrm{X}}=n_{\mathrm{X}}/n_{\mathrm{S}}$. In the above expressions, $\sigma_{g}$ is the geometrical cross section of a grain surface, $n_{g}$ is the number density of grains, $n_{\mathrm{X}}$ is the number density of molecules of species X in the surface layers of the ice mantle, and $n_{\mathrm{S}}$ is the total number of species in the surface layers of the ice mantle. In two-phase models, there is no distinction between the surface coverage factor and the bulk coverage factor. However, these numbers will likely differ in three-phase and multilayer models.

The photodesorption rate is thus given by,

$$ k_{\mathrm{X}}^{\mathrm{pd}} = \theta_{\mathrm{X}} Y_{\mathrm{X}} F_{\mathrm{UV}} \sigma_{g} n_{g}~\mbox{cm}^{-3}\,\mbox{s}^{-1}. $$

(36)

Derivation of the rate beginning instead with the cross-sectional area of an individual grain surface site leads to an identical expression due to the relation between number density of surface sites, grain surface area, and total number of species per monolayer. Note that the experimentally determined yield implicitly includes the effects of ice opacity to FUV radiation (since the yield is given per incident photon) and the fact that photodesorption only occurs from the top few monolayers of the ice. This results in a zeroth-order rate coefficient: the photodesorption rate is independent of ice thickness (or number density of ice molecules) for thick ices. Within this framework, it is also possible to include branching ratios of the expected products, for those species for which photodissociation is the primary process (e.g., H₂O and CO₂). Several of the experimental groups also provide analytical fits to the photodesorption yields as a function of ice thickness (in monolayers) and temperature (Öberg et al. 2009a,c). This allows the inclusion of a more realistic treatment of the photodesorption rate at low surface coverage (important in models which consider evolution from diffuse to dense conditions).

To show the effect of photodesorption, Fig. 8 again shows a simulation of a laminar solar nebula where we have lowered the photodesorption yield from $10^{-2}$ to $10^{-5}$. Locally the abundance can change by several orders of magnitudes, but integrated over $z$ the changes are small, except for gas phase CO₂ which is strongly enhanced in the outer regions.

Fig. 8

The change in molecular abundance of a selection of species in a laminar solar nebula at 1 Myr, when using a photodesorption yield of $1\times 10^{-5}$ instead of $1\times 10^{-2}$. Relative gas phase abundances are in the left panels, the corresponding ice abundances in the panels on the right-hand side

Now that wavelength-dependent data are available, astrochemical models for which the FUV spectrum is included explicitly (e.g., in protoplanetary disks or PDRs), the total photodesorption yield per species can be found by integrating over the experimentally determined yield and the FUV spectrum,

$$ Y_{\mathrm{X}} = \int_{\lambda}Y_{\mathrm{X}}(\lambda) F_{\mathrm{UV}}( \lambda) d\lambda. $$

(37)

These data are now available for pure CO, N₂, O₂, and CO₂ ices (Fayolle et al. 2011, 2013; Fillion et al. 2014). One caveat of this method is that one needs to take into account that the wavelength dependence of the desorption yield will depend upon the dominant molecular component of the ice mantle in the subsurface layers of the ice (Bertin et al. 2013). This is possible to track in three-phase and multilayer models, but non trivial in two-phase models. This is because additional assumptions are required, for instance, whether the ice is sequestered into layers and if so, with which layers (H₂O-dominated or CO-dominated) more minor species are most likely to be associated. In the case of H₂O-dominated ices, wavelength dependence has not been determined explicitly but is likely to follow the photodissociation (or absorption) cross section of amorphous solid water (see, e.g., Cruz-Diaz et al. 2014a). It is also worth bearing in mind that, in the regime where water ice is the dominant ice component, the ice temperature could be such that thermal desorption will dominate, except for those molecules which are comparable in volatility with water (e.g., CH₃OH).

In current gas-grain chemical models, photodissociation is treated as a parallel process to photodesorption. Because interstellar ices are optically thin to FUV radiation, photodissociation can occur throughout the ice mantle (in contrast to photodesorption which occurs only from the top few monolayers). Photodissociation rates for specific species are computed by integrating over the photodissociation cross section multiplied by the FUV spectrum, exactly as is done in the gas phase (e.g., van Dishoeck et al. 2006),

$$ k_{\mathrm{X}}^{\mathrm{ph}} = \int_{\lambda}\sigma_{\mathrm{X}}(\lambda)F_{\mathrm{UV}}( \lambda) d\lambda, $$

(38)

with products and associated branching ratios usually adopted from the equivalent gas-phase measurements or estimates. Photoabsorption cross sections in the VUV range are now available for several molecules in the solid phase and are usually a good proxy for computing the photodissociation rates (e.g., Cruz-Diaz et al. 2014a,b). For other species, the gas-phase data is typically used and should be accurate to within a factor of a few. The derivation of photodissociation cross sections and branching ratios in the solid phase remains a complex process requiring careful modelling of ice experiments which explore the full range of multi-dimensional parameter space (Öberg et al. 2009c).

By treating photodissociation as a parallel process to photodesorption in full gas-grain models, models can implicitly account for the production of reactive radicals within and on the ice mantle. The models are also able to treat the fate of such radicals including loss to the gas phase via desorption (whether photo driven or thermally driven), or reaction with other radicals and atoms present in the ice, as discussed in the introduction to this Section. However, in models in which chemical desorption is also included, one runs into the issue in which photodesorption is accounted for twice: one explicitly as outlined above, and one implicitly, where radicals formed via photodissociation recombine to reform the molecule which then has an assumed efficiency for chemical (or reactive) desorption. This latter process is essentially direct photodesorption. Given the inherent uncertainties which exist in large astrochemical models due mainly to the lack of available experimental data, the effect of this “double counting” is likely minor. A potential solution to the problem is to treat photodissociation products as a separate (and likely) excited species in the ice mantle which have a timescale for thermalization guided by the MD calculations.

9 Precision Considerations in Rate-Based Modeling

A basic, two-phase model of interstellar gas and grain surface chemistry may be solved easily using the standard rate-equation approach that is commonly used for pure gas-phase chemistry. However, as explained earlier, this method can sometimes lead to unphysically high grain surface production rates, under conditions where two reactants on the surface have mean populations per dust grain that are less than unity. In a mathematically exact treatment, if the reaction between two such species is sufficiently rapid once both particles are present on the grain surface, then the surface populations of those species should become anti-correlated. The rate-equation method cannot reproduce such behavior, as it calculates only the mean, individual populations $\langle n_{\mathrm{s}}(A)\rangle$ and $\langle n_{\mathrm{s}}(B)\rangle$, rather than the mean of the product, $\langle n_{\mathrm{s}}(A) n_{\mathrm{s}}(B)\rangle$.

Garrod (2008) proposed a modification scheme (i.e. a so-called “modified rate” approach) to take account of the inaccuracies of the rate-equation method, which Garrod et al. (2009) showed to be effective over a wide range of physical conditions. However, the full implementation of this method involves the calculation of quantities such as $(1-\langle n_{\mathrm{s}}(A)\rangle)$. In a numerical simulation, the value of this differential becomes highly imprecise as $\langle n_{\mathrm{s}}(A)\rangle$ approaches 1, due to the limited fractional precision available in the calculations (for double precision calculations, the fractional precision is on the order of $10^{-16}$). This imprecision can induce the solver to use very small timesteps, significantly slowing the calculations. For this reason, the maximum available precision should be used for modified-rate calculations wherever possible. Several Fortran compilers now include quadruple precision as a standard feature. Converting a chemical code and differential equation solver to quadruple precision is a straight-forward task.

The use of a three-phase gas-grain model can also introduce precision problems. As explained by Hasegawa and Herbst (1993b), the three-phase treatment calculates a net deposition rate onto the dust grains that represents the difference between the total accretion rate of all gas-phase species onto the grains, and the total desorption rate of all surface species back into the gas (as well as any processes such as reaction that may also decrease the total number of surface particles). This net rate is used to determine the rate at which material in the surface layer is pushed down into the bulk ice (in the case of net deposition) or brought up from the bulk into the surface layer (net desorption), in order to maintain the same total population of particles in the surface layer. The net deposition rate therefore strongly couples the calculated abundances of all surface and mantle species with the abundances of whichever gas-phase or surface species have the greatest influence over that net deposition rate. Improving the overall precision of the calculations (i.e. moving to quadruple precision) allows the net rate to be more precisely determined, improving model run-times.

Unfortunately, in the case of three-phase models, simply increasing the overall precision of the calculations is often not enough to avoid a failure of the solver routine. A more careful consideration of the error associated with individual abundances is also required. We may illustrate this point with an example in which the net deposition rate onto dust grains is solely determined by the difference between (i) the rate of deposition of gas phase H₂ and (ii) the rate of desorption of grain surface H₂. If one rate strongly dominates over the other, then the fractional error on the net deposition rate is mainly dependent upon the solver-enforced fractional error (i.e. “EPS”, error per step) on either the gas-phase or surface H₂ abundance. However, if the two rates are nearly equal, then the error on the net rate will depend upon the error on both abundances. Crucially, in this latter case, if the net rate is very much smaller than either the deposition or desorption rates, then the absolute error on the net rate may become much larger even than the rate itself. Because the calculated gas-phase and grain surface H₂ abundances may fluctuate anywhere within their allowed error limits, this situation can produce wild fluctuations in the net deposition rate. This prompts a shift by the solver routine to smaller and smaller timesteps and ultimately a failure to integrate over manageable timescales.

To avoid this problem, we must enforce an error value that is more stringent than EPS for certain species in the model that affect the net deposition rate. The general EPS value used in a model may be, say, $10^{-4}$. The solver compares this EPS with the fractional error on each species resulting from the time-integration of the differential equations. If the EPS exceeds the computed error value for all chemical species, then convergence is assumed, and the timestep is complete. In order to impose stronger error handling on the net deposition rate, we may bias this comparison by some factor, according to the sensitivity of the net deposition rate to changes in the abundance of a particular species. If a small, fractional change in a chemical abundance induces a larger fractional change in the net deposition rate, then the ratio of these values is used to bias the error comparison to enforce smaller errors. The sensitivity of the net deposition rate to the abundance of each species can be determined at each Jacobian call by the solver. Typically only a few chemical abundances are important to the error dependence of the net deposition rate, such as gas-phase and grain surface H, H₂ and CO.

This more stringent error-handling approach is necessary for the use of three-phase models in all cases where the net deposition rate changes sign during the course of the integration, because such conditions induce extreme fractional errors as the net rate approaches zero.

10 Test Case: Hydrogenation of CO Ice Leading to Complex Molecules

Catastrophic freeze-out of CO molecules from the gas phase leads to the formation of an apolar (CO-rich) ice. The surface hydrogenation of CO in inter- and circumstellar ices drives the formation of more complex species. Here we will use the hydrogenation of CO as an example of how a combined experimental and modeling approach can result in an increased understanding of complex molecule formation. In 1992, Charnley et al. proposed that formaldehyde (H₂CO) and methanol (CH₃OH) form in the solid phase by the hydrogenation of CO ice. Formaldehyde is the first organic polyatomic molecule detected in the interstellar medium (ISM) (Snyder et al. 1969). Methanol has been observed abundantly in the ISM (Belloche et al. 2013) and is considered to be a resource for the formation of more complex molecules through surface reactions and after its sublimation in the gas phase (Balucani et al. 2015). Originally, gas-phase complex molecules were postulated to arise from dissociative recombination of protonated methanol, which itself was formed via radiative association of smaller ions and molecules (e.g., ${\mbox{CH}_{3}^{+} + \mbox{H}_{2}\mbox{O} \longrightarrow \mbox{CH}_{3}\mbox{OH}_{2}^{+} + \mbox{h}\nu}$ followed by ${\mbox{CH}_{3}\mbox{OH}_{2}^{+} + \mbox{e}^{-} \longrightarrow \mbox{CH}_{3}\mbox{OH} + \mbox{H}}$). Pioneering measurements with CRYRING, a low-temperature ion-storage ring, showed that dissociative recombination led predominantly to fragmentation of the protonated species, for instance, for the example listed above, only 3 % of dissociative recombinations where found to lead to methanol (Geppert et al. 2006). On the other hand, the surface formation of methanol has been long under debate, as the first sets of laboratory data gave conflicting results. In their study, Hiraoka et al. (2002) reported only the formation of formaldehyde, whereas Watanabe and Kouchi (2002) confirmed an efficient formation of formaldehyde and methanol ice. This controversy was then solved in 2009 by an independent investigation of the same reaction scheme performed for a wider range of H-atom fluxes and ice temperatures (Fuchs et al. 2009), which showed the different results from the two groups were due to the use of different H-atom fluxes and fluences.

The first ($\mbox{CO}+\mbox{H}$) and third ($\mbox{H}_{2}\mbox{CO}+\mbox{H}$) reaction steps in this H-atom addition scheme presents barriers that have to be overcome, see also Fig. 3. The efficiency of the process therefore depends on the barrier height and how much excess energy is available to cross over. Figure 9 shows the time evolution of the surface abundance (in molecules $\mbox{cm}^{-2}$) of CO, H₂CO and CH₃OH during H-atom bombardment of CO ice with an H-atom flux of $5\times 10^{13}~\mbox{cm}^{-2}\,\mbox{s}^{-1}$ at different temperatures. This figure is taken from Fuchs et al. (2009). At all the temperatures, the amount of CO decreases as the abundance of H₂CO increases. After an H-atom fluence of $1\times 10^{17}~\mbox{H}$ atoms $\mbox{cm}^{-2}$, the formation of methanol starts to increase at the expense of the H₂CO abundance. Similar trends are also reported in Watanabe et al. (2006). For temperatures higher than 15 K, a clear drop in the production rate of methanol is observed by Fuchs et al. (2009). This is mainly due to a lower lifetime of H atoms on the surface of the ice that causes the H-surface abundance to drop substantially and therefore reduce the probability of hydrogenation reactions occurring in the ice at temperatures higher than 15 K.

Fig. 9

Time evolution of the surface abundance (in molecules $\mbox{cm}^{-2}$) of CO, H₂CO and CH₃OH during H-atom bombardment of CO ice with a H-atom flux of $5\times 10^{13}~\mbox{cm}^{-2}\,\mbox{s}^{-1}$ at surface temperatures of 12.0 K (a), 13.5 K (b), 15.0 K (c), and 16.5 K (d). Experimental data (symbols) and Monte Carlo simulation results (solid lines) are shown as well. Figure reproduced from Fuchs et al. (2009). Credit: G.W. Fuchs, H.M. Cuppen, S. Ioppolo, C. Romanzin, S.E. Bisschop, S. Andersson, E.F. van Dishoeck, and H. Linnartz: Astron. Astrophys. 505, 629 (2009), reproduced with permission, copyright ESO

A detailed interpretation of laboratory measurements in terms of understanding processes like deposition of H atoms, diffusion, reaction and desorption is only possible if laboratory work is supported by simulations. Fuchs et al. (2009) used a microscopic kinetic Monte Carlo algorithm to fit effective reaction barriers for $\mbox{CO}+\mbox{H}$ and $\mbox{H}_{2}\mbox{CO}+\mbox{H}$ to a large set of H₂CO and CH₃OH production and CO loss data as a function of H-atom fluence and temperature indicated by the solid lines in Fig. 9. The outcome of the model yielded typical values in the range of 400 to 500 K. In general, simulations and experimental data are in very good agreement at 12 K. The agreement for 13.5, 15, and 16.5 K is far less good, probably because of missing mechanisms that promote the penetration into the ice. The reaction rates for $\mbox{H}+\mbox{CO}$ and $\mbox{H}+\mbox{H}_{2}\mbox{CO}$ have only a minor temperature dependence, which suggests that there is a clear tunneling component for the reaction at low temperature. Values from Fuchs et al. (2009) are in good absolute agreement with the barriers found by Awad et al. (2005), who also found a similar temperature behavior. Their values were obtained using a rate equation analysis for $T$ = 10, 15, and 20 K using the data from Watanabe et al. (2006).

Temperature-dependent values are important - within their experimental accuracies - for inclusion in astrochemical models that predict astronomical column densities. The Fuchs et al. Monte Carlo model is further extended to interstellar conditions to compare with observational H₂CO/CH₃OH data. The H₂CO and CH₃OH interstellar abundances obtained in the simulations do not scale directly with H-atom fluence because of the different relative importance of H₂ production and CO hydrogenation in space compared with the laboratory. However, laboratory experiments are required to derive the necessary rates that serve as input to the Monte Carlo program. The obtained H₂CO/CH₃OH ratios for the interstellar simulations are in closer agreement with observational limits than a direct translation of the experimental results. Cuppen et al. (2009) extended the Fuchs et al. (2009) model to simulate in detail the formation of both H₂CO and CH₃OH under interstellar conditions. Unlike the master equation (Stantcheva et al. 2002), rate equation (Ruffle and Herbst 2000), and macroscopic Monte Carlo (Ruffle and Herbst 2000) studies of methanol formation, the continuous-time, random-walk Monte Carlo technique accounts for the layering of the CO. This layering is crucial, because the constant addition of fresh CO tends to protect the underlying layers from hydrogenation and limits the time available for reaction.

The link between solid CO and CH₃OH was also recently confirmed observationally. Using VLT data on the CO band of a large collection of YSOs (Pontoppidan et al. 2003), peak width, peak position and upper limits on the $2152~\mbox{cm}^{-1}$ band were extracted and compared to laboratory spectra. This proved that the astronomical spectra are fully consistent with CO ice intimately mixed with CH₃OH, in agreement with a common chemical history as described above (Fraser et al. 2004; Cuppen et al. 2011). A more recent study expands this proposition by comparing literature laboratory spectra of different ice mixtures containing CO to high-resolution Keck/NIRSPEC spectra, of the massive YSO AFGL 7009S and of the low mass YSO L 1489 IRS (Penteado et al. 2015). Accurate band profile comparisons of both the ${}^{12}\mbox{CO}$ and ${}^{13}\mbox{CO}$ ice bands were carried out. Moreover, Penteado et al. gave evidence that a gradient of different ratios of CO:CH₃OH mixtures might exist in the grain mantles and that the spectrum of high-mass YSOs requires a higher content of methanol than the low-mass YSOs to reproduce the observed CO band profile.

In the past decade, to fully characterize the methanol formation reaction network, deuteration experiments were also performed on CO ices (Hidaka et al. 2007), which resulted in both deuterated formaldehyde and methanol. In this study, it was found that deuteration rate of CO is 12.5 times lower than the rate of CO hydrogenation under the same experimental conditions. A similar CO deuteration rate was also recently confirmed by Fedoseev et al. (2015). In addition to $\mbox{CO}+\mbox{D}$ experiments, deuteration studies of H₂CO and CH₃OH ices have also been performed (Nagaoka et al. 2005; Hidaka et al. 2009; Watanabe and Kouchi 2008). Figure 8 of Hidaka et al. (2009) presents the complete surface reaction network as obtained when CO is exposed to H and D atoms based on several studies (Watanabe et al. 2006; Watanabe and Kouchi 2008; Hidaka et al. 2009). Reaction rate constants are shown relative to those of $\mbox{CO}+\mbox{H}$ at 15 K (Hidaka et al. 2004). The general conclusion here is that although D-atom additions to CO/H₂CO proceed slower than H-atom additions to CO/D₂CO, D-atom induced substitution reactions in H₂CO and CH₂OH take place faster than the corresponding H-atom induced substitutions in D₂CO and CD₂OD.

Recently, Fedoseev et al. (2015) showed the first laboratory evidence that surface hydrogenation of CO ice leads to the formation of molecules with more than one carbon atom along with the formation of H₂CO and CH₃OH. In their experiments, Fedoseev et al. confirmed the surface formation of two larger molecules of astrobiological importance, i.e., glycolaldehyde (HC(O)CH₂OH) and ethylene glycol (H₂C(OH)CH₂OH). The key step in this process is believed to be the recombination of two HCO radicals followed by the formation of a C-C bond as also indicated by modelling work (Woods et al. 2013). The dimerization of HCO is considered by Woods et al. (2013) as one of the formation routes of glycolaldehyde. Their astrophysical model matches the observed estimates in the hot molecular core $\mbox{G}31.41+0.31$ and low-mass binary protostar IRAS 1629-32422. Fedoseev et al. (2015) experimentally determined that the recombination reaction between two HCO radicals yields glyoxal (HC(O)CHO) at 13 K; sequential hydrogenation by 2 or 4 H-atoms turns HC(O)CHO into glycoladehyde and ethylene glycol, respectively. This is an important experimental finding, as so far surface hydrogenation reactions were mainly shown to be effective in the formation of relatively simple interstellar relevant species (e.g., H₂O, NH₃, and CH₃OH), while radical-radical recombination reactions were considered efficient in the ISM only at temperatures higher than 20–30 K. The latter assumption comes from the fact that radicals are commonly assumed to become mobile around 30 K (e.g., Garrod 2013; Butscher et al. 2015). Therefore, diffusion of radicals is assumed to be the limiting factor to the formation of COMs through radical-radical recombination reactions. Fedoseev et al. (2015) gave the first laboratory evidence that diffusion is not necessarily needed for radicals to find each other. In CO-rich ices, HCO radicals can be formed next to one another through the hydrogenation of CO molecules and react without diffusing in the ice.

To further strengthen the conclusions of Fedoseev et al. (2015), experimental results are implemented into the model used by Cuppen et al. (2009) based on the continuous-time, random walk Monte-Carlo method described before. This model allows for the simulation of microscopic grain-surface chemistry for the long timescales typical in interstellar space, including the layering of ice during the CO freeze out. The choice of microscopic simulations is because HCO radicals have to stay in close vicinity to each other in order for a reaction to occur. Therefore, a method that accounts for the position of the species in the ice lattice is required for the best representation of the system. In the simulations reported by Fedoseev et al. (2015), five new species, including glycolaldehyde and ethylene glycol, are introduced in the simulations (i.e., HC(O)CHO, HC(O)CH₂O/HC(O)CHOH, HC(O)CH₂OH, H₂C(O)CH₂OH/HC(OH)CH₂OH, and H₂C(OH)CH₂OH). The glycolaldehyde abundance was found to correlate with the formaldehyde abundance, while the abundance of ethylene glycol correlates with the abundance of methanol. This is not surprising since both CH₃OH and H₂C(OH)CH₂OH are hydrogen saturated species while both H₂CO and HC(O)CH₂OH are not. Finally, the correlation found between the abundances of HC(O)CH₂OH and H₂C(OH)CH₂O with those of formaldehyde and methanol is as expected typically in ratios of the order of a few percent and in good agreement with observations (Jørgensen et al. 2012; Schöier et al. 2002).

More recently, Chuang et al. (2016) presented new laboratory experiments on the low-temperature solid-state formation of three complex molecules, i.e., methyl formate (HC(O)OCH3), glycolaldehyde and ethylene glycol, through recombination of free radicals formed via H-atom addition and abstraction reactions that occur during the hydrogenation of CO ice at 15 K. The experiments include the co-deposition of H atoms and pure H₂CO, H atoms and pure CH₃OH as well as H atoms and mixtures of H₂CO with CO, H₂CO with CH₃OH, and CH₃OH with CO, in fact extending previous CO hydrogenation studies (e.g., Fuchs et al. 2009; Fedoseev et al. 2015). The aim of these experiments is to resemble the physical-chemical conditions typical of the CO freeze-out stage in dark molecular clouds, when H₂CO and CH₃OH form by the recombination of accreting CO molecules and H atoms on ice grains. The occurrence of H-atom abstraction reactions in ice mantles leads to the formation of a higher number of intermediates (HCO, CH₃O, and CH₂OH) than previously thought, when assuming sequential H-atom addition reactions only. This enhances the probability to form COMs through radical-radical recombination without the need of UV photolysis or cosmic rays as external triggers.

Figure 10 shows the resulting COMs formation network based on the investigated reaction routes. From top to bottom, a sequence of H-atom addition and abstraction tunneling reactions is shown to lead to the formation of CH₃OH. The figure shows that H₂CO can undergo an abstraction reaction induced by H atoms to form HCO radicals, which can be further dehydrogenated to form CO. This pathway increases the total number of HCO formation events and its lifetime in the ice mantle. H₂CO also participates in addition reactions with H-atoms. Formation of CH₃O radicals is confirmed by the detection of solid methyl formate, a species formed most likely under the detection limit in Fedoseev et al. (2015). However, formation of CH₂OH radicals is not fully excluded by Chuang et al. (2016). In contrast, H-atom induced abstraction reaction involving CH₃OH likely yield CH₂OH while no proof for CH₃O formation is found by the authors. The barrier-less recombination of HCO intermediates yielding glyoxal followed by consequent hydrogenation (i.e., the mechanism investigated in Fedoseev et al. 2015) yields HC(O)CH₂OH and H₂C(OH)CH₂OH and is presented in the right panel of the reaction scheme in Fig. 10. Alternatively, the intermediate HCO radicals can directly recombine with CH₃O or CH₂OH to form HC(O)OCH₃ and HC(O)CH₂OH, respectively. According to the authors, the recombination of two CH₂OH radicals (dash arrow) seems to contribute less to the formation of H₂C(OH)CH₂OH.

Fig. 10

Extended COM formation network as obtained from the CO, H₂CO, and CH₃OH hydrogenation experiments by Fedoseev et al. (2015) and Chuang et al. (2016). Solid arrows indicate the reaction less efficient pathways. Figure reproduced from Chuang et al. (2016)

Chuang et al. (2016) showed that both H-atom addition and abstraction reactions play an important role in effectively increasing the number of interaction events between reactants as well as the lifetime of radicals in the ice. This way, more COMs are formed through recombination of reactive intermediates than previously assumed at 10 K under interstellar conditions. It is important to note that Fuchs et al. (2009) did not include hydrogen abstraction reactions when calculating reaction rates. Clearly, the final products they observed in their experiments were the result of an equilibrium between addition and abstraction reactions. Therefore, the rates inferred by Fuchs et al. (2009) should be regarded as effective rates. The recent laboratory and modeling advances in our understanding of interstellar surface reactions discussed in this section give a strong indication that the COMs observed in cold dense molecular clouds have a solid phase non-energetic origin and are most likely linked to a combination of hydrogenation addition/abstraction reactions and radical-radical recombination reactions at 10 K. CO-rich ices formed upon the catastrophic freeze-out of CO molecules onto interstellar grains seem to be the right environment for these surface reaction to occur efficiently in the ISM.

11 Outlook

In this review, we have summarized the “state-of-the-art” in terms of the methods used and the data available for treating grain-surface chemistry under astrophysical conditions. Since the first conceptual models of interstellar grain-surface chemistry, now more than 40 years old (e.g., Watson and Salpeter 1972a; Allen and Robinson 1977; Tielens and Hagen 1982), the general structure of a gas-grain code has not significantly changed. What has come on leaps and bounds is the advances in experimental techniques that have enabled the first quantification of many of the parameters required for models including binding energies, reaction barriers, diffusion rates, and desorption yields. Such experiments remain challenging to perform, and as reaction systems become more complex, the various potential chemical pathways and branching ratios remain difficult to disentangle. As discussed here, using sophisticated stochastic astrochemical models to also simulate laboratory data can help extract this information from the experiments which is then input into the larger rate-equation astrochemical models for simulating astrophysical environments.

Despite vast advances in the quantification of interstellar grain-surface chemistry in the past decades, there remains much more to learn. Typical grain-surface networks consist of sometimes hundreds of grain-surface species, the vast majority of which are reactive radicals. Such species are transient by nature and remain difficult to isolate and thus study in the laboratory. We urge the experimental community to continue to push the boundaries of what is possible in the laboratory and develop innovative methods to quantify the necessary data for radicals. In astrochemical models used to simulate interstellar and circumstellar environments, the data adopted for radicals remains the most uncertain.

Thanks to the effort over the past decades, we have a reasonably clear picture of the different processes important for ice chemistry, summarized in each of the Sections in this review. However, one chemical process we have only touched upon is ion-molecule chemistry in ice mantles (e.g., Woon 2011), the gas-phase equivalent of which is arguably the most important for building molecular complexity in the gas phase in space.

A second process we have not covered in detail here, but debate remains on its relative importance, is the direct impact of cosmic-ray particles on ice-covered dust grains. Due the high energies attained by cosmic-ray particles (∼ Mev–GeV), such a process is thought to have catastrophic effect on the composition of the ice due to the creation of an ionized particle track through the grain which can also destroy, ionize, and sublime molecules in an explosive manner from the ice mantle (e.g., D’Hendecourt et al. 1982; Leger et al. 1985). Such a process is also stochastic by nature, and remains challenging to account for in a “mean-field” approach, since each impact can have a different track on each individual grain, and not every dust grain will be impacted by a cosmic-ray particle. This process had fallen out of favor as it became increasingly evident that photodesorption was an efficient (and well understood) mechanism for releasing molecules from the ice mantle at low temperatures. Recent revisitations of this process suggest that the overall effect is of the same order of magnitude as other processes considered to influence the gas-grain balance in dark clouds such as photodesorption (e.g., Ivlev et al. 2015).

It is clear that much work remains in linking the observed gas-phase composition in dark clouds with the known ice composition, and astrochemical models remain critical for constraining the potential processes for determining the gas-grain balance. The outcome of this Lorentz Centre Workshop shows what can be achieved when scientists from various fields and disciplines (modelers, experimentalists, and observers) meet to intensively discuss a naturally interdisciplinary field, such as astrochemistry. We have presented the first detailed review of the data and methods used for interstellar/circumstellar grain-surface chemistry and made recommendations, based on a legacy of past and present experience and research, for the treatment of the reactions rates. We have also compiled in a systematic manner the ever-increasing pool of experimental and theoretical data currently available in the literature. As we move into a new era of observational astronomy, with ALMA and JWST both expected to have a transformational impact on the field of astrochemistry, we anticipate regular future meetings on this topic to review, compile, and update the current “state-of-the-art” for the treatment of grain-surface chemistry in astrochemical models.