Free Access
Issue
A&A
Volume 572, December 2014
Article Number A70
Number of page(s) 12
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201424046
Published online 01 December 2014

© ESO, 2014

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 CO2 are present in molecular clouds at visual extinctions Av ~ 3 mag, whereas solid methanol is only present at Av ≳ 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): 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 H2CO may react with CO to produce HCO. The activation energy was also obtained through an Arrhenius plot between 260340 K using fast Lyman-α spectrophotometry to be about 8.3 ± 2 kJ mol-1 (Wang et al. 1973), whereas energy barriers calculated using quantum chemical methods span the 1521 kJ mol-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), quantum chemical calculations give barriers of 2122 kJ mol-1 (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 H2OCO 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, CH3OH was formed at a lower yield than H2CO. Furthermore, detection of HCO and CH3O 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 H2CO and CH3OH 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 Pr 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 Eb 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 H2O, CO, CO2, H2CO, and CH3OH. The density of atomic H was fixed to 1.2 cm-3, while the abundances of O and CO evolve with time. The initial abundance of CO was set to 5 × 10-5 cm-3, 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 (nH ≳ 105 cm-3), and high visual extinctions (AV ≳ 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 temperaturesTg = Td 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 ad 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 Ed/Eb 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 × 10-5 and 1.2 × 10-4. Most of the atomic O freezing-out onto grains is converted into water, giving abundances of ~ 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, ΔE, and reaction energies, ΔEr) 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 (ΔE+ZPE) are also included. Hereafter, for the sake of clarity, an energetic value calculated at a given theory level is denoted as ΔE/L. For instance, ΔE+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 kJ mol-1 = 120.274 K.

thumbnail 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 H2CO on the C atom, and panel d) hydrogenation of H2CO on the O atom. AS refers to the CO + H zero-energy references for panels a) and b) and to the H2CO + H zero-energy references for panels c) and d). TS refers to the transition-state structures and PROD to the products. Distances in Å.

Open with DEXTER

Table 1

Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO, COH, CH3O, and CH2OH in gas phase.

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; 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 (ΔE/L2 and ΔEr/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 ΔE of the HCO formation; − 80.4 vs. − 81 kJ mol-1 for ΔEr of the HCO formation; 133.1 vs. 137 kJ mol-1 for ΔE of the COH formation; and 88.8 vs. 94 kJ mol-1 for ΔEr of the COH formation.

As mentioned in the Introduction, the addition of a second H atom (i.e., H + HCO H2CO) is a barrierless radical-radical reaction and, accordingly, it was not considered here. The addition of a third H atom (i.e., H + H2CO) can give either CH3O or CH2OH, the calculated structures involved in the processes are presented in Figs. 1c and d. CH2OH is a C-radical centered species with a nonplanar geometry (the C atom adopts a pyramidal structure), whereas the unpaired electron of CH3O is fully localized on the O atom. Because of that, CH2OH is found to be more stable than CH3O by about 35 kJ mol-1, although both processes are largely exoenergetic (ΔE+ZPE/L2 = − 120.1 kJ mol-1 and − 85.7 kJ mol-1). However, the calculated energy barriers go against this trend, as ΔE+ZPE/L2 is higher for CH2OH formation than that for CH3O (44.6 kJ mol-1 and 18.6 kJ mol-1), 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 CH3O into CH2OH 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 CH3O dissociation barrier. Accordingly, and considering the cryogenic temperatures of the interstellar medium at which the processes occur, formation of CH2OH is kinetically hampered in favor of CH3O formation. This conclusion agrees with the experimental exclusive detection of CH3O (Pirim et al. 2010). Remarkably, calculated ΔE/L2 values for CH3O and CH2OH formation (14.5 kJ mol-1 and 41.1 kJ mol-1) agree reasonably well with those obtained by Kamarchik et al. (2012) at the CCSD(T)/aug-cc-pvtz (22.6 kJ mol-1 and 47.1 kJ mol-1). 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 ΔEr values by about 4 kJ mol-1 and 30 kJ mol-1 at most.

thumbnail 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 H2CO 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 H2CO 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 H2CO···W3 + H zero-energy references for panels c) and e). TS refers to the transition-state structures and PROD to the products. Distances in Å.

Open with DEXTER

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.

