H₂ formation on interstellar grains and the fate of reaction energy

STEFANO PANTALEONE,¹, * JOAN ENRIQUE-ROMERO,¹,² CECILIA CECCARELLI,¹ STEFANO FERRERO,² NADIA BALUCANI,¹,³,⁴ ALBERT RIMOLA,² AND PIERO UGLIENGO⁵

¹Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), 38000 Grenoble, France

²Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, 08193, Catalonia, Spain

³Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, Via Elce di Sotto 8, 06123 Perugia, Italy

⁴Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

⁵Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy.

(Received; Revised; Accepted)

Submitted to ApJ

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 H₂ formation on the dust icy mantles, how it is partitioned between the substrate and the newly formed H₂, a process that has a profound impact on the interstellar medium. We carried out state-of-art *ab-initio* molecular dynamics simulations of H₂ 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 ~3.1 eV) is absorbed by the ice in less than 1 ps. The remaining energy (~140 kJ/mol ~1.5 eV) is kept by the newly born H₂. Since it is ten times larger than the H₂ binding energy on the ice, the new H₂ 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 (360–1560 K), to potentially liberate other frozen H₂ 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 H₂ drops from highly excited states (ν = 6) to low (ν ≤ 2) vibrational levels in a timescale of the order of ps.

**Keywords:** H₂ formation, ice mantles, chemical desorption, energy dissipation, molecular dynamics, DFT

1. INTRODUCTION

Molecular hydrogen is the most abundant molecule in the Universe. Its formation is also the first step of the interstellar chemistry and, therefore, a fundamental reaction. In molecular clouds, H₂ is mainly formed

via the H + H association reaction on the interstellar dust grain surfaces, which act as a third body capable to partially absorb the energy (~440 kJ/mol ~4.56 eV) released by the chemical reaction and, consequently, stabilise the newly formed H₂ molecule (e.g. van de Hulst 1946; Hollenbach & Salpeter 1970, 1971; Smoluchowski 1983). Although no doubts exist on the occurrence of this process, many specific details remain poorly known (see e.g. Vidali 2013; Wakelam et al. 2017).

Here, we focus on the fate of the energy released by the reaction, which has been source of debate for decades. Specifically: how much of the reaction energy is absorbed by the dust grain and in what timescale? What fraction of the chemical energy does go into the kinetic

Corresponding authors: Stefano Pantaleone
stefano.pantaleone@univ-grenoble-alpes.fr

<footnote>* present address: Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy.
and
Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, Via Elce di Sotto 8, 06123 Perugia, Italy</footnote>
energy of the newly formed H₂? Is this energy large enough as H₂ breaks the interaction with the surface and leaves into the gas-phase? How much energy does the H₂ molecule possess when it leaves the grain surface and in what ro-vibrational state? And, finally, is the energy transmitted to the dust grain enough for locally warming it up and make adjacent frozen molecules sublimate? These points have a profound impact on several aspects of astrochemistry and, more generally, the physics and chemistry of the interstellar medium (ISM) (Duley & Williams 1993).

Answers to the above questions largely depend on the nature of the substrate, namely the specific dust grain surface or, in practice, the interstellar environment where H₂ forms. Here we focus on cold (~10 K) molecular clouds. In these environments, the grain refractory cores are coated by water-dominated icy mantles either in polycrystalline (PCI) or amorphous (ASW for amorphous solid water) phases, the latter dominating over the former (Smith et al. 1988; Whittet 1993; Boogert et al. 2015).

Unfortunately, laboratory experiments cannot entirely reproduce the interstellar conditions so that, while they can suggest routes and processes, they cannot provide definitive answers to the above-exposed questions (see e.g. Vidali 2013). For example, obtaining experimentally how the H₂ formation nascent energy is partitioned is very difficult (Roser et al. 2003; Hornekær et al. 2003; Hama et al. 2012; Watanabe et al. 2010). Computational chemistry methods can be regarded as a complementary tool and, sometimes, as a unique alternative to laboratory experiments. So far, a limited number of studies devoted to the energy dissipation of the H₂ formation reaction on icy grains has appeared in the literature (Takahashi et al. 1999a; Takahashi & Uehara 2001; Herbst & Cuppen 2006).

