# Combined quantum chemical and modeling study of CO hydrogenation on water ice

Albert Rimola, Vianney Taquet, Piero Ugliengo, Nadia Balucani, Cecilia

Ceccarelli

► To cite this version:

Albert Rimola, Vianney Taquet, Piero Ugliengo, Nadia Balucani, Cecilia Ceccarelli. Combined quantum chemical and modeling study of CO hydrogenation on water ice. Astronomy and Astrophysics - A&A, 2014, 572, 10.1051/0004-6361/201424046. insu-03618155

HAL Id: insu-03618155

https://insu.hal.science/insu-03618155

Submitted on 24 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combined quantum chemical and modeling study of CO hydrogenation on water ice

Albert Rimola<fnref n="1" />, Vianney Taquet<fnref n="2" />,*, Piero Ugliengo<fnref n="3" />, Nadia Balucani<fnref n="4" />, and Cecilia Ceccarelli<fnref n="5" />,<fnref n="6" />

Received 22 April 2014 / Accepted 6 October 2014

## ABSTRACT

*Context.* Successive hydrogenation reactions of CO on interstellar icy grain surfaces are considered one of the most efficient mechanisms in interstellar environments for the formation of H₂CO and CH₃OH, two of the simplest organic molecules detected in space. In the past years, several experimental and theoretical works have been focused on these reactions, providing relevant information both at the macroscopic and atomic scale. However, several questions still remain open, such as the exact role played by water in these processes, a crucial aspect because water is the dominant constituent of the ice mantles around dust grain cores.

*Aims.* We here present a quantum chemical description of the successive H additions to CO both in the gas phase and on the surfaces of several water clusters.

*Methods.* The hydrogenation steps were calculated by means of accurate quantum chemical methods and structural cluster models consisting of 3, 18, and 32 water molecules.

*Results.* Our main result is that the interaction of CO and H₂CO with the water cluster surfaces through H-bonds with the O atoms increases the C-O polarization, thus weakening the C-O bond. Consequently, the C atoms are more prone to receiving H atoms, which in turn lowers the energy barriers for the H additions compared to the gas-phase processes. The calculated energy barriers and transition frequencies associated with the reaction coordinate were adopted as input parameters in our numerical model of the surface chemistry (GRAINOBLE) to simulate the distribution of the H₂CO and CH₃OH ice abundances (with respect to water). Our GRAINOBLE results based on the energy barriers and transition frequencies for the reactions on the 32 water molecule cluster compare well with the observed abundances in low-mass protostars and dark cores.

**Key words.** stars: formation – ISM: clouds – ISM: molecules – ISM: abundances

## 1. Introduction

Methanol and formaldehyde are among the few molecules that have been detected in the solid phase, in the icy mantles that cover interstellar grains. Specifically, observational studies showed that water, CO, and CO₂ are present in molecular clouds at visual extinctions A_v ~ 3 mag, whereas solid methanol is only present at A_v ≥ 15 mag (e.g., Whittet et al. 2011). These observations have been interpreted as a proof that methanol is mostly formed in the coldest and densest phase of the cloud evolution, when CO freezes out onto grains that are already enveloped by layers of water ice. The idea is that methanol is formed by the hydrogenation of iced CO by successive addition of H atoms to the grain surfaces (Tielens & Hagen 1982):

$$CO + H \rightarrow HCO \quad (1)$$

$$HCO + H \rightarrow H_2CO \quad (2)$$

$$H_2CO + H \rightarrow CH_3O \quad (3)$$

$$CH_3O + H \rightarrow CH_3OH. \quad (4)$$

In this widely accepted scenario, formaldehyde forms first, while methanol is the final product of the CO hydrogenation (e.g., Taquet et al. 2012a). In fact, the CO hydrogenation is one of the “simplest” chemical reactions occurring on the interstellar grain surfaces, and as such may be considered as a benchmark for grain surface chemistry. For this reason, the study of the CO hydrogenation has received much attention in the past years, and has triggered theoretical quantum chemistry and laboratory works.

Gas-phase experiments devoted to studying reaction (1) (Reilly et al. 1978; Dane et al. 1988; Rumbles et al. 1990; Sappey et al. 1990) indicate, however, that only relatively energetic H-atoms produced by photolysis or radical-induced dissociation of H₂CO may react with CO to produce HCO. The activation energy was also obtained through an Arrhenius plot between 260–340 K using fast Lyman-α spectrophotometry to be about 8.3 ± 2 kJ mol⁻<fnref n="1" /> (Wang et al. 1973), whereas energy barriers calculated using quantum chemical methods span the 15–21 kJ mol⁻<fnref n="1" /> range, depending on the adopted method and level of theory (Woon 1996, 2002; Marenich & Boggs 2003; Goumans et al. 2007, 2008; Peters et al. 2013). For reaction (3),

<footnote n="1">¹ Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra (Cerdanyola del Vallès), Spain
e-mail: Albert.Rimola@uab.cat</footnote>

<footnote n="2">² Astrochemistry Laboratory and The Goddard Center for Astrobiology, Mailstop 691, NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20770, USA</footnote>

<footnote n="3">³ Dipartimento di Chimica, and NIS - Centre for Nanostructured Interfaces and Surfaces, Università di Torino, via P. Giuria 7, 10125 Torino, Italy</footnote>

<footnote n="4">⁴ Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, via Elce di Sotto 8, 06123 Perugia, Italy</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote>* NASA Postdoctoral Program Fellow.</footnote>
quantum chemical calculations give barriers of 21–22 kJ mol⁻¹ (Woon 1996, 2002; Goumans et al. 2007, 2008), whereas reactions (2) and (4) are considered to be barrierless, since they involve radical-radical reactions.

In the past years, several surface experimental works, based on the H atom bombardment of CO (in the form of pure CO ice or as a H₂O–CO ice mixture), have studied the formation of formaldehyde and methanol (Hiraoka et al. 2002; Watanabe & Kouchi 2002; Watanabe et al. 2003, 2007; Hidaka et al. 2004; Fuchs et al. 2009; Pirim et al. 2010; Pirim & Krim 2011). In all cases, CH₃OH was formed at a lower yield than H₂CO. Furthermore, detection of HCO and CH₃O as intermediate species only occurred when H and CO ices were co-deposited rather than by H-atom bombardment of the CO ice. Awad et al. (2005) used experimental data to calculate the rate constants of these processes by means of chemical kinetic formulations that account for the thickness and morphology of the ices. Excellent reviews focused on surface chemical processes on water ices, including the formation of these two molecules, are available in the literature (Watanabe & Kouchi 2008; Hama & Watanabe 2013).

There is, therefore, a great discrepancy between the results obtained by laboratory experiments and the quantum mechanical gas-phase energy barrier estimates for reactions (1) and (3), since the calculated ones are much too high to proceed at the cryogenic temperatures. In the case of cold interstellar clouds, one possible solution is that the energy barriers can be lowered by the matrix of iced water molecules covering the grain surfaces. In this work, we aim to verify this hypothesis by simulating the H addition to CO using quantum chemical calculations based on density functional theory and post-Hartree Fock wavefunction-based methods. The reactions were simulated both in the gas-phase and on the surfaces of water clusters consisting of 3, 18, and 32 water molecules each as structural models for the interstellar water ice. Moreover, the calculated energy barriers have been used in the GRAINOBLE astrochemical model (Taquet et al. 2012b) to estimate the abundances of H₂CO and CH₃OH ices with respect to water, which in turn have been compared with the observations.

## 2. Methods

### 2.1. Quantum chemical calculations