Table 2

Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO and CH3O for a three-water-molecule cluster (W3).

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 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/L2 value (see Table 2 line a) is only marginally lower than that computed for the same gas-phase reaction (12.6 kJ mol-1 vs. 13.0 kJ mol-1). The same tiny energy barrier lowering was also computed for the H addition to H2CO to give CH3O (ΔE/L2 = 14.1 kJ mol-1 and 14.5 kJ mol-1 for W3 and in the gas phase, see Table 2 line c). The tiny ΔE decreases are due to the polarization of the CO bonds induced by the H-bond interaction with the cluster model. This is shown by the larger CO distances compared with those in the gas phase (1.117 Å and 1.115 Å for CO; 1.194 Å and 1.189 Å for H2CO), and a consequent CO bond weakening, which favors the H additions. ΔEr values for both HCO and H3CO 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/L2 is slightly higher than that obtained in the gas phase (14.0 kJ mol-1 vs. 13.0 kJ mol-1). Here, the H-bond with the cluster model induces a shortening of the CO distance compared to its gas-phase value (1.113 Å vs. 1.115 Å), strengthening the CO 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 CO 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 (hypsochromic) 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 CH2OH isomer (see Fig. 2e and Table 2 line e) are very similar to those computed in the gas phase (i.e., CH2OH is thermodynamically favored, but kinetically hindered compared to CH3O, see above). All attempts to locate activated complexes converting CH3O to CH2OH 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.

thumbnail Fig. 3

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO for a 18-water ice cluster (W18). Formation of HCO a) and CH3O 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) H2CO···W18 + H. Bond lengths in Å.

Open with DEXTER

thumbnail Fig. 4

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO for a 32-water ice cluster (W32). Formation of HCO a) and CH3O 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) H2CO···W32 + H. Bond lengths in Å.

Open with DEXTER

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/L1 values (6.8 kJ mol-1 and 6.6 kJ mol-1 for W18 and W32), which are the lowest values calculated in this work for the first hydrogenation reaction. The corresponding ΔEr/L1 indicates that both H additions are largely exoenergetic and almost independent of the size of the water ice cluster.

The CH3O formation is influenced by the different H-bond interactions between H2CO and the W18 and W32 ice model clusters, as shown in Figs. 3b and 4b. For W18, the calculated ΔE is 8.9 kJ mol-1, 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-1, the lowest value calculated in this work. This decrease is due to the two moderately strong H-bonds holding H2CO in place in the W32 ice model (Fig. 4b), which increases the CO polarization and thereby reduces the energy barrier value.

As mentioned above, hydrogenation of HCO and CH3O to form the final H2CO and CH3OH 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 H2CO and CH3OH 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), the system collapses into the H2CO···W18 (see Fig. 5, singlet state) because of the spontaneous formation of H2CO, indicating that the H coupling to HCO to give H2CO 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.

thumbnail 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 Å.

Open with DEXTER

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 CH3O 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.

thumbnail Fig. 6

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO in the presence of a 4-CO molecule cluster. Formation of HCO a) and CH3O 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 + H2CO···COice. Bond lengths in Å.

Open with DEXTER

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-1). The corresponding reaction to give the COH isomer exhibits a very high kinetic barrier (ΔE+ZPE/L1 = 121.2 kJ mol-1) and a highly endo-energetic character. The formation of H2CO from HCO is obviously barrierless (radical-radical reaction). In contrast, the reaction of H + H2CO brings about either CH2OH or CH3O, both of which are exo-energetic processes with a definite preference for the CH2OH molecule. Nevertheless, the calculated kinetic barrier ΔE associated with the formation of CH3O is lower than that of CH2OH (ΔE+ZPE/L1 = 10.2 kJ mol-1 and 32.9 kJ mol-1) and, therefore, formation of CH2OH 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 CO bond, thus slightly reducing the kinetic barrier for the formation of HCO (ΔE+ZPE/L1 = 9.7 kJ mol-1), whereas for the C-attached complex the CO bond is strengthened with a corresponding increase of the kinetic barrier for the H addition (ΔE+ZPE/L1 = 11.1 kJ mol-1). The polarization of the CO 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-1 and 9.2 kJ mol-1) because of the stronger H-bond interactions and the increased CO polarization. Similarly, the relatively strong H-bond interactions between the surface water molecules of the ice cluster models with the O atom of H2CO is also the reason why the formation of CH3O exhibits lower energy barriers than the gas-phase process (ΔE+ZPE/L1 = 13.9 kJ mol-1, 12.8 kJ mol-1, and 10.8 kJ mol-1 on W3, W18 and W32, respectively).

Despite the effect of water in reducing the energy barriers for the HCO and CH3O 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-1 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-1 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, irrespectively 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 meagre 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 tunneling 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 (). A characteristic parameter is the tunneling cross-over temperature (Tx), below which tunneling becomes dominant and above which tunneling becomes negligible, which can be calculated using the formula of Fermann & Auerbach (2000) (5)where ω = 2πν, ħ= h/ 2π and h is the Planck constant, is the energy barrier including zero-point energy corrections, and kB is the Boltzman constant. From the computed ν and for reactions on W32, we calculated the Tx values. For the formation of HCO, Tx is 142 K, whereas for CH3O 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 relevant. Because of that, future quantum mechanical-based works should certainly focus on using more rigorous methods for reaction 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 atomistic 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 astrochemical 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 barriers 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 consequently 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 relative abundances of CO, H2CO, and CH3OH, we ran three grids of 3125 models each using the GRAINOBLE model by combining the input parameters described in Sect. 2.2; that is, gas temperatures, grain diameters, diffusion-to-binding energy ratios, 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 transmission probabilities: i) those calculated with the W32 cluster model of this work; ii) those calculated by Peters (2013) and Peters et al. (2013) for gas phase conditions; and iii) those measured on laboratory ices by Fuchs et al. (2009) and by introducing the same formalism as in Monte Carlo chemical models 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 reactions (Ea,for and Ea,rev) and the transition frequency of the transition state (ν). The energy barriers of Fuchs et al. (2009), as explained in their paper, are not real energies, but effective energies, 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 1012 s-1.