Here, we present a new theoretical study on the dissipation of the energy released by the H + H → H₂ reaction on water ice using ab-initio molecular dynamics simulations (AIMDs). We simulated the reaction adopting a Langmuir-Hinshelwood (LH) mechanism (i.e. both reactants are adsorbed one the surface in neighboring sites) both on amorphous and crystalline periodic models of interstellar icy mantles.

## 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 (Lip-

pert et al. 1997). The DFT (Density Functional Theory) PBE (Perdew-Burke-Ernzerhof) 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 sec. Ice models), and during the geometry optimization only the reacting H atoms (those forming H₂) were allowed to move, keeping the water molecules fixed at their thermalized positions. All calculations were carried out within the unrestricted formalism as we deal with open-shell systems (see SI for more details). The binding energy (BE) of H₂ was calculated as:

$$BE_{H_2} = E_{CPLX} - (E_{Ice} + E_{H_2}) \quad (1)$$

where E$_{CPLX}$ is the energy of the H₂/Ice system, E$_{Ice}$ that of the bare ice surface, and E$_{H_2}$ the energy of the H₂ alone.

AIMDs were carried out within the NVE (N = number of particles, V = volume, E = energy) ensemble, where the total energy V$_{TOT}$ (i.e. potential + kinetic) is conserved. The evolution of the system was followed for 5 ps for the crystalline model and for 1 ps for the amorphous model, using time-steps of 0.2 fs.

### 2.2. Ice models

For the crystalline ice model, a periodic slab cut from the hexagonal ice bulk structure was used. The periodic cell parameters defining the system are: a = 26.318 Å and b = 28.330 Å, while the slab thickness is ~ 21 Å (corresponding to seven layers), for a total of 576 water molecules. The c parameter of the simulation box, i.e. the non-periodic one, was set to 50 Å to avoid interactions among the fictitious slab replicas. The size of the slab model was chosen to avoid non-physical temperature increase due to the extremely large exothermic reaction (Pantaleone et al. 2020).

The amorphous model was obtained by performing a classical molecular dynamics (MD) simulation on the crystalline structure. The MD was carried out with the TIP3P force-field (Jorgensen et al. 1983) for 200 ps (with a timestep of 0.5 fs) at 300 K within the NVT (N = number of particles, V = volume, T = Temperature) ensemble. A second NVT simulation was performed at 10 K, to cool down the system at the temperature of cold molecular clouds. Finally, a geometry optimization and another NVT-MD simulation at 10 K were run using PBE, to recover the potential energy surface at the same theory level as described in the previous section.

### 2.3. H₂ vibrational state
<figure><img src="image_1.png" /><figcaption>(a) Crystalline ice, top view – ab plane.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>(b) ASW ice, top view – ab plane.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) ASW ice, side view – ac plane.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(d) ASW ice, side view – bc plane.</figcaption></figure>

**Figure 1.** Top view of the crystalline (a), and top (b) and side (c and d) views of the ASW ice models. The yellow circles represent the centers 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 initial H to H distance is about 4 Å (see Supplementary Fig. 1). O atoms are in red and H atoms of the ice are in white.

To calculate the H₂ vibrational state during the AIMDs the anharmonic oscillator model was employed. As a first step, the PES (potential energy surface) of the H₂ isolated molecule was explored by performing a rigid scan of the H–H distance, ranging from 0.3 to 3.0 Å with a step of 0.01 Å. Computed data were fitted with the Morse equation to obtain the force constant of the oscillator, its dissociation energy and, hence, the vibrational levels of H₂. The H₂ turning points calculated with our model were compared with those of the AIMD simulations and the vibrational level was assigned at each H₂ oscillation during the AIMDs.

### 3. RESULTS

#### 3.1. Reactants plus product positions in the ice models

Figure 1 shows the structure of the used crystalline and ASW ice models as well as the starting positions of the center of mass of the two adsorbed H atoms before the reaction. The starting H–H distance, as well as its evolution during the AIMDs is shown in Supplementary Figure 1. The ASW model presents a rugged morphol-

<figure><img src="image_5.png" /><figcaption>(a) R - Crystalline (0.0)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) P - Crystalline (-439.4)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) R - ASW (0.0)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) P - ASW (-436.2)</figcaption></figure>

**Figure 2.** PBE-D3/TZVP optimized geometries of reactants (left panels) and product (right panels) of the H₂ reaction formation on the crystalline (top) and ASW Pos1 (bottom) ice models. The numbers in parenthesis correspond to the relative energy in kJ/mol with respect to the reactants. The H atoms involved in the H + H → H₂ reaction are in yellow, those belonging to surface water molecules are in white, and O atoms are in red.