All calculations were performed using the package program GAUSSIAN03 (Frisch et al. 2004). Methods based on the density functional theory (DFT) have been shown to be cost-effective and have been used to study a wide variety of closed-shell systems with great accuracy (Koch & Holthausen 2001; Sousa et al. 2007; Sholl & Steckel 2009). However, for open-shell systems, calculations carried out by some of us have demonstrated that functionals with a higher percentage of exact exchange, such as BHLYP, provide results in better agreement with those from the highly correlated CCSD(T) method (Poater et al. 2004; Rimola et al. 2006, 2008). This is because functionals based on the generalized gradient approach (GGA) or hybrid approaches with low percentages of exact exchange such as B3LYP (20%) tend to overstabilize electron-delocalized situations as a result of the self-interaction error (Sodupe et al. 1999). Accordingly, the structure of the reactants, products, intermediates, and transition structure analyzed in this work were fully optimized using the hybrid density functional method BHLYP (Lee et al. 1988; Becke 1993) with the standard 6-311++G(d,p) Pople basis set (Krishnan et al. 1980; hereafter referred to as

the L1 theory level) because this functional, which includes the 50% of the exact exchange, has been demonstrated to provide reasonably accurate values for the energy barriers of H addition reactions of astrochemical interest (Andersson & Grüning 2004). Despite this, we also computed single-point energy calculations at the highly-correlated coupled cluster CCSD(T) level (Pople et al. 1987) with the aug-cc-pvtz Dunning basis set (Dunning 1989; hereafter referred to as the L2 theory level) on the BHLYP-optimized geometries for the structures involved in the gas-phase processes and in the presence of the three-water-molecule cluster to evaluate the reliability of L1.

Open-shell calculations were based on an unrestricted formalism. All structures were characterized as minima (reactants, products, and intermediates) and saddle points (transition states) by calculating the analytical harmonic frequencies. In some cases, we also carried out intrinsic reaction coordinate calculations at the same theory level to ensure the nature of the minima connected by a given transition state. Thermochemical corrections to the potential energy values to obtain zero-point energy-corrected values were carried out using the standard rigid rotor/harmonic oscillator formulas (McQuarrie 1986).

### 2.2. GRAINOBLE model

To follow the macroscopic formation of interstellar ices, we used the GRAINOBLE astrochemical model presented in previous studies (Taquet et al. 2012a, 2013, 2014). Briefly, GRAINOBLE couples the gas-phase and grain-surface chemistry with the rate equations approach introduced by Hasegawa et al. (1992) to follow the evolution of chemical abundances with time for constant physical conditions and for assumed initial abundances in the gas phase until the time reaches one million years. The considered gas-grain surface processes are i) accretion of gas-phase species on the surface of spherical grains with a constant size; ii) diffusion of adsorbed species via thermal hopping; iii) reaction between two particles when they meet each other in the same site by the Langmuir-Hinshelwood mechanism, in which the reaction rate is given by the product of the number of times that the two reactants meet each other and the transmission probability P<sub>r</sub> of the reaction; iv) desorption of adsorbed species into the gas phase by thermal desorption, which exponentially depends on the binding energy of each species E<sub>b</sub> relative to the substrate (see Taquet et al. 2014, for a list of binding energies used in the model), and cosmic-ray-induced heating of grains following Hasegawa & Herbst (1993a) and adapted for the binding energies considered in this work.

We used the multilayer approach developed by Hasegawa & Herbst (1993b) to follow the multilayer formation of interstellar ices, which considers three sets of differential equations: one for gas-phase species, one for surface species, and one for bulk species. The equations governing chemical abundances on the surface and in the bulk are linked by an additional term that is proportional to the rate of growth or evaporation of the grain mantle. Surface species are continuously trapped into the bulk because of the accretion of new particles.

Following Taquet et al. (2012a), we only considered the accretion of gaseous H, O, and CO onto interstellar grains to model the formation of the main ice components H₂O, CO, CO₂, H₂CO, and CH₃OH. The density of atomic H was fixed to 1.2 cm⁻³, while the abundances of O and CO evolve with time. The initial abundance of CO was set to 5 × 10⁻⁵ cm⁻³, while the initial abundance of atomic oxygen is a free parameter (see next paragraph). The surface chemical network was based on the one
presented by Taquet et al. (2013) but does not include the formation of ammonia, methane, and deuterated species.

Methanol ice is believed to mainly form at the center of prestellar cores at low temperatures ($T < 12$ K), relatively high densities ($n_H \ge 10^5$ cm$^{-3}$), and high visual extinctions ($A_V \ge 15$ mag; see Whittet et al. 2011; Taquet et al. 2012b). As discussed extensively in Taquet et al. (2012b, 2013), interstellar surface chemistry depends on several physical, chemical, and surface parameters that either evolve with time, are poorly constrained, or show distributions of values depending on the ice morphology. We considered a large grid of models by linearly varying the following input parameters within their ranges of values given by the literature to explore the influence of these parameters on the chemical composition of ices:

- Five values for the dust and gas temperatures $T_g = T_d$ between 8 and 12 K, representative of typical dense cores where interstellar methanol is believed to mainly form.

- Five values for the diameter of interstellar grains $a_d$ between 0.1 and 0.3 µm, corresponding to the uncertainty in the upper limit of the grain size distributions derived in the ISM.

- Five values for the diffusion-to-binding energy ratio $E_d/E_b$ between 0.4 and 0.8, a range given by theoretical and experimental studies.

- Five values for the binding energy of atomic hydrogen between 400 and 600 K given by theoretical and experimental studies.

- Five values for the initial abundance of atomic O between $8 \times 10^{-5}$ and $1.2 \times 10^{-4}$. Most of the atomic O freezing-out onto grains is converted into water, giving abundances of $\sim 10^{-4}$ relative to the total number of H nuclei with some uncertainty (Tielens et al. 1991). The range of values for the initial abundance of atomic O is chosen to take the uncertainty of the observed abundance of water ice into account.

### 3. Results of quantum chemical calculations

This section is organized as follows: first, results focused on the H additions to CO in the gas phase (i.e., in absence of water) are presented. Then, the very same reactions in the presence of a three-water-molecule cluster are shown. Finally, results for the hydrogenation processes occurring on two ice water cluster models consisting of 18 and 32 water molecules are reported.

#### 3.1. H additions to CO in the gas phase

Figure 1 shows the structures involved in the CO hydrogenation in the gas phase, and Table 1 shows the energetics (i.e., energy barriers, $\Delta E^\ddagger$, and reaction energies, $\Delta E_r$) of the calculated processes at the L1 and L2 levels. Moreover, energetic data including the zero-point energy (ZPE) corrections calculated at the L1 level to both L1 and L2 potential energies ($\Delta E^\ddagger + ZPE$) are also included. Hereafter, for the sake of clarity, an energetic value calculated at a given theory level is denoted as $\Delta E/L$. For instance, $\Delta E^\ddagger + ZPE/L2$ means an energy barrier including ZPE at the L2 theory level. Following the International System of Units (ISU), all energy units are given in kJ mol$^{-1}$, whose conversion factor to K is $1 \text{ kJ mol}^{-1} = 120.274 \text{ K}$.

Figures 1a and 1b show the first hydrogenation step (i.e., H + CO) to form either HCO (by an H-C coupling) or COH (by an H-O coupling). The energetic data indicate that only HCO is thermodynamically stable and kinetically feasible compared to COH at either theory level, in agreement with the exclusive detection of the HCO isomer in experiments (Pirim et al. 2010;

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** BHLYP/6-311++G(d,p)-optimized geometries of the stationary points involved in the hydrogenation processes in the gas phase: panel a) hydrogenation of CO on the C atom, panel b) hydrogenation of CO on the O atom, panel c) hydrogenation of $H_2CO$ on the C atom, and panel d) hydrogenation of $H_2CO$ on the O atom. AS refers to the CO + H zero-energy references for panels a) and b) and to the $H_2CO + H$ zero-energy references for panels c) and d). TS refers to the transition-state structures and PROD to the products. Distances in Å.</figcaption></figure>

