© J. Enrique-Romero et al. 2021

1 Introduction

Interstellar dust grains are known to be an essential component of the interstellar medium (ISM) for a large variety of reasons. Among them, dust grains provide the surfaces for chemical reactions that are difficult (or impossible) to take place in the gas phase. An emblematic example is the formation of the most abundant molecule in the Universe, H₂, which largely occurs on grain surfaces (e.g., Hollenbach & Salpeter 1970; Vidali 2013; Wakelam et al. 2017). Other important examples are the formation of water (e.g., Dulieu et al. 2010; He & Vidali 2014; Lamberts & Kästner 2017; Molpeceres et al. 2018) and methanol (e.g., Tielens & Hagen 1982; Watanabe & Kouchi 2002; Rimola et al. 2014), which are also abundant molecules predominantly synthesized on grain surfaces. As a matter of fact, in cold regions, the refractory core of the grains, made up of silicate or carbonaceous material, are coated by icy mantles mostly formed by amorphous water ice synthesized on these surfaces (e.g., Boogert et al. 2015).

It has also been suggested that more complex molecules, the so-called interstellar complex organic molecules (hereinafter iCOMs: Ceccarelli et al. 2017), composed of at least six atoms and containing at least one heavy element other than C, can also be a grain-surface chemistry product (e.g., Garrod & Herbst 2006; Herbst & van Dishoeck 2009; Garrod et al. 2009; Ruaud et al. 2015; Aikawa et al. 2020; Barger & Garrod 2020). One crucial step of this theory is the formation of iCOMs from the combination of two radicals when they meet on the grain icy surfaces. In the majority of the current astrochemical models, the reaction is assumed to proceed barrierless and without competitive channels.

However, previous theoretical works have shown that this is not necessarily the case (Rimola et al. 2018; Enrique-Romero et al. 2019, 2020). For example, theoretical calculations on the energetics showed that the formation of acetaldehyde (CH₃ CHO) on the icy surfaces via HCO + CH₃ is in competition with the formation of CO + CH₄ via direct H-abstraction. In addition, both reactions present barriers, caused by the orientation of the species on the ices, which are governed by the interactions created between the surface water molecules and the two radicals. Indeed, the height of the barriers depends on the site where the reaction occurs, whether the two radicals are on a “plain” ice surface or in a “cavity”, namely on the interactions between the radicals and the ice water molecules.

In the present work, we pursue the above theoretical studies and present new computations to evaluate the efficiency of the radical–radical combination and H-abstraction reactions as a function of the temperature of the icy surfaces with the goal to provide values that can be easily incorporated in astrochemical models. In particular, here we focus on acetaldehyde (CH₃ CHO), one of the most abundant and common iCOMs (e.g., Blake et al. 1987; Cazaux et al. 2003; Vastel et al. 2014; Lefloch et al. 2017; Sakai et al. 2018; Bianchi et al. 2019; Csengeri et al. 2019; Lee et al. 2019; Scibelli & Shirley 2020), as a study case.

Following the works of Enrique-Romero et al. (2019, 2020), we here consider the two competing reactions that can arise from the reactivity between HCO and CH₃: $$\begin{array}{cc}(1)& \text{HCO}\text{+}{\text{CH}}_{\text{3}}\to {\text{CH}}_{\text{3}}\text{HCO}\\ (2)& \text{HCO}\text{+}{\text{CH}}_{\text{3}}\to {\text{CH}}_{\text{4}}\text{+}\text{CO}.\end{array}$$

Our goal is to evaluate the efficiency of each of the two reactions occurring on the interstellar icy surfaces. To this end, we computed the kinetics of the two reactions, using the previous energetic calculations by Enrique-Romero et al. (2019) as a base and the models of amorphous solid water (ASW) for the ice described in Rimola et al. (2018) and Enrique-Romero et al. (2019). In our calculations, we assume that the two radicals are in the most stable energetic configuration prior to reaction, an assumption motivated by the long surface residence timescale (1–10 Myr, the molecular life timescale: e.g. Chevance et al. 2020) that radicals would experience before they become mobile and react with each other.

The article is organized as follows. The definition of the reaction efficiency and our choices for the various assumptions entering in the computations are discussed in Sect. 2. Section 3 describes the adopted methodology. The results are reported in Sect. 4 and we discuss the implications of our new calculations in Sect. 5.

2 Efficiency of radical–radical reaction products on icy surfaces

2.1 Surface-reaction rate definition

Generally, astrochemical models solve the time-dependent equations of the species densities by computing the formation and destruction rates of each species at a given time, both for species in the gas and on the grain surfaces. In particular, the rate R_ij of the formation reaction from two reactant species i and j is expressed as R_ij = k_ijn_in_j, where n_i and n_j are the densities of species i and j, and k_ij is the rate constant at a given temperature. For surface reactions the latter is given by (Hasegawa et al. 1992): $$k_{ij}={\epsilon }_{ij}\times \frac{R_{\text{diff},i}+R_{\text{diff},j}}{n_{\text{d}}},$$(1)

where ε_ij is an efficiency factor which accounts for chemical barriers, n_d is the dust grain density and R_diff,i and R_diff,j are the diffusion rates for species i and j, respectively. These diffusion rates are defined as 1∕t_diff,k, where t_diff,k is the time it takes the species k to scan the whole grain (e.g., Garrod & Herbst 2006). Thus, the sum R_diff,i + R_diff,j gives the rate at which species i and j meet on the surface.

Regarding the efficiency factor ε_ij, different approaches exist in order to derive it. Hasegawa et al. (1992) set it to either 1, in barrierless reactions, or to the tunnelling probability, if the reaction has an activation energy barrier and one of the reactants is light enough to tunnel through it. Later models include also the thermal probability for reaction, if there is an activation energy barrier (e.g., Garrod & Herbst 2006). However, in the presence of an activation energy barrier, reactants need to be close to each other for a certain amount of time for the reaction to occur (Tielens & Hagen 1982). In order to take this into account, Chang et al. (2007) redefined the efficiency taking into account the competition between diffusion and desorption of the most mobile species, as follows: $${\epsilon }_{ij}=\frac{k_{\text{aeb}}(ij)}{k_{\text{aeb}}(ij)+k_{\text{diff}}(i)+k_{\text{des}}(i)},$$(2)

where k_aeb(ij) is the rate constant accounting for the reaction activation energy barrier, which is described by either classical thermal kinetics or quantum tunnelling; (i.e., a frequency times a Boltzmann factor or the tunnelling probability); k_diff (i) is the rate constant for the diffusion of the most mobile species and k_des(i) is its desorption rate constant. Garrod & Pauly (2011) further modified Eq. (2) by removing the desorption term and adding the diffusion of the other reaction partner, j, in the denominator.

For reactions involving radicals, ε_ij is normally assumed equal to 1 (e.g., Garrod & Herbst 2006), as they are considered to react via barrierless exothermic channels.

