Where Does the Energy Go during the Interstellar NH₃ Formation on Water Ice? A Computational Study

A set of static calculations, both in the gas phase and on the solid ice surfaces (crystalline and amorphous water ices; see below for more details), were performed with the CP2K software package, using the Quickstep module (Kühne et al. 2020). We used the DFT Perdew–Burke–Ernzerhof (PBE) method (Perdew et al. 1996), plus the a posteriori Grimmes D3(BJ) correction for dispersion interactions (Grimme et al. 2010, 2011). Core electrons were described by the Goedecker–Teter–Hutter pseudopotentials (Goedecker et al. 1996), whereas valence electrons were represented with a triple zeta MOLOPT basis set (VandeVondele & Hutter 2007) including a single polarization function (TZVP) with mixed Gaussian Plane Waves (GPW) approach (Lippert et al. 1997). The plane-wave cutoff was set to 600 ryd. In geometry optimizations, the reactive species were optimized on the thermalized slab, in which the atoms of the water molecules of the ice were fixed to keep the 10 K thermalized structure.

The Orca software (Neese et al. 2020) was also employed to check the accuracy of PBE-D3(BJ) for the energetics of the reactions in the gas phase. Calculations using the range-separated hybrid functional ωB97x-D3(BJ) (Najibi & Goerigk 2018), employing a def2-TZVP basis set (Weigend & Ahlrichs 2005), have been carried out, as well as single-point energy evaluations at CCSD(T)-F12 level of theory (Hattig et al. 2012), employing a cc-pVTZ-F12 basis set with an auxiliary cc-pVTZ-F12-CABS basis set (Valeev 2004; Peterson et al. 2008), to check the reliability of both DFT methods. All the calculations have been performed within the unrestricted formalism using the broken symmetry approach for biradical systems (Neese 2004). Results indicate that both methods provide similar results to CCSD(T) for reaction energies, hence demonstrating the robustness of the PBE-D3(BJ) method for the subsequent energy dissipation study with AIMD simulations at this theory level.

AIMD simulations were employed to study the fate of the energy released along the various reactions. For these simulations, we only focused on the amorphous water ice surface model, which was previously thermalized at 10 K in the NVT (canonical ensample with fixed number of particles (N), volume (V), and temperature (T)) for 1 ps (using a CSVR thermostat (Bussi et al. 2007) to ensure a sampling at a constant T) with a time step of 1 fs. The energy dissipation process was studied by means of NVE (microcanonincal ensample with fixed number of particles (N), volume (V), and total energy (E)) AIMD simulations following the energy released by every reaction and how this energy is partitioned between the newly formed molecules and the surface, in a similar way to works by some of us (Pantaleone et al. 2020, 2021). To estimate observables, moving averages of 50 MD steps have been calculated. The starting velocities of the water molecules belonging to the ice surface were recovered from those of the NVT thermalization procedure. In the NVE simulations, the time step was lowered to 0.5 fs to ensure a correct energy conservation of the system.

Finally, in order to gain insights into the energy dissipation mechanisms, we calculated the vibrational density of states (VDOS) for the water molecules of our surface model, as the Fourier transform of the velocity autocorrelation function of our trajectories using the TRAVIS analyzer (Brehm & Kirchner 2011; Brehm et al. 2020):

$$\begin{eqnarray}&&P(\omega ,t)={\int }_{{t}_{1}}^{{t}_{2}}\langle {\boldsymbol{v}}(t)\cdot {\boldsymbol{v}}(t+\tau )\rangle {e}^{-i\omega t}d\tau .\end{eqnarray} \tag{ 1 }$$

The VDOS is a powerful tool to assess the energy relaxation because it allows highlighting the surface vibrational modes that are excited by the dissipation of energy from the newly formed molecule. Therefore, it has here been employed to qualitatively understand the vibrational couplings among the components of systems.

In this work, two different models were used to mimic interstellar water ice surfaces. The first (and simplest) model is a crystalline slab constructed by cutting out the proton-ordered P-ice bulk along the (100) plane (Casassa et al. 1997). This periodic ice surface model has 72 water molecules in the unit cell (216 atoms), and the cell parameters are a = 13.156 Å, b = 14.162 Å, and α = β = γ = 90°. The c value, which modulates the interlayer distance among replicated slab images, was fixed to 30 Å, thus ensuring an empty space (about 22 Å) among fictitious slab replicas. This first model was mainly used for static calculations characterizing the potential energy surfaces of the H-addition and H-abstraction reactions. Figure 1(a) shows the structure of the crystalline ice surface slab model.

