© ESO 2020

1 Introduction

About 200 molecules have been identified in interstellar space¹. Approximately 50 of them are classified as complex organic molecules (COMs), that is, they are comprised of six or more atoms. COMs are generally observed in the gas phase, although it is widely accepted that saturated COMs are formed primarily through reactions on the surface of icy dust grains. This is where gas-phase species accrete, meet, and react to form saturated molecules. In the present paper we focus on the formation of glycolaldehyde and ethylene glycol. Glycoaldehyde and ethylene glycol have been detected towards the class 0 protostellar binary IRAS 1629-32422 (Jørgensen et al. 2012) and towards Galactic center source Sagittarius B2 (Hollis et al. 2000, 2002). These COMs are potential prebiotic molecules because they are precursor molecules for sugars, and they serve as an excellent test case for more complex grain surface chemistry networks because their formation involves the creation of a C−C bond.

The precise mechanism of COM formation is unknown, but is subject of many studies. Several different mechanisms have been suggested over the past decades. Charnley (1997) and Charnley & Rodgers (2005, 2009) applied a simple atom addition mechanism where species such as methanol can react with atomic carbon and are subsequently hydrogenated. Another mechanism was proposed by Garrod & Herbst (2006) and is more efficient at elevated temperatures. Cosmic-ray-induced photons can dissociate H₂CO and CH₃OH, creating functional-group radicals such as CH₃ and CH₃O on or within the ice mantle. At low temperatures, these radicals are hydrogenated again, but as the temperature in the core increases to above 20–30 K, the residence time of H atoms on the surface decreases substantially while at the same time radicals become mobile, thereby allowing radical-radical association reactions to become competitive with hydrogenation of radicals. Later work indeed showed that the formation of COMs, including glycolaldehyde and methyl formate, is possible at temperatures above 20 K (Garrod et al. 2008; Vasyunin & Herbst 2013). Additionally, Laas et al. (2011) showed that the branching ratio of radical formation affects the abundance of the recombined species as well.

The detection of acetaldehyde CH₃ CHO, dimethyl ether CH₃ OCH₃, methyl formate CH₃OCHO, and ketene CH₂CO in the prestellar core L1689B suggests that COM synthesis has already started at the prestellar stage and a viable grain surface reaction route exists that is possible at low temperatures and does not require extensive UV irradiation (Bacmann et al. 2012). The scenario described above requires both and therefore cannot be the only mechanism leading to the formation of COMs. Gas-phase routes have been proposed (Skouteris et al. 2018) as well as grain surface formation through cosmic-ray irradiation (Shingledecker & Herbst 2018; Shingledecker et al. 2019). Woods et al. (2013) proposed a mechanism for the formation ofCOMs, where the essential carbon-backbone elongation occurs by recombination of two formyl radicals (HCO). Their astrophysical model matches the observed estimates in the hot molecular core G31.41+0.31 and low-mass binary protostar IRAS 1629-32422. Recently, Fedoseev et al. (2015) confirmed this picture by showing that CO-hydrogenation experiments yield glycoaldehyde and ethylene glycol at 13 K. The authors argued that the formation of the C−C bond must have occurred through a recombination of two HCO radicals although the product glyoxal (HC(O)CHO) was not detected because other proposed mechanisms should have an activation barrier and are thus not feasible at low temperatures. This is an important experimental finding because surface hydrogenation reactions have so far mainly been 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 thought to be only efficient at temperatures higher than 20–30 K, when radicals are commonly assumed to become mobile (e.g., Garrod 2013; Butscher et al. 2015).

To further strengthen the conclusions of Fedoseev et al. (2015), experimental results were implemented into a kinetic Monte Carlo model (Cuppen et al. 2009) that was initially built to study CO hydrogenation up to the formation of methanol. This microscopic model showed that in CO-rich ices, HCO radicals can form next to one another through the hydrogenation of CO molecules and react without diffusing in the ice, even for the low H-atom fluxes typical in interstellar space. The experimental findings were recently complemented by laboratory results on the hydrogen abstraction reactions of H₂ CO and CH₃ OH (Chuang et al. 2016). Abstraction reactions increase the abundance of radicals in the ice with respect to a purely sequential H-atom addition scheme. This enhances the probability that COMs form through radical-radical recombination without the need of UV photolysis or cosmic rays as external triggers.

We here focus on the latter route: surface reactions through radical-radical recombination at low temperatures and without UV irradiation. We apply a microscopic kinetic Monte Carlo (KMC) model to simulate COM formation at low temperatures for a range of different physical conditions. This will give us constraints on the number of COMs that can be formed through cold grain surface chemistry. Gas-phase routes or cosmic-ray-induced chemistry might contribute as well, but these are beyond the scope of this paper. This work aims to aid in the interpretation of COM observations by giving insight into how physical conditions affect the obtained abundance ratios. The main strength of microscopic simulation methods is that the location of individual species is traced throughout the simulation. An earlier study of hydrogenation of CO (Cuppen et al. 2009) showed that this plays an important role and that a gradient mantle builds up during catastrophic CO freeze-out. This is in agreement with ice observations of CO (Penteado et al. 2015). In the past years, several codes have been developed that tackle this problem by introducing multiple phases representing different layers within the ice (Taquet et al. 2012; Vasyunin & Herbst 2013; Furuya et al. 2016). However, the species are still treated in a “mean-field” fashion within a layer, by either rate equations or a macroscopic KMC routine. Crucial to the low-temperature formation mechanism is that two HCO radicals form in close vicinity of each other and can react without the need for diffusion. Obtaining a realistic estimate for the occurrence of this reaction therefore requires knowledge of the exact position of the species with respect to one another. A microscopic KMC model has access to this information and can therefore be used to estimated the relative importance of the different COM formation channels. Moreover, it includes the surface structure, and as our study shows, this has an important effect on the reactivity.