In this work, we include diffusion and desorption rates of the two reactants, which takes into account both the Chang et al. (2007) and Garrod & Pauly (2011) recipes: $${\epsilon }_{ij}=\frac{k_{\text{aeb}}(ij)}{k_{\text{aeb}}(ij)+k_{\text{diff}}(i)+k_{\text{des}}(i)+k_{\text{diff}}(j)+k_{\text{des}}(j)}.$$(3)

In practice, the efficiency for the reaction is equal to unity only when the time scale for the reaction to occur (1∕k_aeb) is shorter than the timescales at which reactants remain on the reaction site (the smallest between 1∕k_diff(i) and 1∕k_diff(j)).

2.2 A novel treatment of surface radical–radical reactions rate constants

The novelty of the present work is the estimate of the k_aeb(ij) coefficient of radical–radical reactions via statistical kinetics calculations based on the Ramsperger–Rice–Kassel–Marcus (RRKM) microcanonical transition state theory. Briefly, RRKM computations provide unimolecular rate constants, namely the rate at which a system A becomes A′ passing through a transition state (only once). In our case, the system A is the ice-water molecules plus the two adsorbed radicals, namely we consider the water-cluster plus the radicals as a super-molecule isolated from its surrounding. The system A′ is the product of the radical–radical reaction on the icy surface, namely the water-cluster plus either the radical–radical recombination (e.g., React. I) or the H-abstraction (e.g., React. II) products.

It is important to note that in order to apply the RRKM theory, we implicitly assume that the intra-molecular energy redistribution of the reaction energy is faster than the reaction itself. This assumption is supported by recent ab initio molecular dynamics (AIMD) computations that show that a large fraction (≥ 50%) of the reaction energy is absorbed by the water ice in less than 1 ps (Pantaleone et al. 2020, 2021). We have checked a posteriori that the timescale of the reactions studied here is indeed longer that 1 ps. Finally, the specific computational details of our proposed RRKM method are reported in Sect. 3.2.

2.3 Desorption and diffusion energies

Equation (3) shows that, in addition to the probability k_aeb (ij) for radicals i and j to react when they meet on a surface site, the efficiency factor ε_ij also depends on $k_{\text{des}}^{-1}$ and on $k_{\text{diff}}^{-1}$, which are related to the residence time of the radicals on the surface and on the diffusion timescale of the radicals on the ice, respectively. The diffusion and desorption timescales t_diff/des are given by the classical Eyring transition state theory (TST), in which t_diff/des is inversely proportional to their rate constant, t_diff/des $\propto k_{\text{diff/des}}^{-1}$. According to TST, the general expression of the rate constant k for a unimolecular reaction (like diffusion and desorption) is: $$k=\frac{k_{\text{B}}T}{h}\frac{Q^{\ne }}{Q_{\text{R}}}\exp (-\Delta V^{\ne }/k_{\text{B}}T),$$(4)

where ΔV^≠ is the zero point energy-corrected energy barrier, Q^≠ and Q_R are the total partition functions of the transition state and the initial state (namely, the reactants), respectively, k_B is Boltzmann’s constant, T is the surface temperature and h is Planck’s constant. We note that we use the classical Eyring equation because, since we are not dealing with light atoms but molecular radicals, tunneling is negligible.

By proper manipulation of Eq. (4), the rate constant becomes expressed as a function of the free energy barrier ΔG^≠ (usually referred to as free energy of activation) at a given temperature: $$k=\frac{k_{\text{B}}T}{h}\exp (-\Delta G^{\ne }/k_{\text{B}}T),$$(5)

in which ΔG^≠ = ΔH^≠ − TΔS^≠ and where ΔH^≠ is the enthalpy of activation and ΔS^≠ the entropy of activation. These terms contain translational, rotational, vibrational and electronic contributions as they arise partly from the total partition functions Q.

With the adopted quantum chemical approach, the application of the Eyring TST allows us to compute desorption-related data (e.g., desorption activation energies and desorption rate constants) for each radical through the outcome of these calculations (electronic energies, vibrational frequencies, partition functions, energy contributions, etc.). It is worth mentioning that, since the radicals are physisorbed on the ice surfaces, the energy barriers of the desorption processes coincide with the desorption energies. In the present case, we only account for the electronic and vibrational contributions to both ΔH^≠ and ΔS^≠ to arrive at the radical desorption energies as follows: $$\Delta H=\Delta E_{\text{electronic}}+\Delta ZPE+\Delta E_{\text{vib}}(T)+\Delta H_{\text{rot}}+\Delta H_{\text{trans}},$$(6)

and $$\Delta S=\Delta S_{\text{vib}}+\Delta S_{\text{rot}}+\Delta S_{\text{trans}},$$(7)

where the terms are the energy difference between the desorbed and the adsorbed states for the total electronic energy (ΔE_electronic), for the zero point vibrational energy corrections (ΔZPE), for the thermal vibrational energy corrections (ΔE_vib(T)), for the vibrational entropy (ΔS_vib), and for the rotational and translational contributions to enthalpy (ΔH_rot and ΔH_trans, respectively) and entropy (ΔS_rot and ΔS_trans, respectively). In this case, since we are dealing with the desorption of the radicals, the translational and rotational contributions arise from only the desorbed (free) radicals. Specific details on the calculation of some of these terms are provided in Appendix B. For the sake of simplicity, we refer to this final desorption energy as E_des.

In contrast to desorption, obtaining diffusion-related data with the present calculations is a daunting task, as it requires localizing a large number of transition states for the radical hoping between the different binding sites. Moreover, the use of a relatively small cluster model dramatically constraints the validity of these results because of its limitation in terms of size and surface morphology. Therefore, to obtain a value for the diffusion energy of each radicals, which by analogy we will refer to as E_diff, we resorted to what is usually done in astrochemical modeling, that is, E_diff is taken to be a fraction f of E_des. However, deriving the value of f has proven to be difficult, both theoretically and experimentally. In the published astrochemical models, one can find a quite wide range of adopted f values, from 0.3 to ~0.8 (e.g., Hasegawa et al. 1992; Ruffle & Herbst 2000). Some authors have taken a middle point by setting this ratio to 0.5 (e.g., Garrod & Herbst 2006; Garrod et al. 2008; Garrod & Pauly 2011; Ruaud et al. 2015; Vasyunin et al. 2017; Jensen et al. 2021).

In the past few years, theoretical and experimental works on the diffusion process of species on ASW surfaces have provided constraints to the f value (see also the more extensive discussion in Sect. 5.3). In a theoretical work, Karssemeijer & Cuppen (2014) showed that the range for the E_diff /E_des ratio can be narrowed down to 0.3–0.4 for molecules like CO and CO₂. Minissale et al. (2016) experimentally found that the f ratio of atomic species like N and O is about 0.55, while He et al. (2018) showed that f is 0.3–0.6, being the lower values more suitable for surface coverage lower than one mono-layer.