Pirim & Krim 2011). These dramatic differences are due to the intrinsic different stability of the C-radical versus O-radical centered species present in HCO and COH isomers, respectively. Our computed data for HCO and COH formation ($\Delta E^\ddagger/L2$ and $\Delta E_r/L2$) are very similar to those obtained by other theoretical works in which the geometries were optimized at the CCSD(T) level using augmented Dunning basis sets (this work vs. literature (van Mourik et al. 2000; Marenich & Boggs 2003): 169.2 vs. 176 kJ mol$^{-1}$ for the relative energy between HCO and COH; 13.0 vs. 15 kJ mol$^{-1}$ for $\Delta E^\ddagger$ of the HCO formation; -80.4 vs. -81 kJ mol$^{-1}$ for $\Delta E_r$ of the HCO formation; 133.1 vs. 137 kJ mol$^{-1}$ for $\Delta E^\ddagger$ of the COH formation; and 88.8 vs. 94 kJ mol$^{-1}$ for $\Delta E_r$ of the COH formation.

As mentioned in the Introduction, the addition of a second H atom (i.e., H + HCO $\rightarrow H_2CO$) is a barrierless radical-radical reaction and, accordingly, it was not considered here. The addition of a third H atom (i.e., H + $H_2CO$) can give either $CH_3O$ or $CH_2OH$, the calculated structures involved in the processes are presented in Figs. 1c and d. $CH_2OH$ is a C-radical centered species with a nonplanar geometry (the C atom adopts a pyramidal structure), whereas the unpaired electron of $CH_3O$ is fully localized on the O atom. Because of that, $CH_2OH$ is found to be more stable than $CH_3O$ by about 35 kJ mol$^{-1}$, although both processes are largely exoenergetic ($\Delta E^\ddagger + ZPE/L2 = -120.1 \text{ kJ mol}^{-1}$ and $-85.7 \text{ kJ mol}^{-1}$).
<table><caption>**Table 1.** Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO, COH, CH₃O, and CH₂OH in gas phase.</caption><thead><tr><th rowspan="2">Reaction</th><th colspan="3">L1</th><th colspan="3">L2</th></tr><tr><th>AS</th><th>TS</th><th>PROD</th><th>AS</th><th>TS</th><th>PROD</th></tr></thead><tbody><tr><td>a) CO + H → HCO</td><td>0.0<br>[0.0]</td><td>8.1<br>[9.9]</td><td>-97.3<br>[-75.6]</td><td>0.0<br>[0.0]</td><td>13.0<br>[14.7]</td><td>-80.4<br>[-58.6]</td></tr><tr><td>b) CO + H → COH</td><td>0.0<br>[0.0]</td><td>121.2<br>[124.7]</td><td>66.2<br>[89.8]</td><td>0.0<br>[0.0]</td><td>133.1<br>[136.6]</td><td>+88.8<br>[+112.4]</td></tr><tr><td>c) H₂CO + H → CH₃O</td><td>0.0<br>[0.0]</td><td>10.2<br>[14.3]</td><td>-147.0<br>[-119.9]</td><td>0.0<br>[0.0]</td><td>14.5<br>[18.6]</td><td>-112.8<br>[-85.7]</td></tr><tr><td>d) H₂CO + H → CH₂OH</td><td>0.0<br>[0.0]</td><td>32.9<br>[36.4]</td><td>-164.5<br>[-136.2]</td><td>0.0<br>[0.0]</td><td>41.1<br>[44.6]</td><td>-148.4<br>[-120.1]</td></tr></tbody></table>

**Notes.** Calculations at BHLYP/6-311++G(d,p) (L1) and at CCSD(T)/aug-cc-pVTZ//BHLYP/6-311++G(d,p) (L2) levels. ZPE corrections were calculated at the L1 level and included to the L1 and L2 potential energies. Values in units of kJ mol⁻¹. AS refers to the zero-energy references, TS to the transition-state structures, and PROD to the products. Geometries and labeling of the considered cases are the same as in Fig. 1.

However, the calculated energy barriers go against this trend, as ΔE⁺⁺+ZPE/L2 is higher for CH₂OH formation than that for CH₃O (44.6 kJ mol⁻¹ and 18.6 kJ mol⁻¹), so that the less stable isomer is expected to form faster than the most stable one. Previous works (Saebo et al. 1983; Geers et al. 1993; Walch 1993; Dertinger et al. 1995; Taketsugu et al. 1996; Marcy et al. 2001; Petraco et al. 2002; Kamarchik et al. 2012) examined the isomerization of CH₃O into CH₂OH as a way to reach the most stable isomer, but both experiment and theory concluded that this isomerization is unlikely because the barrier for isomerization is significantly higher than the CH₃O dissociation barrier. Accordingly, and considering the cryogenic temperatures of the interstellar medium at which the processes occur, formation of CH₂OH is kinetically hampered in favor of CH₃O formation. This conclusion agrees with the experimental exclusive detection of CH₃O (Pirim et al. 2010). Remarkably, calculated ΔE⁺⁺/L2 values for CH₃O and CH₂OH formation (14.5 kJ mol⁻¹ and 41.1 kJ mol⁻¹) agree reasonably well with those obtained by Kamarchik et al. (2012) at the CCSD(T)/aug-cc-pvtz (22.6 kJ mol⁻¹ and 47.1 kJ mol⁻¹). The ZPE corrections do not affect the general trends described above (see values in brackets of Table 1), the main effects are an increase of the ΔE⁺⁺ and ΔEₜ values by about 4 kJ mol⁻¹ and 30 kJ mol⁻¹ at most.

## 3.2. H additions to CO on a three-water-molecule cluster

Figure 2 shows the calculated structures for the CO hydrogenation steps occurring on an ice cluster modeled by three water molecules in the shape of the most stable water trimer (hereafter referred to as W3). Table 2 reports the energetics for these reactions. The very same reactions have been calculated for one water molecule (W1), showing very similar results to those with W3 and, thus, for the sake of brevity, they are reported in Appendix A.

The interaction of CO with water can take place through either ends of the CO molecule. For the 1:1 water-carbon monoxide system, the two possible structures were experimentally identified (Lundell & Rasanen 1995), with the conclusion that the two complexes coexist. Results based on highly-correlated quantum chemical calculations agree with experiments (Reed et al. 1986; Parish et al. 1992; Lundell 1995; Sadlej et al. 1995; Lundell et al. 1997). We here predict that the C-attached complex is the most stable one from the computed potential interaction energies for the W3-carbon monoxide complex at L1. However, including ZPE corrections makes both complexes nearly isoenergetic, the interaction through the C-attached

<figure><img src="image_3.png" /></figure>

