Quantum mechanical modeling of the grain-surface formation of acetaldehyde on H₂O:CO dirty ice surfaces

Jessica Perrero

,^1,2 Piero Ugliengo

,² Cecilia Ceccarelli

,³ and Albert Rimola

¹
¹Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, 08193, Catalonia, Spain
²Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy
³Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France E-mail: piero.ugliengo@unito.itE-mail: albert.rimola@uab.cat

(Accepted XXX. Received YYY; in original form ZZZ)

Abstract

Acetaldehyde (CH₃CHO) is one of the most detected interstellar Complex Organic Molecule (iCOM) in the interstellar medium (ISM). These species have a potential biological relevance, as they can be precursors of more complex species from which life could have emerged. The formation of iCOMs in the ISM is a challenge and a matter of debate, whether gas-phase, grain-surface chemistry or both are needed for their synthesis. In the gas-phase, CH₃CHO can be efficiently synthesized from ethanol and/or ethyl radical. On the grain-surfaces, radical-radical recombinations were traditionally invoked. However, several pitfalls have been recently identified, such as the presence of energy barriers and competitive side reactions (i.e., H abstractions). Here we investigate a new grain-surface reaction pathway for the formation of acetaldehyde, namely the reaction between CH₃ and a CO molecule of a dirty water/CO ice followed by hydrogenation of its product, CH₃CO. To this end, we carried out ab initio computations of the reaction occurring on an ice composed by 75% water and 25% CO molecules. We found that the CH₃ + CO_(ice) reaction exhibits barriers difficult to overcome in the ISM, either adopting a Langmuir-Hinshelwood or an Eley-Rideal mechanism. The subsequent hydrogenation step is found to be barrierless, provided that the two reacting species have the correct orientation. Therefore, this pathway seems unlikely to occur in the ISM.

keywords:

astrochemistry – molecular processes – ISM:molecules

^†^†pubyear: 2023^†^†pagerange: Quantum mechanical modeling of the grain-surface formation of acetaldehyde on H₂O:CO dirty ice surfaces–C

1 Introduction

Acetaldehyde is one of the most common detected interstellar Complex Organic Molecules (iCOMs), which are chemical compounds defined as molecules with six or more atoms that contain at least one carbon atom (Herbst & van Dishoeck, 2009; Herbst, 2017; Ceccarelli et al., 2017). Since their detection in star-forming regions, iCOMs rose the interest of scientists. There is evidence that some iCOMs formed in the Interstellar Medium (ISM) were inherited by the small objects of the Solar System (Caselli & Ceccarelli, 2012; Ceccarelli et al., 2014; Cazaux et al., 2003; Ligterink et al., 2018; Bianchi et al., 2019; Drozdovskaya et al., 2019). These species, after thermal and hydrothermal alterations, can be converted into more complex organic species (Alexander et al., 2014; Rotelli et al., 2016; Yabuta et al., 2007; Callahan et al., 2011), therefore potentially paving the way to the emergence of life on Earth.

Acetaldehyde was first detected in Sagittarius B2 by Gottlieb (1973), Fourikis et al. (1974) and Gilmore et al. (1976). Some years later, it was also detected in cold clouds, TMC-1 and L134N, by Matthews et al. (1985). This molecule has been found in a large number of environments: cold prestellar cores (Bacmann et al., 2012; Scibelli & Shirley, 2020), hot cores (Blake et al., 1987; Law et al., 2021), hot corinos (Cazaux et al., 2003; Chahine et al., 2022), protostellar molecular shocks (Lefloch et al., 2017; De Simone et al., 2020), young discs (Codella et al., 2018; Lee et al., 2019) and also in comets (Crovisier et al., 2004; Biver et al., 2021).

The presence of acetaldehyde (CH₃CHO) and other iCOMs in cold prestellar cores demonstrates that their synthesis cannot be the result of grain-surface reactions involving migration of species (i.e. radicals) other than H and O (e.g. Bacmann et al., 2012; Ceccarelli et al., 2022).

The reaction pathways that lead to iCOMs are still matter of debate, as both gas-phase and grain-surface chemistry are invoked to play a crucial role in their synthesis (e.g. Garrod et al., 2008; Balucani et al., 2015; Ceccarelli et al., 2022). Several paradigms for the formation of iCOMs were proposed in the literature, the most popular being schemes based on gas-phase reactions (e.g. Charnley et al., 1992, 1997; Taquet et al., 2016; Balucani et al., 2015; Skouteris et al., 2018; Vazart et al., 2020) and a network of radical-radical couplings occurring on the surface of grains (e.g. Garrod & Herbst, 2006; Garrod et al., 2008; Jin & Garrod, 2020). Other paradigms include a mechanism based on the condensation of atomic C (Ruaud et al., 2015; Krasnokutski et al., 2020), on the excited O-atom insertion (Bergner et al., 2017, 2019), or on the formation of HCO radical on ice surfaces as a parent precursor of other iCOMs (Fedoseev et al., 2015; Simons et al., 2020).

For what concerns acetaldehyde, several works investigated its formation routes from an experimental and theoretical point of view. In the gas-phase, the reaction between ethyl radical and atomic oxygen was proposed to yield CH₃CHO by Charnley (2004), as well as the insertion of CH into methanol (Vasyunin et al., 2017). In KIDA (Wakelam et al., 2012) and UMIST (McElroy et al., 2013) databases, an ionic route involving protonated acetaldehyde (obtained from dimethyl ether) is present. Finally, in 2018, the idea that acetaldehyde (and other iCOMs) can arise from the chemical transformation of ethanol in the gas-phase was proposed (Skouteris et al., 2018). In this chemical network, called the genealogical tree of ethanol, the latter is the precursor (the parent molecule) that give rise to different iCOMs such as glycolaldehyde, acetic acid, formic acid and acetaldehyde, among others. More recently, Vazart et al. (2020) carried out a systematic study of all gas-phase reactions present in the literature, that lead to the formation of acetaldehyde, and performed new ab initio computations for reactions having only guessed product and rate constants. These authors confirmed the pathways starting from ethyl radical (Charnley, 2004) and ethanol (Skouteris et al., 2018) and discarded the others (as they used incorrect product or rate constants). In their study, Vazart et al. (2020) also showed that the ethanol genealogical tree is currently the most promising explanation for the synthesis of acetaldehyde in warm objects. Perrero et al. (2022a) proposed that ethanol would be formed on the grain-surface by the reaction of CCH with a water molecule of the ice, followed by hydrogenation of the produced vinyl alcohol. Remarkably, the presence of frozen ethanol has been tentatively detected by James Web Space Telescope (JWST) observations (although needs confirmation) (McClure et al., 2023; Yang et al., 2022), supporting the hypothesis of the ethanol being the mother of acetaldehyde.

On the grain surfaces, experimental results are in some cases contradictory. Acetaldehyde was formed in ices containing H₂O, CO, CH₄ and CH₃OH processed by UV-irradiation (Moore & Hudson, 1998; Bennett et al., 2005; Öberg et al., 2009, 2010; Paardekooper et al., 2016; Martín-Doménech et al., 2020). In the experiments by Chuang et al. (2021), a number of iCOMs were obtained by irradiating with 200 keV H⁺ ions C₂H₂:H₂O ices, including acetaldehyde. However, in experiments using a matrix isolation technique of UV-illuminated methanol ices, where the presence of radicals could be monitored, Gutiérrez-Quintanilla et al. (2021) found that several iCOMs were formed, except acetaldehyde. Finally, the experiments by Fedoseev et al. (2022) produced acetaldehyde and its precursor, ketene, via a non-energetic pathway, where CO is co-deposited with C, H and H₂O at 10 K.

From a theoretical point of view, the radical-radical coupling of HCO and CH₃ was proposed by Garrod & Herbst (2006) to yield CH₃CHO on icy grain surfaces, as it was supposed that the reaction is barrierless due to taking place between two radical species. However, successive studies showed that not only this reaction has an energy barrier (because the radicals are physisorbed on the ice surface and have to break the intermolecular forces with the surface to react), but it also presents a competitive channel leading to the formation of CH₄ + CO (Enrique-Romero et al., 2019; Enrique-Romero et al., 2020). The same reaction simulated on a model of CO ice gave similar results, culminating in one of these three outcomes: formation of CH₃CHO, H-abstraction yielding CH₄ + CO or no reaction taking place (Lamberts et al., 2019). A recent study by Enrique-Romero et al. (2021) concluded that the efficiency of acetaldehyde formation via radical coupling is a strong function of the mobility of the radicals on the grain surface (and, consequently, of the grain temperature), in which the easier the diffusion of HCO and CH₃, the less acetaldehyde formation efficiency. Finally, recent kinetic calculations of Ben Chouikha et al. (2022) on the reaction between atomic carbon and methanol proposed by Singh et al. (2019) shows the presence of a slow reaction step preceding the barrierless formation of acetaldehyde from the radicals HCO and CH₃. However, at low temperature, the product can be formed thanks to tunnelling effects. We notice, however, that the simultaneous presence of methanol and atomic carbon on the grain-surfaces is unlike, as methanol is mostly formed by the hydrogenation of frozen CO, i.e. when carbon is prevalently locked into CO.