Given the uncertainty on the E_diff /E_des ratio for CH₃ and HCO, we carried out our calculations for three values: 0.3, 0.4 and 0.5.

2.4 Ice model

Regarding the amorphous solid water (ASW) model, there is still little knowledge that constrains the actual internal structure of interstellar ices. Observations suggest that the interstellar water ice is predominantly in the amorphous form (e.g., Smith et al. 1989; Boogert et al. 2015) (with some exceptions: e.g. Molinari et al. 1999). Many laboratory studies havebeen carried out to characterize the possible porosity of the interstellar ices. Typically, laboratory experiments produce porous ices of different densities by condensation of water vapor, even though they probably do not reproduce the interstellar water ice, in which water is believed to form in situ by hydrogenation reactions of frozen O, O₂ and O₃ (e.g., Dulieu et al. 2010; Hama & Watanabe 2013; He & Vidali 2014; Potapov & McCoustra 2021). In general, porous ices are detected in laboratory via the infrared (IR) signature of dangling OH groups, which are, however, missing in interstellar samples (Bar-Nun et al. 1987; Keane et al. 2001), (see also the discussion in e.g. Hama & Watanabe 2013; Zamirri et al. 2018). Several hypothesis have been suggested to explain the absence of the OH dangling signature (Oba et al. 2009; Palumbo 2006; Palumbo et al. 2010), so that, at the end, there is consensus in the community that interstellar water ices are amorphous and porous in nature, even though many details are missing and we do not have a precise picture of the degree of porosity (e.g., Hama & Watanabe 2013; Isokoski et al. 2014; Potapov & McCoustra 2021).

In order to simulate the interstellar icy surfaces, Enrique-Romero et al. (2019) considered a cluster of 33 water molecules. This ice model possesses two major types of surface with respect to the binding capability: a cavity, where species are in general more strongly bonded to the surface, and an elongated side (Rimola et al. 2014, 2018; Enrique-Romero et al. 2019). In this work, we only report the analysis of the reaction occurring in the cavity for the following reason. In astrochemical models, the vast majority of radical–radical reactions take place inside the bulk of the ice (e.g., Garrod & Herbst 2006). Therefore, the cavity site is a better representation of the sites where radical–radical reactions occur than that on the elongated side, which would at best describe the ice layer exposed to the gas and where just a tiny fraction of the reactions can occur, considering that the ice is constituted by more than 100 layers (e.g., Taquet et al. 2012; Aikawa et al. 2020).

Therefore, in this work, we use the Enrique-Romero et al. (2019) ice model and methodology, but we improve the calculations for a better accuracy of the computed energetics, including dispersion, as described in detail in Sect. 3.

3 Methodology

3.1 Electronic structure calculations

Given the importance of inter-molecular interactions in radical–radical reactions, we recomputed the stationary points of the potential energy surfaces (PES) previously reported by Enrique-Romero et al. (2019) using the Grimme’s D3 dispersion term including the Becke–Johnson damping (D3(BJ)) (Grimme et al. 2010, 2011), this way improving the description of the dispersion forces with respect to the previous work.

All DFT calculations were performed with the GAUSSIAN16 program package (Frisch et al. 2016). A benchmark study showed thatthe BHLYP hybrid density functional method is the best suited DFT method to study these reactions, with an average error of 3%, and a maximum error of 5.0% with respect to benchmark multi-reference CASPT2(2,2) calculations using OPENMOLCAS 18.09 (see Annex). We have chosen CASPT2(2,2) as the minimum level of post-HF theory, aware of the fact that a CASPT2 inclusive of full valence states would have been much better. The latter is prevented, however, by the size of our system.

Thus, stationary points were fully optimized using BHLYP (Becke 1993; Lee et al. 1988) combined with the standard 6-31+G(d,p) Pople basis set alongside the D3(BJ) dispersion term (Grimme et al. 2010, 2011). When needed, intrinsic reaction coordinate (IRC) calculations at the optimization theory level were carried out to ensure that the transition states connect with the corresponding minima. To balance the computational cost and chemical accuracy, reaction energetics were then refined by performing full BHLYP-D3(BJ)/6-311++G(2df,2pd) single-point energy calculations on the BHLYP-D3(BJ)/6-31+G(d,p) optimized stationary points. Improving chemical accuracy is a fundamental aspect when aiming at providing kinetic calculations and rate constants (including tunneling effects) (Álvarez-Barcia et al. 2018), as in the present work. Additionally, as shown in Rimola et al. (2018); Enrique-Romero et al. (2019, 2020), DFT is a cost-effective methodology with which a correct description of biradical systems can be achieved by using the unrestricted broken (spin)-symmetry approach (e.g., Neese 2004).

All optimized stationary points were characterized by the analytical calculation of the harmonic frequencies as minima and saddle points. Thermochemical corrections computed at BHLYP-D3(BJ)/6-31+G(d,p) were included to the single point BHLYP-D3(BJ)/6-311++G(2df,2pd) potential energy values using the standard rigid-rotor and harmonic oscillator formulae in order to obtain the zero-point vibrational energy (ZPE) corrections.

3.2 Kinetic calculations

In order to compute the rate constants for the chemical reactions between the radical pairs, we adapted our in-house kinetic code, based on the RRKM scheme for gas-phase reactions (Skouteris et al. 2018), to the surface plus adsorbed radicals case. First, we obtained the microcanonical rate constant k_aeb (E) at a given energy E as: $$k_{\text{aeb}}(E)=\frac{N(E)}{h\rho (E)},$$(8)

where N(E) is the sum of states for the active degrees of freedom in the transition state, ρ(E) is the density of states for the active degrees of freedom in the reactant, and h is the Planck constant. Since we aim to simulate a reaction taking place on a solid surface, only vibrational degrees of freedom are taken into account. Second, the obtained rate constants were Boltzmann-averaged in order to derive the rate constants as a function of the temperature.

For the H abstraction reaction, we took into account tunneling effects adopting the Eckart scheme via the unsymmetric potential energy barrier approach. In order to have a chemical system of reference to compare with, we applied the same method to the well studied reaction H + CO → HCO. In this case, the initial structures of the reaction were taken from the theoretical study by Rimola et al. (2014), which were re-optimized at the present work computational level. Here, from the optimized transition state, intrinsic reaction coordinate (IRC) calculations were run assuming a Langmuir–Hinshelwood (LH) like reaction, contrarily to the Rimola et al. (2014) original computations. All stationary points were characterized by frequency calculations, obtaining their (harmonic) vibrational modes and their zero-point energies. More details of these computations can be found in Appendix D.

4 Results

4.1 Energetics of the reactions