An important starting point for any chemical model is the reaction network, and the accuracy of the simulated data directly depends on the quality of this network. Tunneling reactions and radical–radical reactions play a crucial role in the CO hydrogenation network. In the past years, many tunneling rate constants have been accurately determined, which makes this study very timely. Unfortunately, branching ratios on radical-radical reactions are still lacking. We therefore start this work by determining them.

2 Method

2.1 Electronic structure calculations to determine the branching ratio

To obtain branching ratios for radical–radical reactions considered in the reaction network, an approach similar to that used in Lamberts (2018) was employed. In particular, we took the approach to run geometry optimizations for structures where the two radicals are placed in a random orientation with respect to each other. That is, the optimized structure of a radical in the gas phase is placed in the origin and the optimized structure of the second radical is automatically placed next to it at a random angle. This was done 60 times in total, with the C-atom of the second radical being placed symmetrically at the vertices of a truncated icosahedron at a distance of approximately 2.5–4.5 Å. The initial geometries were then relaxed with an L-BFGS algorithm (Nocedal 1980; Liu & Nocedal 1989) and the optimized geometry indicates the product channel to which a barrierless reaction has lead. The branching ratio of a particular product channel is here defined as its total occurrence divided by 60. Unrestricted density functional theory (DFT) was used with the PW6B95 functional (Zhao & Truhlar 2005) with additional D3 dispersion correction (Grimme et al. 2010) and the def2-TZVP basis set (Weigend et al. 1998). The ChemShell framework was used for all calculations (Sherwood et al. 2003; Metz et al. 2014), using Turbomole (Ahlrichs et al. 1989; Treutler & Ahlrichs 1995) and DL-Find (Kästner et al. 2009) for the geometry optimizations. For radical–radical reactions it is crucial to use an unrestricted broken-symmetry wave function with overall as many spin-up as spin-down electrons, but finite spin density on both radicals (e.g., Lamberts 2018; Enrique-Romero et al. 2019). No Hessian calculations were performed. To test for the influence of the choice of the functional and basis set on the branching ratios, several full geometry optimization runs were repeated using PBEh-3c/mSVP, which is computationally less demanding but includes dispersion correction as well. The major and minor product channels were found to be the same, although the exact values differ, giving credibility to the results. However, we realize that the DFT approach is approximate, and more accurate calculations could be done, particularly with a multi-reference CASSCF approach.

Although surface reactions on interstellar ices are of interest here, for the current study, a model was employed in which the surface molecules were not explicitly taken into account. The relatively weak interaction between the CO-rich surrounding allow for the assumption that the surface influence on the potential energy landscape of barrierless reactions is small. The D3 dispersion correction can also be seen as mimicking some attractive surface interaction, however. The orientation of the two radicals with respect to each other is important for the outcome (Lamberts et al. 2019), however. We therefore intentionally used 60 randomly generated starting geometries where the two radicals were placed in close proximity to each other, mimicking a situation in which they are nearest neighbors. With this approach the restricted rotation of a radical because of interactions with hydrogen-bonded CH₃ OH, for example, cannot be taken into account. Sampling a variety of binding sites on realistic dirty ices and determining branching ratios per binding site is beyond the scope of this work, however.

The initial separation between the two radicals plays a crucial role in whether or not the calculations can reach SCF convergence. When they are too far apart, the interaction will be too small for the radicals to react, whereas when they are too close together, repulsive forces become important that are not representative of the situation on interstellar grains either. We performed some tests for the distance dependence and found that with increasing distance, (i) fewer configurations lead to a reaction and more stable nonreactive geometries are calculated instead, and (ii) the branching ratios leading to complex organic molecules (C–C or C–O bonds) are lower because COM formation is more sterically hindered in general. For the current study, we mimicked the situation that two radicals are on the surface in close proximity (nearest neighbors), as is the case for the KMC simulations as well, with a closest atom distance of 2.5–4.5 Å between the two radicals at hand. This separation depends on the size of the reactants themselves and is comparable to the site-site distance in the KMC model.

2.2 Kinetic Monte Carlo simulations

Ice evolution was simulated using a continuous-time, random-walk (CTRW) algorithm, which is a type of KMC method. For an in-depth theoretical description of this technique, see Cuppen et al. (2013). The program used in this study is similar to the one used in Cuppen & Herbst (2007) and Cuppen et al. (2009). Briefly, the KMC technique simulates the evolution of a grain surface. This grain surface is represented by a 50 × 50 grid with periodic boundary conditions and is bare at the start of the simulation. A random-number generator in combination with the rates determines the sequence of processes. Similar to Cuppen et al. (2009), the types of processes include deposition onto the surface, reaction of two species, hopping between sites, and desorption from the surface. These processes have been extended to also include photodissociation and cosmic-ray dissociation processes (Cuppen & Herbst 2007). It is assumed that the rate coefficients for hopping and desorption are thermally activated by $$R_{\text{X}}=\nu \exp \left(-\frac{E_{\text{X}}}{T}\right),$$(1)

where ν, the preexponential factor, has been approximated by kT∕h ≈ 2 × 10¹¹ s⁻¹ and E_X is the energy barrier of process X in K. The rate coefficients for the photodissociation and cosmic-ray dissociation processes are determined according to $$R_{\text{photo}}={\alpha }_{\text{photo}}\exp \left(-{\gamma }_{\text{photo}}{A}_{\text{V}}\right)$$(2)

and $$R_{\text{CR}}={\alpha }_{\text{CR}}\zeta ,$$(3)