In summary, there is evidence that, in hot cores/corinos, acetaldehyde is formed in the gas-phase by reactions occurring in the warm gas, but it is still unclear whether this is the only mechanism at work, especially in cold environments. In light of this situation, in this work, we decided to delve into the non-energetic formation of acetaldehyde on the surface of dust grains, when radicals cannot diffuse.

When dealing with surface reactions, several mechanisms can operate: i) Langmuir–Hinshelwood (LH: Langmuir, 1922; Hinshelwood, 1930) reactions, which are efficient in the case that one of the reactants can easily diffuse on the surface of the grains (Hasegawa & Herbst, 1993); ii) Eley–Rideal (ER: Eley & Rideal, 1940) reactions, in which species from the gas-phase directly react with surface molecules, avoiding diffusion, but are efficient only if there are sufficient reactive species on the grain surface and if the reactions do not possess an activation energy (Ruaud et al., 2015); iii) Harris-Kasemo hot atom reactions, in which a high energy species has enough energy to overcome the diffusion barriers and travels on the surface until all excess energy is lost or it reacts (Harris & Kasemo, 1981); iv) reactions of suprathermal species generated by the excitation and/or ionization caused by cosmic ray bombardment (Shingledecker & Herbst, 2018; Shingledecker et al., 2018; Paulive et al., 2020).

In the present work, we investigated the acetaldehyde formation through a two-steps mechanism based on a "radical + ice" reaction, as done for the synthesis of formamide (CN + H₂O_(ice): Rimola et al., 2018)) and of vinyl alcohol/ethanol (CCH + H₂O_(ice): Perrero et al., 2022a)). Here, we propose the reaction of a methyl radical, CH₃, with a CO molecule belonging to the ice mantle of the grain, thus avoiding competitive reactions as H-abstractions, followed by a hydrogenation step:

CH₃ + CO

\to

CH₃CO

(1)

CH₃CO + H

\to

CH₃CHO

(2)

In the first step, which involves the coupling of the methyl radical CH₃ with CO_(ice), we assume that either CH₃ is adsorbed close to the CO and diffuses to react with it (LH mechanism) or that it lands directly on the CO, reacting immediately with it (ER mechanism). The second step consists in the hydrogenation of the acetyl (CH₃CO) radical so formed, which is expected to be almost barrierless as it involves a H atom reacting with a radical, whether it diffuses on the surface (LH) or it comes from the gas-phase (ER).

This work is organized as follows: in Section 2 we report the methodology (including benchmark studies), in Section 3 we present the results and Section 4 is dedicated to their discussion and comparison with other studies. Section 5 concludes the article.

2 Methodology

2.1 Computational details

We employed CRYSTAL17 (Dovesi et al., 2018) and Gaussian16 (Frisch et al., 2016) software packages, to execute periodic and molecular calculations, respectively. CRYSTAL17, at variance with other periodic codes, adopts localized Gaussian functions as basis sets (in a similar approach to that for Gaussian16), and it works with systems that range from zero to three periodic dimensions, avoiding the fake 3D replica of surface models that would arise when working with plane waves basis sets.

To reduce the computational cost of the periodic simulations, all the structures were optimized using the approximated HF-3c method followed by a single point energy calculation at the Density Functional Theory (DFT) level, hereafter referred to as DFT//HF-3c. To check the accuracy of this scheme, we performed a benchmark study, which splits in two phases: i) we compared the performance of DFT//HF-3c (CRYSTAL17) vs CCSD(T)//DFT in the gas-phase (Gaussian16) to select the best performing DFT functional, and ii) we applied the ONIOM2 correction (Dapprich et al., 1999) with the extrapolation scheme of Okoshi et al. (2015) in order to compare DFT//HF-3c against ONIOM2//DFT on the surface.