**Fig. 2.** BHLYP/6-311++G(d,p)-optimized geometries of the stationary points involved in the hydrogenation processes for a three-water-molecule cluster (W3): panel **a**) hydrogenation of CO on the C atom when CO interacts with W3 through the O atom, panel **b**) hydrogenation of CO on the C atom when CO interacts with W3 through the C atom, panel **c**) hydrogenation of H₂CO on the C atom, panel **d**) hydrogenation of CO on the O atom when CO interacts with W3 through the C atom, and panel **e**) hydrogenation of H₂CO on the O atom. AS refers to the CO···W3 + H zero-energy reference for panel **a**), to the OC···W3 + H zero-energy references for panels **b**) and **d**), and to the H₂CO···W3 + H zero-energy references for panels **c**) and **e**). TS refers to the transition-state structures and PROD to the products. Distances in Å.
<table><caption>**Table 2.** Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO and CH<sub>3</sub>O for a three-water-molecule cluster (W3).</caption><thead><tr><th rowspan="2">Reaction</th><th colspan="3">L1</th><th colspan="3">L2</th></tr><tr><th>AS</th><th>TS</th><th>PROD</th><th>AS</th><th>TS</th><th>PROD</th></tr></thead><tbody><tr><td>a) CO···W3 + H → HCO···W3</td><td>0.0<br>[0.0]</td><td>7.8<br>[9.7]</td><td>-111.4<br>[-85.4]</td><td>0.0<br>[0.0]</td><td>12.6<br>[14.4]</td><td>-94.3<br>-68.3</td></tr><tr><td>b) OC···W3 + H → HCO···W3</td><td>0.0<br>[0.0]</td><td>9.4<br>[11.1]</td><td>-109.4<br>[-84.7]</td><td>0.0<br>[0.0]</td><td>14.1<br>[15.7]</td><td>-90.5<br>[-65.8]</td></tr><tr><td>c) H<sub>2</sub>CO···W3 + H → CH<sub>3</sub>O···W3</td><td>0.0<br>[0.0]</td><td>9.8<br>[13.9]</td><td>-146.4<br>[-120.7]</td><td>0.0<br>[0.0]</td><td>10.7<br>[14.7]</td><td>-113.5<br>[-87.8]</td></tr><tr><td>d) OC···W3 + H → HOC···W3</td><td>0.0<br>[0.0]</td><td>118.0<br>[122.1]</td><td>+20.9<br>[+51.3]</td><td>0.0</td><td>129.8</td><td>+42.0</td></tr><tr><td>e) H<sub>2</sub>CO···W3 + H → CH<sub>2</sub>OH···W3</td><td>0.0<br>[0.0]</td><td>32.3<br>[36.4]</td><td>-163.8<br>[-136.0]</td><td>0.0</td><td>39.1</td><td>-147.3</td></tr></tbody></table>

**Notes.** Calculations at BHLYP/6-311++G(d,p) (L1) and at CCSD(T)/aug-cc-pVTZ//BHLYP/6-311++G(d,p) (L2) levels. ZPE corrections were calculated at the L1 level and included to the L1 and L2 potential energies. Values in units of kJ mol<sup>-1</sup>. AS refers to the zero-energy references, TS to the transition-state structures, and PROD to the products. Geometries and labeling of the considered cases are the same as in Fig. 2.

complex being marginally more favorable (data available in Appendix B). Accordingly, for the first H addition reaction, both complexes (O-attached and C-attached) are considered in the following.

Figure 2a shows the H addition to the C atom, driving the reaction through the formation of the HCO isomer via a complex with the model ice interacting through the O atom (CO···W3). The computed ΔE<sup>‡</sup>/L2 value (see Table 2 line a) is only marginally lower than that computed for the same gas-phase reaction (12.6 kJ mol<sup>-1</sup> vs. 13.0 kJ mol<sup>-1</sup>). The same tiny energy barrier lowering was also computed for the H addition to H<sub>2</sub>CO to give CH<sub>3</sub>O (ΔE<sup>‡</sup>/L2 = 14.1 kJ mol<sup>-1</sup> and 14.5 kJ mol<sup>-1</sup> for W3 and in the gas phase, see Table 2 line c). The tiny ΔE<sup>‡</sup> decreases are due to the polarization of the C–O bonds induced by the H-bond interaction with the cluster model. This is shown by the larger C–O distances compared with those in the gas phase (1.117 Å and 1.115 Å for CO; 1.194 Å and 1.189 Å for H<sub>2</sub>CO), and a consequent C–O bond weakening, which favors the H additions. ΔE<sub>r</sub> values for both HCO and H<sub>3</sub>CO formation with W3 were computed to be largely exoenergetic, similarly to the gas-phase process.

Figure 2b shows the same HCO formation reaction, but starting from the C-attached complex (OC···W3 structure). At variance with process 2a), the calculated ΔE<sup>‡</sup>/L2 is slightly higher than that obtained in the gas phase (14.0 kJ mol<sup>-1</sup> vs. 13.0 kJ mol<sup>-1</sup>). Here, the H-bond with the cluster model induces a shortening of the C–O distance compared to its gas-phase value (1.113 Å vs. 1.115 Å), strengthening the C–O bond, which increases the barrier toward the H addition compared to the gas-phase reaction. Similar results were already found by Woon et al. (2002) at the MP2 and QCISD levels. The C–O bond strengthening (weakening) for the O-attached (C-attached) complexes have indeed been observed by FTIR matrix isolation spectroscopy (Lundell & Rasanen 1995), as the CO fundamental frequency is bathochromic (hypochromic) shifted from the value of the gas CO. These effects have been described in the literature (Reed et al. 1986; Lundell 1995; Lundell & Latajka 1997) and are related to the bonding/antibonding character of the CO electron lone pair engaged in the H-bond interaction.

Finally, H addition to the O atom to yield COH from a C-attached complex results in a very unfavorable process, in line with gas-phase results (see Fig. 2d and Table 2 line d). Moreover, computed results for the formation of the CH<sub>2</sub>OH isomer (see

Fig. 2e and Table 2 line e) are very similar to those computed in the gas phase (i.e., CH<sub>2</sub>OH is thermodynamically favored, but kinetically hindered compared to CH<sub>3</sub>O, see above). All attempts to locate activated complexes converting CH<sub>3</sub>O to CH<sub>2</sub>OH through a H relay W3 cluster-assisted mechanism, similarly to what water ice does in proton transfer processes, have failed. The negative result, at variance with what is computed for a similar water ice-assisted proton transfer (Koch et al. 2008; Rimola et al. 2010; Peters et al. 2011) is probably due to the transfer of a neutral H rather than a proton, the latter being much more sensitive to H-bond interaction with water ice than the former.

### 3.3. H additions to CO on large water ice clusters

The results shown in the previous sections are useful to understand the intrinsic nature of the CO hydrogenation reactions and the influence that the water-CO interaction can exert on them. However, we here aim to simulate these reactions in the ISM and, accordingly, large and more realistic water ice models should be used. Two different large cluster models with 18 and 32 water molecules, hereafter referred to as W18 and W32 (see Figs. 3 and 4), were considered. Because their size, all energetic data of this section were calculated at the L1 level.

For the first hydrogenation (potential energy surfaces and structures shown in Figs. 3a and 4a for the W18 and W32 ice clusters), the pre-reactant structures involve O-attached complexes formed by one H-bond interaction (see CO···W18 and CO···W32). The C-attached complexes were also computed and found to be nearly degenerate to the O-attached complexes (data available in Fig. B.1 and Table B.1). Nevertheless, since previous results for the W3 cluster model indicate that the attack through the C atom is less effective than through the O atom, only the latter was considered in the larger ice cluster models. Similarly to for W3, the CO activation due to H-bond with the water ice clusters consistently lowers the calculated ΔE<sup>‡</sup>/L1 values (6.8 kJ mol<sup>-1</sup> and 6.6 kJ mol<sup>-1</sup> for W18 and W32), which are the lowest values calculated in this work for the first hydrogenation reaction. The corresponding ΔE<sub>r</sub>/L1 indicates that both H additions are largely exoenergetic and almost independent of the size of the water ice cluster.

The CH<sub>3</sub>O formation is influenced by the different H-bond interactions between H<sub>2</sub>CO and the W18 and W32 ice model
<figure><img src="image_4.png" /></figure>

**Fig. 3.** BHLYP/6-311++G(d,p)-energy profiles (in kJ mol⁻¹) for the H addition to CO for a 18-water ice cluster (W18). Formation of HCO **a**) and CH₃O **b**) species. Bare values are relative potential energy values; values in brackets with the zero-point energy corrections. The zero-energy references are **a**) CO···W18 + H and **b**) H₂CO···W18 + H. Bond lengths in Å.

clusters, as shown in Figs. 3b and 4b. For W18, the calculated ΔE⁻ is 8.9 kJ mol⁻¹, close to that computed for W3 because of the structural similarities between the two cases. For the largest W32 model, ΔE⁻ decreases to 6.5 kJ mol⁻¹, the lowest value calculated in this work. This decrease is due to the two moderately strong H-bonds holding H₂CO in place in the W32 ice model (Fig. 4b), which increases the C–O polarization and thereby reduces the energy barrier value.