respectively. Here α is the photon rate, γ is the extinction coefficient, A_V is the visual extinction, and ζ = 1.3 × 10⁻¹⁰ yr⁻¹ is the cosmic-ray ionization rate. Standard values from the UMIST RATE2012 database were taken (McElroy et al. 2013).

The program was used to simulate laboratory conditions and interstellar conditions. The former was used as a benchmark, the latter to obtain insight into the conditions under which complex organic species are formed. A full overview of all performed simulations can be found in Table 1.

3 Results and discussion

3.1 Reaction network

The reaction network used in this study is given in Table 2. The same color-coding scheme for the main species is used in all figures throughout the paper. The basis of this reaction network is a hydrogen addition chain from carbon monoxide (yellow) through formalaldehyde (orange) to methanol (red) (Watanabe & Kouchi 2002; Fuchs et al. 2009). When the lifetime of formyl (HCO) radicals is sufficient, two of these radicals could recombine to form glyoxal (HC(O)CHO, GX) (Fedoseev et al. 2015; Chuang et al. 2016). Through continued hydrogenation, glyoxal will consecutively be turned into glycolaldehyde (HC(O)CH₂OH, GA, purple) and ethylene glycol (H₂ C(OH)CH₂OH, EG, green). Finally, recombination of a formyl and a methoxy radical leads to the formation of methyl formate (HCOOCH₃, MF, blue).

Two types of reactions are present in the reaction network. First, recombination reactions involving two radical species. These reactions are assumed to be barrierless and are included with a reaction rate k = kT∕h ≈ 2 × 10¹¹ s⁻¹ at these low temperatures. Second, hydrogenation reactions that possess barriers. The reaction pathways of such reactions have been calculated by Andersson et al. (2011), Song & Kästner (2017), and Álvarez-Barcia et al. (2018) using instanton theory, and they were optimized through a quasi-Newton-Rhapson method. A resulting instanton is the most probable tunneling path for a given reaction. The corresponding tunneling rate can be calculated using harmonic quantum transition state theory (HQTST). The reaction rate is governed by tunneling below the critical temperature T_c, which is $$T_{\text{c}}=\frac{\hslash \Omega }{2\pi k_{\text{B}}},$$(4)

where Ω is the imaginary normal-mode frequency of the saddle point. The crossover temperature is equal to 79 K for a typical reaction H + CO → HCO with a frequency of 2171 cm⁻¹ (Andersson et al. 2011).

An important addition to the reaction network with respect to previous versions (Cuppen et al. 2009; Fedoseev et al. 2015) is the inclusion of abstraction reactions such as H₂ CO + H → HCO + H₂. Penteado et al. (2017) have shown that competition between this abstraction and the two hydrogenation reactions toward CH₃O and CH₂OH has a definitive influence on the yields of the species. Chuang et al. (2016) have already shown that hydrogen-addition to formaldehyde yields significantly less CH₂OH than either CH₃O or HCO, which has been confirmed by the calculations of Song & Kästner (2017).

3.1.1 Branching ratios of radical–radical reactions

The radical–radical branching ratios were estimated using the method described in Sect. 2.1. The branching ratios obtained here for a consistent set of reactions are the first to be reported in the literature and serve as a guide for the major and minor product channels. Product channels with branching ratios below 0.05 (<5%) are listed only in Table 2, but are not included in the KMC model. We did not consider any reactions between radicals of the sort C₂ H_nO₂ and HCO, CH₂OH, or CH₃O. Although we are aware that these reactions can take place in a barrierless fashion as well (see also Butscher et al. 2017), they are beyond the scope of this study.

The branching ratios in Table 2 do not add up to one for all reactions. The reactants HCO + HCO are an example. This is because certain starting geometries did not lead to reaction even though the reaction being barrierless in principle. This is most likely because their relative geometry does not lead to a suitable reaction complex that leads to reaction, and hence no reaction occurs. This effect could also occur on a grain surface at low temperatures where two species meet in the wrong geometry and cannot reorient. The two species then remain frozen in close proximity. The KMC program was changed to allow for nonreactive branches to accommodate for this effect. This was realized by first choosing which reaction was added to the table of events in the KMC routine. This choice was based on the branching ratios, and in principle, also no reaction can be added, in which case the species were only allowed to diffuse or desorb. Simply lowering the overall reaction rate will not help to account for this because it is a stochastic process that cannot be repaired in a mean-field way, that is, smearing out the lowered reactions. For HCO recombination, for instance, there is a chance of 0.40 that the reaction 2HCO → HCOCHO with rate 2 × 10¹¹ s⁻¹ is added tothe table of events, a chance of 0.25 for 2HCO → H₂ CO + CO with rate 2 × 10¹¹ s⁻¹, and 0.35 that no reactions are added. We expect this effect to become increasingly important for larger reactants. We show below that this inclusion does not significantly change the results.

3.1.2 Activated reactions

Álvarez-Barcia et al. (2018) calculated rate constants for GX + H leading to H₂ COCHO and HCOHCHO. They were unable to treat the abstraction reaction because their quantum chemical method could not accurately describe this particular reaction. Gas-phase experiments at room temperature show, however, that this is an important reaction that leads to HCOCO, which can spontaneously decay to CO + HCO (Colberg & Friedrichs 2006; Orlando & Tyndall 2001). This two-step reaction is included in the network as the main destruction route of the C–C bond. We assumed the rate constant to be comparable to its hydrogenation. Room-temperature experimental data show a rate constant of roughly three orders of magnitude higher than abstraction of H₂ CO. For the unimolecular dissociation, we used a rate constant of 10⁴ s⁻¹, which is a factor of 600 lower than the rate constant at 250 K (Orlando & Tyndall 2001). With this rate, the entire HCOCO decays to CO and HCO. This is rather insensitive to the exact value.

