Issue 
A&A
Volume 655, November 2021



Article Number  A9  
Number of page(s)  14  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202141531  
Published online  28 October 2021 
Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains
^{1}
Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG),
38000
Grenoble, France
email: juan.enriqueromero@univgrenoblealpes.fr; cecilia.ceccarelli@univgrenoblealpes.fr
^{2}
Departament de Química, Universitat Autònoma de Barcelona, Bellaterra,
08193,
Catalonia, Spain
email: albert.rimola@uab.cat
^{3}
MasterTech,
06123
Perugia, Italy
^{4}
Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia,
Via Elce di Sotto 8,
06123
Perugia, Italy
^{5}
Osservatorio Astrofisico di Arcetri,
Largo E. Fermi 5,
50125
Firenze, Italy
^{6}
Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino,
via P. Giuria 7,
10125,
Torino, Italy
Received:
11
June
2021
Accepted:
19
August
2021
Context. Interstellar grains are known to be important actors in the formation of interstellar molecules such as H_{2}, water, ammonia, and methanol. It has been suggested that the socalled interstellar complex organic molecules (iCOMs) are also formed on the interstellar grain icy surfaces by the combination of radicals via reactions assumed to have an efficiency equal to unity.
Aims. In this work, we aim to investigate the robustness or weakness of this assumption. In particular, we consider the case of acetaldehyde (CH_{3}CHO), one of the most abundant and commonly identified iCOMs, as a starting study case. In the literature, it has been postulated that acetaldehyde is formed on the icy surfaces via the combination of HCO and CH_{3}. Here we report new theoretical computations on the efficiency of its formation.
Methods. To this end, we coupled quantum chemical calculations of the energetics and kinetics of the reaction CH_{3} + HCO, which can lead to the formation of CH_{3}CHO or CO + CH_{4}. Specifically, we combined reaction kinetics computed with the RiceRamsperger–Kassel–Marcus theory (tunneling included) method with diffusion and desorption competitive channels. We provide the results of our computations in the format used by astrochemical models to facilitate their exploitation.
Results. Our new computations indicate that the efficiency of acetaldehyde formation on the icy surfaces is a complex function of the temperature and, more importantly, of the assumed diffusion over binding energy ratio f of the CH_{3} radical. If the ratio f is ≥0.4, the efficiency is equal to unity in the range where the reaction can occur, namely between 12 and 30 K. However, if f is smaller, the efficiency dramatically crashes: with f = 0.3, it is at most 0.01. In addition, the formation of acetaldehyde is always in competition with that of CO + CH_{4}.
Conclusions. Given the poor understanding of the diffusion over binding energy ratio f and the dramatic effect it has on the formation, or not, of acetaldehyde via the combination of HCO and CH_{3} on icy surfaces, model predictions based on the formation efficiency equal to one should to be taken with precaution. The latest measurements of f suggest f = 0.3 and, if confirmed for CH_{3}, this would rule out the formation of acetaldehyde on the interstellar icy surfaces. We recall the alternative possibility, which was recently reviewed, that acetaldehyde could be synthesized in the gas phase starting from ethanol. Finally, our computations show the paramount importance played by the microphysics involved in the interstellar surface chemistry and call for extensive similar studies on different systems believed to form iCOMs on the interstellar icy surfaces.
Key words: astrochemistry / diffusion / ISM: molecules / molecular processes
© J. EnriqueRomero et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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_{2}, 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 socalled 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 grainsurface 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; EnriqueRomero et al. 2019, 2020). For example, theoretical calculations on the energetics showed that the formation of acetaldehyde (CH_{3} CHO) on the icy surfaces via HCO + CH_{3} is in competition with the formation of CO + CH_{4} via direct Habstraction. 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 Habstraction 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_{3} 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 EnriqueRomero et al. (2019, 2020), we here consider the two competing reactions that can arise from the reactivity between HCO and CH_{3}:
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 EnriqueRomero et al. (2019) as a base and the models of amorphous solid water (ASW) for the ice described in Rimola et al. (2018) and EnriqueRomero 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 Surfacereaction rate definition
Generally, astrochemical models solve the timedependent 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_{ij}n_{i}n_{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): (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: (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: (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 icewater molecules plus the two adsorbed radicals, namely we consider the watercluster plus the radicals as a supermolecule isolated from its surrounding. The system A′ is the product of the radical–radical reaction on the icy surface, namely the watercluster plus either the radical–radical recombination (e.g., React. I) or the Habstraction (e.g., React. II) products.
It is important to note that in order to apply the RRKM theory, we implicitly assume that the intramolecular 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 and on , 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} . According to TST, the general expression of the rate constant k for a unimolecular reaction (like diffusion and desorption) is: (4)
where ΔV^{≠} is the zero point energycorrected 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: (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 desorptionrelated 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: (6)
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 diffusionrelated 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_{2}. 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 monolayer.
Given the uncertainty on the E_{diff} /E_{des} ratio for CH_{3} 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_{2} and O_{3} (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 (BarNun 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, EnriqueRomero 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; EnriqueRomero 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 EnriqueRomero 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 intermolecular interactions in radical–radical reactions, we recomputed the stationary points of the potential energy surfaces (PES) previously reported by EnriqueRomero 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 multireference CASPT2(2,2) calculations using OPENMOLCAS 18.09 (see Annex). We have chosen CASPT2(2,2) as the minimum level of postHF 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 631+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 BHLYPD3(BJ)/6311++G(2df,2pd) singlepoint energy calculations on the BHLYPD3(BJ)/631+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) (ÁlvarezBarcia et al. 2018), as in the present work. Additionally, as shown in Rimola et al. (2018); EnriqueRomero et al. (2019, 2020), DFT is a costeffective 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 BHLYPD3(BJ)/631+G(d,p) were included to the single point BHLYPD3(BJ)/6311++G(2df,2pd) potential energy values using the standard rigidrotor and harmonic oscillator formulae in order to obtain the zeropoint 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 inhouse kinetic code, based on the RRKM scheme for gasphase 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: (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 Boltzmannaveraged 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 reoptimized 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 zeropoint 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_{4} by less than 2.5 kJ mol^{−1} with respect to the values quoted by EnriqueRomero et al. (2019). Inversely, the H + CO → HCO reaction has a higher barrier, 13.5 kJ mol^{−1}, than that quoted by Rimola et al. (2014), 9.2 kJ mol^{−1}, 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.
Energetics and related parameters of the reactions and desorption and diffusion of the radicals.
4.2 Rate constants
Figure 1 shows the rate constants as a function of the temperature of the reactions that form CH_{3}CHO and CO + CH_{4} from the coupling and direct Habstraction of CH_{3} + 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_{2}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_{4} formation at temperatures above ~24 K. This is due to its lower barrier and the almost negligible quantum tunnelling contribution to HCO + CH_{3} → CO + CH_{4} at such temperatures. However, as the temperature decreases the tunnelling probability takes over deviating the rate constant of CO + CH_{4} 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. (9)
Fig. 1
Arrhenius plots, namely rate constants as a function of the inverse of temperature, for the reaction CH_{3} + HCO forming acetaldehyde (black solid line) or CO + CH_{4} (black dashed line), and for the reaction H + CO → HCO (gray dotteddashed 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 halflife of 1 Myr. E_{des} values are obtained at BHLYPD3(BJ)/6311++G(2df,2pd)//BHLYPD3(BJ)/631+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 halflife 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_{3} 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^{1} and UMIST^{2}, often used by modellers. Unfortunately, no experimental data on the CH_{3} and HCO E_{des} desorption energy exist, to our best knowledge.
5 Discussion
5.1 Formation of acetaldehyde versus CO + CH_{4}
We used the results of our new calculations (Sect. 4) of the CH_{3} + 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_{3}CHO or CO + CH_{4}. 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_{3}CHO and CO + CH_{4} formation reactions take place is set by the desorption of CH_{3}, 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_{3} 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_{3} has already sublimated, so that the efficiency of the reaction depends on CH_{3} E_{diff} only. Conversely, for f equal to 0.4 and 0.3, the temperatures at which HCO and CH_{3} can diffuse overlap, so that both species contribute to the denominator of Eq. (3).
Rate constants k_{aeb} (in s^{−1}) and efficiency ε of the two possible reactions between HCO and CH_{3}.
Formation efficiency
For bothreactions, formation of CH_{3}CHO and CO + CH_{4}, 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_{4} 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_{4} 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_{3}, 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_{3} 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_{3} 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_{3} → CO + CH_{4} 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 CH_{3} + HCO leading to either CH_{3}CHO (solid lines) or CO + CH_{4} (dashed lines) as a function of the temperature. The computations were obtained adopting three different E_{diff} /E_{des} ratios: 0.3 (green), 0.4 (blue) and 0.5 (red). We note that, for E_{diff}/E_{des} = 0.5 the CH_{3}CHO and CO + CH_{4} (red) curves overlap, namely they are constant and equal to 1. 
Fig. 3
Branching ratio BR(T) of the formation rate of the CH_{3}CHO over CO + CH_{4} (Eq. (10)) as a function of the temperature in the range where the reactions can occur, namely below 30 K (see text), for E_{diff}∕E_{des} 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_{3}CHO over CO + CH_{4} 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_{3} starts to be mobile T_{0}, a value that depends on the assumed f (see above), to the temperature T. It holds: (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_{3} reaction leads to acetaldehyde and CO + CH_{4} 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_{4} 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_{4} are in approximately equal competition in the range of temperatures where the HCO + CH_{3} 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_{2} O, CO, CH_{3}OH and CH_{4} ice mixtures (e.g. Bennett et al. 2005; Öberg et al. 2010; MartínDoménech et al. 2020). In relation to experimental acetaldehyde formation on grain surfaces, Bennett et al. (2005), after irradiation of a CO:CH_{4} ice mixture, detected acetaldehyde and predicted that the orientation of the CH_{3} and HCO radicals are crucial in the efficiency of the reaction. On the other hand, MartínDoménech et al. (2020) conducted laboratory experiments on the formation of acetaldehyde via the CH_{3} + 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_{4}, 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_{4} drops to 1.5 × 10^{−3} at 30 K.
Therefore, two major messages come out from our calculations: (1) the efficiency of the formation of acetaldehyde from the HCO + CH_{3} reaction on icy surfaces is a complex function of the temperature and of the CH_{3} diffusion energy E_{diff} (Fig. 2) and (2) the acetaldehyde formation receives competition with CO + CH_{4} 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 (EnriqueRomero 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 multiatomic 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 timescales are mandatory. This in practice means to adopt classical force fields, in which the electronic structure of the systems is missing. However, radicals are openshell 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 sitetosite 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_{2}, H_{2} O, CH_{4} and NH_{3} (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_{2} 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_{4} formation. On the contrary, if f is equal to 0.3, then the efficiency of acetaldehyde formation (and CO + CH_{4}) 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_{3} + 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_{2}O), their model overproduces the acetaldehyde column densities by more than a factor 10^{3} 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_{2}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 gasphase (e.g., Charnley 2004; Vastel et al. 2014; De Simone et al. 2020). Recently, Vazart et al. (2020) reviewed the gasphase routes leading to acetaldehyde and found that, very likely, the dominant one is that starting from ethanol, the socalled 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_{3}, which can lead to the formation of either acetaldehyde or CH_{4} + 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_{4} + 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_{3} reaction can only occur when the surface temperature is lower than 30 K, because CH_{3} desorbs at larger temperatures.
 2.
Our computations suggest that acetaldehyde is not the dominant product for the reaction HCO + CH_{3}. 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_{4} channel for f ≥ 0.4, while only the efficiency of CO + CH_{4} 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_{4} 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 gasphase 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 SklodowskaCurie for the project “AstroChemical Origins” (ACO), grant agreement no. 811312. AR is indebted to “Ramón y Cajal” program. MINECO (project CTQ201789132P) 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.univgrenoblealpes.fr), which is partly supported by the Equip@Meso project (reference ANR10EQPX2901) 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 2019A0060810797 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 BHLYPD3(BJ) as an accurate, costeffective method for the reactions studied in this work is shown in this section.
We have taken five hybrid DFT dispersioncorrected methods: MPWB1KD3(BJ), M062XD3, PW6B95D3(BJ), wB97XD3 and BHandHLYPD3(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 BHLYPD3(BJ)/6311++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 6311++G(2df,2pd) basis set, and CASPT2/augccPVTZ 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 HartreeFock wavefunction.
BHLYPD3(BJ) gives the best overall performance with and average unsigned error of 3.0%, and a maximum of 5.0% with respect to CASPT2/augccPVTZ. 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.
DFT method benchmark results. % Unsigned Error from CASPT2/augccPVTZ//BHLYPD3(BJ)/6311++G(d,p). Methods are: (1) BHLYPD3(BJ), (2) MPWB1KD3(BJ), (3) M062XD3, (4) PW6B95D3(BJ) and (5) wB97xD3.
Energetics (kJ mol^{−1}) of the stationary points optimized at BHLYPD3(BJ)/6311++G(d,p), together with the frequencies of the transition states (cm^{−1}).
Raw data from the DFT method benchmark. Energies in kJ mol^{−1}.
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 BHLYPD3(BJ)/6311++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 preexponential factor) and then the exponential with the binding energies. We define the attempt frequency within the Eyring equation as (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): (B.2)
where with c the speed of light and the ith vibrational mode frequency in cm^{−1}.
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: (B.3) (B.4)
where σ_{rot} is the rotational symmetry number (6 for CH_{3} and 1 for HCO), θ_{r,i} = B_{i}h∕k_{B} with B_{i} the ith axis rotational constant (in s^{−1}), m is the mass of the radical and P the gas pressure, which was calculated assuming a density of 10^{4} cm^{−3}.
The thermal corrections follow: H = E_{DFT} + ZPE + E_{vib}(T) + k_{B}T for the surface and the surface + radical complex. Here E_{DFT} is the energy obtained from our DFT calculations, ZPE is the zeropoint 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_{B}T + H_{rot} + H_{trans} where . (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:
Corrections incorporated in the calculation of ΔG of binding. Quantities calculated at 20 K and 30 K. Energy and entropy units are kJ mol − 1 and kJ mol − 1 K^{−1}. See that Boltzmann’s constant is about 0.0083 kJ mol^{−1}, therefore thek_{B}T and 3k_{B}T/2 terms are also very 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^{8} –10^{11} s^{−1}, smaller than the normally used approach of the harmonic oscillator^{3} 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.
Attempt frequencies for desorption and the different cases of diffusion considered in this work. Units in s^{−1}. See that at 20 K k_{B}T∕h ~ 4.2 × 10^{11}, so the effect of entropy at such low temperatures is very small.
Appendix C Rate constant comparison
Figure C.1 compares the rates of the radicalradical 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 CH_{3} + 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 diffusiontodesorption 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 UBHandHLYPD3(BJ)/631+G(d,p) level and DFT energies were refined at UBHandHLYPD3(BJ)/6311++G(2df,2pd) level. Reactants and products were obtained by running intrinsic reaction coordinate calculations. 
H + CO → HCO energetic data, in Hartree, at UBHandHLYPD3(BJ)/631+G(d,p) (double ζ) level. U is the DFT energy, D is the Dispersion energy and ZPE is the zeropoint energy. DFT energies were refined by performing single point calculations on double ζ geometries at UBHandHLYPD3(BJ)/6311++G(2df,2pd) (triple ζ) level. Energy units are Hartree (1.0 Hartree are ~ 2625.5 kJ mol^{−1}).
H⋯ CO transition state data. Units are the usual for a GAUSSIAN16 output: frequencies in cm^{−1}, IR intensities in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYPD3(BJ)/631+G(d,p) level.
Appendix E Radicalradical TS data and PES
CH_{3} + HCO energetic data, in Hartree, at UBHandHLYPD3(BJ)/631+G(d,p) level (double ζ). U is the DFT energy, D is the Dispersion energy and ZPE is the zeropoint energy. DFT energies were refined by performing single point calculations on double ζ geometries at UBHandHLYPD3(BJ)/6311++G(2df,2pd) (triple ζ) level. Energy units are Hartree (1.0 Hartree are ~ 2625.5 kJ mol^{−1}).
Features of the transition states studied in this work. Units are the usual for a GAUSSIAN16 output: frequencies in cm^{−1}, IR intensities in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYPD3(BJ)/631+G(d,p) level.
Fig. E.1
Potential energy surfaces Reacts. I and II. Geometries and ZPE energy correction was obtained at UBHandHLYPD3(BJ)/631+G(d,p) level, DFT energy was refined at UBHandHLYPD3(BJ)/6311++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 greencolored 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
Halflives are calculated from the rate constants (k_{i}, with i being either the diffusion or desorption of radicals) following
and shown in Fig. G.1.
Fig. G.1
Diffusion and diffusion temperatures of CH_{3} and HCO assuming a halflife 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 E_{diff} ∕E_{des} = 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 + CH_{4} formation using E_{diff}∕E_{des} = 0.5, 0.4 and 0.3 (left to right panels). 
References
 Aikawa, Y., Furuya, K., Yamamoto, S., & Sakai, N. 2020, ApJ, 897, 110 [Google Scholar]
 ÁlvarezBarcia, S., Russ, P., Kästner, J., & Lamberts, T. 2018, MNRAS, 479, 2007 [Google Scholar]
 Andersson, S., Goumans, T. P. M., & Arnaldsson, A. 2011, Chem. Phys. Lett., 513, 31 [NASA ADS] [CrossRef] [Google Scholar]
 BarNun, A., Dror, J., Kochavi, E., & Laufer, D. 1987, Phys. Rev. B, 35, 2427 [NASA ADS] [CrossRef] [Google Scholar]
 Barger, C. J., & Garrod, R. T. 2020, ApJ, 888, 38 [Google Scholar]
 Becke, A. D. 1993, J. Chem. Phys., 98, 1372 [Google Scholar]
 Bennett, C. J., Jamieson, C. S., Osamura, Y., & Kaiser, R. I. 2005, ApJ, 624, 1097 [Google Scholar]
 Bianchi, E., Codella, C., Ceccarelli, C., et al. 2019, MNRAS, 483, 1850 [Google Scholar]
 Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621 [Google Scholar]
 Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541 [Google Scholar]
 Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 593, L51 [Google Scholar]
 Ceccarelli, C., Caselli, P., Fontani, F., et al. 2017, ApJ, 850, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, Q., Cuppen, H. M., & Herbst, E. 2007, A&A, 469, 973 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charnley, S. B. 2004, Adv. Space Res., 33, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Chevance, M., Kruijssen, J. M. D., VazquezSemadeni, E., et al. 2020, Space Sci. Rev., 216, 50 [Google Scholar]
 Cooke, I. R., Öberg, K. I., Fayolle, E. C., Peeler, Z., & Bergner, J. B. 2018, ApJ, 852, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Csengeri, T., Belloche, A., Bontemps, S., et al. 2019, A&A, 632, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 [Google Scholar]
 De Simone, M., Codella, C., Ceccarelli, C., et al. 2020, A&A, 640, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, A&A, 512, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 EnriqueRomero, J., Rimola, A., Ceccarelli, C., et al. 2019, ACS Earth Space Chem., 3, 2158 [NASA ADS] [CrossRef] [Google Scholar]
 EnriqueRomero, J., ÁlvarezBarcia, S., Kolb, F. J., et al. 2020, MNRAS, 493, 2523 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian’16 Revision C.01 (Wallingford, CT: Gaussian Inc.) [Google Scholar]
 Garrod, R. T. 2013, ApJ, 778, 158 [NASA ADS] [CrossRef] [Google Scholar]
 Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15 [Google Scholar]
 Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283 [Google Scholar]
 Garrod, R. T., Vasyunin, A. I., Semenov, D. A., Wiebe, D. S., & Henning, T. 2009, ApJ, 700, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Ghesquière, P., Mineva, T., Talbi, D., et al. 2015, Phys. Chem. Chem. Phys., 17, 11455 [CrossRef] [Google Scholar]
 Goerigk, L., Hansen, A., Bauer, C., et al. 2017, Phys. Chem. Chem. Phys., 19, 32184 [CrossRef] [Google Scholar]
 Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, J. Chem. Phys., 132, 154104 [Google Scholar]
 Grimme, S., Ehrlich, S., & Goerigk, L. 2011, J. Comput. Chem., 32, 1456 [Google Scholar]
 Hama, T., & Watanabe, N. 2013, Chem. Rev., 113, 8783 [Google Scholar]
 Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [Google Scholar]
 He, J., & Vidali, G. 2014, ApJ, 788, 50 [Google Scholar]
 He, J., Emtiaz, S. M., & Vidali, G. 2017, ApJ, 837, 65 [Google Scholar]
 He, J., Emtiaz, S., & Vidali, G. 2018, ApJ, 863, 156 [CrossRef] [Google Scholar]
 Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Hollenbach, D., & Salpeter, E. 1970, J. Chem. Phys., 53, 79 [CrossRef] [Google Scholar]
 Hudson, R. L., & Moore, M. H. 1997, Icarus, 126, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Isokoski, K., Bossa, J. B., Triemstra, T., & Linnartz, H. 2014, Phys. Chem. Chem. Phys., 16, 3456 [NASA ADS] [CrossRef] [Google Scholar]
 Jensen, S. S., Jørgensen, J. K., Furuya, K., Haugbølle, T., & Aikawa, Y. 2021, A&A, 649, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, A117 [Google Scholar]
 Karssemeijer, L. J., & Cuppen, H. M. 2014, A&A, 569, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keane, J. V., Tielens, A. G. G. M., Boogert, A. C. A., Schutte, W. A., & Whittet, D. C. B. 2001, A&A, 376, 254 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kouchi, A., Furuya, K., Hama, T., et al. 2020, ApJ, 891, A22 [Google Scholar]
 Lamberts, T., & Kästner, J. 2017, ApJ, 846, 43 [Google Scholar]
 Lamberts, T., Markmeyer, M. N., Kolb, F. J., & Kästner, J. 2019, ACS Earth Space Chem., 3, 958 [Google Scholar]
 Lauck, T., Karssemeijer, L., Shulenberger, K., et al. 2015, ApJ, 801, 118 [Google Scholar]
 Lee, C., Yang, W., & Parr, R. G. 1988, Phys. Rev. B, 37, 785 [Google Scholar]
 Lee, C.F., Codella, C., Li, Z.Y., & Liu, S.Y. 2019, ApJ, 876, 63 [Google Scholar]
 Lefloch, B., Ceccarelli, C., Codella, C., et al. 2017, MNRAS, 469, L73 [Google Scholar]
 MartínDoménech, R., Öberg, K. I., & Rajappan, M. 2020, ApJ, 894, 98 [Google Scholar]
 Maté, B., Cazaux, S., Satorre, M. Á., et al. 2020, A&A, 643, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minissale, M., Congiu, E., & Dulieu, F. 2016, A&A, 585, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mispelaer, F., Theulé, P., Aouididi, H., et al. 2013, A&A, 555, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Molinari, S., Ceccarelli, C., White, G. J., et al. 1999, ApJ, 521, L71 [CrossRef] [Google Scholar]
 Molpeceres, G., Rimola, A., Ceccarelli, C., et al. 2018, MNRAS, 482, 5389 [Google Scholar]
 Neese, F. 2004, J. Phys. Chem. Solids, 65, 781 [NASA ADS] [CrossRef] [Google Scholar]
 Oba, Y., Miyauchi, N., Hidaka, H., et al. 2009, ApJ, 701, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Öberg, K. I., van Dishoeck, E. F., Linnartz, H., & Andersson, S. 2010, ApJ, 718, 832 [Google Scholar]
 Palumbo, M. E. 2006, A&A, 453, 903 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Palumbo, M. E., Baratta, G. A., Leto, G., & Strazzulla, G. 2010, J. Mol. Struct., 972, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Pantaleone, S., EnriqueRomero, J., Ceccarelli, C., et al. 2020, ApJ, 897, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Pantaleone, S., EnriqueRomero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49 [CrossRef] [Google Scholar]
 Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71 [Google Scholar]
 Potapov, A., & McCoustra, M. 2021, Int. Rev. Phys. Chem., 40, 299 [CrossRef] [Google Scholar]
 Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, 572, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rimola, A., Skouteris, D., Balucani, N., et al. 2018, ACS Earth Space Chem., 2, 720 [Google Scholar]
 Ruaud, M., Loison, J. C., Hickson, K. M., et al. 2015, MNRAS, 447, 4004 [Google Scholar]
 Ruffle, D. P., & Herbst, E. 2000, MNRAS, 319, 837 [Google Scholar]
 Sakai, T., Yanagida, T., Furuya, K., et al. 2018, ApJ, 857, 35 [Google Scholar]
 Scibelli, S., & Shirley, Y. 2020, ApJ, 891, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Skouteris, D., Balucani, N., Ceccarelli, C., et al. 2018, ApJ, 854, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1989, ApJ, 344, 413 [CrossRef] [Google Scholar]
 Taquet, V., Ceccarelli, C., & Kahane, C. 2012, A&A, 538, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tielens, A. G.G. M., & Allamandola, L. J. 1987, in Interstellar Processes, eds. D. J. Hollenbach, & H. A. Thronson (Dordrecht: Springer Netherlands), 397 [NASA ADS] [CrossRef] [Google Scholar]
 Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 [Google Scholar]
 Vastel, C., Ceccarelli, C., Lefloch, B., & Bachiller, R. 2014, ApJ, 795, L2 [Google Scholar]
 Vasyunin, A. I., Caselli, P., Dulieu, F., & JiménezSerra, I. 2017, ApJ, 842, 33 [Google Scholar]
 Vazart, F., Ceccarelli, C., Balucani, N., Bianchi, E., & Skouteris, D. 2020, MNRAS, 499, 5547 [Google Scholar]
 Vidali, G. 2013, Chem. Rev., 113, 8752 [Google Scholar]
 Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1 [Google Scholar]
 Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173 [Google Scholar]
 Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427 [NASA ADS] [CrossRef] [Google Scholar]
 Zhitnikov, R. A., & Dmitriev, Y. A. 2002, A&A, 386, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Energetics and related parameters of the reactions and desorption and diffusion of the radicals.
Rate constants k_{aeb} (in s^{−1}) and efficiency ε of the two possible reactions between HCO and CH_{3}.
DFT method benchmark results. % Unsigned Error from CASPT2/augccPVTZ//BHLYPD3(BJ)/6311++G(d,p). Methods are: (1) BHLYPD3(BJ), (2) MPWB1KD3(BJ), (3) M062XD3, (4) PW6B95D3(BJ) and (5) wB97xD3.
Energetics (kJ mol^{−1}) of the stationary points optimized at BHLYPD3(BJ)/6311++G(d,p), together with the frequencies of the transition states (cm^{−1}).
Corrections incorporated in the calculation of ΔG of binding. Quantities calculated at 20 K and 30 K. Energy and entropy units are kJ mol − 1 and kJ mol − 1 K^{−1}. See that Boltzmann’s constant is about 0.0083 kJ mol^{−1}, therefore thek_{B}T and 3k_{B}T/2 terms are also very small.
Attempt frequencies for desorption and the different cases of diffusion considered in this work. Units in s^{−1}. See that at 20 K k_{B}T∕h ~ 4.2 × 10^{11}, so the effect of entropy at such low temperatures is very small.
H + CO → HCO energetic data, in Hartree, at UBHandHLYPD3(BJ)/631+G(d,p) (double ζ) level. U is the DFT energy, D is the Dispersion energy and ZPE is the zeropoint energy. DFT energies were refined by performing single point calculations on double ζ geometries at UBHandHLYPD3(BJ)/6311++G(2df,2pd) (triple ζ) level. Energy units are Hartree (1.0 Hartree are ~ 2625.5 kJ mol^{−1}).
H⋯ CO transition state data. Units are the usual for a GAUSSIAN16 output: frequencies in cm^{−1}, IR intensities in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYPD3(BJ)/631+G(d,p) level.
CH_{3} + HCO energetic data, in Hartree, at UBHandHLYPD3(BJ)/631+G(d,p) level (double ζ). U is the DFT energy, D is the Dispersion energy and ZPE is the zeropoint energy. DFT energies were refined by performing single point calculations on double ζ geometries at UBHandHLYPD3(BJ)/6311++G(2df,2pd) (triple ζ) level. Energy units are Hartree (1.0 Hartree are ~ 2625.5 kJ mol^{−1}).
Features of the transition states studied in this work. Units are the usual for a GAUSSIAN16 output: frequencies in cm^{−1}, IR intensities in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYPD3(BJ)/631+G(d,p) level.
All Figures
Fig. 1
Arrhenius plots, namely rate constants as a function of the inverse of temperature, for the reaction CH_{3} + HCO forming acetaldehyde (black solid line) or CO + CH_{4} (black dashed line), and for the reaction H + CO → HCO (gray dotteddashed line), described in the main text. 

In the text 
Fig. 2
Reaction efficiency ε (Eq. (3)) of the reaction CH_{3} + HCO leading to either CH_{3}CHO (solid lines) or CO + CH_{4} (dashed lines) as a function of the temperature. The computations were obtained adopting three different E_{diff} /E_{des} ratios: 0.3 (green), 0.4 (blue) and 0.5 (red). We note that, for E_{diff}/E_{des} = 0.5 the CH_{3}CHO and CO + CH_{4} (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 CH_{3}CHO over CO + CH_{4} (Eq. (10)) as a function of the temperature in the range where the reactions can occur, namely below 30 K (see text), for E_{diff}∕E_{des} 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 BHLYPD3(BJ)/6311++G(d,p) level. Distances in Å. 

In the text 
Fig. C.1
Comparison of reaction, diffusion and desorption rate constants involved in the CH_{3} + 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 diffusiontodesorption 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 UBHandHLYPD3(BJ)/631+G(d,p) level and DFT energies were refined at UBHandHLYPD3(BJ)/6311++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 UBHandHLYPD3(BJ)/631+G(d,p) level, DFT energy was refined at UBHandHLYPD3(BJ)/6311++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 greencolored 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 CH_{3} and HCO assuming a halflife 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 E_{diff} ∕E_{des} = 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 + CH_{4} formation using E_{diff}∕E_{des} = 0.5, 0.4 and 0.3 (left to right panels). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.