As mentioned above, hydrogenation of HCO and CH₃O to form the final H₂CO and CH₃OH species, was not considered because they are barrierless radical-radical reactions. However, for large water ice clusters, the incoming H atoms can stick on the water ice surface to diffuse toward the adsorbed intermediates, instead of reacting directly with the radical intermediates. When this is the case, the energy barrier for the formation of H₂CO and CH₃OH is no longer associated with the radical-radical coupling, but with the diffusion process. To characterize this process, the adsorption of an H atom on the surface of the HCO/W18 model cluster was computed. When geometry optimization is carried out in a triplet electronic state (i.e., the unpaired electrons of H and HCO with spins of the same sign), the resulting structure presents H and HCO co-adsorbed on W18 distanced by 5.62 Å (see Fig. 5, triplet state). In contrast, when this structure is optimized in a singlet electronic state (i.e., the H and HCO spins couple to define a closed-shell state),

<figure><img src="image_5.png" /></figure>

**Fig. 4.** BHLYP/6-311++G(d,p)-energy profiles (in kJ mol⁻¹) for the H addition to CO for a 32-water ice cluster (W32). Formation of HCO **a**) and CH₃O **b**) species. Bare values are relative potential energy values; values in brackets with the zero-point energy corrections. The zero-energy references are **a**) CO···W32 + H and **b**) H₂CO···W32 + H. Bond lengths in Å.

<figure><img src="image_6.png" /></figure>

**Fig. 5.** BHLYP/6-311++G(d,p)-optimized geometries of an H atom and one HCO molecule co-adsorbed on W18 in the triplet electronic state (**left**) and upon geometry relaxation in the singlet electronic state (**right**). Bond lengths in Å.

the system collapses into the H₂CO···W18 (see Fig. 5, singlet state) because of the spontaneous formation of H₂CO, indicating that the H coupling to HCO to give H₂CO on the W18 ice surface is a barrierless channel. Despite this result, the structures of Fig. 5 cover the situation in which H and HCO are in close proximity; larger ice clusters or periodic water surfaces will allow modeling situations in which H and HCO will be farther away and probably H jumps with sizeable barriers will be mandatory to encounter one another.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** BHLYP/6-311++G(d,p)-energy profiles (in kJ mol⁻¹) for the H addition to CO in the presence of a 4-CO molecule cluster. Formation of HCO **a**) and CH₃O **b**) species. Bare values are relative potential energy values; values in brackets with the zero-point energy corrections. The zero-energy references are **a**) H + COice and **b**) H + H₂CO · · · COice. Bond lengths in Å.</figcaption></figure>

### 3.4. H additions to pure CO ice

Since most of the experiments on the CO hydrogenation are based on the H additions to pure CO ice films, we have also calculated the very same processes using a pure CO ice cluster consisting of four CO molecules. This model is derived from the unit cell of the crystalline CO (Cromer et al. 1983). Since the interaction between the CO molecules is mainly dictated by dispersive forces, which are missing in the definition of the BHLYP method, the geometry optimizations of the stationary points involved in these processes were made keeping the positions of the O atoms fixed. Results for the formation of HCO and CH₃O at L1 are shown in Fig. 6. The calculated energetic values (i.e., energy barriers and reaction energies) are very similar to the gas-phase values (see Table 1, line a). This is consistent with the highly inert chemical behavior of CO, in which activation of CO is not expected.

## 4. Discussion

### 4.1. Effect of water ice on the energy barriers

Gas-phase results indicate that the most favorable channel for the reaction of H + CO is toward the formation of the HCO species through a very exo-energetic reaction with low kinetic barriers ($ΔE^‡+ZPE/L1 = 9.9$ kJ mol⁻¹). The corresponding reaction to give the COH isomer exhibits a very high kinetic barrier

($ΔE^‡+ZPE/L1 = 121.2$ kJ mol⁻¹) and a highly endo-energetic character. The formation of H₂CO from HCO is obviously barrierless (radical-radical reaction). In contrast, the reaction of H + H₂CO brings about either CH₂OH or CH₃O, both of which are exo-energetic processes with a definite preference for the CH₂OH molecule. Nevertheless, the calculated kinetic barrier $ΔE^‡$ associated with the formation of CH₃O is lower than that of CH₂OH ($ΔE^‡+ZPE/L1 = 10.2$ kJ mol⁻¹ and 32.9 kJ mol⁻¹) and, therefore, formation of CH₂OH is kinetically hampered at cryogenic temperatures.

For W3, W18, and W32, CO can be attached by H-bonding either through the C or O ends. Both cases were found to be nearly energetically degenerate. Calculations limited to W3 indicate that the H-bond given in the O-attached complex weakens the C–O bond, thus slightly reducing the kinetic barrier for the formation of HCO ($ΔE^‡+ZPE/L1 = 9.7$ kJ mol⁻¹), whereas for the C-attached complex the C–O bond is strengthened with a corresponding increase of the kinetic barrier for the H addition ($ΔE^‡+ZPE/L1 = 11.1$ kJ mol⁻¹). The polarization of the C–O bond and its subtle activation toward H attack was studied in more detail for the hydrogenation reactions occurring on the larger W18 and W32 water ice cluster models. In these large models, the addition of the first H to CO gives $ΔE^‡+ZPE/L1$ values lower than in the smaller W3 model (9.0 kJ mol⁻¹ and 9.2 kJ mol⁻¹) because of the stronger H-bond interactions and the increased C–O polarization. Similarly, the relatively strong H-bond interactions between the surface water molecules of the ice cluster models with the O atom of H₂CO is also the reason why the formation of CH₃O exhibits lower energy barriers than the gas-phase process ($ΔE^‡+ZPE/L1 = 13.9$ kJ mol⁻¹, 12.8 kJ mol⁻¹, and 10.8 kJ mol⁻¹ on W3, W18 and W32, respectively).

Despite the effect of water in reducing the energy barriers for the HCO and CH₃O formation with respect to the gas-phase processes, the calculated values are exceedingly high considering the very low temperatures of the dense clouds. One possible source of error for the accuracy of the calculations is the adopted theory level, both the functional and the basis set. However, the barrier value for small systems computed with the adopted BHLYP functional is even underestimated with respect to CCSD(T), which in turn are still underestimated with respect to the real values as reported by Peters et al. (2013), so that we are already biasing the barrier toward lower values. Furthermore, we did not expect any dramatic change due to basis set improvement over the adopted L2 level. According to the results obtained using either the L1 or L2 levels for the gas-phase reactions and in the presence of W3, it is expected that calculations based on the larger model clusters using L2 will probably increase the barrier by about 4 kJ mol⁻¹ with respect to the values computed at the L1 level. Another possible limitation that may affect the energy barrier values is the size of the large water clusters, since we modeled water ice particles of micrometer sizes with structural molecular models of about 10 Å (W18) and 15 Å (W32). However, the barrier decrease experienced when the reaction occurs on these two water models is practically the same (the latter being 0.2 kJ mol⁻¹ lower than the former), which presumably means that a comparable barrier lowering effects are expected for more extended water surfaces. A third possible effect on the barrier height may also come from the fact that we modeled the H/CO reactions at the surface of an icy grain, whereas the same reaction can also occur in the core of the ice. In that case, however, we expect an increase in the actual barrier due to the cost of displacing the surrounding water molecules that interact with the target CO. Furthermore, diffusion has to be taken into account
for the H atom to reach a CO molecule hidden in the ice core, and accordingly, the adopted model of surface ice reaction is probably biased toward a lowering of the kinetic barrier value compared to bulky ice models. Taking these considerations into account, we believe, in essence, that the CO hydrogenations, ir- respectively of where they occur (in the gas phase and on water ice surfaces) present high energy barriers to be efficient enough in the cold deep-space. Even though water ice induces a decrease of the entrance channel barrier, such a lowering is only a mea- gre 10%. However, tunneling effects can play a significant role considering the very low temperatures and the lightness of the H atoms involved in the processes. The importance of tunnel- ing effects primarily depends on the curvature of the barrier, which is controlled by the transition frequency (ν$^{‡}$) and, to a lower degree, on the height of barrier (ΔU$^{‡}$). A characteristic parameter is the tunneling cross-over temperature (T$_{x}$), below which tunneling becomes dominant and above which tunneling becomes negligible, which can be calculated using the formula of Fermann & Auerbach (2000)