The HF-3c method was used in the most time-consuming periodic calculations, namely, geometry optimizations and frequency calculations (Sure & Grimme, 2013). HF-3c is based on the Hartree-Fock (HF) method computed with a minimal basis set (MINIX, (Tatewaki & Huzinaga, 1980), in which three empirical corrections (3c) were added: the dispersion energy D3(BJ) (Grimme et al., 2010), the basis set superposition error (BSSE) correction with the geometrical counterpoise (gCP) (Kruse & Grimme, 2012), and the short range bond (SRB) correction to fix overestimated covalent bond lengths for electronegative elements (Grimme et al., 2011; Brandenburg et al., 2013). Subsequently, the energies were refined with single point calculations onto the HF-3c optimized structures at the desired DFT level of theory, following an approach that provided accurate results in several cases, from molecular crystals (Cutini et al., 2016), polypeptides (Cutini et al., 2017), pure-silica zeolites (Cutini et al., 2019) to the computation of binding energies on crystalline and amorphous pure water ice (Ferrero et al., 2020; Perrero et al., 2022b).

In the periodic calculations, all the stationary points of a potential energy surfaces (PES) were characterized by the calculation of the harmonic frequencies at $\Gamma$ point as minima (reactants, products) or first order saddle points (transition states). Each Hessian matrix element was computed numerically by means of a 3-points formula based on two displacements of 0.003 Å from the minimum along each Cartesian coordinate. The zero-point energy correction was computed with the standard rigid rotor/harmonic oscillator formalism (McQuarrie, 1975). Since the systems are open-shell in nature, calculations were performed within the unrestricted formalism. The threshold parameters for the evaluation of the Coulomb and exchange bi-electronic integrals were set equal to 10^-7, 10^-7, 10^-7, 10^-7, and 10^-25. The sampling of the reciprocal space was conducted with a Pack–Monkhorst mesh (Pack & Monkhorst, 1977), with a shrinking factor of 2, which generates 4k points in the first Brillouin zone.

2.2 \chH2O:CO ice surface models

The H₂O:CO surfaces employed in this study were obtained from a bulk model of the pure H₂O crystalline P-ice (Casassa et al., 1997), in which one every four water molecules was replaced by CO (Figure 1). Three surfaces were cut along the planes (001), (010), (100) (see Figure 2) and then fully relaxed at HF-3c level of theory, optimizing both cell parameters and atomic positions, which led to a heavy reorganization of each structure. Although we imposed a lattice to the system through periodic boundary conditions, we modeled large unit cells of H₂O:CO dirty ice (within P1 space symmetry) whose geometry relaxation resulted in an amorphous-like surface model. The cell parameters, the dipole and the number of atoms of each structure can be found in Table 1. The (001) surface has the smallest number of atoms (192) and the smallest dipole moment across the z axis (-0.87 D at BHLYP-D3(BJ) level of theory). The H₂O and the CO molecules broke the symmetric structure of the bulk, creating clathrate-like cages in which CO molecules are surrounded by a network of water molecules engaged in hydrogen bonds (H-bonds). The outer layers of each surfaces are characterized by the presence of faintly interacting CO molecules that, depending on the coverage percentage, form a relatively complete monolayer of CO. Mixed ice models with a similar behaviour were reported in the theoretical work of Zamirri et al. (2018), in which CO adsorption, entrapment and mixture within H₂O ices was studied. There, it emerged that dispersive and quadrupolar forces have a prominent role in determining the structural features of these ice mixtures, and are responsible for the hydrophobic behaviour of CO. Indeed, when CO is entrapped in the ice structure, it causes large rearrangements of the network of hydrogen bonds, while when it is adsorbed on the top of the surfaces, it forms ordered layers. This gives rise to sections of the surface characterized by different electrostatic potential surfaces depending on the local arrangement of the carbon monoxide molecules.

Refer to caption — Figure 1: Top and lateral view of the dirty ice bulk model where one every four water molecules was replaced by CO. Color code: H, white; C, grey; O atoms belonging to water, blue; O atoms belonging to the CO, red.

Table 1: Cell parameters (a and b in Å, and

\gamma

in degrees) of the three ice surfaces optimized at HF-3c level. Dipole moment (

\mu

, in Debye) across the non-periodic z axis of the surfaces calculated at HF-3c and BHLYP-D3(BJ) (DFT) level of theory. The number of water and carbon monoxide molecules of each surface per unit cell (H₂O:CO) is also reported.

Surfaces	a	b	$\gamma$	$\mu$ (HF-3c)	$\mu$ (DFT)	H₂O:CO
(100)	14.718	12.718	101.294	-1.78	-1.72	56:24
(010)	14.289	12.478	107.261	3.13	3.10	64:16
(001)	12.776	14.604	128.338	-1.52	-0.87	48:24

In our model, the quantity of CO is not sufficient to give rise to a complete monolayer of CO, but the bottom face of (010) surface shows a neat arrangement of CO molecules. We notice very few cases in which CO molecules are engaged in short H-bonds with water molecules (< 2.4 Å). This is because the water molecules tend to form a H-bond network between them and minimize the number of dangling H atoms pointing towards CO. Moreover, the majority of H-bond interactions take place through the O-end of the CO molecule and not through the C-end.

This phenomenon is probably driven by the small energy difference between the two interactions at HF-3c level of theory, for which the H-bond through C atom is favoured over the H-bond through the O atom by only 0.3 kJ mol^-1 (see Figure 9). The C $\cdots$ H interaction is also longer (2.390 Å) than the O $\cdots$ H one (2.212 Å), although the results are comparable with those given by B3LYP-D3(BJ)/6-311G(d,p).

3 Results

3.1 Preliminary benchmark studies

3.1.1 Gas-phase benchmark calculations

In order to select the most appropriate DFT level of theory with which executing the periodic simulations, we firstly performed a benchmark analysis of the reactions 1 and 2 in the gas-phase (namely, in the absence of the whole icy surfaces). For the sake of reliability of the benchmark study adopting these gas-phase reactions, geometry optimizations were performed at HF-3c level in CRYSTAL17, while the subsequent single point energies were computed with six different DFT methods, which were corrected with the Grimme’s D3 or, when available, the D3(BJ) terms (Grimme et al., 2010, 2011) to account for dispersion interactions. Thus, the employed methods are: BHLYP-D3(BJ) (Becke, 1993a; Lee et al., 1988), B3LYP-D3(BJ) (Lee et al., 1988; Becke, 1988, 1993b), M062X-D3 (Zhao & Truhlar, 2008), M052X-D3 (Zhao et al., 2006), MPWB1K-D3(BJ) (Zhao & Truhlar, 2004) and $\omega$ B97X-D3 (Chai & Head-Gordon, 2008). They were used in conjunction with the Pople-based 6-311G(d,p) basis set. Diffuse functions were neglected because their small exponents can cause high linear dependencies in the wavefunction of the periodic systems we aim to work with. (Klahn & Bingel, 1977; VandeVondele & Hutter, 2007; Peintinger et al., 2013; Vilela Oliveira et al., 2019)

The DFT//HF-3c results were compared with those obtained by performing optimizations at the same DFT/6-311G(d,p) levels, followed by single point energy calculations at single- and double- electronic excitation coupled-cluster method with an added perturbative description of triple excitations (CCSD(T)) combined with Dunning’s aug-cc-pVTZ basis set in Gaussian16 (Raghavachari et al., 1989). The reaction path is represented in Figure 3. Reaction 1 presents a barrier due to being the coupling of the radical CH₃ with the closed-shell CO. Reaction 2 is a radical-radical coupling and, therefore, is barrierless in the gas-phase.

Table 2: Relative energies (in kJ mol^-1) of the stationary points of the gas-phase reactions 1 and 2 for all the DFT/6-311G(d,p) methods (computed with CRYSTAL17) and the CCSD(T)/aug-cc-pVTZ method (computed with Gaussian16). R stands for isolated reactants, PRE-R for pre-reactant complex, TS for transition state and P for product. The data are not ZPE-corrected.

Reaction	Step	A//H^†	B//H^†	C//H^†	D//H^†	E//H^†	F//H^†	CCSD(T)//A^†
R1: CO + CH₃ $\to$ CH₃CO	R	0.03	-0.03	4.84	6.58	1.32	1.35	2.03
	PRE-R	0.0	0.0	0.0	0.0	0.0	0	0.0
	TS	17.9	0.2	9.4	15.6	4.8	3.5	22.5
	P	-54.3	-66.4	-58.0	-51.1	-74.9	-71.6	-56.1
	err% TS R1	20.6%	99.3%	58.4%	30.9%	78.6%	84.3%
	err% P R1	3.1%	18.5%	3.4%	8.9%	33.5%	27.7%
	mean err% R1	11.9%	58.9%	30.9%	19.9%	56.1%	56.0%
R2: CH₃CO + H $\to$ CH₃CHO	R	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	P	-390.7	-388.0	-394.8	-396.8	-383.0	-402.6	-396.1
	err % P R2	1.4%	2.0%	0.3%	0.2%	3.3%	1.6%

^† A=BHLYP-D3(BJ), B = B3LYP-D3(BJ), C = M052X-D3, D = M062X-D3, E = MPWB1K-D3(BJ), F =

\omega

B97X-D3, H = HF-3c

The results of this benchmark study for the gas-phase reaction models are shown in Table 2. The different CCSD(T)//DFT calculations return comparable energy values, meaning that results are almost unaffected by the small differences in the optimized geometries obtained with the different DFT methods. Thus, we only report here the results computed at CCSD(T)//BHLYP-D3(BJ), which can be taken as reference values. The error percentages on the energy barriers are high when comparing the DFT//HF-3c results with the CCSD(T)//BHLYP-D3(BJ) ones. To have deeper insights into that, we have analysed the resulting optimized structures with each method. Table 3 reports the pivotal bond lengths of each species, where one can observe that HF-3c does not reproduce well the structures for the pre-reactant (PRE-R) and the transition state (TS) structures compared with the DFT methods. In the PRE-R and in the TS, the \chC\bondsingleC length at HF-3c versus any DFT level shows differences going from 0.1 Å to 0.3 Å and, indeed, have a decisive impact on the energies that are computed afterwards at single point DFT.

Table 3: Distances (in Å) of selected bonds in the structures of the pre-reactants (PRE-R), transition state (TS) and products (CH₃CO and CH₃CHO) of the gas-phase reaction models 1 and 2 computed at different theory levels.

Structure	BHLYP-D3(BJ)	B3LYP-D3(BJ)	M052X-D3	M062X-D3	MPWB1K-D3(BJ)	$\omega$ B97X-D3	HF-3c
CO	1.114	1.127	1.120	1.121	1.116	1.125	1.135
\chC\bondsingleC (PRE-R)	3.369	3.190	3.211	3.206	3.431	3.294	3.154
\chC\bondsingleC (TS)	2.146	2.238	2.171	2.147	2.243	2.209	2.015
\chC\bondsingleC (CH₃CO)	1.505	1.515	1.513	1.517	1.498	1.513	1.552
\chC\bondsingleO (CH₃CO)	1.166	1.180	1.174	1.173	1.167	1.177	1.187
\chC\bondsingleH (CH₃CHO)	1.513	1.114	1.106	1.110	1.105	1.113	1.107
\chC\bondsingleO (CH₃CHO)	1.177	1.204	1.200	1.199	1.191	1.200	1.208

Energetic results indicate that B3LYP-D3(BJ) is the least suitable functional (99.3% error on the potential energy barrier for Reaction 1), while BHLYP-D3(BJ) is the best one (20.6% error though, corresponding to an absolute error of -4.6 kJ mol^-1). The BHLYP-D3(BJ) method also gives the smallest error percentage in the reaction energies for the formation of the CH₃CO radical. On the other hand, every functional describes well the reaction of CH₃CO with H.

3.1.2 Grain-surface benchmark calculations

Based on the error percentage on the energy of the gas-phase TS, we chose BHLYP-D3(BJ) to compute the reactions on the surface. However, although BHLYP-D3(BJ)//HF-3c performs globally well compared with CCSD(T)//BHLYP-D3(BJ), one has to pay special care when modeling the transition state of Reaction 1, which is the pivotal step to determine whether the reaction is feasible or not in the ISM. As probably the same discrepancies can affect the periodic calculations, we also benchmarked them to assess the quality of the employed methods. Indeed, it is clear that we need to improve the quality of our data directly on the periodic surface, where the interaction between the reactants and the ice could change again the geometrical features of the interaction between CH₃ and CO, and therefore draw more discrepancies between HF-3c, BHLYP-D3(BJ) and CCSD(T) results.

A possible solution to this problem should be performing ONIOM2 calculations (Dapprich et al., 1999) combining BHLYP-D3(BJ) with CCSD(T), as low and high energy levels. This methodology has been previously applied to the computation of binding energies (Ferrero et al., 2020; Perrero et al., 2022b). However, this would require obtaining optimized structures at BHLYP-D3(BJ) of the full periodic systems, which in our case is not feasible due to the high number of atoms present in the unit cells. Instead, we selected three test model cases onto which simulating the acetaldehyde formation and computing the energies applying the ONIOM2 refinement, in order to compare the performance of BHLYP-D3(BJ)//HF-3c against the ONIOM2-corrected values on the reaction barrier.

The test cases are based on two pure H₂O crystalline periodic ice models, in which one water molecule is replaced by a CO, and a molecular cluster of pure CO ice. They are shown in Figure 8 and are represented by: the 2x1 supercell of the (010) P-ice surface in which i) a water molecule exposing a dangling oxygen atom (\chH2O Ice (a)) and ii) a water molecule exposing a dangling hydrogen atom (\chH2O Ice (b)) was substituted by a CO molecule, and iii) a cluster model made of 20 CO molecules.

Each system was divided in two parts (model and real systems), described by two different levels of theory (high and low). The model system (represented by the CH₃ and the CO) was described by the high level of theory, CCSD(T). The real system (that is, the whole system) was described by the low level of theory, BHLYP-D3(BJ). In this ONIOM correction, CCSD(T) was used in combination with the Dunning’s aug-cc-pVNZ (with N = D,T) basis sets (Dunning, 1989) and, with these data, the OAN(C) extrapolation scheme to the complete basis set (CBS) limit was applied (Okoshi et al., 2015).

The ONIOM2-corrected energy barrier ( $\Delta$ E_TS) was computed as

\Delta E_{TS(ONIOM2)}=\Delta E_{TS}(low,real)+\\ \Delta E_{TS}(high,model)-\Delta E_{TS}(low,model)

(3)

where the $\Delta$ E_corr = $\Delta$ E_TS(high,model) - $\Delta$ E_TS(low,model) represents the correction term to the energy of the real system.

In this work, for the calculation of the ONIOM2-corrected barriers, $\Delta$ E_ONIOM2, equation 3 can be rewritten as:

\Delta E_{TS(ONIOM2)}=\Delta E_{TS}(DFT;all)+\\ \Delta E_{TS}(CCSD(T)/CBS;fragm)-\Delta E_{TS}(DFT;fragm)

(4)

where $\Delta$ E_TS(DFT; all) is the activation barrier computed at DFT//DFT. The $\Delta$ E_TS of the model system (CH₃ + CO; fragm) is computed through single point energy calculations at CCSD(T)/aug-cc-pVNZ with n = D,T and extrapolated to the CBS limit thanks to the OAN(C) equation.

E_{CBS}^{OAN(C)}=\frac{3^{3}E(T)-s^{3}E(D)}{3^{3}-s^{3}}

(5)

In this equation, s=2.091 based on the choice of method and basis set, E(T) is the energy calculated with the aug-cc-pVTZ basis set and E(D) corresponds to that computed with the aug-cc-pVDZ basis set.

We performed Reaction 1 on the three test models, both at BHLYP-D3(BJ)//HF-3c level of theory and at the ONIOM2 scheme, chosen as reference. Data are available in Table 4.

Table 4: Potential energy barrier (in kJ mol^-1) of Reaction 1 computed on the three test cases for the grain-surface benchmark. DFT stands for BHLYP-D3(BJ)/6-311G(d,p)

Structure	real		model		final
	DFT//HF-3c	DFT//DFT	DFT	CBS	ONIOM2
H₂O Ice (a)	27.0	49.5	13.8	19.6	55.2
H₂O Ice (b)	15.1	40.1	9.0	2.3	33.3
CO Ice	11.3	15.0	18.0	20.4	17.4

In the plot of Figure 4, we compare the ONIOM2//DFT results against the DFT//HF-3c ones. The regression line (DFT//HF-3c=0.491 $\cdot$ ONIOM2//DFT, R²=0.932) shows that the method of choice for the periodic calculations, DFT//HF-3c, underestimates the reaction barrier, as already emerging from the gas-phase benchmark study.

Since we have very few cases, we do not aim to adopt the slope as a correction factor. On the other hand, computing the reactions on the dirty ice surface models at full BHLYP-D3(BJ) level is almost unpractical. This benchmark study, however, allows to figure out the error bar associated with the present periodic calculations.

3.2 Grain-surface reactions adopting an LH mechanism

Once we have checked and chosen a reasonably suitable methodology to calculate Reactions 1 and 2, we simulated them on the dirty H₂O:CO icy surface models.

To this end, we first adsorbed the CH₃ on the ice models by manually placing the radical in positions characterized by different local environments. Only the atomic positions were relaxed, while the cell parameters were kept frozen, to be consistent during the successive steps of the reaction and avoid structural deformations of the ice surface models. To calculate the binding energy (BE) of CH₃ on the mixed H₂O:CO surfaces, we followed the same computational scheme as in Ferrero et al. (2020) and Perrero et al. (2022b), that is, by correcting the adsorption energy $\Delta$ E_ads = E_complex - E_ice - E ${}_{CH_{3}}$ for the basis set superposition error (BSSE, which generates from using a finite basis set) through the counterpoise method by Boys & Bernardi (1970).

BE=-\Delta E_{ads}+BSSE

(6)

For each case we computed the deformation energies of the surface and of the radical to verify that heavy structural rearrangements were not affecting our models during the geometry optimizations. We also computed the BE(0) at 0K, by correcting the BE for the $\Delta$ ZPE as in equation 7.

BE(0)=BE-\Delta ZPE

(7)

To identify the transition states, we adopted the distinguished reaction coordinate (DRC) procedure by performing a scan calculation along the \chC\bondsingleC length. The maximum energy structure of the DRC pseudo-PES was used to localize and optimize the actual TS structure (as implemented in the CRYSTAL code, (Rimola et al., 2010)). In the DRC process, we avoided using internal redundant coordinates. Instead, we selected a mixed coordinate system made up of a set of selected valence internal parameters (bond lengths, angles and dihedrals) plus the full set of symmetry adapted fractional displacements and elastic distortions (evoked by the keyword INTLMIXED in CRYSTAL17). This represents an advantage both in reducing the number of valence internal parameters that are automatically generated for the system and in solving the quasi-linear dependencies that arises in systems with a lack of connectivity, which are reflected in a very high condition number of the Wilson B-Matrix. The latter condition usually makes the optimization either to fail or to exhibit an erratic behavior (Dovesi et al., 2018).

To simulate Reaction 2, we selected one newly formed CH₃CO/ice complex, and we manually placed the hydrogen atom to simulate its adsorption. We first set up five starting geometries and optimized the structures as open-shell triplets (the two unpaired electrons with the same sign), which were then optimized as open-shell singlets (the two unpaired electrons with opposite signs). When the orientation of H was in favor of the formation of the \chC\bondsingleH bond, but no acetaldehyde was obtained spontaneously, we performed DRC calculations to characterize the PES of the process.

3.2.1 Adsorption of methyl radical on the dirty ices

For each surface, CH₃ was adsorbed in four different positions (two on the top-face and two on the bottom-face) to sample different binding sites. The adsorption structures at the (100) surface are shown in Figure 5. Table 5 summarizes the computed BEs and their contribution, along with the adsorption enthalpy BE(0). The same Table 5 presents the nomenclature adopted for each binding site.

Table 5: Binding energy (BE) values (in kJ mol^-1) of CH₃ on different surface sites of the dirty H₂O:CO ice model. The contributions from the pure potential energy values (

\Delta E_{ads}

), the BSSE corrections (

\Delta BSSE

), the dispersive (Disp) and the non-dispersive (No Disp) terms of the BE, the zero point energy corrections (

\Delta ZPE

) and the resulting adsorption enthalpy (BE(0)) are shown.

Surface	Binding Site	$\Delta E_{ads}$	$\Delta BSSE$	BE	Disp	No Disp	$\Delta ZPE$	BE(0)
(100)	H0	-15.2	-5.9	9.3	15.0	-5.7	4.0	5.3
	H1	-18.8	-8.9	9.8	17.8	-8.0	4.7	5.1
	H2	-17.4	-6.8	10.7	13.3	-2.6	3.7	7.0
	H3	-22.2	-10.8	11.3	20.2	-8.8	4.7	6.7
(010)	K0	-27.3	-13.4	14.0	19.8	-5.8	4.5	9.5
	K1	-14.5	-6.2	8.3	12.9	-4.6	4.4	3.9
	K2	-31.8	-13.4	18.4	18.8	-0.4	4.7	13.7
	K3	-29.4	-11.0	18.5	16.9	1.5	4.3	14.1
(001)	L0	-13.8	-7.3	6.5	6.7	-0.2	5.5	1.0
	L1	-13.0	-7.1	5.9	27.9	-22.0	6.6	-0.7
	L2	-18.9	-7.6	11.3	14.7	-3.3	3.9	7.4
	L3	-40.3	-6.9	33.4	24.5	8.8	4.3	29.1

The methyl radical possesses an unpaired electron localized on the carbon atom, and its electrostatic potential surface is neutral almost everywhere. Therefore, this species will not form strong electrostatic interactions with the surface, especially with the polar water molecules. Thus, dispersion interactions are the key to explain the behavior of the methyl onto the surface. The BSSE can be as large as 50% of the BE, due to the fact that the 6-311G(d,p) basis set has a rather small number of functions.

The BSSE-non-corrected adsorption energies vary between -13 and -40.3 kJ mol^-1. Considering the BSSE, the resulting BEs range from 5.9 to 33.4 kJ mol^-1. In nine out of the twelve characterized binding sites, the BEs are between 6 and 14 kJ mol^-1, while in two cases a value of 18.5 kJ mol^-1 is found. These values are similar to those obtained by Ferrero et al. (2020), where the BE computed on the crystalline water ice is 18.2 kJ mol^-1, and on the amorphous water span the 9.2 to 13.8 kJ mol^-1 range, indicating a similar interaction of methyl radical with the two ice models. However, in the case of the binding site referred to as L3, the BE overcomes 30 kJ mol^-1. This is the only adsorption complex in which the interaction would be attractive even if we were not accounting for the dispersion. Indeed, while in all other cases the interaction would be repulsive, here the hydrogen atoms of the CH₃ are sufficiently close to the oxygen atoms of both CO and H₂O, therefore determining an advantageous electrostatic interaction of 8.8 kJ mol^-1. The ZPE correction is between 4-5 kJ mol^-1 in all the cases, as we would expect from the limited rearrangement of both the surface and the radical upon adsorption.

3.2.2 Acetyl radical formation

We modeled the formation of acetyl on a selected number of structures (see Figure 10 for an example). For each surface, we chose the three adsorption complexes characterized by the largest differences in the chemical environment of CH₃. The purpose of this choice is to probe the effect of the chemical environment on the potential energy barrier of the reaction. According to the LH mechanism, the CH₃ diffuses towards the closest CO on the surface to form the chemical bond. In all the nine cases analyzed (see Table 6), a potential energy barrier has to be overcome in order to form the acetyl radical. The average barrier is around 10-15 kJ mol^-1, with some exceptions. On the (100) surface, the complex H3 has a barrier of only 5.8 kJ mol^-1. We notice that in this case, in the scan calculation, it is the CO that approaches the CH₃ and not vice versa, as we observed in the other simulations. We suppose that the particular geometry of this surface allows an easy diffusion of the CO, turning into a low potential energy barrier. On the other hand, on the (010) surface, we found two unfavourable mechanisms, K2 and K3, presenting energy barriers up to 36.8 kJ mol^-1. In K2, a H-bond between the carbon of the CO and a water molecule of the surface needs to break to make CO available for the reaction with CH₃. In K3, the large distance (3.9 Å) between the reactants is the responsible of the high barrier.