Table 3

Energetic parameters used in the different astrochemical modeling simulations: energy barriers of the forward (Ea,for) and reverse (Ea,rev) reactions, and transition frequency of the transition state (ν).

thumbnail Fig. 7

Distributions of the ice abundances of H2CO (upper panel) and CH3OH (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.

Open with DEXTER

The abundance distributions (with respect to water) of H2CO and CH3OH 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 H2CO and CH3OH 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 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 (25% compared to the 1020% 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-10. 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 H2CO and CH3OH than the observed values. The grid using the W32 calculations agrees best with the H2CO observed abundances since its distribution peak is close to the range of observed values, whereas for CH3OH 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 H2CO and CH3OH; 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 straighforward. 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 CH3O as radical precursors for H2CO and CH3OH 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 H2CO with water cluster surfaces forming O-attached complexes (i.e., H-bond through O atoms) slightly increases the CO polarization, making the C atoms more prone to react with H atoms to form HCO and CH3O. 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 CH3O formation) indicates that tunneling contributions can operate for these reactions.

Finally, simulations with the GRAINOBLE model to reproduce the abundances of CO, H2CO, and CH3OH using the calculated energy barriers and transition frequencies can reproduce the H2CO and CH3OH 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 CH3OH 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.

Acknowledgments

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

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

thumbnail 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 H2CO 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 H2CO 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 H2CO···W3 + H zero-energy references for panels c) and e). TS refers to the transition-state structures. PROD refers to the products. Distances in Å.

Open with DEXTER

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 CH3O for one water molecule.

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

Table B.1

Interaction energies (ΔEint) and relative energies (ΔErel), and including the zero-point energy (ZPE) corrections (ΔEint + ZPE and ΔErel + ZPE) for the different water-carbon monoxide systems.

thumbnail 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.

Open with DEXTER

All Tables

Table 1

Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO, COH, CH3O, and CH2OH in gas phase.

Table 2

Relative potential energies (bare values) and with the zero-point energy (ZPE) corrections (values in brackets) for the formation of HCO and CH3O for a three-water-molecule cluster (W3).

Table 3

Energetic parameters used in the different astrochemical modeling simulations: energy barriers of the forward (Ea,for) and reverse (Ea,rev) reactions, and transition frequency of the transition state (ν).

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 CH3O for one water molecule.

Table B.1

Interaction energies (ΔEint) and relative energies (ΔErel), and including the zero-point energy (ZPE) corrections (ΔEint + ZPE and ΔErel + ZPE) for the different water-carbon monoxide systems.

All Figures

thumbnail 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 H2CO on the C atom, and panel d) hydrogenation of H2CO on the O atom. AS refers to the CO + H zero-energy references for panels a) and b) and to the H2CO + H zero-energy references for panels c) and d). TS refers to the transition-state structures and PROD to the products. Distances in Å.

Open with DEXTER
In the text
thumbnail 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 H2CO 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 H2CO 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 H2CO···W3 + H zero-energy references for panels c) and e). TS refers to the transition-state structures and PROD to the products. Distances in Å.

Open with DEXTER
In the text
thumbnail Fig. 3

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO for a 18-water ice cluster (W18). Formation of HCO a) and CH3O 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) H2CO···W18 + H. Bond lengths in Å.

Open with DEXTER
In the text
thumbnail Fig. 4

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO for a 32-water ice cluster (W32). Formation of HCO a) and CH3O 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) H2CO···W32 + H. Bond lengths in Å.

Open with DEXTER
In the text
thumbnail 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 Å.

Open with DEXTER
In the text
thumbnail Fig. 6

BHLYP/6-311++G(d,p)-energy profiles (in kJ mol-1) for the H addition to CO in the presence of a 4-CO molecule cluster. Formation of HCO a) and CH3O 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 + H2CO···COice. Bond lengths in Å.

Open with DEXTER
In the text
thumbnail Fig. 7

Distributions of the ice abundances of H2CO (upper panel) and CH3OH (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.

Open with DEXTER
In the text
thumbnail 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 H2CO 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 H2CO 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 H2CO···W3 + H zero-energy references for panels c) and e). TS refers to the transition-state structures. PROD refers to the products. Distances in Å.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.