Table 1 presents the 0 K enthalpies (i.e., potential energies plus ZPE corrections) of the studied reactions. The improvement in the dispersion correction and the refinement of the DFT energy slightly decrease the energy barriers of each one of the reactions to form acetaldehyde and CO + CH₄ by less than 2.5 kJ mol⁻¹ with respect to the values quoted by Enrique-Romero et al. (2019). Inversely, the H + CO → HCO reaction has a higher barrier, 13.5 kJ mol⁻¹, than that quoted by Rimola et al. (2014), 9.2 kJ mol⁻¹, for two reasons: (i) Rimola et al. assumed an Eley–Rideal reaction (namely, the H atom comes from the gas phase and reacts with frozen CO), while here we have considered a LH mechanism (Sect. 3.2), and (ii) they did not consider dispersion corrections.

4.2 Rate constants

Figure 1 shows the rate constants as a function of the temperature of the reactions that form CH₃CHO and CO + CH₄ from the coupling and direct H-abstraction of CH₃ + HCO, respectively.The figure also reports the case of HCO formation from H + CO, for the sake of reference.

The rate constants of the reactions leading to CO + H₂CO and HCO take tunneling into account, which is evidenced by their deviation from linearity. It is also evident the strong temperaturedependence of the radical–radical reactions studied, as compared to HCO formation.

The rate constants of the acetaldehyde formation are larger than those of CO + CH₄ formation at temperatures above ~24 K. This is due to its lower barrier and the almost negligible quantum tunnelling contribution to HCO + CH₃ → CO + CH₄ at such temperatures. However, as the temperature decreases the tunnelling probability takes over deviating the rate constant of CO + CH₄ formation from linearity, becoming faster than the formation of acetaldehyde. On the contrary, HCO formation has a much weaker temperaturedependence and higher rate constants over the considered temperature range. This is the result of the dominant strong quantum tunnelling of the H atom through the reaction barrier, in agreement with the literature results (e.g., Andersson et al. 2011; Rimola et al. 2014).

In order to facilitate the introduction of the new rate constants in astrochemical models, we fit the reactions rate constants with the standard formula in Eq. (9). The values of α, β and γ are listed in Table 2. $$k_{\text{aeb}}(T)=\alpha {\left(\frac{T}{300\text{K}}\right)}^{\beta }\exp (-\gamma /T).$$(9)

Fig. 1 Arrhenius plots, namely rate constants as a function of the inverse of temperature, for the reaction CH3 + HCO forming acetaldehyde (black solid line) or CO + CH4 (black dashed line), and for the reaction H + CO → HCO (gray dotted-dashed line), described in the main text.

4.3 Desorption and diffusion temperatures

Table 1 reports the computed E_des and the temperature for desorption T_des and diffusion T_diff derived assuming a half-life of 1 Myr. E_des values are obtained at BHLYP-D3(BJ)/6-311++G(2df,2pd)//BHLYP-D3(BJ)/6-31+G(d,p) level following the procedure explained in Sect. 2.3, which moreover are corrected for deformation and basis set superposition energy. The T_des and T_diff values are obtained by using the standard equation for the half-life time, t_1∕2 = ln(2)∕k_diff/des(T). These timescales provide an estimation of the characteristic temperatures for desorption and diffusion of the two radicals, CH₃ and HCO, involved in the formation of acetaldehyde on the icy surface.

Finally, we note that our E_des are consistent with those computed by Ferrero et al. (2020) on a substantially larger ASW ice model. Specifically, our E_des in Table 1 lies in the high end of the Ferrero et al. range. On the contrary, and as already discussed in Ferrero et al. (2020), our E_des are different than those reported in the astrochemical databases KIDA¹ and UMIST², often used by modellers. Unfortunately, no experimental data on the CH₃ and HCO E_des desorption energy exist, to our best knowledge.

5 Discussion

5.1 Formation of acetaldehyde versus CO + CH₄

We used the results of our new calculations (Sect. 4) of the CH₃ + HCO reaction kinetics, desorption and diffusion rate constants to compute the efficiency ε (Eq. (3)) of the two channels leading to the formation of either CH₃CHO or CO + CH₄. As discussed in Sect. 2.3, given the uncertainty on its value, we considered three cases for the E_diff/E_des ratio f: 0.3, 0.4 and 0.5. Figure 2 shows the resulting ε as a function of the temperature and Table 2 reports the α, β and γ values obtained by fitting the ε curves with Eq. (9), for the three cases of f.

We note that, although we computed the efficiency of the reactions in the 5–100 K range, they will only take place as long as one of the two radicals can diffuse and scan the ASW sites and both radicals stay on the reaction site, namely they do not desorb (see Table 1). Consequently, the upper limit to the temperature where the CH₃CHO and CO + CH₄ formation reactions take place is set by the desorption of CH₃, as it has a lower desorption energy than HCO (30 and 68 K, respectively, see Table 1). Likewise, the lower limit is also set by the CH₃ diffusion energy only, which is equal to 15, 12 and 9 K for f equal to 0.5, 0.4 and 0.3, respectively. In the case of f equal to 0.5, HCO starts to be mobile when CH₃ has already sublimated, so that the efficiency of the reaction depends on CH₃ E_diff only. Conversely, for f equal to 0.4 and 0.3, the temperatures at which HCO and CH₃ can diffuse overlap, so that both species contribute to the denominator of Eq. (3).

Formation efficiency

For bothreactions, formation of CH₃CHO and CO + CH₄, the efficiency ε is about 1 in the 9–15 K range regardless of the f value (between 0.3 and 0.5) with one exception, acetaldehyde formation with f = 0.3, which starts at very low efficiency values and monotonically increases. For f = 0.5, either reactions have efficiencies of about unity in the whole range of temperatures (up to 30 K). For f = 0.4, formation of CO + CH₄ distances from unity at temperatures above ~22 K, going into lower values so that at 30 K it reaches ε ~0.3, while the efficiency of acetaldehyde formation stays about unity up to 30 K, where it takes a value of ~0.8. On the other hand, for f = 0.3 things are very different. The efficiency of CO + CH₄ crashed at higher temperatures, reaching values of about 0.001 at 30 K, while that of acetaldehyde never goes above ~0.01.

This is because, for relatively large f values (≥0.4), the most mobile radical, CH₃, moves slowly and the two radicals have plenty of time to react when they meet before one of them moves away: ε is, therefore, close to unity. However, when the timescale for diffusion becomes smaller than the reaction timescale (i.e., k_diff ≫ k_aeb), CH₃ moves away before having the time to react and the efficiency drops below unity. In practice, the smaller the E_diff∕E_des ratio, the faster CH₃ moves and the smaller ε. However, since both k_diff and k_aeb have an exponential dependence on the temperature, a change in behavior occurs when the reaction activation energy γ (Eq. (9)) is similar to E_diff and the efficiency ε strongly depends on the temperature. For the formation of acetaldehyde, γ = 663 K (Table 2) and, therefore, the change of behavior occurs when E_diff∕E_des ~ 0.40. In these cases, the lower the temperature, the larger k_diff with respect to k_aeb and the smaller ε, as shown in Fig. 2. Similar arguments hold also for the HCO + CH₃ → CO + CH₄ reaction. The only difference is that, at low temperatures, k_aeb deviates from the exponential law because of the kicking in of the tunneling effect that greatly increases k_aeb (giving a negative γ values: see Table 2). Since the tunneling is more efficient for decreasing temperature, the k_aeb/k_diff ratio decreases at increasing temperatures and, consequently, ε decreases.