When adding the ZPE contributions, according to the equation $\Delta$ H_TS = $\Delta$ E_TS + $\Delta$ ZPE, each barrier increases by about 11 kJ mol^-1. For the gas-phase reaction, we found barriers of $\Delta$ E_TS = 17.9 kJ mol^-1 and $\Delta$ H_TS = 28.7 kJ mol^-1, meaning a $\Delta$ ZPE_TS = 10.8 kJ mol^-1, very close to that obtained on the surface. If we focus solely on the energy barrier, in seven cases out of nine the barriers of the surface reactions are lower than the gas-phase one. However, the grain-surface benchmark calculations warn us that these barriers are underestimated and, therefore, it is highly probable that only the grain-surface H3 case is slightly more favourable than the gas-phase reaction.

Table 6: ZPE-corrected potential energy barriers (

\Delta

H(0)_TS, in kJ mol^-1) for the formation of acetyl on different binding sites of each surface. The internal potential energy values (

\Delta E_{TS}

) and the zero point energy corrections (

\Delta ZPE_{TS}

) are displayed.

Surface	Binding Site	$\Delta E_{TS}$	$\Delta ZPE_{TS}$	$\Delta H(0)_{TS}$
(100)	H1	11.7	10.7	22.4
	H2	13.4	11.3	24.7
	H3	5.8	11.0	16.8
(010)	K1	11.6	10.3	21.9
	K2	23.7	10.0	33.7
	K3	36.8	10.1	46.9
(001)	L1	10.8	10.8	21.6
	L2	14.5	12.3	26.8
	L3	14.4	11.9	26.5
Gas-phase	-	17.9	10.8	28.7

3.2.3 Acetaldehyde formation

Among the different acetyl products obtained, we selected the H1 case to model Reaction 2. The reason of this choice is because, in this structure, the acetyl is bound to a water molecule through a H-bond involving the oxygen of the carbonyl group, namely, the carbon atom of interest is not hindered because it is pointing towards the gas-phase (it is not facing the inner side of the surface), and thus it is prone to react with a hydrogen atom diffusing on the surface (see Figure 11).

The adsorption of atomic hydrogen followed by the optimization of the CH₃CO + H complex in an open-shell triplet spin state resulted, as expected, in no reaction. However, a change of the electronic spin state to an open-shell singlet brought to the spontaneous formation of acetaldehyde in two (out of five) complexes without harming the H-bond interaction between the carbonyl moiety and the H₂O_(ice). We noticed that if the H atom is at a maximum distance of 3.5 Å from the CO moiety and it is correctly oriented (meaning that it approaches the CO from the less hindered side), the radical coupling is barrierless.

In the other cases, either the H atom is hindered by the methyl moiety, or it has to overcome a small diffusion barrier to get close enough to the reactive center. With the methodology chosen for the optimization, we notice that a very small diffusion barrier, of about 1 kJ mol^-1, is limiting the free diffusion of the hydrogen on the surface, as the gradient in the optimization process falls to zero and a minimum is found.

3.3 Grain-surface reactions adopting an ER mechanism

In the ER mechanism, the CH₃ gas-phase species directly reacts with a surface CO icy component. To simulate this, we generated, with a Python script, several geometries in which the distance between the approaching CH₃ from the gas-phase and the reactive CO of the surface progressively decreases. We ran the optimization of these structures without relaxing the geometry of the surface to avoid its thermalization during the approach of the CH₃. Only after the product forms, the entire system was relaxed. Therefore, we provide the energy profiles for Reactions 1 and 2 without characterizing the actual transition states but only giving an estimate. This is because, due to the geometrical constraints applied to the surface, we would obtain high order stationary points with imaginary frequencies associated with the motion of the frozen atoms of the surface. For the same reason, we also do not provide the ZPE-corrected profiles for these processes. We simulated both steps of the reaction assuming an ER mechanism. We chose the (100) surface as the reactive one, focusing on the H1 binding site. We did so in order to compare LH and ER mechanisms for steps 1 and 2, but also because the CO involved in the \chC\bondsingleC bond formation is exposing the C-end towards the gas-phase, therefore being available for this type of reactivity. The energy profile of CH₃ + CO is plotted in Figure 6. Both the profile computed at HF-3c level and the one refined at BHLYP-D3(BJ) level are shown. One can clearly see the disagreement between the two methodologies in terms of energetics. The BHLYP-D3(BJ) profile shows that below 5 Å of distance there is an attractive interaction between the methyl radical and the CO. The presence of a van der Waals complex well in the energy profile appears at a distance of 3.3 Å. When getting close to 2 Å, we find the maximum energy point, whose value is about 8.5 kJ mol^-1 compared with the energy at infinite distance. The $\Delta$ E between the minimum and maximum energy point is 16 kJ mol^-1, which is compatible with the barrier computed for the LH mechanism. According to these results, with both the formation of such complex and the presence of a barrier, the reaction will hardly yield acetyl.

On the other hand, considering that acetyl is already on the surface and a hydrogen atom approaches from the gas-phase, two situations can take place. If the overall spin of the system is a triplet, no coupling takes place and the energy of the system becomes highly repulsive as the distance between the two reactants increase. In contrast, if the system is in a singlet open-shell electronic state, the formation of the \chC\bondsingleH bond occurs spontaneously. We also simulated the competitive reaction, in which H falls onto the surface and reacts with an adjacent surface CO molecule (considering an overall triplet spin state) forming the HCO radical with a small barrier of 4.5 kJ mol^-1. Thus, in case of a H atom approaching the surface, the formation of acetaldehyde via H addition to acetyl is the most probable reaction.

4 Discussion and Astrophysical Implications

In this work, we simulated the formation of acetaldehyde on H₂O:CO dirty ice surfaces models through the reaction of a methyl radical with a CO belonging to the ice surface, in which, subsequently, the newly formed acetyl radical gets hydrogenated to yield acetaldehyde. We performed these two-steps process on the basis of both LH and ER surface mechanisms.

Thanks to this "radical + ice" reaction mechanism, step 1 has no competitive reactions, at variance with the prevailing radical-radical coupling mechanism (Enrique-Romero et al. (2022)). Indeed, in our case, CH₃ cannot extract a hydrogen atom from water to yield CH₄, neither we consider such a high abundance for CH₃ to be able to form ethane (CH₃CH₃). As far as Reaction 2 is concerned, there may be a competitive channel, i.e., H + CO $\rightarrow$ HCO, known to have a barrier of $\sim 5$ kJ mol^-1 and to occur through tunnelling (e.g. Rimola et al., 2014; Pantaleone et al., 2020). In our simulations, we also observe the presence of a barrier. Even though we modeled the process through an ER mechanism and we did not isolate the exact transition state of this reaction, the energy profile reaches its maximum of energy at 4.5 kJ mol^-1. In contrast, the formation of acetaldehyde is barrierless and, therefore, favoured.

The benchmark that we performed for the acetyl formation reaction on the icy surfaces stresses that the BHLYP-D3(BJ)//HF-3c computational scheme seems to underestimate the potential energy barrier of the process. Thus, a direct comparison with the gas-phase reaction (computed at full DFT level) is, therefore, not possible. However, if we roughly assume that the computed $\Delta$ E(0)_TS values on the surfaces have to be doubled to be fairly compared with the gas-phase ones (as emerged from the benchmark study), it turns out that the grain-surface reactions are less favourable than the gas-phase ones in the majority of the cases, namely, just one case out of nine is characterized by a $\Delta$ H(0)_TS lower than the 28.7 kJ mol^-1 found for the gas-phase reaction. However, even in the most favourable case, the barrier is still high so that Reaction 1 is unlikely to take place in the cold molecular clouds. Additionally, in environments where the temperature can be higher (namely, with more thermal energy available to overcome the barrier), the sublimation of CO and CH₃ (above 20 K: e.g. Ferrero et al., 2020) will prevent the reaction from taking place.

Therefore, we can summarize the effect of the ice surface on the acetyl formation through the reaction of CH₃ with iced CO in two points: i) the CO embedded in the surface is not strongly activated towards the reaction by the surrounding icy water molecules, which in fact, if possible, avoid interacting with the CO molecules by creating clathrate-like structures, as it appears during the optimization of the surfaces; and ii) in almost all the surface reactions, the potential energy barrier is fairly larger than that in the gas-phase, meaning that there is an additional, although modest, interaction between the CH₃ radical and the surface that needs to be broken in order to allow the radical to approach the CO and form a \chC\bondsingleC bond.

On the other hand, the second step (formation of acetaldehyde by H-addition to acetyl) is favoured by the presence of the ice, as the surface is well known to act as a reactant concentrator (e.g. Ioppolo et al., 2011; Rimola et al., 2014; Fedoseev et al., 2015; Simons et al., 2020; Ferrero et al., 2023b) as well as third body (i.e., energy dissipator, hence stabilising the newly formed product) (e.g. Pantaleone et al., 2020; Pantaleone et al., 2021; Ferrero et al., 2023a; Molpeceres et al., 2023) in hydrogenation of atoms and small molecules. Moreover, in the event that hydrogen is approaching the acetyl with an unfavourable orientation, thanks to its easy diffusion, it is able to move towards the CO moiety and yield acetaldehyde.