The only radical-neutral reaction with a radical other than the H atom that we considered is H₂ CO + CH₃O → CH₃OH + HCO. Although Butscher et al. (2017) indicated that the reaction between the HCO radical and H₂ CO might also take place, in fact Álvarez-Barcia et al. (2018) showed that reactions where heavy atoms play a role in the bond formation are not significantly enhanced by tunneling. In other words, the rate constants for the reactions H₂ CO + CH₃O → CH₃OCH₂O and H₂CO + HCO → (HCO)CH₂O are below 10⁻⁹ s⁻¹ at low temperatures. Therefore, other reactions of the same type, such as the formation of a C–C or C–O bond through radical-neutral reactions, were discarded as well. The reaction HCO + CO → CO + HCO may lead to transport of H atoms in the CO bulk. This reaction is expected to have a high barrier, however, because a C–H bond needs to be broken and the reaction H + CO → HCO already has a high barrier itself.

In principle, branching ratios for different channels in activated reactions can be determined by simply equating them to the relative rate constants for the individual reactions as determined quantum chemically. However, this is an oversimplification and might lead to ignoring important product channels. For the reaction GA + H, this can result in four product channels of which hydrogen addition to CH₂OHCH₂O (k₁ = 2.8 × 10⁵ s⁻¹) and abstraction leading to HCOHCHO (k₂ = 6.8 × 10⁷ s⁻¹) are the major outcomes. The ratio between the rate constants k₂ ∕k₁ suggests that abstraction is 200 times more likely than addition. On a grain surface this is not the case, however. When a hydrogen atom approaches GA, the possible product channel will depend on the direction of approach, whether this reaction will indeed occur depends on the residence time of the H atom at close vicinity and on possible roaming of the H atom during this time. Considering diffusion rates of hydrogen atoms, which are much slower than both reaction rates, we expect both reactions to occur when meeting in the correct geometry for that particular reaction. The branching ratios of both channels will be much more comparable than based on their relative rates. We therefore chose to treat reactions with multiple channels in a similar way as the nonactivated reactions: a reaction is chosen to be added to the table of events with more or less equal probability (see branching ratio column in Table 2). This reaction is added with its corresponding rate. If GA and H are in adjacent sites, there is a probability of 0.40 that GA + H → CH₂OHCH₂O with rate k₁ is added to the table of events and a probability of 0.40 for GA + H → HCOHCHO + H₂ with rate k₂. This leads to a roughly equal occurrence of both reactions.

3.2 Experimental simulations

The reaction network was benchmarked against the experimental results as published by Chuang et al. (2016, 2017). Table 3 shows an overview of their results as they are relevant to our study. All three experiments are codeposition experiments of either CO and H or H₂ CO and H over 6 h. The surface abundances of CO, H₂CO, and CH₃OH were followed by reflection absorption infrared spectroscopy (RAIRS) in situ during the experiments, whereas relative abundances of MF, GA, and EG were only obtained after heating the sample at the end of the experiment. The abundances of CO, H₂ CO, and CH₃OH in Table 3 are shown in monolayers (MLs). For the COMs with C–C bond, the final ratios for MF:GA:EG are given. In the benchmark study, the simulations ran for an equivalent time of 6 h or 21 600 s, at which time a steady state was reached in most cases, and the overall surface coverage was then compared against the experimental results. We used an effective H-atom flux of f_effective = f_exp × S, where S is the sticking fraction of H on a surface, which we assumed to be 0.20 (Buch & Czerminski 1991; Al-Halabi et al. 2002, 2003, 2004; Batista et al. 2005; Veeraghattam et al. 2014) for a room-temperature atom landing on a cold surface. Reducing the sticking fraction from unity to the more realistic value of 0.2 leads to a slightly better agreement.

As mentioned earlier, the majority of our reaction rate constants are based on Álvarez-Barcia et al. (2018). One of the most critical reactions to obtain agreement with the experiments is H + H₂ CO. This reaction with three possible product channels has been studied quantum-chemically by Song & Kästner (2017). They studied this reaction both in the gas phase and on a water-ice (ASW) surface. The ASW changes the transition state for these reactions through hydrogen bonding, which results in significantly higher rates than in the gas phase. Because we work on a CO surface without hydrogen-bonding capabilities, we expect the current system to be best represented by gas-phase conditions. When we use these rates directly, H₂ CO is predominantly destroyed to form HCO and H₂. Both methanol and formaldehyde are highly underproduced. The actual rate constants calculated by Song & Kästner (2017) are based on DFT calculations. The applied functional is first compared against stationary calculations at CCSD(T)-F12 level, which is more accurate, but also much more computationally demanding. This comparison shows that the reaction barriers are underestimated, and the reaction rate constants are accordingly overestimated. This effect is stronger for the abstraction reaction than for the hydrogenation reaction leading to CH₃O. Decreasing the rates and using equal branching ratios leads to much better results that are within the uncertainty of the method.