Fig. 2 Reaction efficiency ε (Eq. (3)) of the reaction CH3 + HCO leading to either CH3CHO (solid lines) or CO + CH4 (dashed lines) as a function of the temperature. The computations were obtained adopting three different Ediff /Edes ratios: 0.3 (green), 0.4 (blue) and 0.5 (red). We note that, for Ediff/Edes = 0.5 the CH3CHO and CO + CH4 (red) curves overlap, namely they are constant and equal to 1.

Fig. 3 Branching ratio BR(T) of the formation rate of the CH3CHO over CO + CH4 (Eq. (10)) as a function of the temperature in the range where the reactions can occur, namely below 30 K (see text), for Ediff∕Edes equal to 0.3 (green), 0.4 (blue) and 0.5 (red).

Branching ratio

Figure 3 shows the branching ratio BR of the formation rate of CH₃CHO over CO + CH₄ as a function of the temperature, for the three f values (0.5, 0.4 and 0.3). The BR is obtained integrating Eq. (1) from the temperature at which CH₃ starts to be mobile T₀, a value that depends on the assumed f (see above), to the temperature T. It holds: $$\text{BR}(T)=\frac{{\displaystyle {\int }_{T_0}^T\text{d}}{T}^{\prime }{\epsilon }_{{\text{CH}}_3\text{CHO}}\times (R_{\text{diff},{\text{CH}}_3}+R_{\text{diff},\text{HCO}})}{{\displaystyle {\int }_{T_0}^T\text{d}}{T}^{\prime }({\epsilon }_{{\text{CH}}_3\text{CHO}}+{\epsilon }_{\text{CO}\text{+}{\text{CH}}_{\text{4}}})\times (R_{\text{diff},{\text{CH}}_3}+R_{\text{diff},\text{HCO}})}.$$(10)

The different effects commented above can be clearly seen in Fig. 3 and can be summarized as follows. For f = 0.5, the branching ratio BR is constant and equal to 0.5, namely the HCO + CH₃ reaction leads to acetaldehyde and CO + CH₄ in equal quantities. For f = 0.4, BR lies in the range 0.4–0.5 up to 25 K and then it becomes larger, because the tunneling gain in the CO + CH₄ production at low temperatures vanishes. For f = 0.3 (and, in general, ≤0.4), BR is <0.5 at temperatures less than ~25 K and rises to ~0.9 at 30 K.

In other words, for f ≥ 0.4, acetaldehyde and CO + CH₄ are in approximately equal competition in the range of temperatures where the HCO + CH₃ reaction can occur. However, for f < 0.4, acetaldehyde is a very minor product for temperatures lower than about 25 K and flips to be a major product above it.

5.2 The experimental point of view

Experiments studying the formation of acetaldehyde from radical–radical coupling date back to the 1990s (Hudson & Moore 1997). They are mainly based on energetic (UV or particles) irradiation of different H₂ O, CO, CH₃OH and CH₄ ice mixtures (e.g. Bennett et al. 2005; Öberg et al. 2010; Martín-Doménech et al. 2020). In relation to experimental acetaldehyde formation on grain surfaces, Bennett et al. (2005), after irradiation of a CO:CH₄ ice mixture, detected acetaldehyde and predicted that the orientation of the CH₃ and HCO radicals are crucial in the efficiency of the reaction. On the other hand, Martín-Doménech et al. (2020) conducted laboratory experiments on the formation of acetaldehyde via the CH₃ + HCO reaction, concluding that this channel is not efficient enough to reproduce the astronomical observations.

As discussed by various authors, although laboratory experiments are primordial in suggesting possible mechanisms operating in the ISM and, specifically, the possible formation routes of molecules on the interstellar grain surfaces, they cannot provide the exact ISM conditions or a detailed description of the mechanisms at the atomic level. Despite this, the improvement of radical detection methods, such as the electron paramagnetic resonance (EPR) technique, will help to clarify the role of radicals generated ininterstellar ice analogs (e.g., Zhitnikov & Dmitriev 2002). In this respect, therefore, theoretical computations as those reported in this work constitute a complementary, if not unique, tool to understand the interstellar surface chemistry.

5.3 Astrophysical implications

In astrochemical models, it is generally assumed that reactions between radicals on the surface of interstellar ices are barrierless and, consequently, that their efficiency is equal to 1 (Sect. 2). In addition, it is also often assumed that there are no competition channels to the production of iCOMs. At variance with these simpleassumptions, our new calculations presented in Sect. 5.1, indicate that, at low (≤ 15 K) temperatures, the efficiency of the acetaldehyde formation is close to unity, for a E_diff∕E_des ratio f ≥ 0.40. However, there is a competing channel leading to CO + CH₄, for which the efficiency is also equal to 1, so that, at low temperatures and for f ≥ 0.40 the two channels are equally probable. The acetaldehyde formation efficiency remains close to unity in the temperature range where the reaction can occur, namely at ≤30 K, for f ≥ 0.40. However, the situation drastically changes for f < 0.40. Specifically, for f = 0.3, the efficiency of acetaldehyde formation crashes to very low values and increases with temperature to a maximum of 0.01 at 30 K. Similarly, the formation of CO + CH₄ drops to 1.5 × 10⁻³ at 30 K.

Therefore, two major messages come out from our calculations: (1) the efficiency of the formation of acetaldehyde from the HCO + CH₃ reaction on icy surfaces is a complex function of the temperature and of the CH₃ diffusion energy E_diff (Fig. 2) and (2) the acetaldehyde formation receives competition with CO + CH₄ formation, which cannot be neglected and whose efficiency is also a complex function of temperature and E_diff (Fig. 3).

While the dependence on the temperature and the importance of the competition of other products were already recognized (Enrique-Romero et al. 2019), the paramount importance of the diffusion energy E_diff in the radical–radical reactions efficiency was not appreciated, at least not at the extent indicated by this study (because of the assumption of astrochemical models that the efficiency of the radical–radical reactions on grains is 1). Penteado et al. (2017), for example, carried out an extensive study of the surface chemistry on the binding energies (namely, our E_des) showing how critical they are. Our new study suggests that E_diff is as much, or even more, crucial in the reactions involving two radicals on ASW.