The ER mechanism does not bring any improvement over the LH one. For the former, to be effective, the acetyl formation should not have a barrier, this way avoiding the preliminary adsorption and forming directly the product. Ruaud et al. (2015) stresses the importance of this mechanism for reactions between atomic carbon and icy components. However, in our case, both mechanisms present a barrier with similar heights, so that none of the mechanisms dominate over the other (at least from an energetic standpoint). For the hydrogenation of the acetyl radical to form acetaldehyde, both mechanism are feasible to yield the product, although the LH one is probably favoured due to the fact that at a temperature as low as 10 K the majority of atomic hydrogen would probably be adsorbed on the surface. We have seen, moreover, that the relative orientation between H and acetyl is not a hampering factor so that probably, at 10 K, H atoms have enough mobility to jump between different adsorption sites and find the orientation that would favour the formation of acetaldehyde.

Therefore, although the CH₃ + CO_(ice) mechanism seems not to particularly enhance the formation of acetaldehyde, this result is important to constraint the active synthetic paths that drive the formation of acetaldehyde in the ISM.

Other studies have focused on acetaldehyde grain-surface formation through radical recombination between HCO and CH₃ (e.g. Enrique-Romero et al., 2016, 2019, 2021; Lamberts et al., 2019), which already highlighted the difficulty of the synthesis of this iCOM. In Enrique-Romero et al. (2016) and Enrique-Romero et al. (2019), the reactivity between CH₃ and HCO performed on H₂O ice resulted in either formation of acetaldehyde or in CH₄ + CO due to a competitive H abstraction, both reactions having small or no potential energy barriers. In a subsequent work, Enrique-Romero et al. (2021) found that the efficiency of acetaldehyde formation was overall low in comparison with that of CH₄ + CO. In Lamberts et al. (2019), the radical pair reactivity was investigated on CO ices by means of ab initio molecular dynamics. By adopting a sufficient configurational sampling (namely, different trajectories based on different initial guess structures), the authors found that the reactivity results in either no reaction, formation of CH₃CHO or CH₄ + CO, the last two outcomes being barrierless, while in the first case the non-formation of a product was due to the presence of a barrier. Therefore, according to these results, acetaldehyde synthesis on icy surfaces is not a favourable path. This is in agreement with the theoretical study by Simons et al. (2020), who stated that the reaction network obtained by the hydrogenation of CO may cause the production of several iCOMs (glycolaldehyde, ethylene glycol and to a less extent methyl formate), alongside methanol and formaldehyde, but acetaldehyde is not one of them. Likewise, the experimental study by Gutiérrez-Quintanilla et al. (2021) found that acetaldehyde is not formed by the combination of radicals created by the UV illumination of methanol ice.

An alternative way of forming acetaldehyde on the grain-surfaces discussed in the literature is via the reaction of carbon atoms with the CO molecules of the ices. Fedoseev et al. (2022) experimentally studied the reaction of C + CO co-deposited with H₂O and H at 10 K, following the theoretical study of Papakondylis & Mavridis (2019). The authors observed the formation of CCO and its hydrogenated counterpart, the ketene CH₂CO. Successive hydrogenation steps result in CH₃CHO. Likewise, Ferrero et al. (2023b) carried out a theoretical study on the formation of ketene from the reaction of C with CO_(ice) and the potential successive reactions with radicals (e.g. OH and NH₂) that could form other iCOMs and found that only hydrogenation (via H-tunneling), eventually leading to acetaldehyde could occur, because of important energy barriers. However, as Ferrero et al. (2023b) noticed, the simultaneous presence of gaseous atomic carbon landing on the grain-surfaces and a CO-rich ice is very unlikely in the molecular ISM, except in Photo-Dissociation Regions (PDRs), as abundant frozen CO implies an evolved molecular cloud or prestellar core where atomic carbon has an extremely low abundance.

Assuming that CH₃CHO is formed on the grain-surfaces, then a thermal and/or non-thermal mechanism is needed to partially transfer it into the gas-phase, where it is observed. This is particularly problematic when considering the detection of acetaldehyde in several cold prestellar cores (Scibelli & Shirley, 2020). The parameter establishing whether a species stays bound to the icy mantle or enriches the gas-phase is its BE. A recent computational work by Ferrero et al. (2022) highlights how the acetaldehyde BE, ranging from 3000 K to 7000 K on water ice, represent an obstacle to explain its presence in the gas-phase of cold environments. The studies by Corazzi et al. (2021) and Molpeceres et al. (2022) reported a value which is close to the lower end of the distribution outlined in Ferrero et al. (2022), but it is still too large to allow desorption in cold (10 K) astronomical objects. Usually, to explain the desorption process of this and other iCOMs, non-thermal mechanisms such as photo-desorption, reactive desorption, and cosmic-ray desorption are invoked (e.g. Dulieu et al., 2013; Chuang et al., 2018; Dartois et al., 2019), each of these mechanism having their drawbacks (e.g. Bertin et al., 2016; Pantaleone et al., 2020). Clearing this matter is out of the scope of this work, which aims at investigating acetaldehyde formation on icy surfaces.

In summary, all studies so far carried out, theoretical and experimental, tend to agree that acetaldehyde is unlikely to be a grain-surface product. On the contrary, several studies now seem to favor the hypothesis of acetaldehyde formed in the gas phase. We would like to mention in particular the work by Vazart et al. (2020), who showed a very good agreement between the observed and measured abundance of acetaldehyde in hot corinos.

5 Conclusions

In this work, we studied the formation of acetaldehyde (CH₃CHO) on the grain surfaces through a two-steps mechanism consisting of: i) the reaction of a methyl (CH₃) radical with a CO molecule belonging to a periodic surface of H₂O:CO dirty ice to form acetyl (CH₃CO), and ii) the hydrogenation of the newly formed acetyl. We characterized the process via the Langmuir-Hinshelwood and Eley-Rideal mechanisms and compared the results obtained against those simulated for the same reaction in the gas-phase. Our grain-surface benchmark on three test models stressed that the barriers computed at BHLYP-D3(BJ)//HF-3c level are fairly underestimated against CCSD(T)//BHLYP-D3(BJ).

We modeled three periodic dirty ice surfaces by cutting along the planes (001), (010), and (100) a bulk of crystalline H₂O P-ice, where one fourth of the water molecules was substituted by CO. We obtained three amorphous-like surfaces characterized by a variable content in CO onto which we simulated several adsorption structures of CH₃, finding that dispersion interactions are crucial for the adsorption process. We computed the ZPE-corrected potential energy surfaces of acetyl formation in three cases for each H₂O:CO ice surface, results pointing out that in most of the cases the reaction presents high energy barriers insurmountable in the ISM. The formation of acetyl via CH₃ + CO_(ice) on the surface is less favourable than its gas-phase counterpart in the majority of the cases.

We have elucidated that these overall unfavourable reactions are due to that: i) the CO molecule is hardly activated by the surface and the CH₃ does not present an enhanced reactivity (it does not form hemibonded systems), and ii) the barriers for acetyl formation depend on the local environment of the CO, that is, they are high (in most of the cases) when the distance between the reactants is large or when the CO/H₂O interactions have to be broken for the reaction to take place. On the other hand, the hydrogenation of acetyl is barrierless, given that the electronic spin of the system is a singlet open-shell and that the H atom is so mobile that it always finds a proper orientation towards the carbonyl group of CH₃CO. The outcome of the reaction, specifically its barrier, does not improve when adopting the ER mechanism against the LH one.

In summary, several previous studies have challenged the hypothesis that acetaldehyde is formed on the grain-surfaces via the combination of on-surface radicals, specifically CH₃ and HCO (e.g. Enrique-Romero et al., 2016, 2019; Lamberts et al., 2019; Simons et al., 2020; Enrique-Romero et al., 2021; Gutiérrez-Quintanilla et al., 2021). Leveraging on other studies showing that some iCOMs (ethanol and formamide: Perrero et al., 2022a; Rimola et al., 2018) could be formed on the grain surfaces via reactions of radicals with molecules belonging to the ice, in this work, we investigated the reaction of the CH₃ radical with one CO molecule of the ice. Our computations show that also this path is unfavorable to the acetaldehyde formation, reducing the possibilities that the latter is formed on the grain surfaces. Alternatively, gas-phase reactions could be at the origin of the almost ubiquitous presence of acetaldehyde in the ISM (e.g. Vazart et al., 2020).

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 “Quantum Chemistry on Interstellar Grains” (QUANTUMGRAIN), grant agreement No 865657 and “The Dawn of Organic Chemistry” (DOC), grant agreement No 741002. The authors acknowledge funding from the European Union’s Horizon 2020 research and innovation program Marie Sklodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement No 811312. MICINN (project PID2021-126427NB-I00) and Italian MUR (PRIN 2020, Astrochemistry beyond the second period elements, Prot. 2020AFB3FX) are also acknowledged. CSUC supercomputing center is acknowledged for allowance of computer resources. We thank Prof. Gretobape for fruitful and stimulating discussions. J.P. is indebted to Joan Enrique-Romero for his help in python scripting.

