A&A 655, A9 (2021) https://doi.org/10.1051/0004-6361/202141531 Astronomy & © J. Enrique-Romero et al. 2021 Astrophysics Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains Joan Enrique-Romero1,2 , Cecilia Ceccarelli1 , Albert Rimola2 , Dimitrios Skouteris3 , Nadia Balucani1,4,5 , and Piero Ugliengo6 1 Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France e-mail: juan.enrique-romero@univ-grenoble-alpes.fr; cecilia.ceccarelli@univ-grenoble-alpes.fr 2 Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, 08193, Catalonia, Spain e-mail: albert.rimola@uab.cat 3 Master-Tech, 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 ABSTRACT Context. Interstellar grains are known to be important actors in the formation of interstellar molecules such as H2 , water, ammonia, and methanol. It has been suggested that the so-called interstellar complex organic molecules (iCOMs) are also formed on the inter- stellar 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 acetalde- hyde (CH3 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 CH3 . 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 CH3 + HCO, which can lead to the formation of CH3 CHO or CO + CH4 . Specifically, we combined reaction kinetics computed with the Rice–Ramsperger- –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 CH3 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 + CH4 . 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 CH3 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 CH3 , 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 micro-physics 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 1. Introduction Rimola et al. 2014), which are also abundant molecules predom- inantly synthesized on grain surfaces. As a matter of fact, in Interstellar dust grains are known to be an essential component cold regions, the refractory core of the grains, made up of sil- of the interstellar medium (ISM) for a large variety of reasons. icate or carbonaceous material, are coated by icy mantles mostly Among them, dust grains provide the surfaces for chemical reac- formed by amorphous water ice synthesized on these surfaces tions that are difficult (or impossible) to take place in the gas (e.g., Boogert et al. 2015). phase. An emblematic example is the formation of the most It has also been suggested that more complex molecules, abundant molecule in the Universe, H2 , which largely occurs on the so-called interstellar complex organic molecules (hereinafter grain surfaces (e.g., Hollenbach & Salpeter 1970; Vidali 2013; iCOMs: Ceccarelli et al. 2017), composed of at least six atoms Wakelam et al. 2017). Other important examples are the for- and containing at least one heavy element other than C, can mation of water (e.g., Dulieu et al. 2010; He & Vidali 2014; also be a grain-surface chemistry product (e.g., Garrod & Herbst Lamberts & Kästner 2017; Molpeceres et al. 2018) and methanol 2006; Herbst & van Dishoeck 2009; Garrod et al. 2009; Ruaud (e.g., Tielens & Hagen 1982; Watanabe & Kouchi 2002; et al. 2015; Aikawa et al. 2020; Barger & Garrod 2020). One A9, page 1 of 14 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. A&A 655, A9 (2021) crucial step of this theory is the formation of iCOMs from expressed as Ri j = ki j ni n j , where ni and n j are the densities of the combination of two radicals when they meet on the grain species i and j, and ki j is the rate constant at a given tempera- icy surfaces. In the majority of the current astrochemical mod- ture. For surface reactions the latter is given by (Hasegawa et al. els, the reaction is assumed to proceed barrierless and without 1992): competitive channels. However, previous theoretical works have shown that this is Rdiff,i + Rdiff, j ki j = εi j × , (1) not necessarily the case (Rimola et al. 2018; Enrique-Romero nd et al. 2019, 2020). For example, theoretical calculations on the energetics showed that the formation of acetaldehyde (CH3 CHO) where εi j is an efficiency factor which accounts for chemical bar- on the icy surfaces via HCO + CH3 is in competition with the riers, nd is the dust grain density and Rdiff,i and Rdiff, j are the formation of CO + CH4 via direct H-abstraction. In addition, diffusion rates for species i and j, respectively. These diffusion both reactions present barriers, caused by the orientation of the rates are defined as 1/tdiff,k , where tdiff,k is the time it takes the species on the ices, which are governed by the interactions cre- species k to scan the whole grain (e.g., Garrod & Herbst 2006). ated between the surface water molecules and the two radicals. Thus, the sum Rdiff,i + Rdiff, j gives the rate at which species i and Indeed, the height of the barriers depends on the site where the j meet on the surface. reaction occurs, whether the two radicals are on a “plain” ice Regarding the efficiency factor εi j , different approaches exist surface or in a “cavity”, namely on the interactions between the in order to derive it. Hasegawa et al. (1992) set it to either 1, radicals and the ice water molecules. in barrierless reactions, or to the tunnelling probability, if the In the present work, we pursue the above theoretical stud- reaction has an activation energy barrier and one of the reac- ies and present new computations to evaluate the efficiency of tants is light enough to tunnel through it. Later models include the radical–radical combination and H-abstraction reactions as a also the thermal probability for reaction, if there is an activa- function of the temperature of the icy surfaces with the goal to tion energy barrier (e.g., Garrod & Herbst 2006). However, in provide values that can be easily incorporated in astrochemical the presence of an activation energy barrier, reactants need to models. In particular, here we focus on acetaldehyde (CH3 CHO), be close to each other for a certain amount of time for the reac- one of the most abundant and common iCOMs (e.g., Blake et al. tion to occur (Tielens & Hagen 1982). In order to take this into 1987; Cazaux et al. 2003; Vastel et al. 2014; Lefloch et al. 2017; account, Chang et al. (2007) redefined the efficiency taking into Sakai et al. 2018; Bianchi et al. 2019; Csengeri et al. 2019; Lee account the competition between diffusion and desorption of the et al. 2019; Scibelli & Shirley 2020), as a study case. most mobile species, as follows: Following the works of Enrique-Romero et al. (2019, 2020), kaeb (i j) we here consider the two competing reactions that can arise εi j = , (2) from the reactivity between HCO and CH3 : kaeb (i j) + kdiff (i) + kdes (i) (1) HCO + CH3 → CH3 HCO where kaeb (i j) is the rate constant accounting for the reaction (2) HCO + CH3 → CH4 + CO. activation energy barrier, which is described by either classical thermal kinetics or quantum tunnelling; (i.e., a frequency times a Our goal is to evaluate the efficiency of each of the two Boltzmann factor or the tunnelling probability); kdiff (i) is the rate reactions occurring on the interstellar icy surfaces. To this end, constant for the diffusion of the most mobile species and kdes (i) we computed the kinetics of the two reactions, using the previ- is its desorption rate constant. Garrod & Pauly (2011) further ous energetic calculations by Enrique-Romero et al. (2019) as a modified Eq. (2) by removing the desorption term and adding base and the models of amorphous solid water (ASW) for the the diffusion of the other reaction partner, j, in the denominator. ice described in Rimola et al. (2018) and Enrique-Romero et al. For reactions involving radicals, εi j is normally assumed (2019). In our calculations, we assume that the two radicals are equal to 1 (e.g., Garrod & Herbst 2006), as they are considered in the most stable energetic configuration prior to reaction, an to react via barrierless exothermic channels. assumption motivated by the long surface residence timescale In this work, we include diffusion and desorption rates of (1–10 Myr, the molecular life timescale: e.g. Chevance et al. the two reactants, which takes into account both the Chang et al. 2020) that radicals would experience before they become mobile (2007) and Garrod & Pauly (2011) recipes: and react with each other. The article is organized as follows. The definition of the kaeb (i j) εi j = . (3) reaction efficiency and our choices for the various assumptions kaeb (i j) + kdiff (i) + kdes (i) + kdiff ( j) + kdes ( j) entering in the computations are discussed in Sect. 2. Section 3 describes the adopted methodology. The results are reported in In practice, the efficiency for the reaction is equal to unity only Sect. 4 and we discuss the implications of our new calculations when the time scale for the reaction to occur (1/kaeb ) is shorter in Sect. 5. than the timescales at which reactants remain on the reaction site (the smallest between 1/kdiff (i) and 1/kdiff (j)). 2. Efficiency of radical–radical reaction products 2.2. A novel treatment of surface radical–radical reactions on icy surfaces rate constants 2.1. Surface-reaction rate definition The novelty of the present work is the estimate of the kaeb (i j) coefficient of radical–radical reactions via statistical kinet- Generally, astrochemical models solve the time-dependent equa- ics calculations based on the Ramsperger–Rice–Kassel–Marcus tions of the species densities by computing the formation and (RRKM) microcanonical transition state theory. Briefly, RRKM destruction rates of each species at a given time, both for species computations provide unimolecular rate constants, namely the in the gas and on the grain surfaces. In particular, the rate Ri j rate at which a system A becomes A0 passing through a transi- of the formation reaction from two reactant species i and j is tion state (only once). In our case, the system A is the ice-water A9, page 2 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains molecules plus the two adsorbed radicals, namely we con- and sider the water-cluster plus the radicals as a super-molecule isolated from its surrounding. The system A0 is the product ∆S = ∆S vib + ∆S rot + ∆S trans , (7) of the radical–radical reaction on the icy surface, namely the water-cluster plus either the radical–radical recombination (e.g., where the terms are the energy difference between the des- React. I) or the H-abstraction (e.g., React. II) products. orbed and the adsorbed states for the total electronic energy It is important to note that in order to apply the RRKM (∆Eelectronic ), for the zero point vibrational energy correc- theory, we implicitly assume that the intra-molecular energy tions (∆ZPE), for the thermal vibrational energy corrections redistribution of the reaction energy is faster than the reaction (∆Evib (T )), for the vibrational entropy (∆S vib ), and for the rota- itself. This assumption is supported by recent ab initio molecular tional and translational contributions to enthalpy (∆Hrot and dynamics (AIMD) computations that show that a large fraction ∆Htrans , respectively) and entropy (∆S rot and ∆S trans , respec- (≥50%) of the reaction energy is absorbed by the water ice in less tively). In this case, since we are dealing with the desorption of than 1 ps (Pantaleone et al. 2020, 2021). We have checked a pos- the radicals, the translational and rotational contributions arise teriori that the timescale of the reactions studied here is indeed from only the desorbed (free) radicals. Specific details on the longer that 1 ps. Finally, the specific computational details of our calculation of some of these terms are provided in Appendix B. proposed RRKM method are reported in Sect. 3.2. For the sake of simplicity, we refer to this final desorption energy as Edes . 2.3. Desorption and diffusion energies In contrast to desorption, obtaining diffusion-related data with the present calculations is a daunting task, as it requires Equation (3) shows that, in addition to the probability kaeb (i j) for localizing a large number of transition states for the radical hop- radicals i and j to react when they meet on a surface site, the ing between the different binding sites. Moreover, the use of a −1 −1 efficiency factor εi j also depends on kdes and on kdiff , which are relatively small cluster model dramatically constraints the valid- related to the residence time of the radicals on the surface and ity of these results because of its limitation in terms of size and on the diffusion timescale of the radicals on the ice, respectively. surface morphology. Therefore, to obtain a value for the diffu- The diffusion and desorption timescales tdiff/des are given by the sion energy of each radicals, which by analogy we will refer to as classical Eyring transition state theory (TST), in which tdiff/des Ediff , we resorted to what is usually done in astrochemical mod- −1 is inversely proportional to their rate constant, tdiff/des ∝ kdiff/des . eling, that is, Ediff is taken to be a fraction f of Edes . However, According to TST, the general expression of the rate constant k deriving the value of f has proven to be difficult, both theoreti- for a unimolecular reaction (like diffusion and desorption) is: cally and experimentally. In the published astrochemical models, one can find a quite wide range of adopted f values, from 0.3 to kB T Q, ∼0.8 (e.g., Hasegawa et al. 1992; Ruffle & Herbst 2000). Some k= exp(−∆V , /kB T ), (4) h QR 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; where ∆V , is the zero point energy-corrected energy barrier, Ruaud et al. 2015; Vasyunin et al. 2017; Jensen et al. 2021). Q, and QR are the total partition functions of the transition In the past few years, theoretical and experimental works state and the initial state (namely, the reactants), respectively, on the diffusion process of species on ASW surfaces have pro- kB is Boltzmann’s constant, T is the surface temperature and h vided constraints to the f value (see also the more extensive is Planck’s constant. We note that we use the classical Eyring discussion in Sect. 5.3). In a theoretical work, Karssemeijer & equation because, since we are not dealing with light atoms but Cuppen (2014) showed that the range for the Ediff /Edes ratio can molecular radicals, tunneling is negligible. be narrowed down to 0.3–0.4 for molecules like CO and CO2 . By proper manipulation of Eq. (4), the rate constant becomes Minissale et al. (2016) experimentally found that the f ratio of expressed as a function of the free energy barrier ∆G, (usually atomic species like N and O is about 0.55, while He et al. (2018) referred to as free energy of activation) at a given temperature: showed that f is 0.3–0.6, being the lower values more suitable kB T for surface coverage lower than one mono-layer. k= exp(−∆G, /kB T ), (5) Given the uncertainty on the Ediff /Edes ratio for CH3 and h HCO, we carried out our calculations for three values: 0.3, 0.4 in which ∆G, = ∆H , − T∆S , and where ∆H , is the and 0.5. enthalpy of activation and ∆S , the entropy of activation. These terms contain translational, rotational, vibrational and elec- 2.4. Ice model tronic contributions as they arise partly from the total partition functions Q. Regarding the amorphous solid water (ASW) model, there is still With the adopted quantum chemical approach, the applica- little knowledge that constrains the actual internal structure of tion of the Eyring TST allows us to compute desorption-related interstellar ices. Observations suggest that the interstellar water data (e.g., desorption activation energies and desorption rate ice is predominantly in the amorphous form (e.g., Smith et al. constants) for each radical through the outcome of these cal- 1989; Boogert et al. 2015) (with some exceptions: e.g. Molinari culations (electronic energies, vibrational frequencies, partition et al. 1999). Many laboratory studies have been carried out to functions, energy contributions, etc.). It is worth mentioning characterize the possible porosity of the interstellar ices. Typ- that, since the radicals are physisorbed on the ice surfaces, the ically, laboratory experiments produce porous ices of different energy barriers of the desorption processes coincide with the densities by condensation of water vapor, even though they prob- desorption energies. In the present case, we only account for the ably do not reproduce the interstellar water ice, in which water is electronic and vibrational contributions to both ∆H , and ∆S , to believed to form in situ by hydrogenation reactions of frozen O, arrive at the radical desorption energies as follows: O2 and O3 (e.g., Dulieu et al. 2010; Hama & Watanabe 2013; He & Vidali 2014; Potapov & McCoustra 2021). In general, porous ∆H = ∆Eelectronic + ∆ZPE + ∆Evib (T ) + ∆Hrot + ∆Htrans , (6) ices are detected in laboratory via the infrared (IR) signature of A9, page 3 of 14 A&A 655, A9 (2021) dangling OH groups, which are, however, missing in interstellar fundamental aspect when aiming at providing kinetic calcula- samples (Bar-Nun et al. 1987; Keane et al. 2001), (see also the tions and rate constants (including tunneling effects) (Álvarez- discussion in e.g. Hama & Watanabe 2013; Zamirri et al. 2018). Barcia et al. 2018), as in the present work. Additionally, as shown Several hypothesis have been suggested to explain the absence in Rimola et al. (2018); Enrique-Romero et al. (2019, 2020), DFT of the OH dangling signature (Oba et al. 2009; Palumbo 2006; is a cost-effective methodology with which a correct description Palumbo et al. 2010), so that, at the end, there is consensus in the of biradical systems can be achieved by using the unrestricted community that interstellar water ices are amorphous and porous broken (spin)-symmetry approach (e.g., Neese 2004). in nature, even though many details are missing and we do not All optimized stationary points were characterized by the have a precise picture of the degree of porosity (e.g., Hama analytical calculation of the harmonic frequencies as min- & Watanabe 2013; Isokoski et al. 2014; Potapov & McCoustra ima and saddle points. Thermochemical corrections computed 2021). at BHLYP-D3(BJ)/6-31+G(d,p) were included to the single In order to simulate the interstellar icy surfaces, Enrique- point BHLYP-D3(BJ)/6-311++G(2df,2pd) potential energy val- Romero et al. (2019) considered a cluster of 33 water molecules. ues using the standard rigid-rotor and harmonic oscillator for- This ice model possesses two major types of surface with respect mulae in order to obtain the zero-point vibrational energy (ZPE) to the binding capability: a cavity, where species are in gen- corrections. eral more strongly bonded to the surface, and an elongated side (Rimola et al. 2014, 2018; Enrique-Romero et al. 2019). In this 3.2. Kinetic calculations work, we only report the analysis of the reaction occurring in the cavity for the following reason. In astrochemical models, the vast In order to compute the rate constants for the chemical reactions majority of radical–radical reactions take place inside the bulk between the radical pairs, we adapted our in-house kinetic code, of the ice (e.g., Garrod & Herbst 2006). Therefore, the cavity based on the RRKM scheme for gas-phase reactions (Skouteris site is a better representation of the sites where radical–radical et al. 2018), to the surface plus adsorbed radicals case. First, reactions occur than that on the elongated side, which would at we obtained the microcanonical rate constant kaeb (E) at a given best describe the ice layer exposed to the gas and where just a energy E as: tiny fraction of the reactions can occur, considering that the ice is constituted by more than 100 layers (e.g., Taquet et al. 2012; N(E) Aikawa et al. 2020). kaeb (E) = , (8) hρ(E) Therefore, in this work, we use the Enrique-Romero et al. (2019) ice model and methodology, but we improve the calcula- where N(E) is the sum of states for the active degrees of free- tions for a better accuracy of the computed energetics, including dom in the transition state, ρ(E) is the density of states for the dispersion, as described in detail in Sect. 3. 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 3. Methodology account. Second, the obtained rate constants were Boltzmann- 3.1. Electronic structure calculations averaged in order to derive the rate constants as a function of the temperature. Given the importance of inter-molecular interactions in radical– For the H abstraction reaction, we took into account tun- radical reactions, we recomputed the stationary points of the neling effects adopting the Eckart scheme via the unsymmetric potential energy surfaces (PES) previously reported by Enrique- potential energy barrier approach. In order to have a chemi- Romero et al. (2019) using the Grimme’s D3 dispersion term cal system of reference to compare with, we applied the same including the Becke–Johnson damping (D3(BJ)) (Grimme et al. method to the well studied reaction H + CO → HCO. In this case, 2010, 2011), this way improving the description of the dispersion the initial structures of the reaction were taken from the theoret- forces with respect to the previous work. ical study by Rimola et al. (2014), which were re-optimized at All DFT calculations were performed with the GAUSSIAN 16 the present work computational level. Here, from the optimized program package (Frisch et al. 2016). A benchmark study showed transition state, intrinsic reaction coordinate (IRC) calculations that the BHLYP hybrid density functional method is the best were run assuming a Langmuir–Hinshelwood (LH) like reac- suited DFT method to study these reactions, with an average tion, contrarily to the Rimola et al. (2014) original computations. error of 3%, and a maximum error of 5.0% with respect to bench- All stationary points were characterized by frequency calcula- mark multi-reference CASPT2(2,2) calculations using O PEN - tions, obtaining their (harmonic) vibrational modes and their M OLCAS 18.09 (see Annex). We have chosen CASPT2(2,2) as zero-point energies. More details of these computations can be the minimum level of post-HF theory, aware of the fact that a found in Appendix D. 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 4. Results (Becke 1993; Lee et al. 1988) combined with the standard 6- 4.1. Energetics of the reactions 31+G(d,p) Pople basis set alongside the D3(BJ) dispersion term (Grimme et al. 2010, 2011). When needed, intrinsic reaction Table 1 presents the 0 K enthalpies (i.e., potential energies plus coordinate (IRC) calculations at the optimization theory level ZPE corrections) of the studied reactions. The improvement in were carried out to ensure that the transition states connect with the dispersion correction and the refinement of the DFT energy the corresponding minima. To balance the computational cost slightly decrease the energy barriers of each one of the reactions and chemical accuracy, reaction energetics were then refined to form acetaldehyde and CO + CH4 by less than 2.5 kJ mol−1 by performing full BHLYP-D3(BJ)/6-311++G(2df,2pd) single- with respect to the values quoted by Enrique-Romero et al. point energy calculations on the BHLYP-D3(BJ)/6-31+G(d,p) (2019). Inversely, the H + CO → HCO reaction has a higher optimized stationary points. Improving chemical accuracy is a barrier, 13.5 kJ mol−1 , than that quoted by Rimola et al. (2014), A9, page 4 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains Table 1. Energetics and related parameters of the reactions and desorp- tion and diffusion of the radicals. Product ∆H ‡ ∆H RX CH3 CHO 5.5 –324.5 CO + CH4 7.2 –328.9 HCO 13.5 –91.6 Quantity (K) CH3 HCO Edes 1715 3535 T des 30 68 T diff (0.5) 15 32 T diff (0.4) 12 25 T diff (0.3) 9 19 Notes. Top half: activation (∆H ‡ ) and reaction (∆H RX ) enthalpies (in kJ mol−1 ) at 0 K (i.e., sum of electronic energies at BHLYP-D3(BJ)/6- 311+G(2df,2pd)//BHLYP-D3(BJ)/6-31+G(d,p) and ZPE at BHLYP- D3(BJ)/6-31+G(d,p)) for each radical–radical reaction. Values for the H + CO → HCO reference reaction are also shown. Bottom half: des- Fig. 1. Arrhenius plots, namely rate constants as a function of the orption energies (Edes ) and desorption (T des ) and diffusion (T diff ) tem- inverse of temperature, for the reaction CH3 + HCO forming acetalde- peratures (in K) derived using the Edes assuming diffusion-to-desorption hyde (black solid line) or CO + CH4 (black dashed line), and for the energy ratios of 0.5, 0.4 and 0.3, see Sect. 2.3. reaction H + CO → HCO (gray dotted-dashed line), described in the main text. 9.2 kJ mol−1 , for two reasons: (i) Rimola et al. assumed an Eley– 4.3. Desorption and diffusion temperatures Rideal reaction (namely, the H atom comes from the gas phase Table 1 reports the computed Edes and the temperature for and reacts with frozen CO), while here we have considered a LH desorption T des and diffusion T diff derived assuming a half- mechanism (Sect. 3.2), and (ii) they did not consider dispersion life of 1 Myr. Edes values are obtained at BHLYP-D3(BJ)/6- corrections. 311++G(2df,2pd)//BHLYP-D3(BJ)/6-31+G(d,p) level following the procedure explained in Sect. 2.3, which moreover are cor- 4.2. Rate constants rected for deformation and basis set superposition energy. The T des and T diff values are obtained by using the standard equation Figure 1 shows the rate constants as a function of the tempera- for the half-life time, t1/2 = ln(2)/kdiff/des (T ). These timescales ture of the reactions that form CH3 CHO and CO + CH4 from the provide an estimation of the characteristic temperatures for coupling and direct H-abstraction of CH3 + HCO, respectively. desorption and diffusion of the two radicals, CH3 and HCO, The figure also reports the case of HCO formation from H + CO, involved in the formation of acetaldehyde on the icy surface. for the sake of reference. Finally, we note that our Edes are consistent with those com- The rate constants of the reactions leading to CO + H2 CO and puted by Ferrero et al. (2020) on a substantially larger ASW ice HCO take tunneling into account, which is evidenced by their model. Specifically, our Edes in Table 1 lies in the high end of the deviation from linearity. It is also evident the strong temperature Ferrero et al. range. On the contrary, and as already discussed in dependence of the radical–radical reactions studied, as compared Ferrero et al. (2020), our Edes are different than those reported in to HCO formation. the astrochemical databases KIDA1 and UMIST2 , often used by The rate constants of the acetaldehyde formation are larger modellers. Unfortunately, no experimental data on the CH3 and than those of CO + CH4 formation at temperatures above ∼24 K. HCO Edes desorption energy exist, to our best knowledge. This is due to its lower barrier and the almost negligible quantum tunnelling contribution to HCO + CH3 → CO + CH4 at such tem- 5. Discussion peratures. However, as the temperature decreases the tunnelling probability takes over deviating the rate constant of CO + CH4 5.1. Formation of acetaldehyde versus CO + CH4 formation from linearity, becoming faster than the formation of acetaldehyde. On the contrary, HCO formation has a much We used the results of our new calculations (Sect. 4) of the weaker temperature dependence and higher rate constants over CH3 + HCO reaction kinetics, desorption and diffusion rate con- the considered temperature range. This is the result of the dom- stants to compute the efficiency ε (Eq. (3)) of the two channels inant strong quantum tunnelling of the H atom through the leading to the formation of either CH3 CHO or CO + CH4 . As reaction barrier, in agreement with the literature results (e.g., discussed in Sect. 2.3, given the uncertainty on its value, we Andersson et al. 2011; Rimola et al. 2014). considered three cases for the Ediff /Edes ratio f : 0.3, 0.4 and 0.5. In order to facilitate the introduction of the new rate con- Figure 2 shows the resulting ε as a function of the temperature stants in astrochemical models, we fit the reactions rate constants and Table 2 reports the α, β and γ values obtained by fitting the with the standard formula in Eq. (9). The values of α, β and γ are ε curves with Eq. (9), for the three cases of f . listed in Table 2. We note that, although we computed the efficiency of the reactions in the 5–100 K range, they will only take place as long  T β 1 kaeb (T ) = α exp(−γ/T ). (9) http://kida.astrophy.u-bordeaux.fr/ 300 K 2 http://udfa.ajmarkwick.net/ A9, page 5 of 14 A&A 655, A9 (2021) Table 2. Rate constants kaeb (in s−1 ) and efficiency ε of the two possible reactions between HCO and CH3 . Rate Temperature α β γ 9K 20 K 30 K Ediff /Edes constant (K) (s−1 ) (K) Reaction (1): HCO + CH3 → CH3 CHO kaeb 2.7 × 10−21 1.9 × 10−3 160.5 12 9–30 3.1 × 10 0.70 663 ε 1.0 1.0 1.0 9–19 1.0 0.0 0.0 0.5 19–26 0.99 –3.4 × 10−3 0.06 26–30 0.98 –0.01 0.28 ε 0.88 0.81 0.81 0.4 9–13 0.43 –0.21 0.28 13–30 1.0 0.14 −3.4 ε 4.9 × 10−8 7.6 × 10−3 0.01 0.3 9–30 3.3 0.12 161.2 Reaction (2): HCO + CH3 → CH4 + CO kaeb 6.8 × 10−5 0.05 20.1 9 9–15 6.1 × 10 10.4 –39.7 15–30 1.7 × 1023 25.9 –274 ε 1.0 1.0 1.0 0.5 9–24 1.0 0.0 0.0 24–28 0.83 –0.12 2.79 28–30 0.53 –0.44 11.4 ε 1.0 1.0 0.35 9–21 1.0 0.0 0.0 0.4 21–25 7.4 × 10−11 –13.8 278.8 25–30 3.2 × 10−6 –5.8 53.3 ε 1.0 0.02 1.5 × 10−3 9–15 0.96 –0.01 0.02 0.3 15–19 2.4 × 10−76 –84.9 1205.1 19–30 2.4 × 107 19.9 –660.4 Notes. For each reaction, we report the values of α, β and γ of the rate constant kaeb and the efficiency ε calculated assuming Ediff /Edes equal to 0.5, 0.4 and 0.3 (first column). The last three columns report the values of kaeb and ε at 9, 20 and 30 K. as one of the two radicals can diffuse and scan the ASW sites going into lower values so that at 30 K it reaches ε ∼0.3, while and both radicals stay on the reaction site, namely they do not the efficiency of acetaldehyde formation stays about unity up to desorb (see Table 1). Consequently, the upper limit to the tem- 30 K, where it takes a value of ∼0.8. On the other hand, for perature where the CH3 CHO and CO + CH4 formation reactions f = 0.3 things are very different. The efficiency of CO + CH4 take place is set by the desorption of CH3 , as it has a lower crashed at higher temperatures, reaching values of about 0.001 desorption energy than HCO (30 and 68 K, respectively, see at 30 K, while that of acetaldehyde never goes above ∼0.01. Table 1). Likewise, the lower limit is also set by the CH3 diffu- This is because, for relatively large f values (≥0.4), the most sion energy only, which is equal to 15, 12 and 9 K for f equal to mobile radical, CH3 , moves slowly and the two radicals have 0.5, 0.4 and 0.3, respectively. In the case of f equal to 0.5, HCO plenty of time to react when they meet before one of them moves starts to be mobile when CH3 has already sublimated, so that the away: ε is, therefore, close to unity. However, when the timescale efficiency of the reaction depends on CH3 Ediff only. Conversely, for diffusion becomes smaller than the reaction timescale (i.e., for f equal to 0.4 and 0.3, the temperatures at which HCO and kdiff  kaeb ), CH3 moves away before having the time to react CH3 can diffuse overlap, so that both species contribute to the and the efficiency drops below unity. In practice, the smaller the denominator of Eq. (3). Ediff /Edes ratio, the faster CH3 moves and the smaller ε. How- Formation efficiency. For both reactions, formation of ever, since both kdiff and kaeb have an exponential dependence on CH3 CHO and CO + CH4 , the efficiency ε is about 1 in the 9– the temperature, a change in behavior occurs when the reaction 15 K range regardless of the f value (between 0.3 and 0.5) with activation energy γ (Eq. (9)) is similar to Ediff and the effi- one exception, acetaldehyde formation with f = 0.3, which starts ciency ε strongly depends on the temperature. For the formation at very low efficiency values and monotonically increases. For of acetaldehyde, γ = 663 K (Table 2) and, therefore, the change f = 0.5, either reactions have efficiencies of about unity in the of behavior occurs when Ediff /Edes ∼ 0.40. In these cases, the whole range of temperatures (up to 30 K). For f = 0.4, formation lower the temperature, the larger kdiff with respect to kaeb and the of CO + CH4 distances from unity at temperatures above ∼22 K, smaller ε, as shown in Fig. 2. Similar arguments hold also for the A9, page 6 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains f (see above), to the temperature T . It holds: RT T0 dT 0 εCH3 CHO × (Rdiff,CH3 + Rdiff,HCO ) BR(T ) = R T . T0 dT 0 (εCH3 CHO + εCO + CH4 ) × (Rdiff,CH3 + Rdiff,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 branch- ing ratio BR is constant and equal to 0.5, namely the HCO + CH3 reaction leads to acetaldehyde and CO + CH4 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 + CH4 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 + CH4 are Fig. 2. Reaction efficiency ε (Eq. (3)) of the reaction CH3 + HCO lead- in approximately equal competition in the range of temperatures ing to either CH3 CHO (solid lines) or CO + CH4 (dashed lines) as a where the HCO + CH3 reaction can occur. However, for f < 0.4, function of the temperature. The computations were obtained adopting acetaldehyde is a very minor product for temperatures lower than three different Ediff /Edes ratios: 0.3 (green), 0.4 (blue) and 0.5 (red). We about 25 K and flips to be a major product above it. note that, for Ediff /Edes = 0.5 the CH3 CHO and CO + CH4 (red) curves overlap, namely they are constant and equal to 1. 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 par- ticles) irradiation of different H2 O, CO, CH3 OH and CH4 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 irra- diation of a CO:CH4 ice mixture, detected acetaldehyde and predicted that the orientation of the CH3 and HCO radicals are crucial in the efficiency of the reaction. On the other hand, Martín-Doménech et al. (2020) conducted laboratory experi- ments on the formation of acetaldehyde via the CH3 + HCO reaction, concluding that this channel is not efficient enough to reproduce the astronomical observations. As discussed by various authors, although laboratory exper- iments are primordial in suggesting possible mechanisms oper- ating in the ISM and, specifically, the possible formation routes of molecules on the interstellar grain surfaces, they cannot pro- vide the exact ISM conditions or a detailed description of the mechanisms at the atomic level. Despite this, the improvement of Fig. 3. Branching ratio BR(T ) of the formation rate of the CH3 CHO radical detection methods, such as the electron paramagnetic res- over CO + CH4 (Eq. (10)) as a function of the temperature in the range onance (EPR) technique, will help to clarify the role of radicals where the reactions can occur, namely below 30 K (see text), for generated in interstellar ice analogs (e.g., Zhitnikov & Dmitriev Ediff /Edes equal to 0.3 (green), 0.4 (blue) and 0.5 (red). 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. HCO + CH3 → CO + CH4 reaction. The only difference is that, at low temperatures, kaeb deviates from the exponential law because 5.3. Astrophysical implications of the kicking in of the tunneling effect that greatly increases kaeb (giving a negative γ values: see Table 2). Since the tun- In astrochemical models, it is generally assumed that reactions neling is more efficient for decreasing temperature, the kaeb /kdiff between radicals on the surface of interstellar ices are barrierless ratio decreases at increasing temperatures and, consequently, ε and, consequently, that their efficiency is equal to 1 (Sect. 2). In decreases. addition, it is also often assumed that there are no competition channels to the production of iCOMs. At variance with these Branching ratio. Figure 3 shows the branching ratio BR of simple assumptions, our new calculations presented in Sect. 5.1, the formation rate of CH3 CHO over CO + CH4 as a function of indicate that, at low (≤15 K) temperatures, the efficiency of the the temperature, for the three f values (0.5, 0.4 and 0.3). The acetaldehyde formation is close to unity, for a Ediff /Edes ratio BR is obtained integrating Eq. (1) from the temperature at which f ≥ 0.40. However, there is a competing channel leading to CH3 starts to be mobile T 0 , a value that depends on the assumed CO + CH4 , for which the efficiency is also equal to 1, so that, at A9, page 7 of 14 A&A 655, A9 (2021) low temperatures and for f ≥ 0.40 the two channels are equally and experiments on species such as CO, CO2 , H2 O, CH4 and probable. The acetaldehyde formation efficiency remains close NH3 (e.g., Mispelaer et al. 2013; Karssemeijer & Cuppen 2014; to unity in the temperature range where the reaction can occur, Lauck et al. 2015; Ghesquière et al. 2015; He et al. 2017, 2018; namely at ≤30 K, for f ≥ 0.40. However, the situation drastically Cooke et al. 2018; Maté et al. 2020; Kouchi et al. 2020). These changes for f < 0.40. Specifically, for f = 0.3, the efficiency of studies give a value for f between 0.3 and 0.6, as mentioned in acetaldehyde formation crashes to very low values and increases Sect. 2.3. However, rigorously speaking, the experiments do not with temperature to a maximum of 0.01 at 30 K. Similarly, the necessarily measure the same diffusion processes as in interstel- formation of CO + CH4 drops to 1.5 × 10−3 at 30 K. lar conditions, for at least the reasons of the surface coverage Therefore, two major messages come out from our calcula- (He et al. 2018) and its dependence of the nature of the ice, tions: (1) the efficiency of the formation of acetaldehyde from specifically its degree of porosity (Maté et al. 2020), which is the HCO + CH3 reaction on icy surfaces is a complex func- poorly known in the case of interstellar ices. As a matter of fact, tion of the temperature and of the CH3 diffusion energy Ediff using experimental and a theoretical Monte Carlo code, Maté (Fig. 2) and (2) the acetaldehyde formation receives competi- et al. (2020) found that “the microscopic diffusion is many times tion with CO + CH4 formation, which cannot be neglected and faster than the macroscopic diffusion measured experimentally". whose efficiency is also a complex function of temperature and Most recently, Kouchi et al. (2020) obtained a direct measure- Ediff (Fig. 3). ment of the diffusion energy of CO and CO2 on ASW, using While the dependence on the temperature and the importance the transmission electron microscopy (TEM) technique, which of the competition of other products were already recognized allows a direct measurement of the surface diffusion coefficients (Enrique-Romero et al. 2019), the paramount importance of the (against the often used technique of IR spectroscopy, which only diffusion energy Ediff in the radical–radical reactions efficiency indirectly estimates the diffusion energy). Kouchi and coworkers was not appreciated, at least not at the extent indicated by this found an f ratio equal to 0.3. study (because of the assumption of astrochemical models that We have seen that, in the case of acetaldehyde, this uncer- the efficiency of the radical–radical reactions on grains is 1). tainty on f has a dramatic effect. If f is >0.4, the efficiency of Penteado et al. (2017), for example, carried out an extensive acetaldehyde formation is equal to 1 and it is about equal to that study of the surface chemistry on the binding energies (namely, of the CO + CH4 formation. On the contrary, if f is equal to 0.3, our Edes ) showing how critical they are. Our new study sug- then the efficiency of acetaldehyde formation (and CO + CH4 ) gests that Ediff is as much, or even more, crucial in the reactions crashes, to a maximum value of 0.02. The most recent measure- involving two radicals on ASW. ments by Kouchi et al. (2020) point out the latter case as the most What makes the situation actually critical is that, while stud- probable. If the value f = 0.3 is confirmed, then acetaldehyde is ies of the binding energy of radicals can and have been estimated unlikely to be formed on the interstellar icy grain surfaces. in experimental and theoretical works (e.g., see the recent works One could be tempted to use the astronomical observations by Penteado et al. 2017; Ferrero et al. 2020), evaluating the dif- against the astrochemical model predictions to add constraints to fusion energy of multi-atomic radicals on cold icy surfaces has the f value in (real) interstellar ices. Of course, given the large proven to be extremely complicated and, to the best of our knowl- number of parameters associated with the astrochemical models edge, no experimental or theoretical studies exist in the literature it could be a dangerous exercise. Nonetheless, we can analyze (see also e.g., Cuppen et al. 2017; Potapov & McCoustra 2021). two cases, as illustrative examples. Barger & Garrod (2020) com- Indeed, as mentioned above, obtaining Ediff experimentally is pared the predictions of their model, where the formation of hitherto hampered by technical limitations on the instrumenta- acetaldehyde is dominated by the reaction CH3 + HCO assumed tion used to detect the radicals (i.e., EPR measurements). In to have ε = 1, with the observations toward various hot cores relation to theoretical investigations, this lacking in bibliogra- and found that in two of them, NGC 7538 IRS 1 and W3(H2 O), phy is due to the convergence of methodological difficulties their model overproduces the acetaldehyde column densities by that make the study of diffusion with computational simula- more than a factor 103 with respect to the observed ones. If f is tions intrinsically complex (but also compelling). Diffusion can equal to 0.3, introducing our new values for ε could possibly cure currently be studied by means of molecular dynamics (MD) or this mismatch. On the contrary, Jørgensen et al. (2016) found a kinetic Monte Carlo (kMC) simulations. With the first, to obtain good agreement between the observed abundances of acetalde- a sufficient representativeness of the species diffusion, long hyde in IRAS16293B and SgrB2(N) and those predicted by the simulation time-scales are mandatory. This in practice means Garrod (2013) model. In this case, the agreement would point to adopt classical force fields, in which the electronic struc- to f ≥ 0.4. In other words, our new computations might solve ture of the systems is missing. However, radicals are open-shell the mismatch observed toward NGC 7538 IRS 1 and W3(H2 O) species (with at least one unpaired electron) and accordingly if f = 0.3, but they would create a mismatch on the observations electrons have to be accounted for. Thus MD simulations should toward IRAS16293B and SgrB2(N), or vice versa. Alternatively, be grounded within the quantum mechanics realm, which are it is possible that f varies in different sources, belonging to dif- much more expensive than the classical ones, making the MD ferent environments. For example, one could think that sources simulations unfeasible. The alternative would be the adoption in cold quiescent regions have ices different, more or less porous, of kMC simulations. However, these simulations require build- from those in warm and chaotic ones. The two examples dis- ing a complete network of the site-to-site radical hopping, in cussed above, unfortunately, do not lead to a coherent behavior, which for each hopping the corresponding rate constant has to as, for example, IRAS16293B and SgrB2(N) could not belong to be known a priori. This actually means to localize for each hop- more different environments. ping the corresponding transition state structure (at a quantum In conclusion, the acetaldehyde formation by radical–radical chemical level), in which by using a realistic ASW model (i.e., recombination on the ices is such a strong function of the dif- large, amorphous and accordingly plenty of binding sites) makes fusion energy, likely linked to the nature of the ice, that a the problem unpractical. little variation of the Ediff /Edes (by 0.1) value can shift the effi- Usually, astrochemical models assume that the radical Ediff is ciency from 1 to less than 0.01. The most recent estimates of a fraction f of Edes and the value f is derived from computations Ediff /Edes suggest a value of 0.3 (Kouchi et al. 2020), which A9, page 8 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains would make the formation of acetaldehyde on the grain surfaces Finally, this new study calls for specific similar computations unlikely. Anyway, the important message here is that astrochem- on the radical–radical reactions assumed to form iCOMs in astro- ical model predictions should be taken with a certain precaution. chemical models as assuming that they have efficiency ε equal to On the contrary, our new computations clearly show the huge 1 and are the only reaction product could be highly misleading. importance of better knowing the microprocesses involved in the radical–radical chemistry on the icy interstellar grains and Acknowledgements. This project has received funding within the European the urgent need of extensive studies, similar to the one pre- Union’s Horizon 2020 research and innovation programme from the Euro- sented here, on different systems believed to form iCOMs on the pean Research Council (ERC) for the projects “The Dawn of Organic Chemistry” (DOC), grant agreement no. 741002 and “Quantum Chemistry interstellar icy surfaces. on Interstellar Grains” (QUANTUMGRAIN), grant agreement no. 865 657, For the sake of completeness, the formation of acetalde- and from the Marie Sklodowska-Curie for the project “Astro-Chemical Ori- hyde via radical–radical reactions on surfaces with lower binding gins” (ACO), grant agreement no. 811312. AR is indebted to “Ramón y energies such as solid CO could have a higher efficiency, due Cajal” program. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged. Finally, we thank Prof. Gretobape for fruit- to the low radical–surface interactions (Lamberts et al. 2019). ful and stimulating discussions. Most of the quantum chemistry calculations Finally, it is worth reminding that acetaldehyde can alternatively presented in this paper were performed using the GRICAD infrastructure be synthesized in the gas-phase (e.g., Charnley 2004; Vastel et al. (https://gricad.univ-grenoble-alpes.fr), which is partly supported by 2014; De Simone et al. 2020). Recently, Vazart et al. (2020) the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme reviewed the gas-phase routes leading to acetaldehyde and found Investissements d’Avenir supervized by the Agence Nationale pour la Recherche. Additionally this work was granted access to the HPC resources of IDRIS under that, very likely, the dominant one is that starting from ethanol, the allocation 2019-A0060810797 attributed by GENCI (Grand Equipement the so-called ethanol tree (Skouteris et al. 2018). In particular, the National de Calcul Intensif). We thank prof. Gretobape for stimulating discus- ethanol tree route reproduces quite well both the acetaldehyde sions, and J. Perrero for useful contributions. and glycolaldehyde abundances in the sources where ethanol was also observed, including IRAS16293B (Vazart et al. 2020). References Aikawa, Y., Furuya, K., Yamamoto, S., & Sakai, N. 2020, ApJ, 897, 110 6. Conclusions Álvarez-Barcia, S., Russ, P., Kästner, J., & Lamberts, T. 2018, MNRAS, 479, 2007 In this work, we report new computations on the energetics Andersson, S., Goumans, T. P. M., & Arnaldsson, A. 2011, Chem. Phys. Lett., and kinetics of the reaction HCO + CH3 , which can lead to 513, 31 the formation of either acetaldehyde or CH4 + CO. Specifically, Bar-Nun, A., Dror, J., Kochavi, E., & Laufer, D. 1987, Phys. Rev. B, 35, 2427 Barger, C. J., & Garrod, R. T. 2020, ApJ, 888, 38 we compute the rate constants of both reactions as a function Becke, A. D. 1993, J. Chem. Phys., 98, 1372 of temperature as well as the efficiency of the formation of Bennett, C. J., Jamieson, C. S., Osamura, Y., & Kaiser, R. I. 2005, ApJ, 624, acetaldehyde and CH4 + CO, respectively, combining reaction 1097 kinetics at RRKM (tunneling included) with diffusion and des- Bianchi, E., Codella, C., Ceccarelli, C., et al. 2019, MNRAS, 483, 1850 orption competitive channels. We provide analytical formulae so Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621 that the computed rate constants and efficiency can be easily Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, introduced in astrochemical models. 541 The main conclusions of our study are the following: Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 593, L51 1. The HCO + CH3 reaction can only occur when the surface Ceccarelli, C., Caselli, P., Fontani, F., et al. 2017, ApJ, 850, 176 Chang, Q., Cuppen, H. M., & Herbst, E. 2007, A&A, 469, 973 temperature is lower than 30 K, because CH3 desorbs at larger Charnley, S. B. 2004, Adv. Space Res., 33, 23 temperatures. Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E., et al. 2020, 2. Our computations suggest that acetaldehyde is not the Space Sci. Rev., 216, 50 dominant product for the reaction HCO + CH3 . The efficiency ε Cooke, I. R., Öberg, K. I., Fayolle, E. C., Peeler, Z., & Bergner, J. B. 2018, ApJ, of its formation strongly depends on the Ediff /Edes ratio, provid- 852, 75 Csengeri, T., Belloche, A., Bontemps, S., et al. 2019, A&A, 632, A57 ing dramatic variations between 0.3–0.5 values, the most usually Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 used values in astrochemical models. De Simone, M., Codella, C., Ceccarelli, C., et al. 2020, A&A, 640, A75 3. At low (≤15 K) temperatures, ε is close to unity for both Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, A&A, 512, A30 the formation of acetaldehyde and its competing CO + CH4 chan- Enrique-Romero, J., Rimola, A., Ceccarelli, C., et al. 2019, ACS Earth Space nel for f ≥ 0.4, while only the efficiency of CO + CH4 is unity Chem., 3, 2158 Enrique-Romero, J., Álvarez-Barcia, S., Kolb, F. J., et al. 2020, MNRAS, 493, at these temperatures for f = 0.3 thanks to quantum tunnelling. 2523 The efficiency of acetaldehyde formation remains unity in the Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11 range of temperatures where the reaction can occur (≤30 K) Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian’16 Revision for Ediff /Edes ≥ 0.40. For lower Ediff /Edes ratios, ε becomes C.01 (Wallingford, CT: Gaussian Inc.) Garrod, R. T. 2013, ApJ, 778, 158 1 and increases with increasing temperature: in the case of Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 Ediff /Edes = 0.3, it reaches a maximum of ∼0.01 at 30 K. Con- Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15 versely, the efficiency of the formation of CO + CH4 increases Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283 with decreasing temperature because of the tunneling. Garrod, R. T., Vasyunin, A. I., Semenov, D. A., Wiebe, D. S., & Henning, T. 4. These variant ε values as a function of Ediff /Edes go against 2009, ApJ, 700, L43 Ghesquière, P., Mineva, T., Talbi, D., et al. 2015, Phys. Chem. Chem. Phys., 17, the assumption made in many astrochemical models, in which 11455 ε is equal to 1. This might have a substantial impact on the Goerigk, L., Hansen, A., Bauer, C., et al. 2017, Phys. Chem. Chem. Phys., 19, acetaldehyde abundance predicted by these models, which may 32184 overestimate it by a few orders of magnitude. Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, J. Chem. Phys., 132, 154104 5. We discussed the example of IRAS16293B and suggested Grimme, S., Ehrlich, S., & Goerigk, L. 2011, J. Comput. Chem., 32, 1456 that, in this object, acetaldehyde is likely synthesized by a gas- Hama, T., & Watanabe, N. 2013, Chem. Rev., 113, 8783 phase reaction route that starts from ethanol. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 A9, page 9 of 14 A&A 655, A9 (2021) He, J., & Vidali, G. 2014, ApJ, 788, 50 Öberg, K. I., van Dishoeck, E. F., Linnartz, H., & Andersson, S. 2010, ApJ, 718, He, J., Emtiaz, S. M., & Vidali, G. 2017, ApJ, 837, 65 832 He, J., Emtiaz, S., & Vidali, G. 2018, ApJ, 863, 156 Palumbo, M. E. 2006, A&A, 453, 903 Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427 Palumbo, M. E., Baratta, G. A., Leto, G., & Strazzulla, G. 2010, J. Mol. Struct., Hollenbach, D., & Salpeter, E. 1970, J. Chem. Phys., 53, 79 972, 64 Hudson, R. L., & Moore, M. H. 1997, Icarus, 126, 233 Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, ApJ, 897, 56 Isokoski, K., Bossa, J. B., Triemstra, T., & Linnartz, H. 2014, Phys. Chem. Chem. Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49 Phys., 16, 3456 Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71 Jensen, S. S., Jørgensen, J. K., Furuya, K., Haugbølle, T., & Aikawa, Y. 2021, Potapov, A., & McCoustra, M. 2021, Int. Rev. Phys. Chem., 40, 299 A&A, 649, A66 Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, 572, A70 A117 Rimola, A., Skouteris, D., Balucani, N., et al. 2018, ACS Earth Space Chem., 2, Karssemeijer, L. J., & Cuppen, H. M. 2014, A&A, 569, A107 720 Keane, J. V., Tielens, A. G. G. M., Boogert, A. C. A., Schutte, W. A., & Whittet, Ruaud, M., Loison, J. C., Hickson, K. M., et al. 2015, MNRAS, 447, 4004 D. C. B. 2001, A&A, 376, 254 Ruffle, D. P., & Herbst, E. 2000, MNRAS, 319, 837 Kouchi, A., Furuya, K., Hama, T., et al. 2020, ApJ, 891, A22 Sakai, T., Yanagida, T., Furuya, K., et al. 2018, ApJ, 857, 35 Lamberts, T., & Kästner, J. 2017, ApJ, 846, 43 Scibelli, S., & Shirley, Y. 2020, ApJ, 891, 73 Lamberts, T., Markmeyer, M. N., Kolb, F. J., & Kästner, J. 2019, ACS Earth Skouteris, D., Balucani, N., Ceccarelli, C., et al. 2018, ApJ, 854, 135 Space Chem., 3, 958 Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1989, ApJ, 344, 413 Lauck, T., Karssemeijer, L., Shulenberger, K., et al. 2015, ApJ, 801, 118 Taquet, V., Ceccarelli, C., & Kahane, C. 2012, A&A, 538, A42 Lee, C., Yang, W., & Parr, R. G. 1988, Phys. Rev. B, 37, 785 Tielens, A. G. G. M., & Allamandola, L. J. 1987, in Interstellar Processes, eds. Lee, C.-F., Codella, C., Li, Z.-Y., & Liu, S.-Y. 2019, ApJ, 876, 63 D. J. Hollenbach, & H. A. Thronson (Dordrecht: Springer Netherlands), 397 Lefloch, B., Ceccarelli, C., Codella, C., et al. 2017, MNRAS, 469, L73 Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 Martín-Doménech, R., Öberg, K. I., & Rajappan, M. 2020, ApJ, 894, 98 Vastel, C., Ceccarelli, C., Lefloch, B., & Bachiller, R. 2014, ApJ, 795, L2 Maté, B., Cazaux, S., Satorre, M. Á., et al. 2020, A&A, 643, A163 Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33 Minissale, M., Congiu, E., & Dulieu, F. 2016, A&A, 585, A146 Vazart, F., Ceccarelli, C., Balucani, N., Bianchi, E., & Skouteris, D. 2020, Mispelaer, F., Theulé, P., Aouididi, H., et al. 2013, A&A, 555, A13 MNRAS, 499, 5547 Molinari, S., Ceccarelli, C., White, G. J., et al. 1999, ApJ, 521, L71 Vidali, G. 2013, Chem. Rev., 113, 8752 Molpeceres, G., Rimola, A., Ceccarelli, C., et al. 2018, MNRAS, 482, Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1 5389 Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173 Neese, F. 2004, J. Phys. Chem. Solids, 65, 781 Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427 Oba, Y., Miyauchi, N., Hidaka, H., et al. 2009, ApJ, 701, 464 Zhitnikov, R. A., & Dmitriev, Y. A. 2002, A&A, 386, 1129 A9, page 10 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains Appendix A: Benchmark study Table A.3. Raw data from the DFT method benchmark. Energies in kJ mol−1 . The quality of BHLYP-D3(BJ) as an accurate, cost-effective DFT/6-311++G(2df,2pd) method for the reactions studied in this work is shown in this BHLYP-D3(BJ) MPWB1K-D3(BJ) M062X-D3 PW6B95-D3(BJ) wB97X-D3 React. I -306.5278673 -306.5199277 -306.5443283 -307.0275371 -306.6077421 section. TS I -306.5274906 -306.5195273 -306.5429782 -307.0270829 -306.6074261 We have taken five hybrid DFT dispersion-corrected meth- Prod. I -306.6651096 -306.6673671 -306.6888795 -307.1714705 -306.7523061 React. II -306.5278665 -306.5199185 -306.5443285 -307.0275286 -306.6077403 ods: MPWB1K-D3(BJ), M062X-D3, PW6B95-D3(BJ), wB97X- TS II -306.5273893 -306.5195091 -306.5435295 -307.027106 -306.6074694 D3 and BHandHLYP-D3(BJ), recommended in Goerigk et al. Prod. II -306.6561571 -306.6555907 -306.6825872 -307.1623168 -307.0257242 (2017) for their good overall performance. Multi-reference methods We have then compared their performance with respect to React. I CASSCF(2,2)/aug-cc-PVTZ CASPT2/aug-cc-PVTZ -304.99849 -306.092193 CASPT2 by studying reactions I and II on 2 water molecules. TS I -304.99793 -306.0918077 We proceeded in two steps: (i) we performed geometry opti- Prod. I React. II -305.14355 -304.99849 -306.2278841 -306.0921923 mized and run frequency calculations at BHLYP-D3(BJ)/6- TS II -304.9976 -306.0916902 311++G(d,p) for each reaction channel, finding and checking Prod. II -305.14411 -306.2251986 each stationary point (i.e., reactants, transition state and prod- ucts); (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 ref- erence. The unrestricted broken symmetry scheme was adopted for all DFT calculations, and the CASPT2 guess wave func- tion 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 Fig. A.1. Geometries of the stationary points (reactants, transition state respect to CASPT2/aug-cc-PVTZ. The rest of DFT methods and products) of reactions a) I and b) II on top of two water molecules have average errors between 10 and 80% (Table A.1). The raw optimized at BHLYP-D3(BJ)/6-311++G(d,p) level. Distances in Å. energy values are shown in Table A.3. Table A.1. DFT method benchmark results. % Unsigned Error Appendix B: Calculation of desorption and from CASPT2/aug-cc-PVTZ//BHLYP-D3(BJ)/6-311++G(d,p). diffusion rate constants Methods are: (1) BHLYP-D3(BJ), (2) MPWB1K-D3(BJ), (3) In order to calculate the rates of diffusion and desorption we have M062X-D3, (4) PW6B95-D3(BJ) and (5) wB97x-D3. used Eyring’s equation (Eq. 5), where entropy and thermal cor- System / Method 1 2 3 4 5 rections to the enthalpy are accounted for since energy term in TS dHa on W2 5.0 18.5 59.1 15.8 46.0 the exponential should be the Gibbs free energy (G = H − T S ). HCO + CH4 on W2 3.5 2.0 3.9 1.3 214.3 In the astrochemistry community, desorption and diffusion rates TS Rc on W2 2.2 3.9 250.4 17.9 18.0 consist of two parts, first the attempt frequency (i.e., the pre- CH3CHO on W2 1.1 8.7 6.5 6.1 6.5 exponential factor) and then the exponential with the binding Average 3.0 8.3 80.0 10.3 71.2 energies. We define the attempt frequency within the Eyring equation as Table A.2. Energetics (kJ mol−1 ) of the stationary points opti- kB T mized at BHLYP-D3(BJ)/6-311++G(d,p), together with the fre- ν= × exp(∆S /kB ), (B.1) h quencies of the transition states (cm−1 ). where h, kB are the Planck and Boltzmann constants, T the tem- OPT: BHLYP-D3(BJ)/6-311++G(d,p) perature and S the entropy. For the ice model and the complex E iv surface+radical the entropy contributions are S = S vib , where the React. I -306.5101581 vibrational counterpart is (Eq. B.2): TS I -306.5101000 -26.1495 Prod. I -306.6472250 X " θv,i /T # React. II -306.5106910 −θv,i /T  S vib = kB − ln 1 − e (B.2) TS II -306.5098973 -32.8207 i eθv,i /T − 1 Prod. II -306.6377604 where θv,i = hcνi /kB with c the speed of light and νi 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: π1/2 T 3/2 3 " ! # S rot = kB ln + , (B.3) σrot (θr,x θr,y θr,z )1/2 2 A9, page 11 of 14 A&A 655, A9 (2021)   2πmkB 3/2 kB T      5  ! S trans = kB ln   + , (B.4) h2 P  2 where σrot is the rotational symmetry number (6 for CH3 and 1 for HCO), θr,i = Bi h/kB with Bi 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 104 cm−3 . The thermal corrections follow: H = E DFT + ZPE + Evib (T ) + kB T 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 Evib (T ) is the thermal vibrational energy calculated with Eq. B.5. In analogy to to the entropy, the rota- tional and translational contributions are also included for the free radicals: H = E DFT + ZPE + Evib (T ) + kB T + Hrot + Htrans where Hrot = Htrans = 32 kB T . X θv,i Evib (T ) = kB (B.5) i eθv,i /T − 1 Fig. C.1. Comparison of reaction, diffusion and desorption rate con- stants involved in the CH3 + HCO system. Notice that the desorption The magnitude of these contributions is shown in Table B.1, rate of HCO is not seen as it appears at very low rate constant values. where it can be seen that the most important contribution to the Numbers in brackets indicate the diffusion-to-desorption energy barrier enthalpy comes from the vibrational modes, while for entropy all ratio. contributions are rather small: Table B.1. Corrections incorporated in the calculation of ∆G of Appendix D: H + CO PES and data 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 Boltz- mann’s constant is about 0.0083 kJ mol−1 , therefore the kB T and 3kB T /2 terms are also very small. T = 20 K T = 30 K rad rad rad rad ∆Evib ∆S vib S rot S trans ∆Evib ∆S vib S rot S trans CH3 0.124 -0.007 0.010 0.005 0.293 -0.003 0.015 0.011 HCO 0.077 -0.012 0.034 0.011 0.228 -0.009 0.039 0.016 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 multi- ply the ∆G of desorption times these fractions (assumed to be 0.3, 0.4 and 0.5). With this, we get different attempt frequencies Fig. D.1. Potential energy surface of the H + CO → HCO reaction, for diffusion and desorption (see Table B.2), all around 108 –1011 in kJ mol−1 . Energies are corrected for dispersion and ZPE. Geome- s−1 , smaller than the normally used approach of the harmonic tries and ZPE energies were obtained at UBHandHLYP-D3(BJ)/6- oscillator3 in astrochemical models (e.g., Tielens & Allaman- 31+G(d,p) level and DFT energies were refined at UBHandHLYP- dola (1987); Hasegawa et al. (1992)) due to the use of Eyring’s D3(BJ)/6-311++G(2df,2pd) level. Reactants and products were obtained relation and the inclusion of entropy. by running intrinsic reaction coordinate calculations. Table B.2. Attempt frequencies for desorption and the different Table D.1. H + CO → HCO energetic data, in Hartree, at cases of diffusion considered in this work. Units in s−1 . See that UBHandHLYP-D3(BJ)/6-31+G(d,p) (double ζ) level. U is the at 20 K kB T/h ∼ 4.2 × 1011 , so the effect of entropy at such low DFT energy, D is the Dispersion energy and ZPE is the zero- temperatures is very small. point energy. DFT energies were refined by performing single Eyring equation (computed at 20 K) Harmonic point calculations on double ζ geometries at UBHandHLYP- Desorption Diffusion oscillator D3(BJ)/6-311++G(2df,2pd) (triple ζ) level. Energy units are Diff-to-Des – 0.5 0.4 0.3 – Hartree (1.0 Hartree are ∼ 2625.5 kJ mol−1 ). ν(CH3 ) 4.7 × 108 1.1 × 1011 1.4 × 1011 1.9 × 1011 1.5 × 1012 ν(HCO) 2.9 × 1010 1.4 × 1010 2.8 × 1010 5.4 × 1010 1.5 × 1012 H + CO Reactant TS Product U (double ζ) -2558.84099479 -2558.83665194 -2558.88783666 D (double ζ) -0.08428582 -0.08266012 -0.08407096 ZPE (double ζ) 0.86306000 0.86239000 0.872578000 Appendix C: Rate constant comparison U (triple ζ) -2559.79323262 -2559.78741062 -2559.83765264 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. ν = 2N s Ebind /π2 m, with NS the density of sites, ∼1015 and m the p 3 mass of the particle. A9, page 12 of 14 J. Enrique-Romero et al.: Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains Table D.2. H· · · CO transition state data. Units are the usual Appendix F: Efficiency figures, separated by for a GAUSSIAN 16 output: frequencies in cm−1 , IR intensities E di f f /E des ratios in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYP-D3(BJ)/6-31+G(d,p) level. iν red mass F ctn IR int H + CO 605.7354 1.1042 0.2387 14.8487 Appendix E: Radical-radical TS data and PES: Table E.1. CH3 + HCO energetic data, in Hartree, at UBHandHLYP-D3(BJ)/6-31+G(d,p) level (double ζ). U is the DFT energy, D is the Dispersion energy and ZPE is the zero- point energy. DFT energies were refined by performing single point calculations on double ζ geometries at UBHandHLYP- D3(BJ)/6-311++G(2df,2pd) (triple ζ) level. Energy units are Hartree (1.0 Hartree are ∼ 2625.5 kJ mol−1 ). React. I Reactant TS Product U (double ζ) -2675.130087 -2675.127479 -2675.266272 D (double ζ) -0.091309097 -0.091229233 -0.088094692 ZPE (double ζ) 0.933864 0.933528 0.94406 U (triple ζ) -2676.122287 -2676.119848 -2676.253986 React. II Reactant TS Product U (double ζ) -2675.130087 -2675.12475 -2675.255587 D (double ζ) -0.091309097 -0.092369778 -0.091838262 ZPE (double ζ) 0.933864 0.933255 0.938711 U (triple ζ) -2676.122287 -2676.118937 -2676.250317 Table E.2. Features of the transition states studied in this work. Units are the usual for a GAUSSIAN 16 output: frequencies in cm−1 , IR intensities in KM/Mole, reduced masses in AMU and force constants in mDyne/Å. At UBHandHLYP-D3(BJ)/6- Fig. F.1. CH3 + HCO reaction efficiencies ε (Eq. 3), assuming diffusion 31+G(d,p) level. barriers 0.3, 0.4 and 0.5 times those of desorption (panels from top to reaction I reaction II bottom). The green-colored regions indicate the diffusion and desorp- iν red mass F ctn IR int iν red mass F ctn IR int tion temperatures limits of CH3, while the red ones are the same for CH3 + HCO 91.0685 4.7645 0.0233 0.5366 168.6488 1.8100 0.0303 16.6868 HCO. 1.174 1.209 1.486 3.207 1.877 Appendix G: Calculation of diffusion and desorption temperatures .I Re act Half-lives are calculated from the rate constants (ki , with i being 2.477 2.242 either the diffusion or desorption of radicals) following ΔH‡: 5.5 ΔHrx: -324.5 2.391 ln(2) = ln(2) ,, 2.346 1.174 N1/2 N =N 1/2 = N00 exp(−k exp(−kiit) t) =⇒ =⇒ tt1/2 1/2 = Re kkii 2.075 1.084 act .I I H 2.387 1.123 O C and shown and shown in in Fig. Fig. G.1. G.1. (HO)C---C(H3): 4.077 (CO)H---C(H3): 4.655 1.169 ΔH‡: 7.2 ΔHrx: -328.9 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 . Fig. G.1. Diffusion and diffusion temperatures of CH3 and HCO assuming a half-life of 1 Myrs for desorption. A9, page 13 of 14 A&A 655, A9 (2021) 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 Edi f f /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 Edi f f /Edes = 0.5, 0.4 and 0.3 (left to right panels). A9, page 14 of 14