What makes the situation actually critical is that, while studies of the binding energy of radicals can and have been estimated in experimental and theoretical works (e.g., see the recent works by Penteado et al. 2017; Ferrero et al. 2020), evaluating the diffusion energy of multi-atomic radicals on cold icy surfaces has proven to be extremely complicated and, to the best of our knowledge, no experimental or theoretical studies exist in the literature (see also e.g., Cuppen et al. 2017; Potapov & McCoustra 2021). Indeed, as mentioned above, obtaining E_diff experimentally is hitherto hampered by technical limitations on the instrumentation used to detect the radicals (i.e., EPR measurements). In relation to theoretical investigations, this lacking in bibliography is due to the convergence of methodological difficulties that make the study of diffusion with computational simulations intrinsically complex (but also compelling). Diffusion can currently be studied by means of molecular dynamics (MD) or kinetic Monte Carlo (kMC) simulations. With the first, to obtain a sufficient representativeness of the species diffusion, long simulation time-scales are mandatory. This in practice means to adopt classical force fields, in which the electronic structure of the systems is missing. However, radicals are open-shell species (with at least one unpaired electron) and accordingly electrons have to be accounted for. Thus MD simulations should be grounded within the quantum mechanics realm, which are much more expensive than the classical ones, making the MD simulations unfeasible. The alternative would be the adoption of kMC simulations. However, these simulations require building a complete network of the site-to-site radical hopping, in which for each hopping the corresponding rate constant has to be known a priori. This actually means to localize for each hopping the corresponding transition state structure (at a quantum chemical level), in which by using a realistic ASW model (i.e., large, amorphous and accordingly plenty of binding sites) makes the problem unpractical.

Usually, astrochemical models assume that the radical E_diff is a fraction f of E_des and the value f is derived from computations and experiments on species such as CO, CO₂, H₂ O, CH₄ and NH₃ (e.g., Mispelaer et al. 2013; Karssemeijer & Cuppen 2014; Lauck et al. 2015; Ghesquière et al. 2015; He et al. 2017, 2018; Cooke et al. 2018; Maté et al. 2020; Kouchi et al. 2020). These studies give a value for f between 0.3 and 0.6, as mentioned in Sect. 2.3. However, rigorously speaking, the experiments do not necessarily measure the same diffusion processes as in interstellar conditions, for at least the reasons of the surface coverage (He et al. 2018) and its dependence of the nature of the ice, specifically its degree of porosity (Maté et al. 2020), which is poorly known in the case of interstellar ices. As a matter of fact, using experimental and a theoretical Monte Carlo code, Maté et al. (2020) found that “the microscopic diffusion is many times faster than the macroscopic diffusion measured experimentally". Most recently, Kouchi et al. (2020) obtained a direct measurement of the diffusion energy of CO and CO₂ on ASW, using the transmission electron microscopy (TEM) technique, which allows a direct measurement of the surface diffusion coefficients (against the often used technique of IR spectroscopy, which only indirectly estimates the diffusion energy). Kouchi and coworkers found an f ratio equal to 0.3.

We have seen that, in the case of acetaldehyde, this uncertainty on f has a dramatic effect. If f is >0.4, the efficiency of acetaldehyde formation is equal to 1 and it is about equal to that of the CO + CH₄ formation. On the contrary, if f is equal to 0.3, then the efficiency of acetaldehyde formation (and CO + CH₄) crashes, to a maximum value of 0.02. The most recent measurements by Kouchi et al. (2020) point out the latter case as the most probable. If the value f = 0.3 is confirmed, then acetaldehyde is unlikely to be formed on the interstellar icy grain surfaces.

One could be tempted to use the astronomical observations against the astrochemical model predictions to add constraints to the f value in (real) interstellar ices. Of course, given the large number of parameters associated with the astrochemical models it could be a dangerous exercise. Nonetheless, we can analyze two cases, as illustrative examples. Barger & Garrod (2020) compared the predictions of their model, where the formation of acetaldehyde is dominated by the reaction CH₃ + HCO assumed to have ε = 1, with the observations toward various hot cores and found that in two of them, NGC 7538 IRS 1 and W3(H₂O), their model overproduces the acetaldehyde column densities by more than a factor 10³ with respect to the observed ones. If f is equal to 0.3, introducing our new values for ε could possibly cure this mismatch. On the contrary, Jørgensen et al. (2016) found a good agreement between the observed abundances of acetaldehyde in IRAS16293B and SgrB2(N) and those predicted by the Garrod (2013) model. In this case, the agreement would point to f ≥ 0.4. In other words, our new computations might solve the mismatch observed toward NGC 7538 IRS 1 and W3(H₂O) if f = 0.3, but they would create a mismatch on the observations toward IRAS16293B and SgrB2(N), or vice versa. Alternatively, it is possible that f varies in different sources, belonging to different environments. For example, one could think that sources in cold quiescent regions have ices different, more or less porous, from those in warm and chaotic ones. The two examples discussed above, unfortunately, do not lead to a coherent behavior, as, for example, IRAS16293B and SgrB2(N) could not belong to more different environments.

In conclusion, the acetaldehyde formation by radical–radical recombination on the ices is such a strong function of the diffusion energy, likely linked to the nature of the ice, that a little variation of the E_diff∕E_des (by 0.1) value can shift the efficiency from 1 to less than 0.01. The most recent estimates of E_diff∕E_des suggest a value of 0.3 (Kouchi et al. 2020), which would make the formation of acetaldehyde on the grain surfaces unlikely. Anyway, the important message here is that astrochemical model predictions should be taken with a certain precaution. On the contrary, our new computations clearly show the huge importance of better knowing the microprocesses involved in the radical–radical chemistry on the icy interstellar grains and the urgent need of extensive studies, similar to the one presented here, on different systems believed to form iCOMs on the interstellar icy surfaces.

For the sake of completeness, the formation of acetaldehyde via radical–radical reactions on surfaces with lower binding energies such as solid CO could have a higher efficiency, due to the low radical–surface interactions (Lamberts et al. 2019). Finally, it is worth reminding that acetaldehyde can alternatively be synthesized in the gas-phase (e.g., Charnley 2004; Vastel et al. 2014; De Simone et al. 2020). Recently, Vazart et al. (2020) reviewed the gas-phase routes leading to acetaldehyde and found that, very likely, the dominant one is that starting from ethanol, the so-called ethanol tree (Skouteris et al. 2018). In particular, the ethanol tree route reproduces quite well both the acetaldehyde and glycolaldehyde abundances in the sources where ethanol was also observed, including IRAS16293B (Vazart et al. 2020).

6 Conclusions

In this work, we report new computations on the energetics and kinetics of the reaction HCO + CH₃, which can lead to the formation of either acetaldehyde or CH₄ + CO. Specifically, we compute the rate constants of both reactions as a function of temperature as well as the efficiency of the formation of acetaldehyde and CH₄ + CO, respectively, combining reaction kinetics at RRKM (tunneling included) with diffusion and desorption competitive channels. We provide analytical formulae so that the computed rate constants and efficiency can be easily introduced in astrochemical models.

The main conclusions of our study are the following:

• 1.

  The HCO + CH₃ reaction can only occur when the surface temperature is lower than 30 K, because CH₃ desorbs at larger temperatures.

• 2.

  Our computations suggest that acetaldehyde is not the dominant product for the reaction HCO + CH₃. The efficiency ε of its formation strongly depends on the E_diff∕E_des ratio, providing dramatic variations between 0.3–0.5 values, the most usually used values in astrochemical models.

