H₂ Formation on Interstellar Grains and the Fate of Reaction Energy

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 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 the Appendix for more details). The binding energy (BE) of H₂ was calculated as:

$$\begin{eqnarray}&&{\mathrm{BE}}_{{{\rm{H}}}_{2}}={E}_{\mathrm{CPLX}}-({E}_{\mathrm{Ice}}+{E}_{{{\rm{H}}}_{2}}),\end{eqnarray} \tag{ 1 }$$

where E_CPLX is the energy of the H₂/ice system, E_Ice is the energy of the bare ice surface, and ${E}_{{{\rm{H}}}_{2}}$ is 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.

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 nonperiodic one, was set to 50 Å to avoid interactions among the fictitious slab replicas. The size of the slab model was chosen to avoid nonphysical 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 time step 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.

To calculate the H₂ vibrational state during the AIMDs the anharmonic oscillator model was employed. As a first step, the potential energy surface (PES) 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.

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 Figure A1. The ASW model presents a rugged morphology (e.g., holes, cavities, and 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. The position marked as Pos1 is the deepest one, within a cavity, while positions Pos2 and Pos3 are in the outermost parts of the surface.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

AIMDs results are summarized in Table 1 and shown in Figure 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 no depend much on the surface structural details of the ice.

Zoom In Zoom Out Reset image size

Table 1. Results of the AIMD for the H₂ Formation on the Crystalline (First Row) and the ASW Three Position (Figure 1; Bottom Three Rows) Ices

Model	$\tfrac{{K}_{\mathrm{ice}}}{{K}_{\mathrm{ice}}+{K}_{{{\rm{H}}}_{2}}}$	$\tfrac{{K}_{{{\rm{H}}}_{2}}}{{K}_{\mathrm{ice}}+{K}_{{{\rm{H}}}_{2}}}$	BE${}_{{{\rm{H}}}_{2}}$	ν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	⋯	43.0
Pos2	0.46	0.54	10.6	6	⋯	7.0
Pos3	0.47	0.53	12.4	5	⋯	7.0

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

Download table as:  ASCIITypeset image

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 being similar (Figure 3, left panels) the H₂ diffusion (Figure 3, right panels) is different.

Zoom In Zoom Out Reset image size

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 Figure A2, which constrains the H₂ diffusion to the perpendicular a-axis, within a channel-like structure created in between the two alternated potential regions (panels (a) and (c)). In contrast, on the ASW ice model, H₂ 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, H₂ achieves a maximum height of ∼10 Å over the surface in a timescale of 0.4 ps and does not come back (Figure 3(d)), thus definitely leaving the surface. On the contrary, in Pos2 H₂ slides and hops on the surface within the 1 ps of the simulation (Figure 3(f)). A similar behavior is observed for Pos3, but with wider hops and jumps (Figure 3(h)). 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.

Figure 4 shows how the energy absorbed by the ice is first transmitted from the 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 toward 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).

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 Figure A3).

While there are 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 of 10, so that H₂ will very likely be injected into the gas phase. Even in 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 onto an imperfection of the ice, and thus its trajectory will be deviated and the molecule will 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 is equal to 0.1 μm and each icy site size is equal to 3 Å the number of sites to scan is about 10⁵. 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 than 1 ps; therefore, to scan the 10⁵ it will take less than 10 ns. However, this is strictly true if the H₂ molecule does not loose 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:

$$\begin{eqnarray}&&{t}_{\mathrm{scan}}={N}_{s}{\nu }_{0}^{-1}{e}^{\tfrac{{E}_{d}}{{k}_{B}{T}_{\mathrm{dust}}}},\end{eqnarray} \tag{ 2 }$$

where ν₀ is about 10¹² s⁻¹ 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 ≤ 1000 yr. In laboratory analogs 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 nanoseconds 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% (Table 1). This energy propagates like a wave from the point where the H + H reaction occurs (see Figure 4 and also Pantaleone et al. 2020). In ASW ice, a water molecule close to the 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–1680 K, average value counting all the water molecules within the first shell). 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 nonthermal 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 timescale 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 (1993) and 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 1999; Takahashi et al. 1999b; 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 ≤ 40 K for a time too short (∼10 fs) to allow CO to sublimate. However, despite the fact that Takahashi's simulations are impressive considering that they were carried out more than 20 yr 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 the 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 (Figure 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 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=\tfrac{v({{\rm{H}}}_{2})}{v(\mathrm{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 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³–10⁴ cm⁻³, temperature 10 K, a gaseous (undepleted) CO abundance with respect to H₂ equal to 2 × 10⁻⁴, and a cosmic-ray ionization rate of 3 × 10⁻¹⁷ 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⁴ cm⁻³) molecular clouds and solve a decades-long mystery. However, to put this affirmation on 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.

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) and Takahashi & Uehara (2001) predicted that H₂ formed on amorphous ice would be vibrationally excited to ν = 7–8. In contrast, in some laboratory experiments (Hornekær et al. 2003; Roser 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) ν (Table 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 "thermalized" 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); and
• (ii)  
  the possibility to detect newly formed H₂ with near-future James Webb Space Telescope (JWST) observations and, consequently, measure the rate of H₂ formation in molecular clouds and put constraints to theories.