Zoom In Zoom Out Reset image size

To study the dissipation and partition of the nascent reaction energies with (NVE) AIMD simulations, we used a larger and amorphous slab model—larger with the purpose of avoiding nonphysical temperature increases caused by size limitations of the surface ice model, and amorphous with the purpose of simulating as realistic as possible an interstellar water icy mantle. To have a robust ice surface model accounting for the first aspect, we estimated the temperature increase of the entire system caused by the occurrence of the chemical reaction through the equipartition theorem, in which such an increase of temperature (ΔT) is given by

$$\begin{eqnarray}&&{\rm{\Delta }}T=\displaystyle \frac{2{E}_{\mathrm{nasc}}}{{{RN}}_{\mathrm{at}}},\end{eqnarray} \tag{ 2 }$$

where E_nasc is the nascent energy due to the reaction leading to the chemical bond formation (in this work a N–H for each hydrogenation step) and N_at is the number of atoms of the system. According to our calculations, the reaction energies lay between −340 and −440 kJ mol⁻¹. These values are similar to the reaction energy of the H₂ formation through the coupling of two H atoms. Accordingly, the expected ΔT is the same (30 K at most), and thus we used the same amorphous water ice surface model adopted by some of us to investigate the fate of the reaction energy in the interstellar H₂ formation (Pantaleone et al. 2021). This slab model was initially built up as a 2 × 2 supercell of the crystalline slab (a = 26.32 Å, b = 28.33 Å, and α = β = γ = 90°), in which more layers of water molecules along the z-axis (the nonperiodic direction) were added, resulting in 576 water molecules in the unit cell (1728 atoms). Because of these additions, the c value was enlarged to 50 Å (about 30 Å of real vacuum space). Such a crystalline supercell was then "amorphized" by a classical MD annealing process using the rigid TIP3P potential (Jorgensen et al. 1983). That is, a first run of 200 ps was performed in the NVT ensemble at 300 K to break the crystallinity of the initial structure, followed by a second NVT run at 10 K for 200 ps to cool down the ice at this interstellar temperature. The final disordered ice structure was optimized at the PBE-D3(BJ) level and then thermalized as described above. Figure 1(b) shows the structure of the amorphous ice surface slab model.

Previous to the AIMD simulations, we first carried out static calculations in order to elucidate the energetics of the processes under study. The reaction steps toward NH₃ formation and the competitive channels are shown in Figure 2.

Zoom In Zoom Out Reset image size

Hydrogenation of N to form NH₃ involves three steps (see Figure 2(a)), leading to the formation of NH, NH₂, and NH₃. Figure 2(a) also shows the electronic ground states of the involved species in the gas phase, i.e., N(⁴S), NH(³Σ⁻¹), NH₂(²B₁), and NH₃(¹A₁). On the water ice surface models, these electronic ground states are kept as such (i.e., no electronic intersystem crossings take place upon adsorption). Hydrogenation reactions, however, have competitive processes, i.e., H-abstraction, in which an incoming H atom abstracts an H atom belonging to the nitrogenated species, reverting the hydrogenation reaction and forming H₂. The possible H-abstraction reactions are shown in Figures 2(b) and (c). Figure 2(b) represents the H-abstraction on NH to form N and H₂. The resulting N can be either in the ⁴S ground electronic state or in the ²D first excited electronic state. Figure 2(c) represents the H-abstraction of NH₂ to form NH and H₂, in which the formed NH can be either in the ³Σ⁻¹ ground electronic state or in the ¹Δ first excited electronic state (Figure 2(c)). All these possible situations have been considered in this work.

Since the amount of energy to dissipate is, for the scope of this work, a fundamental quantity, we computed the reaction energies of the three hydrogenation steps in the gas phase at the PBE-D3(BJ) level, using CP2K and Orca. Results were compared with the reference CCSD(T)-F12 values and, to get deeper insights into the comparison, also with those obtained at ωB97x-D3(BJ). Additionally, energy barriers (if existent) of these reactions, as well as the energetics of the H-abstraction processes (see Figure 2), were also computed.