• 3.

  At low (≤15 K) temperatures, ε is close to unity for both the formation of acetaldehyde and its competing CO + CH₄ channel for f ≥ 0.4, while only the efficiency of CO + CH₄ is unity at these temperatures for f = 0.3 thanks to quantum tunnelling. The efficiency of acetaldehyde formation remains unity in the range of temperatures where the reaction can occur (≤30 K) for E_diff∕E_des ≥ 0.40. For lower E_diff∕E_des ratios, ε becomes ≪1 and increases with increasing temperature: in the case of E_diff∕E_des = 0.3, it reaches a maximum of ~0.01 at 30 K. Conversely, the efficiency of the formation of CO + CH₄ increases with decreasing temperature because of the tunneling.

• 4.

  Thesevariant ε values as a function of E_diff∕E_des go against the assumption made in many astrochemical models, in which ε is equal to 1. This might have a substantial impact on the acetaldehyde abundance predicted by these models, which may overestimate it by a few orders of magnitude.

• 5.

  We discussed the example of IRAS16293B and suggested that, in this object, acetaldehyde is likely synthesized by a gas-phase reaction route that starts from ethanol.

Finally,this new study calls for specific similar computations on the radical–radical reactions assumed to form iCOMs in astrochemical models as assuming that they have efficiency ε equal to 1 and are the only reaction product could be highly misleading.

Acknowledgements

This project has received funding within the European Union’s Horizon 2020 research and innovation programme from the European Research Council (ERC) for the projects “The Dawn of Organic Chemistry” (DOC), grant agreement no. 741002 and “Quantum Chemistry on Interstellar Grains” (QUANTUMGRAIN), grant agreement no. 865 657, and from the Marie Sklodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement no. 811312. AR is indebted to “Ramón y Cajal” program. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged. Finally, we thank Prof. Gretobape for fruitful and stimulating discussions. Most of the quantum chemistry calculations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervized by the Agence Nationale pour la Recherche. Additionally this work was granted access to the HPC resources of IDRIS under the allocation 2019-A0060810797 attributed by GENCI (Grand Equipement National de Calcul Intensif). We thank prof. Gretobape for stimulating discussions, and J. Perrero for useful contributions.

Appendix A Benchmark study

The quality of BHLYP-D3(BJ) as an accurate, cost-effective method for the reactions studied in this work is shown in this section.

We have taken five hybrid DFT dispersion-corrected methods: MPWB1K-D3(BJ), M062X-D3, PW6B95-D3(BJ), wB97X-D3 and BHandHLYP-D3(BJ), recommended in Goerigk et al. (2017) for their good overall performance.

We have then compared their performance with respect to CASPT2 by studying reactions I and II on 2 water molecules. We proceeded in two steps: (i) we performed geometry optimized and run frequency calculations at BHLYP-D3(BJ)/6-311++G(d,p) for each reaction channel, finding and checking each stationary point (i.e., reactants, transition state and products); (ii) we then run single point calculations on top of these geometries at each DFT method combined with the 6-311++G(2df,2pd) basis set, and CASPT2/aug-cc-PVTZ for reference. The unrestricted broken symmetry scheme was adopted for all DFT calculations, and the CASPT2 guess wave function was generated using a CASSCF(2,2) calculation where the active space is formed by the two unpaired electrons and their molecular orbitals, starting from a triplet Hartree-Fock wave-function.

BHLYP-D3(BJ) gives the best overall performance with and average unsigned error of 3.0%, and a maximum of 5.0% with respect to CASPT2/aug-cc-PVTZ. The rest of DFT methods have average errors between 10 and 80% (Table A.1). The raw energy values are shown in Table A.3.

Fig. A.1 Geometries of the stationary points (reactants, transition state and products) of reactions a) I and b) II on top of two water molecules optimized at BHLYP-D3(BJ)/6-311++G(d,p) level. Distances in Å.

Appendix B Calculation of desorption and diffusion rate constants

In order to calculate the rates of diffusion and desorption we have used Eyring’s equation (Eq. 5), where entropy and thermal corrections to the enthalpy are accounted for since energy term in the exponential should be the Gibbs free energy (G = H − TS). In the astrochemistry community, desorption and diffusion rates consist of two parts, first the attempt frequency (i.e., the pre-exponential factor) and then the exponential with the binding energies. We define the attempt frequency within the Eyring equation as $$\nu =\frac{k_{B}T}{h}\times \exp (\Delta S/k_B),$$(B.1)

where h, k_B are the Planck and Boltzmann constants, T the temperature and S the entropy. For the ice model and the complex surface+radical the entropy contributions are S = S_vib, where the vibrational counterpart is (Eq. B.2): $$S_{vib}=k_B{\sum }_i\left[\frac{{\theta }_{v,i}/T}{e^{{\theta }_{v,i}/T}-1}-\ln \left(1-e^{-{\theta }_{v,i}/T}\right)\right]$$(B.2)

where ${\theta }_{v,i}=hc{\overline{\nu }}_i/k_B$ with c the speed of light and ${\overline{\nu }}_i$ the ith vibrational mode frequency in cm⁻¹.

For the free radicals, additionally, we also take into account the rotational and translational contributions: S = S_vib + S_rot + S_trans, where the rotational and translational counterparts are given by eqs. B.3 and B.4 respectively: $$S_{rot}=k_B\left[\ln \left(\frac{{\pi }^{1/2}}{{\sigma }_{rot}}\frac{T^{3/2}}{{({\theta }_{r,x}{\theta }_{r,y}{\theta }_{r,z})}^{1/2}}\right)+\frac{3}{2}\right],$$(B.3) $$S_{trans}=k_B\left(\ln \left[{\left(\frac{2\pi m{k}_B}{h^2}\right)}^{3/2}\frac{k_{B}T}{P}\right]+\frac{5}{2}\right),$$(B.4)

where σ_rot is the rotational symmetry number (6 for CH₃ and 1 for HCO), θ_r,i = B_ih∕k_B with B_i the ith axis rotational constant (in s⁻¹), m is the mass of the radical and P the gas pressure, which was calculated assuming a density of 10⁴ cm⁻³.

The thermal corrections follow: H = E_DFT + ZPE + E_vib(T) + k_BT for the surface and the surface + radical complex. Here E_DFT is the energy obtained from our DFT calculations, ZPE is the zero-point energy and E_vib(T) is the thermal vibrational energy calculated with Eq. B.5. In analogy to to the entropy, the rotational and translational contributions are also included for the free radicals: H = E_DFT + ZPE + E_vib(T) + k_BT + H_rot + H_trans where $H_{rot}=H_{trans}=\frac{3}{2}{k}_{B}T$. $$E_{vib}(T)=k_B{\sum }_i\frac{{\theta }_{v,i}}{e^{{\theta }_{v,i}/T}-1}$$(B.5)