Figure 1 shows the simulated build-up during codeposition of H and CO or H₂ CO. The H input flux of simulation 1 (Fig. 1a) is equal to 1.6 × 10¹² cm⁻² and the flux of the codeposited species equals 6.5 × 10¹¹ cm⁻² s⁻¹. For simulations 2 and 3 (Fig. 1b,c), these fluxes are equal to 1.4 × 10¹³ and 7.0 × 10¹¹ cm⁻² s⁻¹, respectively. In all simulations, a molecular hydrogen flux was added equal to that of the atomic hydrogen flux to account for recombination in the cooling nozzle of the particle beam. Along the y-axis, the graphs show the deposition of several species relative to the CO fluence, for example, the total number of monolayers of deposited CO. The colored translucent bands indicate the results obtained by Chuang et al. (2016), maintaining an uncertainty factor of 1.5. The three different rows correspond to three different sets of H-atom binding energy. We used a fixed diffusion barrier for hydrogen on a flat CO surface of 200 K. Diffusion barriers are site-specific and depend on the binding energy before and after the hop, ΔE: $E_{\text{diff}}=E_{\text{diff}}^{\text{flat}}+\max (\Delta E,0)$ K. The value of 200 K is in line with experimental results by Kimura et al. (2018) of H diffusion on CO ice and agrees well with calculations of long-range H diffusion on ASW by Ásgeirsson et al. (2017) and Senevirathne et al. (2017). The binding energy of H to a CO surface is less well constrained, however. To test the sensitivity of the results to the H-binding energy, we varied this. Panels I correspond to a binding energy of 250 K (flat surface), panels II to 420 K, and panels III to 670 K. Small differences can be observed between the various plots, but the final abundance is rather independent of the choice of binding energy. Experiment 2 is quite well reproduced by the model, whereas the model slightly overproduces CH₃OH for experiment 3 and slightly underproduces it for experiment 1. Even when we varied the rates, we were not able to arrive at a model that could reproduce all experimental results simultaneously. In the remainder of the paper, the intermediate value of E_bind(H₂) = 420 K is used.

Table 4 summarizes the abundances at the end of simulations of Fig. 1 II and can be directly compared to Table 3. It is quite a triumph for the computational chemistry that the model that uses purely calculated rate constants and branching ratios can reproduce the experimental results this well. The main discrepancy is the formation of MF. The formation trends are the same between experiment and simulation, but MF is much more abundant in the simulations than in the experiments. In our model, MF is formed through HCO + CH₃O → HCOOCH₃, and once it is formed, there is no destruction reaction. The overabundance in our model is hence either due to an overestimation of the branching ratio of HCO + CH₃O → HCOOCH₃ or due to missing reactions. From quantum chemistry calculations, we know that abstraction reactions such as HCOOCH₃ + H → COOCH₃ + H₂ have very low rates (Álvarez-Barcia et al. 2018). Other destruction reactions are not obvious based on the species available in our system.

The removal of nonreactive branches has only a slight effect on the final results, as can be seen from Fig. 2 and comparing it to Fig. 1 II for the same conditions. When all radical-radical encounters lead to reaction, slightly more glycolaldehyde is formed at the expense of CH₃OH, but other species remain rather unaffected. Whether two radicals can indeed remain frozen next to each other in the ice depends on their reorientation timescale, which is unknown. The present system is a CO-rich ice where hydrogen bonding plays a minor role and species might reorient. On a water-ice surface, we expect species with hydrogen-bonding capabilities to be more easily frozen into a specific geometry, and nonreactive branches can then become important in the grain surface chemistry, especially for larger species. How nonreactive branches should be considered is a topic for further investigation. To present a conservativelower limit of the COM abundances, we included nonreactive branches in our simulations of ISM conditions.

Fig. 1 Simulated codeposition of CO+H (a and b) and H2 CO+H (c) or simulations 1-3 of Table 1. The three different rows correspond to different choices of H-atom binding energy (see main text). I: 250 K, II: 420 K, and III: 670 K.

Fig. 2 Same as for Fig. 1II, but without nonreactive branches.

3.3 Comparison with UMIST rates

The UMIST RATE2012 database² (McElroy et al. 2013) is generally used as a source of rate constants for astrochemical models. Our new network, which is based on a consistent set of calculated tunneling rates and is benchmarked against experiments, is a substantial improvement over the UMIST values for this network. The large number of missing reactions in the latter network is immediately apparent when the new rates are compared with the UMIST rates in Table 2. This involves mostly abstraction reactions toward less complex species and the HCO + HCO reactions and the subsequent hydrogenation reactions toward GA and EG. The second point that stands out is that the rates for H + CO and H + H₂ CO are two orders of magnitude higher in the new network. The low-rate hydrogenation reactions in the UMIST database are hence an effective rate for the full hydrogenation chain. We would also like to point out the reaction constant of hydrogenation of MF, which is 36 orders of magnitude lower according to Álvarez-Barcia et al. (2018) than the constant included in UMIST. Rate constants in UMIST are typically calculated based on a barrier height, assumed barrier width, and an estimated effective mass using a crude formula for tunneling. The effective mass estimation completely fails for this reaction because it does not only involve movement of H, but also breaking a C–O bond.

We performed the same set of three simulations using rates from the UMIST network (Table 2). The results of these simulations are shown in Fig. 3. Whereas the agreement for the experiment a is quite reasonable, methanol is heavily overproduced for b and c. The reason for this is that the UMIST network uses a lower effective rate instead of the full hydrogenation or abstraction network. As long as sufficient H atoms are available, CO is fully hydrogenated toward CH₃OH. When methanol is formed, no destruction pathway leadis away from this. The dependence of the simulated coverages on the chosen sticking fraction of the H atoms to the surface is also much larger than for the new network. The need of an abstraction route is further apparent in the H₂ CO-hydrogenation simulation (panel c), which shows no CO at all. Only MF and GA are formed of the two-carbon COMs. Their abundance is very similar, whereas experiments show that MF is much less produced than GA, and their ratio also differs between the different experimental conditions.

Based on these simulations, we therefore expect the UMIST network to predict a too high methanol and MF abundance in dark clouds and a too low abundance of GA and EG. This would automatically favor the photodissociation heat-up mechanism by Garrod & Herbst (2006).

Fig. 3 Same as Fig. 1 II using UMIST as a source for reaction rate constants instead of Álvarez-Barcia et al. (2018) (Table 2).

3.4 Interstellar conditions

