The Astrophysical Journal, 917:49 (11pp), 2021 August 10 https://doi.org/10.3847/1538-4357/ac0142 © 2021. The American Astronomical Society. All rights reserved. H2 Formation on Interstellar Grains and the Fate of Reaction Energy Stefano Pantaleone1,6 , Joan Enrique-Romero1,2 , Cecilia Ceccarelli1 , Stefano Ferrero2 , Nadia Balucani1,3,4 , Albert Rimola2 , and Piero Ugliengo5 1 Université Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), F-38000 Grenoble, France stefano.pantaleone@univ-grenoble-alpes.fr 2 Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, E-08193, Catalonia, Spain 3 Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, Via Elce di Sotto 8, I-06123 Perugia, Italy 4 Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy 5 Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, I-10125, Torino, Italy Received 2021 March 26; revised 2021 April 27; accepted 2021 May 12; published 2021 August 16 Abstract Molecular hydrogen is the most abundant molecular species in the universe. While no doubts exist that it is mainly formed on the interstellar dust grain surfaces, many details of this process remain poorly known. In this work, we focus on the fate of the energy released by the H2 formation on the dust icy mantles: how it is partitioned between the substrate and the newly formed H2, a process that has a profound impact on the interstellar medium. We carried out state-of-the-art ab initio molecular dynamics simulations of H2 formation on periodic crystalline and amorphous ice surface models. Our calculations show that up to two-thirds of the energy liberated in the reaction (∼300 kJ mol−1 ∼3.1 eV) is absorbed by the ice in less than 1 ps. The remaining energy (∼140 kJ mol−1 ∼1.5 eV) is kept by the newly born H2. Since it is 10 times larger than the H2 binding energy on the ice, the new H2 molecule will eventually be released into the gas phase. The ice water molecules within ∼4 Å from the reaction site acquire enough energy, between 3 and 14 kJ mol−1 (360–1560 K), to potentially liberate other frozen H2 and, perhaps, frozen CO molecules. If confirmed, the latter process would solve the long standing conundrum of the presence of gaseous CO in molecular clouds. Finally, the vibrational state of the newly formed H2 drops from highly excited states (ν = 6) to low (ν „ 2) vibrational levels in a timescale of the order of picoseconds. Unified Astronomy Thesaurus concepts: Astrochemistry (75) 1. Introduction astrochemistry and, more generally, the physics and chemistry of the interstellar medium (ISM; Duley & Williams 1993). Molecular hydrogen is the most abundant molecule in the Answers to the above questions largely depend on the nature universe. Its formation is also the first step of the interstellar of the substrate, namely the specific dust grain surface or, in chemistry and, therefore, a fundamental reaction. In molecular practice, the interstellar environment where H2 forms. Here we clouds, H2 is mainly formed via the H + H association reaction focus on cold (∼10 K) molecular clouds. In these environ- on the interstellar dust grain surfaces, which act as a third body ments, the grain refractory cores are coated by water-dominated capable of partially absorbing the energy (∼440 kJ mol−1 icy mantles either in polycrystalline or amorphous (ASW for ∼4.56 eV) released by the chemical reaction and, consequently, amorphous solid water) phases, the latter dominating over the stabilize the newly formed H2 molecule (e.g., van de Hulst former (Smith et al. 1988; Whittet 1993; Boogert et al. 2015). 1946; Hollenbach & Salpeter 1970, 1971; Smoluchowski Unfortunately, laboratory experiments cannot entirely repro- 1983). Although no doubts exist regarding the occurrence of duce the interstellar conditions so that, while they can suggest this process, many specific details remain poorly known (see, routes and processes, they cannot provide definitive answers to e.g., Vidali 2013; Wakelam et al. 2017). the above-exposed questions (see, e.g., Vidali 2013). For Here, we focus on the fate of the energy released by the example, obtaining experimentally how the H2 formation reaction, which has been a source of debate for decades. nascent energy is partitioned is very difficult (Hornekær et al. Specifically, how much of the reaction energy is absorbed by 2003; Roser et al. 2003; Watanabe et al. 2010; Hama et al. the dust grain and in what timescale? What fraction of the 2012). Computational chemistry methods can be regarded as a chemical energy does go into the kinetic energy of the newly complementary tool and, sometimes, as a unique alternative to formed H2? Is this energy large enough as H2 breaks the laboratory experiments. So far, a limited number of studies interaction with the surface and proceeds into the gas phase? devoted to the energy dissipation of the H2 formation reaction How much energy does the H2 molecule possess when it leaves on icy grains has appeared in the literature (Takahashi et al. the grain surface and in what rovibrational state? And, finally, 1999a; Takahashi & Uehara 2001; Herbst & Cuppen 2006). is the energy transmitted to the dust grain enough to locally Here, we present a new theoretical study on the dissipation warm it up and make adjacent frozen molecules sublimate? of the energy released by the H + H → H2 reaction on water ice These points have a profound impact on several aspects of using ab initio molecular dynamics simulations (AIMDs). We 6 simulated the reaction adopting a Langmuir–Hinshelwood Present addresses: Dipartimento di Chimica and Nanostructured Interfaces mechanism (i.e., both reactants are adsorbed on the surface in and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, I-10125 Torino, Italy and Dipartimento di Chimica, Biologia e Biotecnologie, neighboring sites) both on amorphous and crystalline periodic Università di Perugia, Via Elce di Sotto 8, I-06123 Perugia, Italy. models of interstellar icy mantles. 1 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. 2. Computational Details 2.1. Methods All our calculations were carried out with the CP2K package (Hutter et al. 2014). Core electrons were described with the Goedecker–Teter–Hutter pseudopotentials (Goedecker et al. 1996) and valence ones with a mixed Gaussian and plane-wave (GPW) approach (Lippert et al. 1997). The density functional theory (DFT) Perdew–Burke–Ernzerhof (PBE) method was adopted (Perdew et al. 1996), combined with a triple-ζ basis set for valence electrons with a single polarization function (TZVP), and the a posteriori D3 Grimme correction to account for dispersion forces (Grimme et al. 2010, 2011). The plane- wave cutoff was set to 600 Ry. Ice surfaces were thermalized (see Section 2.2), and during the geometry optimization only the reacting H atoms (those forming H2) were allowed to move, keeping the water molecules fixed at their thermalized positions. All calculations were carried out within the unrest- ricted formalism as we deal with open-shell systems (see the Appendix for more details). The binding energy (BE) of H2 was calculated as: BE H2 = E CPLX - (EIce + E H2) , (1 ) where ECPLX is the energy of the H2/ice system, EIce is the energy of the bare ice surface, and EH2 is the energy of the H2 alone. AIMDs were carried out within the NVE (N = number of Figure 1. Top view of the crystalline ice model (a), and top (b) and side ((c) particles, V = volume, E = energy) ensemble, where the total and (d)) views of the ASW ice models. The yellow circles represent the centers energy VTOT (i.e., potential + kinetic) is conserved. The of mass of the two H atoms’ starting positions. In the case of ASW, the numbers mark the positions (Pos1, Pos2, and Pos3) discussed in the text. The evolution of the system was followed for 5 ps for the crystalline initial H to H distance is about 4 Å (see Figure A1). O atoms are in red and H model and for 1 ps for the amorphous model, using time steps atoms of the ice are in white. of 0.2 fs. ranging from 0.3 to 3.0 Å with a step of 0.01 Å . Computed data 2.2. Ice Models were fitted with the Morse equation to obtain the force constant For the crystalline ice model, a periodic slab cut from the of the oscillator, its dissociation energy and, hence, the hexagonal ice bulk structure was used. The periodic cell vibrational levels of H2. The H2 turning points calculated with parameters defining the system are a = 26.318 Å and our model were compared with those of the AIMD simulations b = 28.330 Å , while the slab thickness is ∼21 Å (corresp- and the vibrational level was assigned at each H2 oscillation onding to seven layers), for a total of 576 water molecules. The during the AIMDs. c parameter of the simulation box, i.e., the nonperiodic one, was set to 50 Å to avoid interactions among the fictitious slab 3. Results replicas. The size of the slab model was chosen to avoid nonphysical temperature increase due to the extremely large 3.1. Reactants Plus Product Positions in the Ice Models exothermic reaction (Pantaleone et al. 2020). Figure 1 shows the structure of the used crystalline and ASW The amorphous model was obtained by performing a ice models as well as the starting positions of the center of mass classical molecular dynamics (MD) simulation on the crystal- of the two adsorbed H atoms before the reaction. The starting line structure. The MD was carried out with the TIP3P force H–H distance, as well as its evolution during the AIMDs is field (Jorgensen et al. 1983) for 200 ps (with a time step of shown in Figure A1. The ASW model presents a rugged 0.5 fs) at 300 K within the NVT (N = number of particles, morphology (e.g., holes, cavities, and channels) in comparison V = volume, T = Temperature) ensemble. A second NVT to the crystalline one and, therefore, we carried out simulations simulation was performed at 10 K to cool down the system at in three positions, roughly representing different possible the temperature of cold molecular clouds. Finally, a geometry situations in terms of energetics and surface morphology. The optimization and another NVT-MD simulation at 10 K were position marked as Pos1 is the deepest one, within a cavity, run using PBE to recover the potential energy surface at the while positions Pos2 and Pos3 are in the outermost parts of the same theory level as described in the previous section. surface. As a first step, we optimized the geometries of reactants (the 2.3. H2 Vibrational State two H atoms) and product (the H2 molecule) on the water ice To calculate the H2 vibrational state during the AIMDs the surface in order to obtain the potential energy surface of the anharmonic oscillator model was employed. As a first step, the reaction. No search for transition state structures has been potential energy surface (PES) of the H2 isolated molecule was performed because AIMDs indicated this reaction is barrierless explored by performing a rigid scan of the H–H distance, at 10 K. That is, the starting kinetic energy of the two H atoms 2 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Table 1 Results of the AIMD for the H2 Formation on the Crystalline (First Row) and the ASW Three Position (Figure 1; Bottom Three Rows) Ices Kice KH2 Model Kice + KH2 Kice + KH2 BE H2 ν0.2−1 ν1−5 Epeak (kJ mol−1) (kJ mol−1) Crys. 0.45 0.55 9.0 4–5 1–2 12.0 Pos1 0.65 0.35 9.7 1–2 L 43.0 Pos2 0.46 0.54 10.6 6 L 7.0 Pos3 0.47 0.53 12.4 5 L 7.0 Note. Second and third columns report the fraction of kinetic energy of the ice and the newly formed H2, respectively, at the end of the simulation (5 and 1 ps for crystalline and ASW ices, respectively). Fourth column reports the H2 binding energy (in kilojoules per mole) at the reaction site. Fifth and sixth columns report the vibrational state averaged over 0.2–1 ps and 1–5 ps, respectively (the latter can be computed only for the crystalline case). Last column reports the peak energy (in kilojoules per mole) of an ice-molecule neighbor (i.e., at ∼3 Å) to where the reaction occurs. contrast, on the ASW ice model, H2 diffuses over all three directions, also along the c-axis, which corresponds to direction perpendicular to the ice surface. However, the movement is not the same in the three studied positions. In Pos1, H2 achieves a maximum height of ∼10 Å over the surface in a timescale of Figure 2. PBE-D3/TZVP optimized geometries of reactants (left panels) and 0.4 ps and does not come back (Figure 3(d)), thus definitely product (right panels) of the H2 reaction formation on the crystalline (top) and leaving the surface. On the contrary, in Pos2 H2 slides and hops ASW Pos1 (bottom) ice models. The numbers in parentheses correspond to the on the surface within the 1 ps of the simulation (Figure 3(f)). A relative energy in kJ mol−1 with respect to the reactants. The H atoms involved similar behavior is observed for Pos3, but with wider hops and in the H+H → H2 reaction are in yellow, those belonging to surface water jumps (Figure 3(h)). This is the key point of the simulations: molecules are in white, and O atoms are in red. the kinetic energy remains much larger (∼150 kJ mol−1) than the binding energy (∼10 kJ mol−1) in all cases; therefore, the provided by the 10 K is high enough to overcome the nascent H2 is likely doomed to leave the ice surface. eventually small potential energy barrier. The total energy to be dissipated is around 435–440 kJ mol−1, depending on the surface and starting position. 3.3. Energy Dissipation toward the Ice Water Molecules Figure 2 provides a view of the two H atoms’ starting positions and H2 position immediately after the reaction, on Figure 4 shows how the energy absorbed by the ice is first the crystalline and ASW (Pos1 as an example) models, transmitted from the reaction site to a neighbor water (i.e., at respectively. ∼4 Å) in ∼200–800 fs and, with time, to shells of water molecules with increasing distances. The neighbor water molecule acquires from ∼3 to ∼14 kJ mol−1 and stays with 3.2. Molecular Hydrogen Desorption that energy for more than about 100–200 fs before the energy is AIMDs results are summarized in Table 1 and shown in dissipated toward more external water molecules. These Figure 3. First, at least half of the kinetic energy is absorbed by timescales are in excellent agreement with previous studies of the ice while the remaining one is kept by the newly formed H2 the sound speed in ices (see, e.g., Ruocco et al. 1996). molecule (∼50%–35%). These numbers are similar for both for the crystalline and ASW ice positions, suggesting they do no 3.4. H2 Vibrational State depend much on the surface structural details of the ice. Second, after the H–H bond formation, the newly formed H2 The initial state of the H2 formed on the crystalline ice is molecule diffuses over the surface, as it keeps a large kinetic ν = 4–5 (obtained averaging over 0.2–1 ps) and then it energy, in both crystalline and ASW ice models. However, decreases to ν = 1–2 (1–5 ps average) because of the H2 despite the energetics of the two processes being similar energy dissipation. Conversely, at Pos1 of the ASW model, H2 (Figure 3, left panels) the H2 diffusion (Figure 3, right panels) is formed with a low vibrational state (ν = 1–2, 0.2–1 ps is different. average). On the crystalline surface, the diffusion of the newly formed For the other two positions (Pos2 and Pos3), the vibrational H2 over the ice surface is constrained in a specific direction, state of the newly formed H2 molecule lies in highly excited along the a-axis, parallel to the ice surface, whereas in the b- states, ν = 6 and ν = 5, respectively. However, our simulations axis and c-axis directions the starting position does not change. stop at 1 ps and, likely, as on the crystalline ice, the vibrational This is because the crystalline model has alternated and state of the newly formed H2 will lower to ν = 1–2. The opposite electrostatic potentials along the b-axis direction (see difference in the behavior between the Pos1 and Pos2/3 is Figure A2, which constrains the H2 diffusion to the perpend- probably due to a transient chemical bond between one of the icular a-axis, within a channel-like structure created in between reacting H and one O atom of a neighboring water molecule the two alternated potential regions (panels (a) and (c)). In (see Figure A3). 3 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Figure 3. Results of the AIMDs for the crystalline (top panels) and ASW Pos1, Pos2, and Pos3 (panels on the second, third, and bottom rows) models, respectively. Left panels: evolution over time (in fs) of the most relevant energetic components (in kilojoules per mole) of the H2/ice system for the crystalline (panel (a)) and ASW Pos1, Pos2, and Pos3 (panels (c), (e), and (g)) models, respectively. VTOT is the total potential energy (green lines), K H2 and KIce are the kinetic energies of H2 (red lines) and ice (blue lines), respectively, while BE H2 (10 kJ mol−1 = 1200 K) is the binding energy of the H2 (black lines). Right panels: diffusion of the center of mass of the H2 molecule split into the three Cartesian components—c-axis is the direction perpendicular to the ice surface, while the a- and b-axes are along the ice surface (Section 2). 4. Discussion a factor of 10, so that H2 will very likely be injected into the gas phase. Even in the worst case of crystalline ice, where the 4.1. H2 and Other Molecule Desorption simulations show H2 sliding over the perfect ice surface, the While there are no doubts that, once formed on the icy newly formed H2 will sooner or later stumble onto an surfaces, H2 molecules will be released into the gas phase, our imperfection of the ice, and thus its trajectory will be deviated new computations provide a quantitative and atomistic-view and the molecule will escape into the gas. To have an support to this theory. The newly formed H2 molecules possess estimation of the timescale of this phenomenon one could an energy much larger than their binding energy, by more than consider the timescale for the newly H2 to scan the entire 4 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Figure 4. Kinetic energy (in kilojoules per mole) acquired by the water ice molecules as a function of time for the crystalline (top first row panels) and ASW Pos1-3 (lower panels) ices, as marked by the left-hand labels. Left first column panels: energy acquired by a single neighboring water ice molecule (i.e., at ∼3 Å) from the H2 reaction site. Other column panels: energy acquired by the ice surface divided in concentric shells (normalized by the number of water molecules per shell), centered at the reaction site. The labels on the top mark the shells’ radii. surface of an interstellar grain. Assuming its radius is equal to nanoseconds to a max of about 1000 yr, which is still a very 0.1 μm and each icy site size is equal to 3 Å the number of sites short lifetime with respect to the cloud life. to scan is about 105. A first estimate of the time to scan all the We want to highlight that the choice of using a proton order sites can be directly obtained by our simulations by considering ice is just a matter of convenience to test our simulations before that about 15 sites of our crystalline ice are covered in less than going to the more realistic case of the amorphous surface. On a 1 ps; therefore, to scan the 105 it will take less than 10 ns. more realistic proton disordered crystalline ice we expect a However, this is strictly true if the H2 molecule does not loose more anisotropic behavior of the H2 molecule, and, as a energy in other minor impacts, so ∼10 ns can be considered a consequence, a faster H2 desorption. lower limit. The upper limit can be computed by assuming that Our computations also show that the ice absorbs a significant no residual energy is left to the H2 molecule except the thermal fraction of the reaction energy, 45%–65% (Table 1). This one and, in this case, the velocity to scan is given by the energy propagates like a wave from the point where the H + H hopping rate and the timescale for scanning the entire grain reaction occurs (see Figure 4 and also Pantaleone et al. 2020). surface becomes: In ASW ice, a water molecule close to the reaction site acquires 7–43 kJ mol−1 (840–5160 K), for more than 100-200 fs. Within tscan = Ns n - 1 k T Ed the first shell of radius 4 Å, water molecules acquire energies 0 e B dust , (2 ) from 3 to 14 kJ mol−1 (360–1680 K, average value counting all the water molecules within the first shell). Farther away the where ν0 is about 1012 s−1 and Ed is assumed to be 0.3 times energy acquired by the ice-water molecules decreases to less the binding energy (about 400 K). Inserting the numbers, a than 1.6 kJ mol−1. The energy acquired by the water ice rough estimate of the time that H2 takes to leave a typical molecules within a radius of 4 Å from the reaction site could, interstellar grain is „ 1000 yr. In laboratory analogs of potentially, be enough to release into the gas phase any interstellar grains, the timescale would be even larger, and, molecule whose binding energy is lower than 3–14 kJ mol−1. consequently, not observable. In conclusion, the newly formed This could be the case of another H2 molecule frozen on the H2 will leave the surface in a timescale between a few ice, a possibility predicted by some astrochemical models (see, 5 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. e.g., Hincelin et al. 2015). Indeed, since the binding energy of and R ∼ 26–2.6. That is, the H2 formation rate dominates over H2 is less than ∼5 kJ mol−1 (e.g., Vidali 2013; Ferrero et al. the CO freezing one. Therefore, frozen CO molecules can 2020; Molpeceres & Kästner 2020), several H2 molecules potentially be desorbed by the energy released by the H2 could be kicked into the gas phase. formation on the icy grain surfaces. This process, if confirmed, More interesting is the case of CO, the most abundant might naturally explain the presence of gaseous CO in not too molecule after H2 in galactic cold molecular clouds. It has a dense (nH  104 cm−3) molecular clouds and solve a decades- binding energy between 7 and 16 kJ mol−1 (e.g., He et al. long mystery. However, to put this affirmation on solid ground, 2016; Ferrero et al. 2020) so that frozen CO could potentially specific AIMDs showing the effective sublimation of the CO be liberated into the gas-phase by the formation of H2. It has molecule as well as larger ASW ice models and dedicated been long recognized that, in absence of a nonthermal astrochemical modeling simulations that include this effect are desorption process, CO in molecular clouds should be entirely mandatory. They will be the focus of forthcoming works. frozen onto the interstellar ices in ∼3 × 105 yr (see, e.g., Léger 1983). Several mechanisms have been proposed to 4.2. H2 Vibrational State explain its presence in the gas (see, e.g., Leger et al. 1985; Schutte & Greenberg 1991; Duley & Williams 1993; Shen Previous experimental and computational works on graphite et al. 2004; Ivlev et al. 2015). They are all based on the idea surfaces show that after its formation, the H2 molecule that the grain is locally heated and that the frozen molecules are populates vibrational levels around ν = 3–4 (Latimer et al. statistically desorbed on a timescale of a few thousands seconds 2008; Islam et al. 2010; Casolo et al. 2013). On crystalline and (see, e.g., Hasegawa & Herbst 1993; Herbst & Cuppen 2006), amorphous ice surfaces only a few studies were carried out. except when faster photodesorption processes are involved. In Based on classical MD simulations, Takahashi et al. (1999a) particular, Duley & Williams (1993) and Takahashi & and Takahashi & Uehara (2001) predicted that H2 formed on Williams (2000) focused on the CO desorption caused by the amorphous ice would be vibrationally excited to ν = 7–8. In H2 formation. Using the results of classical MD simulations by contrast, in some laboratory experiments (Hornekær et al. Takahashi and colleagues (Takahashi 1999; Takahashi et al. 2003; Roser et al. 2003; Hama et al. 2012), the authors did not 1999b; Takahashi & Williams 2000) and a statistical approach, detect such excited states. Our new simulations confirm both these authors found that the surface within a radius of findings: depending on the site where the formation occurs, H2 ∼4 Å from the reaction site is heated up to „ 40 K for a time can have large (up to 6) or low (1–2) ν (Table 1). For example, too short (∼10 fs) to allow CO to sublimate. However, despite the case of ASW Pos1, which is the one in a cavity, shows the the fact that Takahashi’s simulations are impressive consider- lowest ν. In the ISM, the vibrational state of the newly formed ing that they were carried out more than 20 yr ago, old force H2 molecules will depend on the probability of H2 leaving the fields (like the TIP3P used by Takahashi et al.) have limitations grain surface before being “thermalized” by collisions with ice in describing the actual intermolecular interactions and the water molecules, as also found and discussed in Watanabe et al. energy transfer from one molecule to another, compared to our (2010). simulations treated at the quantum mechanical level. The vibrational state of the newly formed H2 molecule, i.e., Our AIMDs show that the water molecules within a radius of how much the molecule is vibrationally excited, has important consequences on two major astrophysical aspects: 4 Å from the reaction site can be excited from 3 to 14 kJ mol−1 (depending on the position) for more than 100 fs (Figure 4). (i) to help gas-phase reactions with high activation energy Using the range of values in the literature for the CO binding barriers, some of which are the starting points of the energies (7–16 kJ mol−1: He et al. 2016; Ferrero et al. 2020), entire chemistry in molecular clouds (Agúndez et al. and the usual Eyring equation to estimate the half-life time of (2010); and CO desorption, we obtain that a CO molecule within a radius of (ii) the possibility to detect newly formed H2 with near-future 4 Å from Pos1 would desorb in 33–62 fs, and 69–185 fs from James Webb Space Telescope (JWST) observations and, Pos2. On the contrary, CO molecules close to Pos3 and on consequently, measure the rate of H2 formation in crystalline ice would not desorb. Therefore, based on this rough molecular clouds and put constraints to theories. argument, the formation of H2 can potentially desorb frozen CO molecules. Whether this hypothesis is realistic or not depends on (i) how many sites of the AWS ice water molecules 5. Summary are excited as in Pos1 and Pos2, (ii) the ratio of H2 formation We studied the energy dissipation of the H2 formation v (H ) rate with respect to the CO freezing rate (R = v (CO2 ) ), and (iii) reaction on both crystalline and amorphous (ASW) ice models the probability that CO and H2 are adjacent. by means of state-of-the-art ab initio molecular dynamics While a more realistic model for the ASW ice is needed to simulations (AIMDs). In the ASW ice, we explored three estimate the first quantity (i), one can estimate the second one formation sites meant to represent different situations in terms (ii) as follows. At steady state, assuming that the H + H of energetics and surface morphology: one is in a cavity (Pos1) reaction has efficiency 1 (e.g., Vidali 2013), the ratio R of the and the other two are in the outermost parts of the surface (Pos2 H2 formation rate with respect to the CO freezing rate is equal and Pos3). to the ratio between the rate of H atoms over CO molecules In all simulations, we found that about at least 30% of the landing on the grain surface, divided by 2 (because two H reaction energy (∼440 kJ mol−1) is acquired by the newly atoms are needed for the H2 formation). Considering an formed H2, namely more than 10 times the H2 binding energy average Milky Way molecular cloud with H2 density of (∼10 kJ mol−1), so that H2 is likely doomed to leave the ice 103–104 cm−3, temperature 10 K, a gaseous (undepleted) CO and to be injected into the gas. The remaining two-thirds of the abundance with respect to H2 equal to 2 × 10−4, and a cosmic- reaction energy are absorbed by the water ice in less than 1 ps. ray ionization rate of 3 × 10−17 s−1, one obtains nH ∼ 2 cm−3 The water molecules nearby the reaction site have energy peaks 6 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. −1 of 7–43 kJ mol for more than 100–200 fs, while those within Appendix a 4.0 Å radius of 3–14 kJ mol−1. Computational Details We showed that it is in principle possible that frozen CO Proper reactant spin localization is achieved by relying on molecules close to the H2 formation site sublimate. If the unrestricted-DFT broken (spin) symmetry approach. This is confirmed, this will be a simple explanation to the decades- mandatory when dealing with diradical systems in singlet long conundrum of why gaseous CO is present in cold electronic state. In the H + H specific case, despite the total molecular clouds. In order to quantify this effect, new focused spin moment on the global system being 0, the unpaired AIMDs, adopting larger ASW ice models and dedicated electrons on the two H atoms are expected to have spin α astrochemical modeling, will be necessary. (i.e. +1) and β (i.e. −1) each. Finally, the nascent H2 molecules have large (ν = 6) In order to check the quality of our results, in Table A1 a vibrational states in the first picosecond and later decay to comparison between the H2 formation energy in the gas phase 1–2. This high vibrational state could help reactions with an calculated at PBE and CCSD(T) levels is presented. Moreover, activation barrier involving H2 to occur also in cold gas and be the interaction between H2 and one H2O molecule was also observable by JWST. The authors acknowledge funding from the European Table A1 Adsorption and Reaction (H + H → H2) Energy Data Unionʼs Horizon 2020 research and innovation program, the European Research Council (ERC) Project “the Dawn of Reaction Energy (kJ mol−1) Organic Chemistry,” grant agreement No. 741002, the Crystalline −438.9 European Research Council (ERC) Project “Quantum Chem- Pos1 −436.2 istry on Interstellar Grains,” grant agreement No 865657, and Pos2 −435.9 the Marie Skodowska-Curie project “Astro-Chemical Origins” Pos3 −440.3 (ACO), grant agreement No. 811312. A.R. is indebted to the H2 CP2Ka −440.6 H2 Gaussianb −457.4 “Ramón y Cajal” program. MINECO (project CTQ2017- 89132-P) and DIUE (project 2017SGR1323) are acknowl- Binding Energy (kJ mol−1) edged. BSC-MN and OCCIGEN HPCs are kindly acknowl- CP2Ka 4.1 edged for the generous allowance of supercomputing time Gaussianb 2.6 through the QS-2019-2-0028 and 2019-A0060810797 projects, respectively. S.P., N.B., and P.U. acknowledge the Italian Notes. Space Agency for cofunding the Life in Space Project (ASI N. a H2 gas-phase reaction energy calculated at PBE-D3/TZVP level. b 2019-3-U.O) H2 gas-phase reaction energy calculated at CCSD(T)/aug-cc-pv5z level. 7 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Figure A1. Evolution of the H–H bond distance during the AIMD simulations. evaluated at the same levels of theory, giving results in good step of 1 fs) only for the bare ice surface. This ensures an initial agreement. Additionally, the interaction between H2 and one thermally equilibrated ice. Then the velocities of the equilibrated H2O molecule was also evaluated at the same levels of theory, ice surfaces are used as input for the NVE production runs. Those giving results in good agreement. The difference between the velocities correspond to the initial kinetic energy of the ice, as two methodologies is acceptable, considering the energy presented in Figure 3. The H velocities, instead, were manually released by the reaction. set to favor the H–H bond formation respecting the 10 K Regarding the AIMD simulations, we run an equilibration limitation. In Figures A1–A3 the evolution of the H–H distance, AIMD in the NVT ensemble (using the CSVR thermostat, with a the electrostatic potential map of the bare surface, and the time constant of 20 femtoseconds) at 10 K for 1 ps (with a time snapshot of the adduct (H3O+–H−) are reported, respectively. 8 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Figure A2. Top view of the electrostatic potential map of the crystalline (a) and amorphous (b) ice models. Lateral view of the crystalline ice model (c). White arrows represent the H2 diffusion direction. The light blue circle in the bottom panel represents the cavity inside of which the H2 diffuses over the crystalline ice model. Red and blue zones of the electrostatic maps correspond to positive and negative potentials. O atoms in red, H atoms in white. 9 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Figure A3. Evolution of the kinetic energy of H2 (panel (a)) and a neighbor water molecule (panel (b)) for the amorphous Pos1 case. Starting position of the AIMD simulation (panel (c)) and formation of the H−/H3O+ adduct after 0.2 ps of simulation (panel (d)). O atoms are in red, and H atoms are in white. The H atoms highlighted in green are the reactants. ORCID iDs References Stefano Pantaleone https://orcid.org/0000-0002- Agúndez, M., Goicoechea, J., Cernicharo, J., Faure, A., & Roueff, E. 2010, 2457-1065 ApJ, 713, 662 Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, Joan Enrique-Romero https://orcid.org/0000-0002- 53, 541 2147-7735 Casolo, S., Tantardini, G. F., & Martinazzo, R. 2013, PNAS, 110, 6674 Cecilia Ceccarelli https://orcid.org/0000-0001-9664-6292 Duley, W. W., & Williams, D. A. 1993, MNRAS, 260, 37 Stefano Ferrero https://orcid.org/0000-0001-7819-7657 Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11 Goedecker, S., Teter, M., & Hutter, J. 1996, PhRvB, 54, 1703 Nadia Balucani https://orcid.org/0000-0001-5121-5683 Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, JChPh, 132, 154104 Albert Rimola https://orcid.org/0000-0002-9637-4554 Grimme, S., Ehrlich, S., & Goerigk, L. 2011, JCoCh, 32, 1456 Piero Ugliengo https://orcid.org/0000-0001-8886-9832 Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, ApJ, 757, 185 10 The Astrophysical Journal, 917:49 (11pp), 2021 August 10 Pantaleone et al. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 Perdew, J., Burke, K., & Ernzerhof, M. 1996, PhRvL, 77, 3865 He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 825, 89 Roser, J., Swords, S., Vidali, G., Manico, G., & Pirronello, V. 2003, ApJL, Herbst, E., & Cuppen, H. M. 2006, PNAS, 103, 12257 596, L55 Hincelin, U., Chang, Q., & Herbst, E. 2015, A&A, 574, A24 Ruocco, G., Sette, F., Bergmann, U., et al. 1996, Natur, 379, 521 Hollenbach, D., & Salpeter, E. 1970, JChPh, 53, 79 Schutte, W., & Greenberg, J. 1991, A&A, 244, 190 Hollenbach, D., & Salpeter, E. 1971, ApJ, 163, 155 Shen, C., Greenberg, J., Schutte, W., & Van Dishoeck, E. 2004, A&A, Hornekær, L., Baurichter, A., Petrunin, V., Field, D., & Luntz, A. 2003, Sci, 415, 203 302, 1943 Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1988, ApJ, 334, 209 Hutter, J., Iannuzzi, M., Schiffmann, F., & VandeVondele, J. 2014, WIREs Smoluchowski, R. 1983, JPhCh, 87, 4229 Comput. Mol. Sci, 4, 15 Takahashi, J. 1999, EP&S, 51, 1215 Islam, F., Cecchi-Pestellini, C., Viti, S., & Casu, S. 2010, ApJ, 725, 1111 Takahashi, J., Masuda, K., & Nagaoka, M. 1999a, ApJ, 520, 724 Ivlev, A., Röcker, T., Vasyunin, A., & Caselli, P. 2015, ApJ, 805, 59 Takahashi, J., Masuda, K., & Nagaoka, M. 1999b, MNRAS, 306, 22 Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., & Takahashi, J., & Uehara, H. 2001, ApJ, 561, 843 Klein, M. L. 1983, JChPh, 79, 926 Takahashi, J., & Williams, D. A. 2000, MNRAS, 314, 273 Latimer, E. R., Islam, F., & Price, S. D. 2008, CPL, 455, 174 van de Hulst, H. C. 1946, RAOU, 11, 2 Léger, A. 1983, A&A, 123, 271 Vidali, G. 2013, ChRv, 113, 8762 Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147 Wakelam, V., Bron, E., Cazaux, S., et al. 2017, MolAs, 9, 1 Lippert, G., Hutter, J., & Parrinello, M. 1997, MolPh, 92, 477 Watanabe, N., Kimura, Y., Kouchi, A., et al. 2010, ApJL, 714, L233 Molpeceres, G., & Kästner, J. 2020, PCCP, 22, 7552 Whittet, D. C. B. 1993, in Dust and Chemistry in Astronomy, ed. T. J. Millar & Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, ApJ, 897, 56 D. A. Williams (Boca Raton, FL: Routledge), 9 11