The magnitude of these contributions is shown in Table B.1, where it can be seen that the most important contribution to the enthalpy comes from the vibrational modes, while for entropy all contributions are rather small:

As it was explained in the introduction, we have adopted the assumption that the barrier for diffusion can be expressed as a fraction of that of desorption, therefore in this model we multiply the ΔG of desorption times these fractions (assumed to be 0.3, 0.4 and 0.5). With this, we get different attempt frequencies for diffusion and desorption (see Table B.2), all around 10⁸ –10¹¹ s⁻¹, smaller than the normally used approach of the harmonic oscillator³ in astrochemical models (e.g., Tielens & Allamandola (1987); Hasegawa et al. (1992)) due to the use of Eyring’s relation and the inclusion of entropy.

Appendix C Rate constant comparison

Figure C.1 compares the rates of the radical-radical reactions with the hopping and desorption rates for each radical species, using three different criteria for the diffusion barrier, namely making it 0.3, 0.4 and 0.5 times those of desorption.

Fig. C.1 Comparison of reaction, diffusion and desorption rate constants involved in the CH3 + HCO system. Notice that the desorption rate of HCO is not seen as it appears at very low rate constant values. Numbers in brackets indicate the diffusion-to-desorption energy barrier ratio.

Appendix D H + CO PES and data

Fig. D.1 Potential energy surface of the H + CO → HCO reaction, in kJ mol−1. Energies are corrected for dispersion and ZPE. Geometries and ZPE energies were obtained at UBHandHLYP-D3(BJ)/6-31+G(d,p) level and DFT energies were refined at UBHandHLYP-D3(BJ)/6-311++G(2df,2pd) level. Reactants and products were obtained by running intrinsic reaction coordinate calculations.

Appendix E Radical-radical TS data and PES

Fig. E.1 Potential energy surfaces Reacts. I and II. Geometries and ZPE energy correction was obtained at UBHandHLYP-D3(BJ)/6-31+G(d,p) level, DFT energy was refined at UBHandHLYP-D3(BJ)/6-311++G(2df,2pd) level. Energy units are in kJ mol−1.

Appendix F Efficiency figures, separated by E_diff∕E_des ratios

Fig. F.1 CH3 + HCO reaction efficiencies ε (Eq. 3), assuming diffusion barriers 0.3, 0.4 and 0.5 times those of desorption (panels from top to bottom). The green-colored regions indicate the diffusion and desorption temperatures limits of CH3, while the red ones are the same for HCO.

Appendix G Calculation of diffusion and desorption temperatures

Half-lives are calculated from the rate constants (k_i, with i being either the diffusion or desorption of radicals) following $$N_{1/2}=N_0\exp (-k_{i}t)⟺t_{1/2}=\frac{\ln (2)}{k_i},$$

and shown in Fig. G.1.

Fig. G.1 Diffusion and diffusion temperatures of CH3 and HCO assuming a half-life of 1 Myrs for desorption.

Appendix H Fittings to k_aeb and ε

Fig. H.1 Fittings to the computed rate constants (Figure 1) with Eq. 9 for reactions I (left hand side panel) and II (left hand side panel).

Fig. H.2 Fittings (solid lines) to the computed efficiency factors (points) using Eq. 9 for acetaldehyde using Ediff ∕Edes = 0.5, 0.4 and 0.3 (left to right panels).

Fig. H.3 Fittings (solid lines) to the computed efficiency factors (points) using Eq. 9 for CO + CH4 formation using Ediff∕Edes = 0.5, 0.4 and 0.3 (left to right panels).

References

All Tables

All Figures

	Fig. 1 Arrhenius plots, namely rate constants as a function of the inverse of temperature, for the reaction CH3 + HCO forming acetaldehyde (black solid line) or CO + CH4 (black dashed line), and for the reaction H + CO → HCO (gray dotted-dashed line), described in the main text.
In the text

Fig. 2 Reaction efficiency ε (Eq. (3)) of the reaction CH3 + HCO leading to either CH3CHO (solid lines) or CO + CH4 (dashed lines) as a function of the temperature. The computations were obtained adopting three different Ediff /Edes ratios: 0.3 (green), 0.4 (blue) and 0.5 (red). We note that, for Ediff/Edes = 0.5 the CH3CHO and CO + CH4 (red) curves overlap, namely they are constant and equal to 1.

In the text

	Fig. 3 Branching ratio BR(T) of the formation rate of the CH3CHO over CO + CH4 (Eq. (10)) as a function of the temperature in the range where the reactions can occur, namely below 30 K (see text), for Ediff∕Edes equal to 0.3 (green), 0.4 (blue) and 0.5 (red).
In the text

	Fig. A.1 Geometries of the stationary points (reactants, transition state and products) of reactions a) I and b) II on top of two water molecules optimized at BHLYP-D3(BJ)/6-311++G(d,p) level. Distances in Å.
In the text

	Fig. C.1 Comparison of reaction, diffusion and desorption rate constants involved in the CH3 + HCO system. Notice that the desorption rate of HCO is not seen as it appears at very low rate constant values. Numbers in brackets indicate the diffusion-to-desorption energy barrier ratio.
In the text

	Fig. D.1 Potential energy surface of the H + CO → HCO reaction, in kJ mol−1. Energies are corrected for dispersion and ZPE. Geometries and ZPE energies were obtained at UBHandHLYP-D3(BJ)/6-31+G(d,p) level and DFT energies were refined at UBHandHLYP-D3(BJ)/6-311++G(2df,2pd) level. Reactants and products were obtained by running intrinsic reaction coordinate calculations.
In the text

	Fig. E.1 Potential energy surfaces Reacts. I and II. Geometries and ZPE energy correction was obtained at UBHandHLYP-D3(BJ)/6-31+G(d,p) level, DFT energy was refined at UBHandHLYP-D3(BJ)/6-311++G(2df,2pd) level. Energy units are in kJ mol−1.
In the text

	Fig. F.1 CH3 + HCO reaction efficiencies ε (Eq. 3), assuming diffusion barriers 0.3, 0.4 and 0.5 times those of desorption (panels from top to bottom). The green-colored regions indicate the diffusion and desorption temperatures limits of CH3, while the red ones are the same for HCO.
In the text

	Fig. G.1 Diffusion and diffusion temperatures of CH3 and HCO assuming a half-life of 1 Myrs for desorption.
In the text

	Fig. H.1 Fittings to the computed rate constants (Figure 1) with Eq. 9 for reactions I (left hand side panel) and II (left hand side panel).
In the text

	Fig. H.2 Fittings (solid lines) to the computed efficiency factors (points) using Eq. 9 for acetaldehyde using Ediff ∕Edes = 0.5, 0.4 and 0.3 (left to right panels).
In the text

	Fig. H.3 Fittings (solid lines) to the computed efficiency factors (points) using Eq. 9 for CO + CH4 formation using Ediff∕Edes = 0.5, 0.4 and 0.3 (left to right panels).
In the text