The second part of this study aims to use the benchmarked network to simulate the evolution of a grain-surface mantle in dark molecular cloud conditions. In the ISM simulations, the flux of CO and H onto the grain surface depends on the gas-phase abundance of both species. For H atoms, we adopted a constant gas-phase abundance throughout the simulation. Because the CO abundance gradually decreases as more CO is depleted, the CO gas-phase abundance was recalculated during the simulation, taking the current surface CO and its derivatives into account. Because molecular hydrogen is nearly inert in our reaction network (see Sect. 3.1), we kept it at a constant low abundance of 10 cm⁻³ for all simulations to save computational time. A subset of the simulations was repeated with a higher H₂ abundance of 100 cm⁻³, and we found no significant difference.

Figure 4 shows the evolution of a grain surface for four different hydrogen abundances in dark molecular cloud conditions. These simulated abundances are n(H) = 1.0, 2.5, 4.0 and 6.0 cm⁻³ in panels a through d, respectively, which span an extreme range in accordance with Goldsmith & Li (2005), who found n(H) ≈ 2–6 cm⁻³. The initial CO abundance is n_initial(CO) = 10 cm⁻³. This relates to a density of n_H = 10⁵ cm⁻³ using the typical relative abundance of 10⁻⁴ (Bergin & Tafalla 2007), and the grain-surface temperature is equal to 10 K in all cases. The left-hand vertical axis (red) denotes the grain-surface build-up in number of particles per surface area, whereas the right-hand axis (blue) denotes the grain-surface abundance with respect to n_H. The total simulated time is 2 × 10⁵ yr for all ISM simulations (Table 1).

Carbon monoxide freeze-out is faster than hydrogenation at the initial stage, as shown in Fig. 4. Some of this carbon monoxide is hydrogenated into formyl radicals (HCO) and subsequently into formaldehyde (H₂ CO). CO remains dominant for the simulated time in all four simulations. A small decrease in mantle build-up can be seen toward the end of the simulations; this is a result of CO depletion from the gas phase. In all simulations hydrogenation and covering of the surface by incoming CO molecules compete.

Panel a of Fig. 4 is the result of a simulation with a low H:CO ratio where the covering by CO is fast with respect to hydrogenation. The mantle is mainly composed of CO, with H₂ CO being the second most abundant species. A cross section of the formed grain mantle can be found in Fig. 5, which shows the formation of compact CO ice with some minority species, predominantly H₂ CO, and covered with a layer of H₂. In these conditions the formation of HCO is relatively rare, and two HCO radicals are hardly formed in close proximity. The low amount of glycolaldehyde that forms indicates that few HCO recombinations have occurred. Enough hydrogen is present to further efficiently hydrogenate HCO to H₂CO. In the latter stages of the simulation, the efficiency of glycolaldehyde formation is sufficient to form ethylene glycol because the CO deposition rate decreases as a result of progressed depletion in the gas phase.

Results of simulations with a slightly higher gas-phase hydrogen are shown in Fig. 4b and in the upper right panel of Fig. 5. The grain mantle now features less CO in favor of H₂ CO and more complex species. COMs are present on the grain early in the simulation. Because the H:CO ratio is higher than in the previously described simulation, the formation of HCO has increased and alongside the probability that two formyl radicals are formed in close proximity. The formation of glycolaldehyde and ethylene glycol, which are the result of subsequent hydrogenation, competes with direct hydrogenation of HCO toward methanol. These effects increase when the abundance of gas-phase hydrogen is raised to n(H) = 4.0 cm⁻³ (Fig. 4c).

At the highest simulated hydrogen abundance n(H) = 6.0 cm⁻³, the amount of methanol equals the amount of formaldehyde, and EG is more abundant than GA. Moreover, methyl formate formation generally increases with hydrogen abundance. This is to be expected because it is dependent on both HCO and CH₃O, whose formation also increases with the presence of hydrogen.

To study the effect of temperature, the simulations were repeated for 8, 12, 14, 16, 18, and 20 K. Figure 6 shows the mantle composition for a selection of these simulations: n(H) = 2.5 cm⁻³ and n_initial(CO) = 10 cm⁻³ and four different temperatures: (a) 8 K, (b) 12 K, (c) 16 K, and (d) 20 K. The surface abundance of the main species is rather independent of the surface temperature in the range 8–16 K. Because we took thegas and grain temperature to be equal, the deposition rate increased with temperature ($\propto \sqrt{T}$). This has some effect on the H:CO ratio at late times because CO depletes faster at 16 K than at 8 K. H₂ CO/CH₃OH and GA/EG decrease with temperature as a consequence. Figure 7 shows the mantle cross sections corresponding to the results in Fig. 6. The build-up of ice at 20 K is initially significantlyreduced due to desorption of CO at these elevated temperatures, but the formation of CH₃OH helps to retain some CO at later times. Comparison between 8, 10 (Figs. 6), 12, and 16 K shows that the surface structure of the ice is smoother for higher temperatures because of the increased mobility of CO. Where the surface is covered with H₂ for 8 to 12 K, this has all disappeared at 16 K. A full overview of all the obtained results of dark cloud simulations can be found in Table 5.

Fig. 4 Simulated time evolution of the grain surface composition for four different levels of hydrogen abundances: 1.0, 2.5, 4.0, and 6.0 cm−3 (a–d respectively). The abundance of gas-phase carbon monoxide is 10.0 cm−3. The temperature of the grain is 10 K. The relative surface abundance is given with respect to nH.

Fig. 5 Cross sections of the formed ice layer for four different levels of hydrogen abundances: 1.0, 2.5, 4.0, and 6.0 cm−3. The abundance of gas-phase carbon monoxide is 10.0 cm−3. The temperature of the grain is 10 K.