The three H-addition reactions to form NH₃ are barrierless in the gas phase, namely, the reactants spontaneously evolve toward the products during the optimization, because they involve radical−radical coupling processes and, thus, are driven by spin couplings. Computed reaction energies (reported in Table 1) are all negative and very large (ranging from roughly −360 to −460 kJ mol⁻¹), in which, moreover, the more hydrogenated the product, the more favorable the reaction energy. There is a good agreement among the values computed with all the quantum chemical methods employed and the two codes, the unsigned standard errors being 1%–11% (using CCSD(T)-F12 as the reference method). Results indicate that the PBE-D3(BJ) level of theory is accurate enough to describe the energetics of the H-additions.

In relation to the competitive H-abstraction reactions (results available as supplementary material; Ferrero 2022), it is found that those processes leading to products in their electronic ground state (i.e., N(S) and NH(³Σ⁻¹)) are exoergic because of the favorable energy of formation of the H–H bond. Nevertheless, this is not the case for processes forming excited state products (i.e., N(²D) and NH(¹Δ)), which counterbalance the formation energy of the H–H bond, resulting in endoergic reactions, not possible at the cold ISM temperatures. Exoergic reactions present energy barriers of about 10–11 kJ mol⁻¹ and 26–29 kJ mol⁻¹ for the formation of N(⁴S) and NH(³Σ⁻¹), respectively (computed at ωB97x-D3(BJ)//ωB97x-D3(BJ) and CCSD(T)-F12//ωB97x-D3(BJ) levels of theory). Thus, we can conclude that H-abstraction reactions cannot be considered as efficient competitive channels to the H-addition ones, since they either are endoergic processes (while H-additions are all largely exoergic) or present high-energy barriers for interstellar conditions (while H-additions are all barrierless). Therefore, to deal with the reaction energy dissipation, we will solely consider the H-additions as the main-occurring reaction channels.

On the water ice surfaces, we studied the reactions according to a Langmuir–Hinshelwood (LH) mechanism, namely the reactive species are on the surface and diffuse until they encounter each other. To study the energetics of the reactions, we make here an important assumption: the diffusion of the species on the ice surface already happened. Accordingly, the initial states involve the two reactants in close proximity to favor the reaction. Thus, we performed static optimizations on the ice surface models, placing the reactants in two adjacent adsorption sites (i.e., distances of 4.0–4.3 Å on the crystalline surface and of 3.6–4.0 Å on the amorphous one), as depicted in Figure 3. We optimized first these pre-reactant species in the electronic states that prevent the occurrence of the H-addition reactions due to the Pauli/exclusion repulsion principle (namely, quintet, quadruplet, and triplet for the first, second, and third H-addition, respectively), this way obtaining an optimized geometry of these complexes based on only species/surface interactions. From these pre-reactant adducts, then optimizations at the reactive electronic states (namely, triplet, doublet, and singlet, respectively) were performed. It is worth mentioning that, at variance to the gas phase, on the surface, the reactant/surface interactions must be partly broken as the reaction takes place. Despite this, our results indicate that, also on the icy surfaces, all the H-addition reactions occur spontaneously during the optimization processes and, accordingly, they are also barrierless processes. To ensure that this chemical behavior is kept also when dynamic effects are accounted for, for each reaction we run a short AIMD simulation at 10 K, and we observed that reactants also spontaneously evolve to products in the first femtoseconds.

Zoom In Zoom Out Reset image size

Figure 3 also shows the calculated reaction energies of the processes. On the surfaces, similarly to the gas phase, the reactions are exoergic (between −340 and −440 kJ mol⁻¹; values obtained taking as reference the pre-reactant complexes in their nonreactive electronic states). Those reactions occurring on the amorphous model are slightly less negative than on the crystalline model (see Figure 3) owing to the different number of hydrogen bonds established by the adsorbed species with the different surfaces. The trend in which at each H-addition step the reaction energy is more favorable is also kept. In view of these results, we therefore can conclude that placing the reactants at nearly 4 Å leads to very favorable H-addition reactions (i.e., barrierless and with large and negative reaction energies) so that the hydrogenation reactions are not energetically hampered. However, it is worth bearing in mind that, depending on the distance between the two reactants, as well as the surface morphology, the reactions are limited by diffusion.

Here we describe the results concerning the energy dissipation in the three H-additions to the atomic N that result from the AIMD simulations.

The initial guess structures were the optimized geometries of the pre-reactant complexes in their nonreactive electronic states, in which the reactive species are roughly 3.6–4.0 Å apart (reactant structure of Figures 3(d), (e), and (f)). This intermolecular distance is a suitable choice because it is short enough to neglect the diffusion of the species but long enough to include at least 94% of the reaction energy in the AIMD simulations. The kinetic energy liberated by the reactions is studied by following the evolution of the velocities of the system particles throughout the entire trajectories. This analysis allowed us to monitor the kinetic energy partitioning between the newly formed species and the ice surface, providing unique information on (i) the kinetic energy fraction retained by the product, (ii) the amount of energy dissipated through the ice surface phonon modes, (iii) the temperature increase of the ice due to this energy dissipation, and (iv) the possibility of desorption of the product due to retaining a large fraction of the kinetic energy. Table 2 summarizes part of these data for the studied processes.

Table 2. Data Obtained from the (NVE) AIMD Simulations for the Reactions Investigated

Reaction	$\tfrac{{T}_{\mathrm{ice}}}{{T}_{\mathrm{ice}}+{T}_{S}}$	$\tfrac{{T}_{S}}{{T}_{\mathrm{ice}}+{T}_{S}}$	Temperature Variation	BE	TKE	TKE-z
N+H → NH	0.76	0.24	10 → 21.6	17.5	2.7	1.4
NH + H → NH2	0.71	0.29	10 → 23	18.8	2.9	1.3
NH2 + H → NH3 (Pos1)	0.73	0.27	10 → 24.1	22.4	5.2	1.4
NH2 + H → NH3 (Pos2)	0.58	0.42	10 → 23.8	34.4	3.5	2.2
NH2 + H → NH3 (Pos3)	0.90	0.10	10 → 27.6	36.3	2.9	2.2

Note. The first column represents the reactions studied. The second and third columns report the fraction of kinetic energy transferred to the ice and retained by the newly born species (S), respectively. The fourth column shows the variation of the temperature (in K) undergone by the ice at the end of the simulation. The fifth column reports the computed binding energies (in kJ mol⁻¹) of the species at the reaction sites. The sixth and seventh columns summarize the translational kinetic energy (TKE) and the TKE along the z-axis (TKE-z) of the newly formed species (in kJ mol⁻¹), respectively.

Download table as:  ASCIITypeset image

According to the static calculations, the first H-addition is the less exothermic step, with a reaction energy of −343.5 kJ mol⁻¹ (see Figure 3(d)), a considerable amount of energy, whose distribution has been analyzed by following the kinetic energy as depicted in Figure 4.

Zoom In Zoom Out Reset image size

Panel (a) shows the initial positions of N and H on the amorphous surface, panel (b) the evolution of the kinetic energy of the newly formed NH molecule, panel (c) the evolution of the NH distance (starting at roughly 3.6 Å), panel (d) the evolution of the kinetic energy of the ice surface, and panel (e) the temperature variation of the ice along the simulation. Figure 4(b) shows initial peaks reaching ≈350 kJ mol⁻¹, which is very close to the NH reaction formation energy on the ice surface, indicating that this is the actual amount of energy released by the reaction. At the beginning of the simulation, the NH species forms in a highly excited vibrational state, as indicated by the high and narrow energy spikes in Figure 4(b). Such transient excited states induce large oscillations of the forming N–H bond (see Figure 4(c)). Nevertheless, already in the first femtoseconds, the water ice molecules close to the reaction site partly absorb the reaction energy, a total of about 157 kJ mol⁻¹ during the whole simulation (see Figure 4(d)). This energy transfer induces a temperature increase of the ice from 10 K to an average of 22 K (see Figure 4(e)). Despite this thermal relaxation of the NH molecule by exchanging energy with the surface, it also retains a certain amount of kinetic energy, about 50 kJ mol⁻¹ (see average line of Figure 4(b)). At the end of the simulation, the kinetic energy liberated by the reaction is distributed as follows: 76% is transferred to the ice, and 24% is retained by the product (see Table 2). In the supplementary material (Ferrero 2022), the evolutions along the AIMD simulation of the most relevant energetic components, i.e., the total energy (E_TOT, potential + kinetic), the potential energy (V_TOT), and the kinetic energy (T_TOT), are compared, showing the energy conservation of our AIMD simulations in consistency with the NVE ensemble.

Results of the second H-addition, i.e., NH + H → NH₂, are shown in Figure 5. The first peak of the kinetic energy of the newly formed NH₂ (Figure 5(b)) is in good agreement with the reaction energy of the process (−382 kJ mol⁻¹). This second hydrogenation shows a faster relaxation in comparison with the first one, as deduced from the narrower oscillations in the kinetic energy of NH₂ (Figure 5(b)) and from the lower oscillations of the formed N–H bond distance (Figure 5(c)). The reason for this faster relaxation is that NH₂ has more internal degrees of freedom than NH, which can strongly couple with the newly formed N–H bond, as well as with the vibrations of the icy water molecules of the surface. As the released energy for this reaction is larger than the first H-addition, more kinetic energy is dissipated through the ice surface, about 180 kJ mol⁻¹ (see Figure 5(d)), leading to a temperature increase of the ice of up to 23 K (see Figure 5(e)), and the internal kinetic energy retained by the newly formed NH₂ molecule is of about 73 kJ mol⁻¹ (see average line of Figure 5(b)), resulting in a distribution of 71% and 29%, respectively (see Table 2).

Zoom In Zoom Out Reset image size

The third H-addition reaction (i.e., NH₂ + H → NH₃) is that releasing the largest amount of energy (≈−440 kJ mol⁻¹). Because of that, and also due to the fact that this is the last step to form the final NH₃ product, for this reaction we calculated the trajectories from three different starting positions of the reactants, with the aim to consider different possible situations as far as the surface morphology is concerned. The chosen initial positions are shown in Figure 6 and present the following features: (i) Pos1 (the same initial positions as the first and second H-additions), in which NH₂ and H are placed in an outermost position of the surface (i.e., not inside a cavity) and NH₂ acts as an H-bond donor; (ii) Pos2, in which NH₂ and H are placed in an outermost position of the surface and NH₂ acts as an H-bond acceptor; and (iii) Pos3, in which NH₂ and H are placed within a cavity and NH₂ acts as an H-bond donor. By proceeding this way, we roughly sample different sites of the amorphous surface that can likely exhibit specific trends in the energy dissipation, as observed for the H₂ formation case (Pantaleone et al. 2021).

Zoom In Zoom Out Reset image size

For the AIMD simulation from Pos2, the trends in the energy partition are similar to those observed for the other H-additions. Graphs on the evolution of the kinetic energies, oscillation of the newly formed N–H bond, and temperature variation of the ice are shown in the supplementary material (Ferrero 2022). The NH₃ molecule remains highly excited during the whole simulation, which is reflected by large oscillations of both the kinetic energy and the N–H bond length. The amount of the reaction energy that is transferred to the surface as kinetic energy is about 159 kJ mol⁻¹, resulting in a temperature increase of the ice from 10 to 24 K along the simulation. In contrast, in the simulated trajectories starting from Pos1 and Pos3, one can distinguish a different energy dissipation behavior, as provided by the graphs shown in Figures 7(a) and (e) (for Pos1 and Pos3, respectively): the newly formed NH₃ molecule exchanges more energy with the surface and in shorter timescales. This can be seen from the abrupt change in the oscillation amplitude of the kinetic energy of NH₃ and from the ice temperature variation, while there is a sudden increase in the evolution of the kinetic energy of the ice at these simulation times. For Pos1, the kinetic energy of NH₃ is about 78 kJ mol⁻¹, while that transferred to the surface is 210 kJ mol⁻¹, and for Pos3 they are 34 and 321 kJ mol⁻¹, indicating that the amount of the nascent energy transferred to the surface is larger in Pos1 and Pos3 than in Pos2. Interestingly, in Pos1 and Pos3, the formation of an ion pair transient species of the kind H₃O⁺/NH${}_{2}^{-}$ at around 500 and 200 fs, respectively, takes place, which is crucial to efficiently dissipate the large amount of reaction energy. The occurrence of this event (i.e., formation of a transient species) and its relevance in the energy dissipation are addressed in the next section.

Zoom In Zoom Out Reset image size

An important aspect to consider here is the possibility that CD of the newly formed products can occur after the reaction. A plausible CD mechanism is that the newly formed species desorbs because it conserves part of the reaction energy in the form of translational energy, which is used to desorb toward the gas phase. To assess this point, we extracted the translational kinetic energy (TKE) of the newly formed species as the kinetic energy associated with the center of mass (Takahashi et al. 1999). The TKE component along the axis normal to the surface (in this case the z-axis) is the component that would bring the newly formed species toward the gas phase. Therefore, the z component of the TKE (TKE-z) was extracted and compared to the binding energies (BEs) of the newly formed species (see Table 2). We find out that, in all cases, TKE-z contribution accounts for just a small fraction of the BE (between 4.3% and 8.0%) and, accordingly, CD is not expected. However, it could also be possible that the newly formed species, while diffusing on the ice surface once formed, can stumble over a surface irregularity (e.g., dangling H atoms), and if it was so, the total TKE could be canalized into the desorption direction. This was indeed the case of H₂ formation (Pantaleone et al. 2021), which was ejected into the gas phase upon stumbling over an outermost water molecule of the surface. Accordingly, to assess this point, we compared the total TKE retained by the born species with their BEs. Like in TKE-z, none of the total TKE values are larger than the computed BEs (total TKEs are between 8% and 16% of the corresponding BEs), so we can conclude, again, that CD is not expected to take place through this mechanism. Therefore, the formed species are doomed to remain attached to the ice surface.

On the other hand, it is worth mentioning that the timescale in which most of the energy dissipates throughout the surface takes place in about 1 ps. However, the H-addition reactions, as they are barrierless, are controlled by the diffusion of the H atoms. It is thus interesting to assess whether the H diffusion can compete with the energy dissipation because, in the case of similar timescales, the intermediate H-addition reactions (i.e., the second and third hydrogenations) can take place with the newly formed species in excited vibrational states. To check qualitatively on this point, we have proceeded as follows. It is estimated in the literature that the H diffusion on water ice can be 5.8 × 10⁻¹¹ cm² s⁻¹ or 3.3 × 10⁻⁷ cm² s⁻¹ at 25 K (Ásgeirsson et al. 2017), depending on the presence/absence of inert species blocking the most stable adsorption sites, and considering that these numbers are upper limits, as they assume that there is always an H atom ready to diffuse, omitting the probability that an H atom lands on the surface (Watanabe et al. 2010; Hama et al. 2012). According to these diffusion coefficients, an H atom scans a surface area of 745.6 Ų (the area of our unit cell surfaces) in 1.28 ms and 0.89 ns, respectively, which are significantly longer than our energy dissipation timescale (1 ps). This supports the idea that the formation of NH₃ by hydrogenation of N is a true stepwise H-addition process, in which the NH formation energy is efficiently dissipated throughout the grain while the newly formed species get relaxed from any vibrationally excited states before entering the next hydrogenation step.

The conclusions described above, i.e., the nonoccurrence of CD and that the energy relaxation of the newly formed species is faster than the H diffusion, are very important aspects for the interstellar chemistry. That is, the NH, NH₂, and NH₃ species are not ejected into the gas phase once formed, thereby ruling out a direct CD mechanism. After each hydrogenation step, the associated reaction energies are promptly dissipated, and the newly formed species remain in low vibrational states, ready for a new hydrogenation (unlike if they become desorbed, which, to be the case, would be in excited vibrational states, presenting a different reactivity compared to their fundamental states). Remarkably, this picture is consistent with the fact that NH₃ is a component of the interstellar ices: after the completion of the full H-addition sequence, NH₃ remains on the water ice surface, which, upon subsequent H₂O accretion and/or in situ formation on the grains, can become a bulky component of the ice mantles, in agreement with the observational IR features of interstellar ices (Boogert et al. 2015). However, it is worth mentioning that our simulations do not account for additional surface phenomena that can induce CD. This would be, for instance, the occurrence of other large exothermic reactions (e.g., the H₂ formation) nearby the newly formed species, which can produce a local heating in their surroundings, triggering the desorption into the gas. In fact, this type of indirect CD is usually included in the astrochemical models and allows us to explain the gas-phase abundances of relevant interstellar molecules (e.g., H₂O, CO, H₂CO, CH₃OH, and several iCOMs; Garrod et al. 2007; Vasyunin & Herbst 2013; Hocuk & Cazaux 2015; Ruaud et al. 2015; Minissale et al. 2016a), by considering that a few percent of these components, formed on the grain surfaces, are ejected into the gas phase through this nonthermal mechanism. Interestingly, direct CD is also invoked to occur through H-abstraction reactions. This has experimentally been investigated for the CO and H₂CO cases (Minissale et al. 2016b), in which authors found that for the [CO + H] and [H₂CO + H] reactive systems CD of CO and H₂CO was detected, presumably due to the ease of reaction and large exothermicity of the HCO + H → CO + H₂ and CH₃O/H₂COH + H → H₂CO + H₂ H-abstraction processes. Nonetheless, similar CD is not expected for the N-bearing species dealt with in this work because the H-abstraction reactions, as shown above, present a limited exothermicity, rendering unlikely the CD of these species.

Due to the large amount of energy released by every H-addition step, each product is formed in a highly excited vibrational state, which relaxes to lower vibrational states at different timescales. This is indicative that different dissipation energy mechanisms are involved in the relaxation of the newly born species.

Upon formation, the vibrational modes of the formed species (stretching and bending for NH₂ and NH₃, and only stretching for NH) are highly excited, while this is not the case for the 10 K thermalized internal modes of the icy water molecules (consisting of stretching and bending vibrations of individual water molecules and the low-energy libration modes arising from the lattice nature of the water ice surface). Accordingly, at first instance, there is a coupling between the excited stretching and bending modes of the newly formed species and the libration modes of the underneath water ice surface. To have evidence on that, we calculated the trajectory VDOS power spectrum (which highlights the active vibrational modes during the simulations) for the water ice surface when the reactions occur and for the isolated water ice surface thermalized at 10 K, and we plot the difference of the computed VDOS power spectra (ΔVDOS) to highlight the difference in the population of the vibrational modes. Figures 8(a), (b), and (c) show the ΔVDOS plot during the reactions involving the formation of NH, NH₂, and NH₃ (Pos2), whereas Figures 8(e) and (e) show the same but for NH₃ formation from Pos1 and Pos3. In all the cases, the bands that become more populated correspond to the libration frequencies. Interestingly, for the formation of NH, NH₂, and NH₃ (Pos2), the water ice bending modes become populated, whereas the stretching ones remain unpopulated, indicating a vibrational coupling between the bending modes of the newly formed species and the ones of the water molecules in the ice surface. These results thus clearly indicate that a main mechanism responsible for the dissipation of the large excess of the reaction energy is the intermolecular vibrational coupling between the vibrational modes of the born species and the libration modes (with a minor contribution to the bending modes) of the water ice.

Zoom In Zoom Out Reset image size

For reactions represented in Figures 8(d) and (e), in addition to the libration modes, there is a significant population of the icy water stretching modes and a larger population of the bending modes compared to the ΔVDOS plots of Figures 8(a), (b), and (c). Thus, for these cases an additional energy mechanism operates, which in this case is based on the formation of the transient H₃O⁺/NH${}_{2}^{-}$ ion pair species (see Figure 9) at 400 and 200 fs during the NH₃ formation (Pos1 and Pos3 trajectories, respectively). The ion pair nature of these species is confirmed by their calculated charges: the N charge is lowering from −0.34 to −0.74, and the charge of the oxygen atom involved in the ion pair is changing from −0.50 to −0.38, approximately. These species are formed by a proton transfer from NH₃ to a nearby surface H₂O molecule. This proton transfer could seem surprising, as in water solution the expected ion pair is NH4⁺/OH⁻. Here the opposite reaction to form H₃O⁺/NH${}_{2}^{-}$ occurs as a result of the large internal energy of the newly formed NH₃ in the early stages of its formation, in which the energy of the highly excited N–H stretching modes is enough to transfer the proton toward a nearby water molecule. As the ion pair is formed, a sudden energy transfer occurs toward the water ice surface (see Figure 7). Formation of these species permits us to populate the high‐frequency bending and stretching modes of nearby water molecules, but in this case the population is not exclusively due to a vibrational coupling between vibrational modes of the two partners but also through a "by-side" chemical reaction (a proton transfer), which injects kinetic energy from the newly formed NH₃ into the stretching modes of the icy water molecules. This "chemical" energy transfer is highlighted in the ΔVDOS plots shown in Figures 8(d) and (e), where the stretching peak around 3600–3700 cm⁻¹ is populated, while completely absent in the spectra of Figures 8(a), (b), and (c), in which the ion pair is not formed. Obviously, the population of these peaks is rather low, as the ion pair has a short lifetime (50–100 fs), while the VDOS spectra is for the entire trajectory (1200 fs). Interestingly, the fate of the transient H₃O⁺/NH${}_{2}^{-}$ ion pair in Pos1 and Pos3 is different, which affects the efficiency of the relaxation mechanism. While in Pos1 the transferred proton reverts directly to NH₃ (see Figure 9(a)), in Pos3 the proton does not rebind to NH${}_{2}^{-}$ but remains on H₃O⁺, which starts a proton transfer to the NH${}_{2}^{-}$ species (recovering NH₃) via a proton shuttle mechanism involving three water molecules (see Figure 9(b)).

Zoom In Zoom Out Reset image size

The proton relay mechanism ensures an even more efficient relaxation of the newly formed NH₃: the transferred proton, engaging more than one water molecule, allows for a higher degree of energy dissipation toward the icy surface. It is worth mentioning that, due to the high computational cost of these simulations, collecting enough trajectories to arrive at a robust statistical picture on the occurrence of this "chemical" energy dissipation mechanisms would be overwhelming. Therefore, with the simulated trajectories, we can have hints on the process: (i) the formation of the transient ion pairs is geometrically constrained and is only possible if the proton belonging to the N–H vibrationally excited bond is close enough to the lone pair of an O atom of the ice surface acting as a proton acceptor, and (ii) the formation of the ion pair is more probable from NH₃ than from NH₂ and NH, in agreement with their deprotonation energy (NH₃ < NH₂ < NH). Despite this, we cannot exclude the occurrence of this mechanism also during the NH and NH₂ formation, since it was also observed in the H₂ formation (Pantaleone et al. 2021).

Finally, we would like to stress that we are not including nuclear quantum effects (NQEs; Ceriotti et al. 2016; Markland & Ceriotti 2018; Huppert et al. 2022), which can be of great significance for hydrogen bonded systems at these very low temperatures. NQEs can, indeed, play a crucial role in both the energy transmission and energy dissipation processes. However, considering the size of the systems investigated, inclusion of NQEs in our simulations is computationally almost out of reach and unpractical. Along this line, efforts dedicated to alleviating the zero-point energy (ZPE) leakage problem in classical trajectory calculations are ongoing (Habershon & Manolopoulos 2009; Brieuc et al. 2016; Buchholz et al. 2018; Mangaud et al. 2019; Mukherjee & Barbatti 2022).

Ammonia is a very important molecule in astrochemistry, as it is a major reservoir of nitrogen in the molecular regions of the ISM. Its chemistry has been studied by various authors, and a simple scheme has been shown in Tinacci et al. (2022) and De Simone et al. (2022), from which the present Figure 10, describing the ammonia formation in cold molecular clouds, has been adapted for helping the discussion here. Briefly, ammonia is prevalently formed at the beginning of the evolution from a translucent to a molecular cloud (e.g., Ceccarelli et al. 2022) by the addition of H atoms to the N atoms frozen onto the interstellar grain surfaces, the process studied in this work (Figure 10). The hydrogenation of frozen nitrogen atoms continues also during the molecular cloud phase, but with smaller efficiency, as the remaining gaseous nitrogen becomes incorporated into molecular nitrogen and some of it goes into ammonia via gas-phase reactions. In cold molecular clouds, the only mechanism to eject the frozen ammonia into the gas phase is via the CD during the N hydrogenation process because the temperature of cold molecular clouds (∼10 K) is too low with respect to the ammonia BE for thermal sublimation to be efficient (Tinacci et al. 2022).

Zoom In Zoom Out Reset image size

In practice, the bulk of ammonia observed in the ISM, when considering the frozen and gaseous form altogether, is built on the grain surfaces at early times, and it is mainly due to the hydrogenation of nitrogen atoms on the grain surfaces. Please note that frozen water coats very quickly the bare silicate/carbonaceous surface (Tinacci et al. 2022), so that the nitrogen hydrogenation mostly occurs on icy surfaces. If the hydrogenation process were counterbalanced by the H-extraction competition or the dissipation of the chemical energy were inefficient, making possible a large degree of CD, the amount of formed ammonia would be reduced. Our computations show, on the contrary, that the H-abstraction processes are completely negligible and that direct CD is an extremely improbable process, supporting the assumptions of the present astrochemical models that are based on indirect CD.