$$
T_x = \frac{\hbar\omega^\dagger\Delta U_0^\dagger/k_B}{2\pi\Delta U_0^\dagger - \hbar\omega^\dagger\ln 2}, \quad (5)
$$

where ω‡ = 2πν‡, ℏ = h/2π and h is the Planck constant, ΔU‡ 0 is the energy barrier including zero-point energy corrections, and k$_{B}$ is the Boltzman constant. From the computed ν‡ and ΔU‡ 0 for reactions on W32, we calculated the T$_{x}$ values. For the formation of HCO, T$_{x}$ is 142 K, whereas for CH$_{3}$O it is 175 K. Therefore, at the very low temperatures of the ISM (as well as at those of the experiments) tunneling contributions are expected to be rel- evant. Because of that, future quantum mechanical-based works should certainly focus on using more rigorous methods for reac- tion dynamics that incorporate tunneling such as the Feynman’s instanton theory implemented as harmonic quantum transition state theory (HQTST).

4.2. Astrochemical modeling

The present quantum chemical results provides a detailed atom-
istic description of the reaction channels for H additions to CO
on water ices. However, these results are limited because of their
atomic-scale nature. Thus, if macroscopic information is desired
to be compared with molecular spectroscopic observations such
as the molecular abundances, one must resort to numerical astro-
chemical models. The predictions and information derived from
these models, however, are strongly dependent on the parameters
used as input data, among which are energetic parameters such
as energy barriers. In most cases, the introduced energy barri-
ers arise from gas-phase experimental results, a rather dubious
procedure when simulating gas-grain processes, as the gas-grain
interaction can deeply alter the value of this parameter, and con-
sequently the predictions. In other cases, the energy barriers are
derived from numerical or experimental fittings giving rise to a
widespread range of energy barriers.

We here proceeded in a different way. To estimate the rela-
tive abundances of CO, H₂CO, and CH₃OH, we ran three grids
of 3125 models each using the GRAINOBLE model by com-
bining the input parameters described in Sect. 2.2; that is, gas
temperatures, grain diameters, diffusion-to-binding energy ra-
tios, binding energies of atomic H, and the initial abundance
of atomic O, and using three sets of energetic data, which in
essence are used to calculate rate constants and quantum trans-
mission probabilities: i) those calculated with the W32 cluster
model of this work; ii) those calculated by Peters (2013) and

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Distributions of the ice abundances of H₂CO (upper panel) and CH₃OH (lower panel) with respect to water predicted by GRAINOBLE. Standard blue lines show the distributions using the water cluster calculations of this work, thin green lines correspond to the gas-phase calculations by Peters et al. (2013), and red thick lines correspond to the ice measurements by Fuchs et al. (2009). Gray boxes show the observed abundances summarized by Oberg et al. (2011) toward low-mass protostars and dark cores.</figcaption></figure>

Peters et al. (2013) for gas phase conditions; and iii) those mea-
sured on laboratory ices by Fuchs et al. (2009) and by intro-
ducing the same formalism as in Monte Carlo chemical mod-
els to compute the probability of surface reactions; that is, the
probability of reaction is not directly given by the probability of
transmission but is a competition between the reaction, diffusion,
and evaporation. The introduced values are reported in Table 3
and concern the energy barriers of the forward and reverse reac-
tions (E<sub>a,for</sub> and E<sub>a,rev</sub>) and the transition frequency of the tran-
sition state (ν<sup>‡</sup>). The energy barriers of Fuchs et al. (2009), as
explained in their paper, are not real energies, but effective en-
ergies, which reflect the energy the reactants need to thermally
react. Accordingly, the reaction rates were calculated using the
Arrhenius equation employing a pre-exponential factor estimate
of 10<sup>12</sup> s<sup>-1</sup>.

The abundance distributions (with respect to water) of H₂CO and CH₃OH induced by the variation of the input parameters are shown in Fig. 7. For the three sets of activation barriers, the chemical abundances of formaldehyde and methanol show a wide range of values since their formation rates exponentially depend on the poorly constrained diffusion and binding energies of light particles and on the dust temperature. The abundance distributions of H₂CO and CH₃OH using the W32 cluster calculations and the experimental values derived by Fuchs et al. (2009) are similar, and both show a peak at abundances higher
<table><caption>**Table 3.** Energetic parameters used in the different astrochemical modeling simulations: energy barriers of the forward (E<sub>a,for</sub>) and reverse (E<sub>a,rev</sub>) reactions, and transition frequency of the transition state (ν<sup>‡</sup>).</caption><thead><tr><th>System</th><th>Reaction</th><th>E<sub>a,for</sub> (K)</th><th>E<sub>a,rev</sub> (K)</th><th>ν<sup>‡</sup> (cm<sup>-1</sup>)</th><th>Reference</th></tr></thead><tbody><tr><td rowspan="2">W32</td><td>CO + H → HCO</td><td>1106</td><td>11 290</td><td>569.7</td><td>this work</td></tr><tr><td>H<sub>2</sub>CO + H → CH<sub>3</sub>O</td><td>1299</td><td>15 880</td><td>701.3</td><td>this work</td></tr><tr><td rowspan="2">Gas-phase</td><td>CO + H → HCO</td><td>1979</td><td>8909</td><td>793.6</td><td>Peters et al. (2013)</td></tr><tr><td>H<sub>2</sub>CO + H → CH<sub>3</sub>O</td><td>2396</td><td>27 290</td><td>1037.0</td><td>Peters (2013)</td></tr><tr><td rowspan="2">laboratory ices</td><td>CO + H → HCO</td><td>390–520</td><td></td><td></td><td>Fuchs et al. (2009)</td></tr><tr><td>H<sub>2</sub>CO + H → CH<sub>3</sub>O</td><td>415–490</td><td></td><td></td><td>Fuchs et al. (2009)</td></tr></tbody></table>

than 1%, although the experimental values tend to give higher formaldehyde abundances and lower methanol abundances than the W32 calculations. Models using the higher gas-phase activation barriers tend to give lower abundances (≈2–5% compared to the ≈10–20% of the two models using ice data). Moreover, the peak of the abundance distribution of the gas-phase model at >1% is less pronounced, while the abundances are more spread over the entire range of abundances, down to 10<sup>-10</sup>. Comparisons between the modeled abundance distributions and interstellar ice observations (Oberg et al. 2011), depicted by the gray boxes in Fig. 7, were carried out. Models using the values derived by Fuchs et al. (2009) tend to predict slightly higher abundances of H<sub>2</sub>CO and CH<sub>3</sub>OH than the observed values. The grid using the W32 calculations agrees best with the H<sub>2</sub>CO observed abundances since its distribution peak is close to the range of observed values, whereas for CH<sub>3</sub>OH it predicts higher abundances on the same order using the values of Fuchs et al. (2009). In contrast, the peaks of the abundance distribution using the gas-phase energetic data agree well with the observations for both H<sub>2</sub>CO and CH<sub>3</sub>OH; however, they are more spread over the entire range of abundances and accordingly remain weaker than the models using ice data.