Fig. 6 Simulated time evolution of the grain surface composition for four different temperatures: 8, 16, 18, and 20 K (a–d resp.). The abundance of gas-phase carbon monoxide is 10.0 cm−3 and of hydrogen atoms is 2.5 cm−3.

Fig. 7 Simulated mantle cross sections for four different temperatures: 8, 12, 16, and 20 K. The abundance of gas-phase carbon monoxide is 10.0 cm−3 and of hydrogen atoms is 2.5 cm−3.

3.5 Formation routes

Figure 8 shows the flux distribution of the network for a standard simulation of n(CO) = 10.0 cm⁻³, n(H) = 2.5 cm⁻³, and T = 10 K. The thickness of the arrows indicates the relative frequency of the reaction. The figure clearly shows that HCO constantly reacts and reforms. In this process of hydrogenation and abstraction, H₂ is formed, which becomes the dominant route to form H₂. This holds for the ISM and experimental simulations. Hydrogen atoms react away before a new hydrogen atom can land on the surface, and the lifetime of a single H atom on the surface is too short to form H₂ through H + H → H₂. Hydrogenation reactions are hence the main loss route for H atoms. This also explains the weak dependence of the results in Fig. 1 on the binding energy of H to the surface.

Furthermore, the recycling of EG is immediately obvious. The rate constant for the abstraction reaction of EG is high, much higher than for GA and CH₃OH. The H₂COHCHOH radical that is formed in this way has its unpaired electron on a carbon. Reaction of this species with HCO could further elongate toward glycerol.

The formation routes of methanol and glycoaldehyde both require at least four hydrogenation steps. Methanol is the more abundantly formed species of the two for all conditions. The majority of methanol is formed from the reaction between CH₃O and H₂ CO, as we show in Fig. 8. This is further quantified in Table 6 for the ISM and experimental conditions. The relative contributions of the different formation routes of H₂ CO and CH₃OH was found to be relatively insensitive to the temperature, but dependent on the H:CO ratio. The table therefore only shows the contribution per simulated H:CO ratio for 10 K. In all cases, CH₃O + H₂ CO → CH₃OH + HCO is the dominant formation route for methanol and not CH₃O + H → CH₃OH, which might be intuitively expected to be the dominant formation route. Penteado et al. (2017) have shown for a full gas-grain model with diffusive surface reactions that methanol is predominantly formed from destruction reactions of more complex species rather than through hydrogenation of CH₃O. The reaction CH₃O + H₂CO → CH₃OH + HCO was underestimated in this study because radical–radical recombination without diffusion was not taken into account. In both cases, the reason is the low surface coverage by H atoms. In the present case, H atoms mostly react with CO, HCO, H₂ CO, and EG and its radical to produce H₂. If two formaldehyde molecules are formed close to each other, only one hydrogen atom is needed for hydrogenation to CH₃O to convert this into CH₃OH and HCO. The latter will be quickly hydrogenated again. CH₃O radicals that are not formed in close proximity to a H₂CO molecule will be quickly hydrogenated by one of the H₂ molecules that covers the surface. Only at higher temperatures, when the surface coverage of H₂ is low, does the reaction CH₃O + H → CH₃OH start to play a role in methanol formation. The other product branch CH₃O + H → H₂ CO + H₂ leading to formaldehyde is similarly only important for elevated temperatures. H₂ CO is predominantly formed through hydrogenation of HCO, with HCO + HCO → H₂ CO + CO of increasing importance with lower H:CO ratio.

Fig. 8 Flux distribution of the network for a standard simulation of n(CO) = 10.0 cm−3, n(H) = 2.5 cm−3 and T = 10 K. The thickness of the arrows indicates the relative frequency of the reaction. Some arrows become too thin and are indicated in gray. Reactions that contribute less than 0.1% with respect to H + CO → HCO are dashed.

3.6 Comparison with observations

To place our simulations into astrophysical context, we would like to compare our results to astrophysical observations. Unfortunately, this cannot be done directly to ice observations. Complex species such as EG, GA, and MF are hard to detect with IR observations because of sensitivity and selectivity issues. They have been observed in the gas phase toward various astrophysical objects. Table 7 shows a selection of observed abundances of H₂ CO, GA, EG, and MF with respectto CH₃OH. For cometary observations, this corresponds to the coma observations where mantle material desorbs by sunlight. For hot corinos and hot cores, the temperature of the grains in the object is sufficiently high to have lost all its mantle material. Finally, for the cold and quiescent sources, the ice mantles are mostly intact and the observed gas-phase species have desorbed from the grain through some nonthermal desorption mechanism such as photodesorption, reactive desorption, and/or cosmic-ray desorption (either through heating or sputtering). In all cases, we assume that the gas-phase composition reflects the composition of the original ice mantle. This only holds if gas-phase chemistry has not sufficiently altered their abundance. This is especially relevant for H₂CO, which also has gas-phase formation routes. In cometary coma, gas-phase chemistry is unlikely to occur. Because for cold cloud observationsnot all ice has desorbed, we need to take into account that the nonthermal desorption efficiencies are not the same for all species. The observed composition is therefore not a direct reflection of the mantle composition. The observations will give us some guidance on the accuracy of our model, however, and might indicate the physical conditions under which the ice mantle has been formed.