Data Availability

The data underlying this article are freely available in Zenodo at https://doi.org/10.5281/zenodo.7937759.

References

Alexander et al. (2014) Alexander C. M. O., Cody G. D., Kebukawa Y., Bowden R., Fogel M. L., Kilcoyne A. L. D., Nittler L. R., Herd C. D. K., 2014, Meteorit. Planet. Sci., 49, 503
Bacmann et al. (2012) Bacmann A., Taquet V., Faure A., Kahane C., Ceccarelli C., 2012, A&A, 541, L12
Balucani et al. (2015) Balucani N., Ceccarelli C., Taquet V., 2015, MNRAS, 449, L16
Becke (1988) Becke A. D., 1988, Phys. Rev. A, 38, 3098
Becke (1993a) Becke A. D., 1993a, J. Chem. Phys., 98, 1372
Becke (1993b) Becke A. D., 1993b, J. Chem. Phys., 98, 5648
Ben Chouikha et al. (2022) Ben Chouikha I., Kerkeni B., Ouerfelli G., Makroni L., Nyman G., 2022, RSC Adv., 12, 18994
Bennett et al. (2005) Bennett C. J., Osamura Y., Lebar M. D., Kaiser R. I., 2005, ApJ, 634, 698
Bergner et al. (2017) Bergner J. B., Öberg K. I., Rajappan M., 2017, ApJ, 845, 29
Bergner et al. (2019) Bergner J. B., Öberg K. I., Rajappan M., 2019, ApJ, 874, 115
Bertin et al. (2016) Bertin M., et al., 2016, ApJ, 817, L12
Bianchi et al. (2019) Bianchi E., et al., 2019, MNRAS, 483, 1850
Biver et al. (2021) Biver N., et al., 2021, A&A, 648, A49
Blake et al. (1987) Blake G. A., Sutton E. C., Masson C. R., Phillips T. G., 1987, ApJ, 315, 621
Boys & Bernardi (1970) Boys S., Bernardi F., 1970, Mol. Phys., 19, 553
Brandenburg et al. (2013) Brandenburg J. G., Alessio M., Civalleri B., Peintinger M. F., Bredow T., Grimme S., 2013, J. Phys. Chem. A, 117, 9282
Callahan et al. (2011) Callahan M. P., Smith K. E., Cleaves H. J., Ruzicka J., Stern J. C., Glavin D. P., House C. H., Dworkin J. P., 2011, Proc. Natl. Acad. Sci., 108, 13995
Casassa et al. (1997) Casassa S., Ugliengo P., Pisani C., 1997, J. Chem. Phys., 106, 8030
Caselli & Ceccarelli (2012) Caselli P., Ceccarelli C., 2012, A&A Rev., 20, 56
Cazaux et al. (2003) Cazaux S., Tielens A. G. G. M., Ceccarelli C., Castets A., Wakelam V., Caux E., Parise B., Teyssier D., 2003, ApJ, 593, L51
Ceccarelli et al. (2014) Ceccarelli C., Caselli P., Bockelée-Morvan D., Mousis O., Pizzarello S., Robert F., Semenov D., 2014, Protostars and Planets VI
Ceccarelli et al. (2017) Ceccarelli C., et al., 2017, ApJ, 850, 176
Ceccarelli et al. (2022) Ceccarelli C., et al., 2022, Protostars and Planets VII, p. arXiv:2206.13270
Chahine et al. (2022) Chahine L., et al., 2022, A&A, 657, A78
Chai & Head-Gordon (2008) Chai J.-D., Head-Gordon M., 2008, J. Chem. Phys., 128, 084106
Charnley (2004) Charnley S., 2004, Advances in Space Research, 33, 23
Charnley et al. (1992) Charnley S. B., Tielens A. G. G. M., Millar T. J., 1992, ApJ, 399, L71
Charnley et al. (1997) Charnley S. B., Tielens A. G. G. M., Rodgers S. D., 1997, ApJ, 482, L203
Chuang et al. (2018) Chuang K.-J., Fedoseev G., Qasim D., Ioppolo S., van Dishoeck E. F., Linnartz H., 2018, ApJ, 853, 102
Chuang et al. (2021) Chuang K.-J., Fedoseev G., Scirè C., Baratta G. A., Jäger C., Henning T., Linnartz H., Palumbo M. E., 2021, A&A, 650, A85
Codella et al. (2018) Codella C., et al., 2018, A&A, 617, A10
Corazzi et al. (2021) Corazzi M. A., Brucato J. R., Poggiali G., Podio L., Fedele D., Codella C., 2021, ApJ, 913, 128
Crovisier et al. (2004) Crovisier J., Bockelée-Morvan D., Colom P., Biver N., Despois D., Lis D. C., the Team for target-of-opportunity radio observations of comets 2004, A&A, 418, 1141
Cutini et al. (2016) Cutini M., Civalleri B., Corno M., Orlando R., Brandenburg J. G., Maschio L., Ugliengo P., 2016, J. Chem. Theory Comput., 12, 3340
Cutini et al. (2017) Cutini M., Corno M., Ugliengo P., 2017, J. Chem. Theory Comput., 13, 370
Cutini et al. (2019) Cutini M., Civalleri B., Ugliengo P., 2019, ACS Omega, 4, 1838
Dapprich et al. (1999) Dapprich S., Komáromi I., Byun K., Morokuma K., Frisch M., 1999, Comput. Theor. Chem., 461-462, 1
Dartois et al. (2019) Dartois E., Chabot M., Id Barkach T., Rothard H., Augé B., Agnihotri A. N., Domaracka A., Boduch P., 2019, A&A, 627, A55
De Simone et al. (2020) De Simone M., et al., 2020, A&A, 640, A75
Dovesi et al. (2018) Dovesi R., et al., 2018, WIREs Comp. Mol. Sci., 8, e1360
Drozdovskaya et al. (2019) Drozdovskaya M. N., van Dishoeck E. F., Rubin M., Jørgensen J. K., Altwegg K., 2019, MNRAS, 490, 50
Dulieu et al. (2013) Dulieu F., Congiu E., Noble J., Baouche S., Chaabouni H., Moudens A., Minissale M., Cazaux S., 2013, Sci. Rep., 3, 1338
Dunning (1989) Dunning T. H., 1989, J. Chem. Phys., 90, 1007
Eley & Rideal (1940) Eley D., Rideal E., 1940, Nature, 146, 401
Enrique-Romero et al. (2016) Enrique-Romero J., Rimola A., Ceccarelli C., Balucani N., 2016, MNRAS, 459, L6
Enrique-Romero et al. (2019) Enrique-Romero J., Rimola A., Ceccarelli C., Ugliengo P., Balucani N., Skouteris D., 2019, ACS Earth Space Chem., 3, 2158
Enrique-Romero et al. (2020) Enrique-Romero J., et al., 2020, MNRAS, 493, 2523
Enrique-Romero et al. (2021) Enrique-Romero J., Ceccarelli C., Rimola A., Skouteris D., Balucani N., Ugliengo P., 2021, A&A, 655, A9
Enrique-Romero et al. (2022) Enrique-Romero J., Rimola A., Ceccarelli C., Ugliengo P., Balucani N., Skouteris D., 2022, ApJS, 259, 39
Fedoseev et al. (2015) Fedoseev G., Cuppen H. M., Ioppolo S., Lamberts T., Linnartz H., 2015, MNRAS, 448, 1288
Fedoseev et al. (2022) Fedoseev G., Qasim D., Chuang K.-J., Ioppolo S., Lamberts T., van Dishoeck E. F., Linnartz H., 2022, ApJ, 924, 110
Ferrero et al. (2020) Ferrero S., Zamirri L., Ceccarelli C., Witzel A., Rimola A., Ugliengo P., 2020, ApJ, 904, 11
Ferrero et al. (2022) Ferrero S., et al., 2022, MNRAS, 516, 2586
Ferrero et al. (2023a) Ferrero S., Pantaleone S., Ceccarelli C., Ugliengo P., Sodupe M., Rimola A., 2023a, The Astrophysical Journal, 944, 142
Ferrero et al. (2023b) Ferrero S., Ceccarelli C., Ugliengo P., Sodupe M., Rimola A., 2023b, ApJ, 951, 150
Fourikis et al. (1974) Fourikis N., Sinclair M. W., Robinson B. J., Godfrey P. D., Brown R. D., 1974, Aust. J. Phys., 27, 425
Frisch et al. (2016) Frisch M. J., et al., 2016, Gaussian 16 Revision C.01
Garrod & Herbst (2006) Garrod R. T., Herbst E., 2006, A&A, 457, 927
Garrod et al. (2008) Garrod R. T., Weaver S. L. W., Herbst E., 2008, ApJ, 682, 283
Gilmore et al. (1976) Gilmore W., Morris M., Johnson D. R., Lovas F. J., Zuckerman B., Turner B. E., Palmer P., 1976, ApJ, 204, 43
Gottlieb (1973) Gottlieb C. A., 1973, in Gordon M. A., Snyder L. E., eds, Molecules in the Galactic Environment. p. 181
Grimme et al. (2010) Grimme S., Antony J., Ehrlich S., Krieg H., 2010, J. Chem. Phys., 132, 154104
Grimme et al. (2011) Grimme S., Ehrlich S., Goerigk L., 2011, J. Comp. Chem., 32, 1456
Gutiérrez-Quintanilla et al. (2021) Gutiérrez-Quintanilla A., et al., 2021, MNRAS, 506, 3734
Harris & Kasemo (1981) Harris J., Kasemo B., 1981, Surf. Sci., 105, L281
Hasegawa & Herbst (1993) Hasegawa T. I., Herbst E., 1993, MNRAS, 263, 589
Herbst (2017) Herbst E., 2017, Int. Rev. Phys. Chem., 36, 287
Herbst & van Dishoeck (2009) Herbst E., van Dishoeck E. F., 2009, Annu. Rev. A&A, 47, 427
Hinshelwood (1930) Hinshelwood C., 1930, Annu. Rep. Prog. Chem., 27, 11
Ioppolo et al. (2011) Ioppolo S., Cuppen H., Linnartz H., 2011, Rend. Fis. Acc. Lincei, 22, 211
Jin & Garrod (2020) Jin M., Garrod R. T., 2020, ApJS, 249, 26
Klahn & Bingel (1977) Klahn B., Bingel W. A., 1977, Int. J. Quantum Chem., 11, 943
Krasnokutski et al. (2020) Krasnokutski S. A., Jäger C., Henning T., 2020, ApJ, 889, 67
Kruse & Grimme (2012) Kruse H., Grimme S., 2012, J. Chem. Phys., 136, 154101
Lamberts et al. (2019) Lamberts T., Markmeyer M. N., Kolb F. J., Kästner J., 2019, ACS Earth Space Chem., 3, 958
Langmuir (1922) Langmuir I., 1922, Trans. Faraday Soc., 17, 621
Law et al. (2021) Law C. J., Zhang Q., Öberg K. I., Galván-Madrid R., Keto E., Liu H. B., Ho P. T. P., 2021, ApJ, 909, 214
Lee et al. (1988) Lee C., Yang W., Parr R. G., 1988, Phys. Rev. B, 37, 785
Lee et al. (2019) Lee C.-F., Codella C., Li Z.-Y., Liu S.-Y., 2019, ApJ, 876, 63
Lefloch et al. (2017) Lefloch B., Ceccarelli C., Codella C., Favre C., Podio L., Vastel C., Viti S., Bachiller R., 2017, MNRAS Letters, 469, L73
Ligterink et al. (2018) Ligterink N. F. W., et al., 2018, A&A, 619, A28
Martín-Doménech et al. (2020) Martín-Doménech R., Öberg K. I., Rajappan M., 2020, ApJ, 894, 98
Matthews et al. (1985) Matthews H. E., Friberg P., Irvine W. M., 1985, ApJ, 290, 609
McClure et al. (2023) McClure M. K., et al., 2023, Nature astronomy, 7, 1
McElroy et al. (2013) McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36
McQuarrie (1975) McQuarrie D. A., 1975, Statistical mechanics. Harper & Row New York
Molpeceres et al. (2022) Molpeceres G., Kästner J., Herrero V. J., Peláez R. J., Maté B., 2022, A&A, 664, A169
Molpeceres et al. (2023) Molpeceres G., Zaverkin V., Furuya K., Aikawa Y., Kästner J., 2023, A&A, 673, A51
Moore & Hudson (1998) Moore M. H., Hudson R. L., 1998, Icarus, 135, 518
Öberg et al. (2009) Öberg K. I., Garrod R. T., Van Dishoeck E. F., Linnartz H., 2009, A&A, 504, 891
Okoshi et al. (2015) Okoshi M., Atsumi T., Nakai H., 2015, J. Comp. Chem., 36, 1075
Paardekooper et al. (2016) Paardekooper D., Bossa J.-B., Linnartz H., 2016, A&A, 592, A67
Pack & Monkhorst (1977) Pack J. D., Monkhorst H. J., 1977, Phys. Rev. B, 16, 1748
Pantaleone et al. (2020) Pantaleone S., Enrique-Romero J., Ceccarelli C., Ugliengo P., Balucani N., Rimola A., 2020, ApJ, 897, 56
Pantaleone et al. (2021) Pantaleone S., Enrique-Romero J., Ceccarelli C., Ferrero S., Balucani N., Rimola A., Ugliengo P., 2021, ApJ, 917, 49
Papakondylis & Mavridis (2019) Papakondylis A., Mavridis A., 2019, J. Phys. Chem. A, 123, 10290
Paulive et al. (2020) Paulive A., Shingledecker C. N., Herbst E., 2020, MNRAS, 500, 3414
Peintinger et al. (2013) Peintinger M. F., Oliveira D. V., Bredow T., 2013, J. Comp. Chem., 34, 451
Perrero et al. (2022a) Perrero J., Enrique-Romero J., Martínez-Bachs B., Ceccarelli C., Balucani N., Ugliengo P., Rimola A., 2022a, ACS Earth and Space Chemistry, 6, 496
Perrero et al. (2022b) Perrero J., Enrique-Romero J., Ferrero S., Ceccarelli C., Podio L., Codella C., Rimola A., Ugliengo P., 2022b, ApJ, 938, 158
Raghavachari et al. (1989) Raghavachari K., Trucks G. W., Pople J. A., Head-Gordon M., 1989, Chem. Phys. Lett., 157, 479
Rimola et al. (2010) Rimola A., Zicovich-Wilson C. M., Dovesi R., Ugliengo P., 2010, J. Chem. Theory Comput., 6, 1341
Rimola et al. (2014) Rimola A., Taquet V., Ugliengo P., Balucani N., Ceccarelli C., 2014, A&A, 572, A70
Rimola et al. (2018) Rimola A., Skouteris D., Balucani N., Ceccarelli C., Enrique-Romero J., Taquet V., Ugliengo P., 2018, ACS Earth Space Chem., 2, 720
Rotelli et al. (2016) Rotelli L., Trigo-Rodríguez J. M., Moyano-Cambero C. E., Carota E., Botta L., di Mauro E., Saladino R., 2016, Sci. Rep., 6, 38888
Ruaud et al. (2015) Ruaud M., Loison J. C., Hickson K. M., Gratier P., Hersant F., Wakelam V., 2015, MNRAS, 447, 4004
Scibelli & Shirley (2020) Scibelli S., Shirley Y., 2020, ApJ, 891, 73
Shingledecker & Herbst (2018) Shingledecker C. N., Herbst E., 2018, Phys. Chem. Chem. Phys., 20, 5359
Shingledecker et al. (2018) Shingledecker C. N., Tennis J., Gal R. L., Herbst E., 2018, ApJ, 861, 20
Simons et al. (2020) Simons M. A. J., Lamberts T., Cuppen H. M., 2020, A&A, 634, A52
Singh et al. (2019) Singh K. K., Tandon P., Misra A., 2019, in Singh D. K., Das S., Materny A., eds, Advances in Spectroscopy: Molecules to Materials. Springer Singapore, Singapore, pp 415–422
Skouteris et al. (2018) Skouteris D., Balucani N., Ceccarelli C., Vazart F., Puzzarini C., Barone V., Codella C., Lefloch B., 2018, ApJ, 854, 135
Sure & Grimme (2013) Sure R., Grimme S., 2013, J. Comp. Chem., 34, 1672
Taquet et al. (2016) Taquet V., Wirström E. S., Charnley S. B., 2016, ApJ, 821, 46
Tatewaki & Huzinaga (1980) Tatewaki H., Huzinaga S., 1980, J. Comp. Chem., 1, 205
VandeVondele & Hutter (2007) VandeVondele J., Hutter J., 2007, J. Chem. Phys., 127, 114105
Vasyunin et al. (2017) Vasyunin A. I., Caselli P., Dulieu F., Jiménez-Serra I., 2017, ApJ, 842, 33
Vazart et al. (2020) Vazart F., Ceccarelli C., Balucani N., Bianchi E., Skouteris D., 2020, MNRAS, 499, 5547
Vilela Oliveira et al. (2019) Vilela Oliveira D., Laun J., Peintinger M. F., Bredow T., 2019, J. Comp. Chem., 40, 2364
Wakelam et al. (2012) Wakelam V., et al., 2012, ApJS, 199, 21
Yabuta et al. (2007) Yabuta H., Williams L. B., Cody G. D., Alexander C. M. O. D., Pizzarello S., 2007, Meteorit. Planet. Sci., 42, 37
Yang et al. (2022) Yang Y.-L., et al., 2022, ApJ, 941, L13
Zamirri et al. (2018) Zamirri L., Casassa S., Rimola A., Segado-Centellas M., Ceccarelli C., Ugliengo P., 2018, MNRAS, 480, 1427
Zhao & Truhlar (2004) Zhao Y., Truhlar D. G., 2004, J. Phys. Chem. A, 108, 6908
Zhao & Truhlar (2008) Zhao Y., Truhlar D. G., 2008, Theor. Chem. Acc., 119
Zhao et al. (2006) Zhao Y., Schultz N. E., Truhlar D. G., 2006, J. Chem. Theory Comput., 2, 364
Öberg et al. (2010) Öberg K. I., van Dishoeck E. F., Linnartz H., Andersson S., 2010, ApJ, 718, 832

Quantum mechanical modeling of the grain-surface formation of acetaldehyde on H₂O:CO dirty ice surfaces

Abstract

keywords:

1 Introduction