In principle, one could consider that the very same astrochemical modeling procedure could be used to reproduce and compare the experimental results obtained in laboratory ices. We tried to adapt our code to model ice experiments, but the procedure is not straightforward. First of all, in our models we consider the temperatures of the ice and the gaseous components to be equal, whereas in the experiments this is probably not the case. Indeed, in the experiments H atoms were remotely produced through either a microwave discharge or thermal cracking source at 2000 K in an external chamber and then cooled down through cold pipers; therefore, the actual temperature at which the reactions occur is poorly constrained. Moreover, the model needs to be adapted to other laboratory conditions such as the beam fluxes and ice coverages, which are significantly higher than in the ISM, and the multilayer regime considered in the experiments also requires a careful treatment.

## 5. Conclusions

Reactions by a successive H addition to CO toward the formation of HCO and CH<sub>3</sub>O as radical precursors for H<sub>2</sub>CO and CH<sub>3</sub>OH molecules were modeled using the BHLYP hybrid density functional method and the highly correlated CCSD(T) method. The reactions were simulated both in the gas phase and on the surfaces of relatively large water ice clusters consisting of 3, 18, and 32 molecules (W3, W18, and W32, respectively), mimicking the water ice present in interstellar dense molecular clouds.

One of the main conclusions is that the interaction of CO and H<sub>2</sub>CO with water cluster surfaces forming O-attached

complexes (i.e., H-bond through O atoms) slightly increases the C–O polarization, making the C atoms more prone to react with H atoms to form HCO and CH<sub>3</sub>O. This activation is shown by a tiny lowering of the energy barriers compared to the gas-phase processes. The energy barrier lowering is more pronounced with the H-bond strength, giving rise the trend of (from higher to lower energy barriers) gas phase >W3 > W18 ≈ W32. Despite this lowering, the calculated energy barriers are still exceedingly high to be overcome in the ISM conditions. Estimation of the cross-over temperature (142 for HCO and 175 K for CH<sub>3</sub>O formation) indicates that tunneling contributions can operate for these reactions.

Finally, simulations with the GRAINOBLE model to reproduce the abundances of CO, H<sub>2</sub>CO, and CH<sub>3</sub>OH using the calculated energy barriers and transition frequencies can reproduce the H<sub>2</sub>CO and CH<sub>3</sub>OH abundances with respect to water observed in low-mass protostars and dark cores. In contrast, when gas-phase data are adopted, a significant part of the distribution underpredicts the observations by several orders of magnitude.

Further work is certainly needed, since understanding the atomistic details of this chemical step is a key point in improving the actual phenomenological models that describe the CH<sub>3</sub>OH formation in the ISM conditions, such as application of the Feynman path integral method, which would mean a more rigorous description of the processes taking quantum effects (i.e., tunneling contributions) fully into account.

*Acknowledgements.* A.R. is indebted to MINECO of the Spanish Government for a Juan de la Cierva contract. Financial support from MINECO (projects CTQ2011-24847/BQU and CTQ2013-40347-ERC) and the use of the Catalonia Supercomputer Centre (CESCA) are gratefully acknowledged. PU acknowledge Progetti di Ricerca di Ateneo-Compagnia di San Paolo-2011-Linea 1A, progetto ORTO11RRT5 for funding. C.C. and V.T. acknowledge the funding by the French agency ANR, project FORCOMS ANR-08-BLAN-0225.

## References

Andersson, S., & Grüning, M. 2004, J. Phys. Chem. A, 108, 7621

Awad, Z., Chigai, T., Kimura, Y., Shalabiea, O. M., & Yamamoto, T. 2005, ApJ, 626, 262

Becke, A. D. 1993, J. Chem. Phys., 98, 1372

Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133

Cromer, D. T., Schiferl, D., leSar, R., & Mills, R. 1983, Acta Crystallograph. C, 1143

Dane, C. B., Lander, D. R., Curl, R. F., et al. 1988, J. Chem. Phys., 88, 2121

Dertinger, S., Geers, A., Kappert, J., Wiebrecht, J., & Temps, F. 1995, Faraday Discuss., 102, 31

Dunning, T. H. 1989, J. Chem. Phys., 90, 1007

Fermann, J. T., & Auerbach, S. 2000, J. Chem. Phys., 112, 6787

Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2004, Gaussian 03 (Wallingford CT: Gaussian Inc.)

Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629

Geers, A., Kappert, J., Temps, F., & Wiebrecht, J. W. 1993, J. Chem. Phys., 99, 2271

Goumans, T. P. M., Wander, A., Catlow, C. R. A., & Brown, W. A. 2007, MNRAS, 382, 1829
Goumans, T. P. M., Catlow, C. R. A., & Brown, W. A. 2008, J. Chem. Phys., 128, 4709

Hama, T., & Watanabe, N. 2013, Chem. Rev., 113, 8783

Hidaka, H., Watanabe, N., Shiraki, T., Nagaoka, A., & Kouchi, A. 2004, ApJ, 614, 1124

Hiraoka, K., Sato, T., Sato, S., et al. 2002, ApJ, 577, 265

Kamarchik, E., Rodrigo, C., Bowman, J. M., Reisler, H., & Krylov, A. I. 2012, J. Chem. Phys., 136, 4304

Koch, W., & Holthausen, M. C. 2001, A Chemist's Guide to Density Functional Theory (Weinheim: Wiley-VCH)

Koch, D. M., Toubin, C., Peslherbe, G. H., & Hynes, J. T. 2008, J. Phys. Chem. C, 112, 2972

Krishnan, R., Binkley, J. S., Seeger, R., & Pople, J. A. 1980, J. Chem. Phys., 72, 650

Lee, C., Yang, W., & Parr, R. G. 1988, Phys. Rev. B, 37, 785

Lundell, J. 1995, J. Phys. Chem., 99, 14290

Lundell, J., & Latajka, Z. 1997, J. Phys. Chem. A, 101, 5004

Lundell, J., & Rasanen, M. 1995, J. Phys. Chem., 99, 14301

Marcy, T. P., Díaz, R. R., Heard, D., et al. 2001, J. Phys. Chem. A, 105, 8361

Marenich, A. V., & Boggs, J. E. 2003, J. Phys. Chem. A, 107, 2343

McQuarrie, D. 1986, Statistical Mechanics (New York: Harper and Row)

Oberg, K. I., Van der Marel, N., Kristensen, L. E., & Van Dishoeck. E. F. 2011, ApJ, 740, 14

Parish, C. A., Augspurger, J. D., & Dykstra, C. E. 1992, J. Phys. Chem., 96, 2069

Peters, P. S. 2013, Ph.D. Thesis, University of Lille 1

Peters, P. S., Duflot, D., Faure, A., et al. 2011, J. Phys. Chem. A, 115, 8983

Peters, P. S., Duflot, D., Wiesenfeld, L., & Toubin, C. 2013, J. Chem. Phys., 139, 4310

Petraco, N. D. K., Allen, W. D., & Schaefer III, H. F. 2002, J. Chem. Phys., 116, 10229

Pirim, C., & Krim, L. 2011, Phys. Chem. Chem. Phys., 13, 19454

Pirim, C., Krim, L., Laffon, C., et al. 2010, J. Phys. Chem. A, 114, 3320

Poater, J., Solàà, M., Rimola, A., Rodríguez-Santiago, L., & Sodupe, M. 2004, J. Phys. Chem. A, 108, 6072

Pople, J. A., Head-Gordon, M., & Raghavachari, K. 1987, J. Chem. Phys., 87, 5968

Reed, A. E., Weinhold, F., Curtiss, L. A., & Pochatko, D. J. 1986, J. Chem. Phys., 84, 5687

Reilly, J. P., Clark, J. H., Moore, C. B., & Pimentel, G. C. 1978, J. Chem. Phys., 69, 4381

Rimola, A., Rodríguez-Santiago, L., & Sodupe, M. 2006, J. Phys. Chem. B, 110, 24189

Rimola, A., Constantino, E., Rodríguez-Santiago, L., & Sodupe, M. 2008, J. Phys. Chem. A, 112, 3444

Rimola, A., Sodupe, M., & Ugliengo, P. 2010, Phys. Chem. Chem. Phys., 12, 5285