ogy (e.g. holes, cavities, channels) in comparison to the crystalline one and, therefore, we carried out simulations in three positions, roughly representing different possible situations in terms of energetics and surface morphology. Position marked as Pos1 is the deepest one, within a cavity, while positions Pos2 and Pos3 are in the outermost parts of the surface.

As a first step, we optimized the geometries of reactants (the two H atoms) and product (the H₂ molecule) on the water ice surface in order to obtain the potential energy surface of the reaction. No search for transition state structures has been performed because AIMDs indicated this reaction is barrierless at 10 K. That is, the starting kinetic energy of the two H atoms provided by the 10 K is high enough to overcome the eventually small potential energy barrier. The total energy to be dissipated is around 435-440 kJ/mol, depending on the surface and starting position.

Figure 2 provides a view of the two H atoms starting positions and H₂ position immediately after the reaction, on the crystalline and ASW (Pos1 as an example) models, respectively.
<table><caption>**Table 1.** Results of the AIMD for the H₂ formation on the crystalline (first row) and the ASW three positions (Fig. 1) (bottom three rows) ices. Second and third columns report the fraction of kinetic energy of the ice and the newly formed H₂, respectively, at the end of the simulation (5 and 1 ps for crystalline and ASW ices, respectively). Fourth column reports the H₂ binding energy (in kJ/mol) 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 kJ/mol) of an ice-molecule neighbor (i.e. at ~3 Å) to where the reaction occurs.</caption><thead><tr><th rowspan="2">Model</th><th>K<sub>ice</sub></th><th>K<sub>H2</sub></th><th rowspan="2">BE<sub>H2</sub></th><th rowspan="2">ν<sup>0.2–1</sup></th><th rowspan="2">ν<sup>1–5</sup></th><th rowspan="2">E<sub>peak</sub></th></tr><tr><th>K<sub>ice</sub>+K<sub>H2</sub></th><th>K<sub>ice</sub>+K<sub>H2</sub></th></tr><tr><th></th><th></th><th></th><th colspan="3">[kJ/mol]</th><th colspan="2">[kJ/mol]</th></tr></thead><tbody><tr><td>Crys.</td><td>0.45</td><td>0.55</td><td>9.0</td><td>4–5</td><td>1–2</td><td>12.0</td><td></td></tr><tr><td>Pos1</td><td>0.65</td><td>0.35</td><td>9.7</td><td>1–2</td><td>–</td><td>43.0</td><td></td></tr><tr><td>Pos2</td><td>0.46</td><td>0.54</td><td>10.6</td><td>6</td><td>–</td><td>7.0</td><td></td></tr><tr><td>Pos3</td><td>0.47</td><td>0.53</td><td>12.4</td><td>5</td><td>–</td><td>7.0</td><td></td></tr></tbody></table>

## 3.2. Molecular hydrogen desorption

AIMDs results are summarised in Table 1 and shown in Fig. 3. First, at least half of the kinetic energy is absorbed by the ice while the remaining one is kept by the newly formed H₂ molecule (~50–35%). These numbers are similar for both for the crystalline and ASW ice positions, suggesting they do not depend much on the surface structural details of the ice.

Second, after the H–H bond formation, the newly formed H₂ molecule diffuses over the surface, as it keeps a large kinetic energy, in both crystalline and ASW ice models. However, despite the energetics of the two processes are similar (Figs. 3 left panels) the H₂ diffusion (Figs. 3 right panels) is different.