The H₂ CO/CH₃OH ratio in Table 7 shows large variation, but remains below unity for all cases. Comparison with Table 5 shows that this means that the H:CO ratio should be at least 0.3. Quiescent clouds toward the Galactic center show an even lower H₂ CO/CH₃OH. This is hardto reconcile directly with our results, but might be caused by a higher reactive desorption efficiency for CH₃OH than for H₂CO. Methanol is predominantly formed in a reaction with two products, whereas an addition reaction is the most dominant route for H₂ CO formation. Bond-breaking reactions are more likely to lead to translational excitation than to an addition reaction. Translational excitation can result in very efficient desorption (Fredon et al. 2017; Fredon & Cuppen 2018). Missing CH₃OH formation routes could be another reason, for instance, CH₃O + HCOOH → CH₃OH + COOH.

The GA/EG ratio also varies strongly. In our simulations, H₂ CO/CH₃OH and GA/EG ratios are directly correlated because both CH₃OH and EG formation at the expense of H₂CO and GA, respectively, requires a high H abundance. In the observational data, no clear correlation between these two ratios can be observed. The reason may be the limited data set with an inhomogeneity of sources. The simulated EG and GA abundance is higher than the observed abundance. This could be due to missing EG destruction routes, as discussed in the previous section, or to missing CH₃OH formation routes through CH₃ + OH.

Methylformate has abundances of a few percent with respect to methanol for all sources except for the hot core source, where MF is more abundant than CH₃OH. The low abundances are consistent with our simulations, where a rather constant MF abundance with respect to CH₃OH can be observed. The high abundance of MF observed in the hot core is probably due to a different formation mechanisms. Our simulations are representative for a cold and dark grain chemistry, whereas during high-mass star formation, the cloud is probably warmer and more exposed to UV light. This leads to more radicals and a higher mobility of these radicals. This favors radical-radical recombination to form MF, which is more in line with the simulations by Garrod & Herbst (2006). The MF/CH₃OH might serve as a probe for cold dark grain chemistry.

4 Conclusions

In summary, the present paper presented Monte Carlo simulations to constrain methyl formate, glycolaldehyde, and ethylene glycol formation on the surface of interstellar dust grains in cold and dark molecular clouds. For this purpose, a new CO-hydrogenation reaction network was presented. Using this network, we arrived at the following conclusions.

• 1.

  The rate constants in this network are based on literature values. The constants were calculated using harmonic quantum transition state theory. This takes hydrogen tunneling correctly into account. The branching ratios for radical–radical reactions in this new network were determined quantum-chemically as well. We present a large consistent set of these values.

• 2.

  This new network was benchmarked against experimental results of Chuang et al. (2016, 2017) on CO and H₂CO hydrogenation. A good agreement is reached with this new model.

• 3.

  A comparison between experiments and simulations using rate constants and branching ratios from UMIST showed an overproduction of CH₃OH, which stresses the importance of considering both hydrogenation and abstraction reactions instead of effective hydrogenation rates.

• 4.

  The new model was used to simulate CO freeze-out in dense molecular clouds. A grid of simulations was performed with varying dust temperatures and H-to-CO ratios. The final COM abundances depend more on the H/CO and less on temperature. Only above 14 K, where CO build-up is less efficient, does temperature start to play a role.

• 5.

  Production of GA and EG through HCO recombination occurs at all H/CO and temperatures we studied, even at temperatures as low as 8 K. The reason is that the HCO + HCO reaction can occur when HCO radicals are formed close to each other and does not require diffusion of the radicals. Relatively low abundances of MF have been formed.

• 6.

  Molecular hydrogen is predominantly formed through abstraction reactions on the surface. The most important reaction leading to methanol is H₂CO + CH₃O → HCO + CH₃OH.

• 7.

  Our simulations agree reasonably well with observed COM ratios for mantles that have been formed at low temperatures.

Acknowledgements

M.S. is grateful for financial support from the Dutch Astrochemistry Network II (DANII) financed by the Dutch Organisation for Scientific Research (NWO). T.L. is grateful for support from NWO via a VENI fellowship (722.017.008).

References

All Tables

All Figures

	Fig. 1 Simulated codeposition of CO+H (a and b) and H2 CO+H (c) or simulations 1-3 of Table 1. The three different rows correspond to different choices of H-atom binding energy (see main text). I: 250 K, II: 420 K, and III: 670 K.
In the text

	Fig. 2 Same as for Fig. 1II, but without nonreactive branches.
In the text

	Fig. 3 Same as Fig. 1 II using UMIST as a source for reaction rate constants instead of Álvarez-Barcia et al. (2018) (Table 2).
In the text

	Fig. 4 Simulated time evolution of the grain surface composition for four different levels of hydrogen abundances: 1.0, 2.5, 4.0, and 6.0 cm−3 (a–d respectively). The abundance of gas-phase carbon monoxide is 10.0 cm−3. The temperature of the grain is 10 K. The relative surface abundance is given with respect to nH.
In the text

	Fig. 5 Cross sections of the formed ice layer for four different levels of hydrogen abundances: 1.0, 2.5, 4.0, and 6.0 cm−3. The abundance of gas-phase carbon monoxide is 10.0 cm−3. The temperature of the grain is 10 K.
In the text

	Fig. 6 Simulated time evolution of the grain surface composition for four different temperatures: 8, 16, 18, and 20 K (a–d resp.). The abundance of gas-phase carbon monoxide is 10.0 cm−3 and of hydrogen atoms is 2.5 cm−3.
In the text

	Fig. 7 Simulated mantle cross sections for four different temperatures: 8, 12, 16, and 20 K. The abundance of gas-phase carbon monoxide is 10.0 cm−3 and of hydrogen atoms is 2.5 cm−3.
In the text

	Fig. 8 Flux distribution of the network for a standard simulation of n(CO) = 10.0 cm−3, n(H) = 2.5 cm−3 and T = 10 K. The thickness of the arrows indicates the relative frequency of the reaction. Some arrows become too thin and are indicated in gray. Reactions that contribute less than 0.1% with respect to H + CO → HCO are dashed.
In the text