Rumbles, G., Lee, E. K. C., & Valentini, J. J. 1990, J. Chem Soc. Faraday Trans., 86, 3837

Sadlej, J., Rowland, B., Devlin, J. P., & Buch, V. 1995, J. Chem. Phys., 102, 4804

Saebo, S., Radom, L., & Schaefer III, H. F. 1983, J. Chem. Phys., 78, 845

Sappey, A. D., & Crosley, D. R. 1990, J. Chem. Phys., 93, 7601

Sholl, D., & Steckel, J. A. 2009, Density Functional Theory: A Practical Introduction (New Jersey: Wiley, Hoboken)

Sodupe, M., Bertran, J., Rodríguez-Santiago, L., & Baerends, E. J. 1999, J. Phys. Chem. A, 103, 166

Sousa, S. F., Fernandes, P. A., & Ramos, M. J. 2007, J. Phys. Chem. A, 111, 10439

Taketsugu, T., Tajima, N., & Hirao, K. 1996, J. Chem. Phys., 105, 1933

Taquet, V., Ceccarelli, C., & Kahane, C. 2012a, A&A, 538, A42

Taquet, V., Ceccarelli, C., & Kahane, C. 2012b, ApJ, 748, L3

Taquet, V., Peters, P. S., Kahane, C., et al. 2013, A&A, 550, A127

Taquet, V., Charnley, S., & Sipilä 2014, ApJ, 791, 1

Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245

Tielens, A. G. G. M., Tokunaga, A. T., Geballe, T. R., & Baas, F. 1991, ApJ, 381, 181

van Mourik, T., Dunning, T. H., & Peterson, K. A. 2000, J. Phys. Chem. A, 104, 2287

Walch, S. P. 1993, J. Chem. Phys., 98, 3076

Wang, H. Y., Eyre, J. A., & Dorfman, L. M. 1973, J. Chem. Phys., 59, 5199

Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173

Watanabe, N., & Kouchi, A. 2008, Progress Surf. Sci., 83, 439

Watanabe, N., Shiraki, T., & Kouchi, A. 2003, ApJ, 588, L121

Watanabe, N., Mouri, O., Nagaoka, A., et al. 2007, ApJ, 668, 1001

Whittet, D. C. B., Cook, A. M., Herbst, E., Chiar, J. E., & Shenoy, S. S. 2011, ApJ, 742, 28

Woon, D. E. 1996, J. Chem. Phys., 105, 9921

Woon, D. E. 2002, ApJ, 569, 541

## Appendix A: Geometries and energetics for CO hydrogenation for one water molecule

<figure><img src="image_9.png" /><figcaption>**Fig. A.1.** BHLYP/6-311++G(d,p)-optimized geometries of the stationary points involved in the H additions to CO for one water molecule (W1): panel **a**) hydrogenation of CO on the C atom when CO interacts with W1 through the O atom; panel **b**) hydrogenation of CO on the C atom when CO interacts with W1 through the C atom; panel **c**) hydrogenation of H₂CO on the C atom; panel **d**) hydrogenation of CO on the O atom when CO interacts with W1 through the C atom; and panel **e**) hydrogenation of H₂CO on the O atom. AS refers to the CO···W1 + H zero-energy reference for panel **a**; to the OC···W1 + H zero-energy references for panels **b** and **d**; and to the H₂CO···W3 + H zero-energy references for panels **c** and **e**. TS refers to the transition-state structures. PROD refers to the products. Distances in Å.</figcaption></figure>
<table><caption>**Table A.1.** Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO and CH₃O for one water molecule.</caption><thead><tr><th rowspan="2">Reaction</th><th colspan="3">L1</th><th colspan="3">L2</th></tr><tr><th>AS</th><th>TS</th><th>PROD</th><th>AS</th><th>TS</th><th>PROD</th></tr></thead><tbody><tr><td>a) CO···W1 + H → HCO···W1</td><td>0.0<br>[0.0]</td><td>7.7<br>[10.7]</td><td>-108.3<br>[-81.6]</td><td>0.0</td><td>12.5</td><td>-90.4</td></tr><tr><td>b) OC···W1 + H → HCO···W1</td><td>0.0<br>[0.0]</td><td>9.5<br>[11.5]</td><td>-106.4<br>[-81.6]</td><td>0.0</td><td>14.0</td><td>-86.8</td></tr><tr><td>c) H₂CO···W1 + H → CH₃O···W1</td><td>0.0<br>[0.0]</td><td>9.6<br>[13.6]</td><td>-144.9<br>[-117.9]</td><td>0.0</td><td>13.6</td><td>-111.0</td></tr><tr><td>d) OC···W1 + H → HOC···W1</td><td>0.0<br>[0.0]</td><td>118.4<br>[122.4]</td><td>+57.3<br>[+83.5]</td><td>0.0</td><td>130.4</td><td>+81.8</td></tr><tr><td>e) H₂CO···W1 + H → CH₂OH···W1</td><td>0.0<br>[0.0]</td><td>35.7<br>[39.1]</td><td>-161.9<br>[-134.8]</td><td>0.0</td><td>43.5</td><td>-143.5</td></tr></tbody></table>

**Notes.** Calculations at BHLYP/6-311++G(d,p) (L1) and at CCSD(T)/aug-cc-pVTZ//BHLYP/6-311++G(d,p) (L2) levels. ZPE corrections were calculated at the L1 level and included to the L1 and L2 potential energies. Values in units of kJ mol⁻¹. AS refers to the zero-energy references, TS to the transition-state structures and, PROD to the products. Geometries and labeling of the considered cases are the same as for Fig. A.1.

## Appendix B: Geometries and energetics for the interaction of CO with the water clusters

<table><caption>**Table B.1.** Interaction energies (ΔE<sub>int</sub>) and relative energies (ΔE<sub>rel</sub>), and including the zero-point energy (ZPE) corrections (ΔE<sub>int</sub> + ZPE and ΔE<sub>rel</sub> + ZPE) for the different water-carbon monoxide systems.</caption><thead><tr><th>Complex</th><th>ΔE<sub>int</sub></th><th>ΔE<sub>rel</sub></th><th>ΔE<sub>int</sub> + ZPE</th><th>ΔE<sub>rel</sub> + ZPE</th></tr></thead><tbody><tr><td>CO···W1</td><td>-4.7</td><td>1.9</td><td>-2.5</td><td>0.0</td></tr><tr><td>OC···W1</td><td>-6.6</td><td>0.0</td><td>-2.5</td><td>0.0</td></tr><tr><td>CO···W3</td><td>-4.6</td><td>2.0</td><td>-2.7</td><td>0.7</td></tr><tr><td>OC···W3</td><td>-6.6</td><td>0.0</td><td>-3.4</td><td>0.0</td></tr><tr><td>CO···W18</td><td>-5.6</td><td>2.9</td><td>-3.2</td><td>1.6</td></tr><tr><td>OC···W18</td><td>-8.5</td><td>0.0</td><td>-4.8</td><td>0.0</td></tr><tr><td>CO···W32</td><td>-14.9</td><td>-1.1</td><td>-10.8</td><td>-1.6</td></tr><tr><td>OC···W32</td><td>-13.8</td><td>0.0</td><td>-9.2</td><td>0.0</td></tr><tr><td>Experimental</td><td></td><td></td><td>-9.5</td><td></td></tr></tbody></table>

**Notes.** In kJ mol⁻¹. Calculated at BHLYP/6-311++G(d,p) (L1). The experimental binding energy of CO on water ice measured by Collings et al. (2004) is also included. Geometries and labeling of the considered cases are the same as in Fig. B.1.
<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** BHLYP/6-311++G(d,p)-optimized geometries of the C-attached and the O-attached complexes for the different water-carbon monoxide systems. Distances in Å. Relative energies are provided in Table B.1.</figcaption></figure>