On the crystalline surface, the diffusion of the newly formed H₂ over the ice surface is constrained in a specific direction, along the a-axis, parallel to the ice surface, whereas in the b-axis and c-axis directions the starting position does not change. This is because the crystalline model has alternated and opposite electrostatic potentials along the b-axis direction (see Supplementary Fig. 2, which constrains the H₂ diffusion to the perpendicular a-axis, within a channel-like structure created in between the two alternated potential regions (see Figs. S2a and S2c). In contrast, on the ASW ice model, H₂ diffuses over all the 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, H₂ achieves a maximum height of ~10 Å over the surface in a timescale of 0.4 ps and does not come back (Fig. 3-d), this way leaving definitely the surface. In Pos2, on the contrary, H₂ slides and hops on the surface within the 1 ps of the simulation (Fig. 3-f). A similar behavior is observed for Pos3, but with wider hops and jumps (Fig. 3-h). However, and this is the key point of the simulations, the kinetic energy remains much larger (~150 kJ/mol) than the binding energy (~10 kJ/mol) in all cases: therefore, the nascent H₂ is likely doomed to leave the ice surface.

## 3.3. Energy dissipation towards the ice water molecules

Fig. 4 shows how the energy absorbed by the ice is first transmitted from reaction site to a neighbor water (i.e. at ~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 and stays with that energy for more than about 100–200 fs before the energy is dissipated towards more external water molecules. These timescales are in excellent agreement with previous studies of the sound speed in ices (see e.g. Ruocco et al. 1996).

## 3.4. H₂ vibrational state

The initial state of the H₂ formed on the crystalline ice is $ν = 4–5$ (obtained averaging over 0.2–1 ps) and then it decreases to $ν = 1–2$ (1–5 ps average) because of the H₂ energy dissipation. Conversely, at Pos1 of the ASW model, H₂ is formed with a low vibrational state ($ν = 1–2, 0.2–1$ ps average).

For the other two positions (Pos2 and Pos3), the vibrational state of the newly formed H₂ molecule lies in highly excited states, $ν = 6$ and $ν = 5$, respectively. However, our simulations stop at 1 ps and, likely, as on the crystalline ice, the vibrational state of the newly formed H₂ will lower to $ν = 1–2$. The difference in the behavior between the Pos1 and Pos2/3 is probably due to a transient chemical bond between one of the reacting H and one O atom of a neighboring water molecule (see Supplementary Figure 3).

# 4. DISCUSSION

## 4.1. H₂ and other molecules desorption

While there were no doubts that, once formed on the icy surfaces, H₂ molecules will be released into the gas-phase, our new computations provide a quantitative and atomistic-view support to this theory. The newly formed H₂ molecules possess an energy much larger than their binding energy, by more than a factor ten, so that very likely H₂ will be injected into the gas-phase. Even in
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Results of the AIMDs for the crystalline (top panels) and ASW Pos1, Pos2, 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 kJ/mol) of the H₂/ice system for the crystalline (panel a) and ASW Pos1, Pos2 and Pos3 (panel c, e and g) models, respectively. V_TOT is the total potential energy (green lines), K_H₂ and K_Ice are the kinetic energies of H₂ (red lines) and ice (blue lines), respectively, while BE_H₂ (10 kJ/mol = 1200 K) is the binding energy of the H₂ (black lines). Right panels: Diffusion of the center of mass of the H₂ molecule split into the three cartesian components: c-axis is the direction perpendicular to the ice surface, while a-axis and b-axis are along the ice surface (§ 2).</figcaption></figure>

the worst case of crystalline ice, where the simulations show H₂ sliding over the perfect ice surface, the newly formed H₂ will sooner or later stumble on an imperfection of the ice so that its trajectory would be deviated and the molecule escape into the gas. To have an estimation of the timescale of this phenomenon one could consider the timescale for the newly H₂ to scan the entire surface of an interstellar grain. Assuming its radius

equal to 0.1 µm and each icy site size equal to 3 Å, the number of sites to scan are about $10^5$. A first estimate of the time to scan all the sites can be directly obtained by our simulations by considering that about 15 sites of our crystalline ice are covered in less 1 ps: therefore, to scan the $10^5$ it will take less than 10 ns. However, this is strictly true if the H₂ molecule does not lose energy in other minor impacts, so ~ 10 ns can be considered
a lower limit. The upper limit can be computed by assuming that no residual energy is left to the H₂ molecule except the thermal one and, in this case, the velocity to scan is given by the hopping rate and the timescale for scanning the entire grain surface becomes:

$$t_{\text{scan}} = N_s \nu_0^{-1} e^{\frac{E_d}{k_B T_{\text{dust}}}} \quad (2)$$

where $\nu_0$ is about $10^{12} \text{ s}^{-1}$ and $E_d$ is assumed to be 0.3 times the binding energy (about 400 K). Inserting the numbers, a rough estimate of the time that H₂ takes to leave a typical interstellar grain is $\le 1000$ yr. In laboratory analogues of interstellar grains, the timescale would be even larger, and, consequently, not observable. In conclusion, the newly formed H₂ will leave the surface in a timescale between a few ns to a max of about 1000 yr, which is still a very short lifetime with respect to the cloud life.

We want to highlight that the choice of using a proton order ice is just a matter of convenience to test our simulations, before going to the more realistic case of the amorphous surface. On a more realistic proton disordered crystalline ice we expect a more anisotropic behavior of the H₂ molecule, and, as a consequence, a faster H₂ desorption.

Our computations also show that the ice absorbs a significant fraction of the reaction energy, 45–65% (Tab. 1). This energy propagates like a wave from the point where the H + H reaction occurs (see Fig. 4 and also Pantaleone et al. 2020). In ASW ice, a water molecule close to reaction site acquires 7–43 kJ/mol, (840–5160 K), for more than 100-200 fs. Within the first shell of radius 4 Å, water molecules acquire energies from 3 to 14 kJ/mol (360 to 1680 K, average value counting all the water molecules within the first shell). and farther away the energy acquired by the ice-water molecules decreases to less than 1.6 kJ/mol. The energy acquired by the water ice molecules within a radius of 4 Å from the reaction site could, potentially, be enough to release into the gas-phase any molecule whose binding energy is lower than 3–14 kJ/mol. This could be the case of another H₂ molecule frozen on the ice, a possibility predicted by some astrochemical models (see e.g. Hincelin et al. 2015). Indeed, since the binding energy of H₂ is less than ~5 kJ/mol (e.g. Vidali 2013; Ferrero et al. 2020; Molpeceres & Kästner 2020), several H₂ molecules could be kicked into the gas-phase.

More interesting is the case of CO, the most abundant molecule after H₂ in galactic cold molecular clouds. It has a binding energy between 7 and 16 kJ/mol (e.g. He et al. 2016; Ferrero et al. 2020) so that frozen CO could potentially be liberated into the gas-phase by the formation of H₂. It has been long recognized that, in absence

of a non-thermal desorption process, CO in molecular clouds should be entirely frozen onto the interstellar ices in ~3×10⁵ yr (see e.g. Léger 1983). Several mechanisms have been proposed to explain its presence in the gas (see e.g. Leger et al. 1985; Schutte & Greenberg 1991; Duley & Williams 1993; Shen et al. 2004; Ivlev et al. 2015). They are all based on the idea that the grain is locally heated and that the frozen molecules are statistically desorbed on a time scale of a few thousands seconds (see e.g. Hasegawa & Herbst 1993; Herbst & Cuppen 2006), except when faster photodesorption processes are involved. In particular, Duley & Williams (Duley & Williams 1993) and Takahashi & Williams (Takahashi & Williams 2000) focused on the CO desorption caused by the H₂ formation. Using the results of classical MD simulations by Takahashi and colleagues (Takahashi et al. 1999b; Takahashi 1999; Takahashi & Williams 2000) and a statistical approach, these authors found that the surface within a radius of ~4 Å from the reaction site is heated up to $\le 40$ K for a too short time (~10 fs) to allow CO to sublimate. However, despite the Takahashi's simulations are impressive considering that they were carried out more than 20 years ago, old force fields (like the TIP3P used by Takahashi et al.) have limitations in describing the actual intermolecular interactions and the energy transfer from one molecule to another, compared to our simulations treated at quantum mechanical level.

Our AIMDs show that the water molecules within a radius of 4 Å from the reaction site can be excited from 3 to 14 kJ/mol (depending on the position) for more than 100 fs (Fig. 4). Using the range of values in the literature for the CO binding energies (7–16 kJ/mol: He et al. 2016; Ferrero et al. 2020), and the usual Eyring equation to estimate the half-life time of CO desorption, we obtain that a CO molecule within a radius of 4 Å from Pos1 would desorb in 33–62 fs, and 69–185 fs from Pos2. On the contrary, CO molecules close to Pos3 and on crystalline ice would not desorb. Therefore, based on this rough argument, the formation of H₂ can potentially desorb frozen CO molecules. Whether this hypothesis is realistic or not, it depends on (i) how many sites of the AWS ice water molecules are excited as in Pos1 and Pos2, (ii) the ratio of H₂ formation rate with respect to the CO freezing rate ($R = \frac{v(H_2)}{v(CO)}$), and (iii) the probability that CO and H₂ are adjacent.

While a more realistic model for the ASW ice is needed to estimate the first quantity (i), one can estimate the second one (ii) as follows. At steady state, assuming that the H + H reaction has efficiency 1 (e.g. Vidali 2013), the ratio R of the H₂ formation rate with respect to the CO freezing rate is equal to the ratio between
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Kinetic energy (in kJ/mol) 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 H₂ 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.</figcaption></figure>

the rate of H atoms over CO molecules landing on the grain surface, divided by 2 (because two H atoms are needed for the H₂ formation). Considering an average Milky Way molecular cloud with H₂ density of $10^3–10^4$ cm⁻³, temperature 10 K, a gaseous (undepleted) CO abundance with respect to H₂ equal to $2 × 10^{-4}$, and a cosmic ray ionisation rate of $3 × 10^{-17}$ s⁻¹, one obtains $n_H ∼ 2$ cm⁻³ and $R ∼ 26–2.6$. That is, the H₂ formation rate dominates over the CO freezing one. Therefore, frozen CO molecules can potentially be desorbed by the energy released by the H₂ formation on the icy grain surfaces. This process, if confirmed, might naturally explain the presence of gaseous CO in not too dense ($n_H ≲ 10^4$ cm⁻³) molecular clouds and solve a decades long mystery. However, to put this affirmation on a solid ground, specific AIMDs showing the effective sublimation of the CO molecule as well as larger ASW ice models and dedicated astrochemical modeling simulations that include this effect are mandatory. They will be the focus of forthcoming works.

## 4.2. H₂ vibrational state

Previous experimental and computational works on graphite surfaces show that after its formation, the H₂ molecule populates vibrational levels around $ν = 3–4$ (Latimer et al. (2008); Islam et al. (2010); Casolo et al. (2013)). On crystalline and amorphous ice surfaces only a few studies were carried out. Based on classical MD simulations, Takahashi et al. (1999a); Takahashi & Uehara (2001) predicted that H₂ formed on amorphous ice would be vibrationally excited to $ν = 7–8$. In contrast, in some laboratory experiments Roser et al. (2003); Hornekær et al. (2003); Hama et al. (2012), the authors did not detect such excited states. Our new simulations confirm both findings: depending on the site where the formation occurs, H₂ can have large (up to 6) or low (1–2) $ν$ (Tab. 1). For example, the case of ASW Pos1, which is the one in a cavity, shows the lowest $ν$. In the ISM, the vibrational state of the newly formed H₂ molecules will depend on the probability of H₂ leaving the grain surface before being "thermalised"
by collisions with ice water molecules, as also found and discussed in Watanabe et al. (2010).

The vibrational state of the newly formed H₂ molecule, i.e. how much the molecule is vibrationally excited, has important consequences on two major astrophysical aspects: (i) to help gas-phase reactions with high activation energy barriers, some of which are the starting points of the entire chemistry in molecular clouds Agúndez et al. (2010); (ii) the possibility to detect newly formed H₂ with near-future JWST (James Webb Space Telescope) observations and, consequently, measure the rate of H₂ formation in molecular clouds and put constraints to theories.

## SUMMARY

We studied the energy dissipation of the H₂ formation reaction on both crystalline and amorphous (ASW) ice models by means of state-of-art ab-initio molecular dynamics simulations (AIMDs). In the ASW ice, we explored three formation sites meant to represent different situations in terms of energetics and surface morphology: one is in a cavity (Pos1) and two in the outermost parts of the surface (Pos2 and Pos3).

In all simulations, we found that about at least 30% of the reaction energy (~440 kJ/mol) is acquired by the newly formed H₂, namely more than ten times the H₂ binding energy (~10 kJ/mol), so that H₂ is likely doomed to leave the ice and to be injected into the gas. The remaining two thirds of the reaction energy are absorbed by the water ice in less than 1 ps. The water molecules nearby to the reaction site have energy peaks

of 7–43 kJ/mol for more than 100–200 fs, while those within a 4.0 Å radius of 3–14 kJ/mol.

We showed that it is in principle possible that frozen CO molecules close to the H₂ formation site sublimate. If confirmed, this will be a simple explanation to the decades-long conundrum of why gaseous CO is present in cold molecular clouds. In order to quantify this effect new focused AIMDs adopting larger ASW ice models and dedicated astrochemical modeling will be necessary.

Finally, the nascent H₂ molecules have large ($\nu$ = 6) vibrational states in the first ps and later decay to 1–2. This high vibrational state could help reactions with an activation barrier involving H₂ to occur also in cold gas and be observable by JWST.

## ACKNOWLEDGMENTS

The authors acknowledge funding from European Union's Horizon 2020 research and innovation program, the European Research Council (ERC) Project "the Dawn of Organic Chemistry", grant agreement No 741002, the European Research Council (ERC) Project "Quantum Chemistry on Interstellar Grains", grant agreement No 865657, and the Marie Skłodowska-Curie project "Astro-Chemical Origins" (ACO), grant agreement No 811312. AR is indebted to the "Ramón y Cajal" program. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged. BSC-MN and OCCIGEN HPCs are kindly acknowledged for the generous allowance of supercomputing time through the QS-2019-2-0028 and 2019-A0060810797 projects, respectively. SP, NB and PU acknowledge the Italian Space Agency for co-funding the Life in Space Project (ASI N. 2019-3-U.O

## REFERENCES

Agúndez, M., Goicoechea, J., Cernicharo, J., Faure, A., & Roueff, E. 2010, The Astrophysical Journal, 713, 662

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, Annual Review of Astronomy and Astrophysics, 53, 541, doi: 10.1146/annurev-astro-082214-122348

Casolo, S., Tantardini, G. F., & Martinazzo, R. 2013, Proceedings of the National Academy of Sciences, 110, 6674

Duley, W. W., & Williams, D. A. 1993, Mon. Not. R. Astron. Soc., 260, 37, doi: 10.1093/mnras/260.1.37

Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, The Astrophysical Journal, 904, 11

Goedecker, S., Teter, M., & Hutter, J. 1996, Phys. Rev. B, 54 (3), 1703, doi: 10.1103/PhysRevB.54.1703

Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, J. Chem. Phys., 132 (15), 154104, doi: 10.1063/1.3382344

Grimme, S., Ehrlich, S., & Goerigk, L. 2011, J. Comput. Chem., 32, 1456, doi: 10.1002/jcc.21759

Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, The Astrophysical Journal, 757, 185

Hasegawa, T. I., & Herbst, E. 1993, Mon. Not. R. Astron. Soc., 261, 83, doi: 10.1093/mnras/261.1.83

He, J., Acharyya, K., & Vidali, G. 2016, The Astrophysical Journal, 825, 89, doi: 10.3847/0004-637X/825/2/89

Herbst, E., & Cuppen, H. M. 2006, Proceedings of the National Academy of Sciences, 103, 12257

Hincelin, U., Chang, Q., & Herbst, E. 2015, Astronomy & Astrophysics, 574, A24
Hollenbach, D., & Salpeter, E. 1970, The Journal of Chemical Physics, 53, 79

—. 1971, The Astrophysical Journal, 163, 155

Hornekær, L., Baurichter, A., Petrunin, V., Field, D., & Luntz, A. 2003, Science, 302, 1943

Hutter, J., Iannuzzi, M., Schiffmann, F., & VandeVondele, J. 2014, WIREs Comput. Mol. Sci., 4, 15, doi: 10.1002/wcms.1159

Islam, F., Cecchi-Pestellini, C., Viti, S., & Casu, S. 2010, The Astrophysical Journal, 725, 1111

Ivlev, A., Röcker, T., Vasyunin, A., & Caselli, P. 2015, The Astrophysical Journal, 805, 59

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., & Klein, M. L. 1983, J. Chem. Phys., 79, 926

Latimer, E. R., Islam, F., & Price, S. D. 2008, Chemical Physics Letters, 455, 174

Léger, A. 1983, Astronomy and Astrophysics, 123, 271

Leger, A., Jura, M., & Omont, A. 1985, Astronomy and Astrophysics, 144, 147

Lippert, G., Hutter, J., & Parrinello, M. 1997, Mol. Phys., 92 (3), 477, doi: 10.1080/002689797170220

Molpeceres, G., & Kästner, J. 2020, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 22, 7552, doi: 10.1039/D0CP00250J

Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, The Astrophysical Journal, 897, 56

Perdew, J., Burke, K., & Ernzerhof, M. 1996, Phys. Rev. Lett., 77 (18), 3865, doi: 10.1103/PhysRevLett.77.3865

Roser, J., Swords, S., Vidali, G., Manico, G., & Pirronello, V. 2003, The Astrophysical Journal Letters, 596, L55

Ruocco, G., Sette, F., Bergmann, U., et al. 1996, Nature, 379, 521, doi: 10.1038/379521a0

Schutte, W., & Greenberg, J. 1991, Astronomy and Astrophysics, 244, 190

Shen, C., Greenberg, J., Schutte, W., & Van Dishoeck, E. 2004, Astronomy & Astrophysics, 415, 203

Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1988, The Astrophysical Journal, 334, 209

Smoluchowski, R. 1983, J. Phys. Chem., 87, 4229

Takahashi, J. 1999, Earth, planets and space, 51, 1215

Takahashi, J., Masuda, K., & Nagaoka, M. 1999a, The Astrophysical Journal, 520, 724

—. 1999b, Monthly Notices of the Royal Astronomical Society, 306, 22

Takahashi, J., & Uehara, H. 2001, The Astrophysical Journal, 561, 843

Takahashi, J., & Williams, D. A. 2000, Monthly Notices of the Royal Astronomical Society, 314, 273

van de Hulst, H. C. 1946, Recherches Astronomiques de l'observatoire d'Utrecht, 11, 2

Vidali, G. 2013, Chemical Reviews, 113, 8762, doi: 10.1021/cr400156b

Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1

Watanabe, N., Kimura, Y., Kouchi, A., et al. 2010, The Astrophysical Journal Letters, 714, L233

Whittet, D. C. B. 1993, Observations of Molecular Ices, ed. T. J. Millar & D. A. Williams, 9
# APPENDIX

## A. COMPUTATIONAL DETAILS

Proper reactants spin localisation is achieved by relying on the unrestricted-DFT broken (spin) symmetry approach. This is mandatory when dealing with diradical systems in singlet electronic state. In the H + H specific case, despite the total spin moment on the global system is 0, the unpaired electrons on the two H atoms are expected to have spin $\alpha$ (i.e. +1) and $\beta$ (i.e. -1) each.

In order to check the quality of our results, in Table 1 a comparison between the $H_2$ formation energy in the gas phase calculated at PBE and CCSD(T) levels is presented. Moreover, also the interaction between $H_2$ and one $H_2O$ molecule was evaluated at the same levels of theory, giving results in good agreement. Moreover, also the interaction between $H_2$ and one $H_2O$ molecule was evaluated at the same levels of theory, giving results in good agreement. The difference between the two methodologies is acceptable, considering the energy released by the reaction.

As regards the AIMD simulations, we run an equilibration AIMD in the NVT ensemble (using the CSVR thermostat, with a time constant of 20 femtoseconds) at 10 K for 1 ps (with a time step of 1 fs) only for the bare ice surface. This ensures an initial thermally equilibrated ice. Then the velocities of the equilibrated ice surfaces are used as input for the NVE production runs. Those velocities, correspond to the initial kinetic energy of the ice, as presented in Figure 3. The H velocities, instead, were manually set to favor the H-H bond formation respecting the 10 K limitation.

<table><caption>**Supplementary Table 1.** Adsorption and reaction (H + H → $H_2$) energy data.</caption><thead><tr><th></th><th>Reaction energy (kJ mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Crystalline</td><td>-438.9</td></tr><tr><td>Pos1</td><td>-436.2</td></tr><tr><td>Pos2</td><td>-435.9</td></tr><tr><td>Pos3</td><td>-440.3</td></tr><tr><td>H<sub>2</sub> CP2K<sup>a</sup></td><td>-440.6</td></tr><tr><td>H<sub>2</sub> Gaussian<sup>b</sup></td><td>-457.4</td></tr></tbody><thead><tr><th></th><th>Binding energy (kJ mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>CP2K<sup>a</sup></td><td>4.1</td></tr><tr><td>Gaussian<sup>b</sup></td><td>2.6</td></tr></tbody></table>

<footnote><sup>a</sup> $H_2$ gas phase reaction energy calculated at PBE-D3/TZVP level</footnote>

<footnote><sup>b</sup> $H_2$ gas phase reaction energy calculated at CCSD(T)/aug-cc-pv5z level</footnote>
<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 1.** Evolution of the H–H bond distance during the AIMD simulations.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 2.** 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 H₂ diffusion direction. The light blue circle in the bottom panel represents the cavity inside of which the H₂ 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.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 3.** Evolution of the kinetic energy of H₂ (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⁻/H₃O⁺ adduct after 0.2 picoseconds of simulation (panel (d)). O atoms in red, H atoms in white. The H atoms highlighted in green are the reactants.</figcaption></figure>