<figure><img src="image_1.png" /></figure>

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

Stefano Ferrero<fnref n="1" />, Stefano Pantaleone<fnref n="2" />, Cecilia Ceccarelli<fnref n="3" />, Piero Ugliengo<fnref n="2" />, Mariona Sodupe<fnref n="1" />, and Albert Rimola<fnref n="1" />

Received 2022 September 27; revised 2022 December 21; accepted 2022 December 22; published 2023 February 21

## Abstract

In the coldest (10–20 K) regions of the interstellar medium, the icy surfaces of interstellar grains serve as solid-state supports for chemical reactions. Among their plausible roles, that of third body is advocated, in which the reaction energies of surface reactions dissipate throughout the grain, stabilizing the product. This energy dissipation process is poorly understood at the atomic scale, although it can have a high impact on astrochemistry. Here we study, by means of quantum mechanical simulations, the formation of NH₃ via successive H-additions to atomic N on water ice surfaces, paying special attention to the third-body role. We first characterize the hydrogenation reactions and the possible competitive processes (i.e., H-abstractions), in which the H-additions are more favorable than the H-abstractions. Subsequently, we study the fate of the hydrogenation reaction energies by means of ab initio molecular dynamics simulations. Results show that around 58%–90% of the released energy is quickly absorbed by the ice surface, inducing a temporary increase of the ice temperature. Different energy dissipation mechanisms are distinguished. One mechanism, more general, is based on the coupling of the highly excited vibrational modes of the newly formed species and the libration modes of the icy water molecules. A second mechanism, exclusive during the NH₃ formation, is based on the formation of a transient H₃O⁺/NH⁻₂ ion pair, which significantly accelerates the energy transfer to the surface. Finally, the astrophysical implications of our findings relative to the interstellar synthesis of NH₃ and its chemical desorption into the gas are discussed.

*Unified Astronomy Thesaurus concepts: Surface ices (2117); Interstellar dust (836); Interstellar molecules (849); Dense interstellar clouds (371); Interstellar medium (847); Solid matter physics (2090); Interstellar dust processes (838); Computational methods (1965)*

## 1. Introduction

Dust grains are the primarily solid-state particles of the interstellar medium (ISM). Their chemical composition, as can be inferred by observationally IR spectroscopic measurements, is based on a core of silicates or carbonaceous materials (Henning 2010; Jones et al. 2013, 2017), coated by ice mantles predominantly made by water. These ice mantles are usually referred to as “dirty” ices because, besides water, they also contain several other volatile molecules, e.g., CO, CO₂, CH₄, NH₃, and CH₃OH, among others (Boogert et al. 2015). Spectroscopic comparisons between spectra from laboratory ice analogs and observational spectroscopic measurements suggest that interstellar ices have an amorphous structure and are relatively porous (Fraser et al. 2004; Watanabe & Kouchi 2008; Potapov et al. 2020). It has long been recognized that interstellar dust grains are of extreme importance for the chemistry of the ISM since, among other roles (e.g., chemical catalysts), they can act as third bodies (or energy sinks), in which the energy released by the reactions occurring on the grain surfaces is dissipated throughout the grains, this way stabilizing the newly formed products. The third-body effect of the interstellar grains on chemical reactions is a subject of interest in astrochemistry because it has direct implications for the occurrence of chemical desorption (CD) of species synthesized on the grain surfaces. CD is a desorption

mechanism driven by chemistry, in which an adsorbed molecule on the grain surface is ejected into the gas phase because, once formed on the surface, it desorbs owing to the local heating caused by the exothermicity of the reaction. Since the interstellar temperatures are very cold (usually below the sublimation of the icy species, making standard thermal desorption inoperative), CD is an important and frequently resorted mechanism accounted for in the astrochemical models to justify the abundances of gas-phase species that otherwise would not be possible without being considered.

Hydrogenation reactions are the most common interstellar reactions taking place on the surfaces of dust grains (Watanabe & Kouchi 2008; Hama & Watanabe 2013), because atomic hydrogen is the most abundant element in space and the most mobile atom that can diffuse over the grain surfaces at the ultracold temperatures of the ISM (10–20 K). An emblematic grain surface reaction is the formation of H₂, the most abundant molecule in the universe, which is formed by the coupling of two H atoms (Vidali 2013; Wakelam et al. 2017). In the gas phase, the coupling of two H atoms in their electronic ground states does not take place because, in the absence of a third body (e.g., a grain), the reaction energy associated with the H–H chemical bond formation cannot be released, consequently leading to its re-dissociation. The partitioning and dissipation of reaction energies are key surface phenomena in the synthesis of compounds taking place in the presence of interstellar grains, such as the formation of simple but relevant molecules (e.g., H₂, H₂O, CH₃OH; Watanabe & Kouchi 2002; Watanabe et al. 2004; Fuchs et al. 2009; Dulieu et al. 2010; Van Dishoeck et al. 2013; Vidali 2013; Rimola et al. 2014; Wakelam et al.

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, E-08193, Catalonia, Spain; mariona.sodupe@uab.cat, albert.rimola@uab.cat</footnote>

<footnote n="2">²Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, I-10125, Torino, Italy</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), F-38000 Grenoble, France</footnote>
2017; Qasim et al. 2018) and of the so-called interstellar complex organic molecules (iCOMs; Ceccarelli et al. 2017; Enrique-Romero et al. 2019; Zamirri et al. 2019; Gutiérrez-Quintanilla et al. 2021; Enrique-Romero et al. 2022; Perrero et al. 2022), these later ones inferring a primogenital organic chemistry to the ISM (Ceccarelli et al. 2022). However, energy partitioning and dissipation are hitherto poorly understood at an atomistic level. Nevertheless, ab initio molecular dynamics (AIMD) based on density-functional theory (DFT) are very powerful computational simulations that allow us to gain insights into these processes because they are based on electronic structure calculations (and hence can deal with the breaking/formation of chemical bonds) and naturally account for anharmonic effects, including the coupling among the vibrational modes of the molecule/surface system. For large atomistic models, this methodology becomes computationally very expensive, and therefore true statistics of the processes under study cannot be provided. This can only be possible by means of classical molecular dynamics (MD) simulations, in which hundreds of simulations with long simulation timescales can be executed, this way allowing a statistical treatment of the results (Fredon et al. 2017; Fredon & Cuppen 2018). However, simulations based on classical force fields cannot cope with the breaking/formation of chemical bonds, and accordingly they are not suitable to investigate the dissipation of energies arising from chemical reactions of molecules synthesis. To cope with these problems, in addition to the AIMD simulations (Pantaleone et al. 2020, 2021), other sophisticated approaches have been applied to the study of energy dissipation at surfaces, such as reactive classical MD (Pezzella & Meuwly 2019; Upadhyay & Meuwly 2021), embedding schemes (Meyer & Reuter 2014; Rittmeyer et al. 2018), nonequilibrium MD (Melani et al. 2021), and also MD based on neural network potentials (Shakouri et al. 2018).

In the present work, we have studied the formation of interstellar NH₃, one of the components forming part of the ice mantles, from the successive hydrogenation (i.e., H-addition) reactions of atomic N on water ice mantles, with the particular focus to analyze, at a molecular scale, how the reaction energies of each elementary step partition and dissipate and to identify the mechanisms through which these events take place, all in all by means of AIMD simulations.

## 2. Methodology

### 2.1. Computational Details

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 & Goergik 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 CVSR 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 (microcanonical 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):

$$P(\omega, t) = \int_{t_1}^{t_2} \langle \mathbf{v}(t) \cdot \mathbf{v}(t + \tau) \rangle e^{-i\omega t} d\tau. \quad (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.

### 2.2. Water Ice Surface Models

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
<figure><img src="image_3.png" /><figcaption>Figure 1. Water ice surface models used in this work. (a) Crystalline model of the (100) water ice cut out from the proton-ordered P-ice bulk. (b) Amorphous surface model derived from the amorphization of the crystalline surface (see text for more details).</figcaption></figure>

cell (216 atoms), and the cell parameters are $a = 13.156$ Å, $b = 14.162$ Å, and $\alpha = \beta = \gamma = 90^{\circ}$. 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.

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 ($\Delta T$) is given by

$$ \Delta T = \frac{2E_{\text{nasc}}}{RN_{\text{at}}}, \qquad (2) $$

where $E_{\text{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_{\text{at}}$ is the number of atoms of the system. According to our calculations, the reaction energies lay between -340 and -440 kJ mol$^{-1}$. These values are similar to the reaction energy of the H$_2$ formation through the coupling of two H atoms. Accordingly, the expected $\Delta 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$_2$ 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 $\alpha = \beta = \gamma = 90^{\circ}$), 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.

## 3. Results

### 3.1. Energetics of the Reactions

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$_3$ formation and the competitive channels are shown in Figure 2.

Hydrogenation of N to form NH$_3$ involves three steps (see Figure 2(a)), leading to the formation of NH, NH$_2$, and NH$_3$. Figure 2(a) also shows the electronic ground states of the involved species in the gas phase, i.e., N($^4$S), NH($^3\Sigma^{-1}$), NH$_2(^2B_1)$, and NH$_3(^1A_1)$. 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$_2$. 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$_2$. The resulting N can be either in the $^4$S ground electronic state or in the $^2$D first excited electronic state. Figure 2(c) represents the H-abstraction of NH$_2$ to form NH and H$_2$, in which the formed NH can be either in the $^3\Sigma^{-1}$ ground electronic state or in the $^1\Delta$ 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 $\omega$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$_3$ 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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Scheme of the sequential H-addition (panel (a)) and of the H-abstraction (panels (b) and (c)) reactions.</figcaption></figure>

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.

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.
<table><caption>**Table 1**

Computed Reaction Energies (ΔE<sub>rx</sub>) for the Sequential H-addition Reactions in the Gas Phase Using Different Quantum Chemical Methods (PBE-D3(BJ), ωB97x-D3 (BJ), and Single Points at CCSD(T)-F12) and Codes (CP2K and Orca)</caption><thead><tr><th rowspan="2">ΔE<sub>rx</sub></th><th colspan="1">CP2K</th><th colspan="4">Orca</th></tr><tr><th>PBE-D3(BJ)</th><th>PBE-D3(BJ)</th><th>CCSD(T)-F12//PBE-D3(BJ)</th><th>ωB97x-D3(BJ)</th><th>CCSD(T)-F12//ωB97x-D3(BJ)</th></tr></thead><tbody><tr><td>N+H → NH</td><td>-359.3</td><td>-367.2</td><td>-330.8</td><td>-351.5</td><td>-331.2</td></tr><tr><td>NH + H → NH<sub>2</sub></td><td>-406.5</td><td>-415.4</td><td>-399.7</td><td>-411.5</td><td>-399.9</td></tr><tr><td>NH<sub>2</sub> + H → NH<sub>3</sub></td><td>-463.6</td><td>-472.3</td><td>-466.2</td><td>-473.4</td><td>-466.0</td></tr></tbody></table>

**Note.** All the H-addition reactions have been found to be barrierless, namely, the products are spontaneously formed during the optimization process. Units are in kJ mol<sup>-1</sup>.

## 3.2. Energy Dissipation in the H-addition Reactions

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.

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

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<sup>-1</sup>, 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<sup>-1</sup> 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<sup>-1</sup> (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<sub>TOT</sub>, potential + kinetic), the potential energy (V<sub>TOT</sub>), and the kinetic energy (T<sub>TOT</sub>), 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<sub>2</sub>, are shown in Figure 5. The first peak of the kinetic energy of the newly formed NH<sub>2</sub> (Figure 5(b)) is in good agreement with the reaction energy of the process (-382 kJ mol<sup>-1</sup>). 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<sub>2</sub> (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<sub>2</sub> 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<sup>-1</sup> (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<sub>2</sub> molecule is of about 73 kJ mol<sup>-1</sup> (see average line of Figure 5(b)), resulting in a distribution of 71% and 29%, respectively (see Table 2).

The third H-addition reaction (i.e., NH<sub>2</sub> + H → NH<sub>3</sub>) is that releasing the largest amount of energy (≈-440 kJ mol<sup>-1</sup>). Because of that, and also due to the fact that this is the last step to form the final NH<sub>3</sub> 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<sub>2</sub> and H are placed in an outermost position of the surface (i.e., not inside a cavity) and NH<sub>2</sub> acts as an H-bond donor; (ii) Pos2, in which NH<sub>2</sub> and H are placed in an outermost position of the surface and NH<sub>2</sub> acts as an H-bond acceptor; and (iii) Pos3, in which NH<sub>2</sub> and H are placed within a cavity and NH<sub>2</sub> acts as an H-bond donor. By proceeding this way, we roughly sample different sites of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** H-addition reactions on the crystalline (panels (a)–(c)) and amorphous (panels (d)–(f)) surfaces. PBE-D3(BJ)-computed reaction energies are also shown, in kJ mol⁻¹.</figcaption></figure>

the amorphous surface that can likely exhibit specific trends in the energy dissipation, as observed for the H₂ formation case (Pantaleone et al. 2021).

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₂⁻ 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.

# 4. Discussion

## 4.1. *Energy Dissipation, Chemical Desorption, and Diffusion*

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
<table><caption>**Table 2**
Data Obtained from the (NVE) AIMD Simulations for the Reactions Investigated</caption><thead><tr><th>Reaction</th><th><i>t</i><sub>ice</sub>/<i>t</i><sub>ice</sub> + <i>t</i><sub>S</sub></th><th><i>t</i><sub>S</sub>/<i>t</i><sub>ice</sub> + <i>t</i><sub>S</sub></th><th>Temperature Variation</th><th>BE</th><th>TKE</th><th>TKE-z</th></tr></thead><tbody><tr><td>N + H → NH</td><td>0.76</td><td>0.24</td><td>10 → 21.6</td><td>17.5</td><td>2.7</td><td>1.4</td></tr><tr><td>NH + H → NH<sub>2</sub></td><td>0.71</td><td>0.29</td><td>10 → 23</td><td>18.8</td><td>2.9</td><td>1.3</td></tr><tr><td>NH<sub>2</sub> + H → NH<sub>3</sub> (Pos1)</td><td>0.73</td><td>0.27</td><td>10 → 24.1</td><td>22.4</td><td>5.2</td><td>1.4</td></tr><tr><td>NH<sub>2</sub> + H → NH<sub>3</sub> (Pos2)</td><td>0.58</td><td>0.42</td><td>10 → 23.8</td><td>34.4</td><td>3.5</td><td>2.2</td></tr><tr><td>NH<sub>2</sub> + H → NH<sub>3</sub> (Pos3)</td><td>0.90</td><td>0.10</td><td>10 → 27.6</td><td>36.3</td><td>2.9</td><td>2.2</td></tr></tbody></table>

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

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results of the NVE AIMD simulations for the first H-addition N + H → NH: (a) initial structure of the simulation; (b) evolution of the kinetic energy of NH; (c) evolution of the N-H distance; (d) evolution of the kinetic energy of the ice; and (e) evolution of the temperature of the ice. Instantaneous and averaged values (in yellow) are reported.</figcaption></figure>

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 \times 10^{-11}$ cm² s⁻¹ or $3.3 \times 10^{-7}$ cm² s⁻¹ at 25 K (Ásgeirsson et al. 2017), depending on the presence/absence of

<caption>**Figure 5.** Results of the NVE AIMD simulations for the second H-addition NH + H → NH₂: (a) initial structure of the simulation; (b) evolution of the kinetic energy of NH₂; (c) evolution of the newly formed N-H bond distance; (d) evolution of the kinetic energy of the ice; and (e) evolution of the temperature of the ice. Instantaneous and averaged values (in yellow) are reported.</caption>

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 \AA^2$ (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
Pos1

<figure><img src="image_7.png" /></figure>

Pos2

<figure><img src="image_8.png" /></figure>

Pos3

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Starting positions of the reactive species for the NVE AIMD simulations relative to the third H-addition $NH_2 + H \rightarrow NH_3$.</figcaption></figure>

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_2$, and $NH_3$ species

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 7.** Results of the NVE AIMD simulations for the third H-addition $NH_2 + H \rightarrow NH_3$ starting from positions Pos1 and Pos3: evolution of the kinetic energy of $NH_3$ (panels (a) and (e)); evolution of the newly formed N-H bond distance (panels (b) and (f)); evolution of the kinetic energy of the ice (panels (c) and (g)); and evolution of the temperature of the ice (panels (d) and (h)). Instantaneous and averaged values (in yellow) are reported.</figcaption></figure>

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_3$ is a component of the interstellar ices: after the completion of the full H-addition sequence, $NH_3$ remains on the water ice surface, which, upon subsequent $H_2O$ 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.

## 4.2. Energy Dissipation Mechanisms

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

<figure><img src="image_18.png" /><figcaption>**Figure 8.** Difference of the simulated VDOS power spectrum (ΔVDOS) between the computed trajectories for the surfaces when the reaction takes place and for the isolated surface.</figcaption></figure>

libration modes (with a minor contribution to the bending modes) of the water ice.

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₂⁻ 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 NH₄⁺/OH⁻. Here the opposite reaction to form H₃O⁺/NH₂⁻ 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
<figure><img src="image_19.png" /><figcaption>**Figure 9.** Snapshots of the trajectories relative to the NH₃ formation (panel (a): Pos1; panel (b): Pos3), in which the transient H₃O⁺/NH₂ ion pair species forms. In panel (a), the ion pair is directly formed. In panel (b), a proton transfer relay mechanism takes place.</figcaption></figure>

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₂ 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₂⁻ but remains on H₃O⁺, which starts a proton transfer to the NH₂⁻ species (recovering NH₃) via a proton shuttle mechanism involving three water molecules (see Figure 9(b)).

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

<figure><img src="image_20.png" /><figcaption>**Figure 10.** Scheme of the ammonia formation in cold molecular clouds (adapted from De Simone et al. 2022; Tinacci et al. 2022). The major process forming ammonia is the hydrogenation of frozen N atoms into NH, followed by NH₂ and finally NH₃. Ammonia formed in this way, on the grain surfaces, remains frozen onto the grain mantles because its BE is too large for the thermal sublimation to be efficient. Only the fraction released by the CD would be gaseous ammonia.</figcaption></figure>

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

### 4.3. Astrophysical Implications

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).
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 hydro- genation 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.

5. Conclusions

In this paper, the formation of interstellar NH₃ through the successive hydrogenation of atomic N, i.e., N(⁴S) + H(²S) → NH(³Σ⁻) + H(²S) → NH₂(²B₁) + H(²S) → NH₃(¹A₁), on water ice surfaces has been investigated at a quantum mechanical DFT PBE-D3(BJ) level, paying special attention to the third-body effect exerted by the ice as energy sink and the mechanisms through which the reaction energies dissipate throughout the water ice surface. Static calculations indicate that the H-addition reactions leading to the final NH₃ are barrierless, on both crystalline and amorphous surfaces, and exhibit large and very negative reaction energies (between −340 and −440 kJ mol⁻¹). The alternative H-abstraction reactions, i.e., NH + H → N + H₂ and NH₂ + H → NH + H₂, do not compete with the NH₃ formation. AIMD simulations within the microcanonical NVE ensemble have been performed to investigate the fate of the energy released in each hydrogenation step and the mechanisms through which this energy flows among the different parts of the system. This has been done by monitoring the evolution over time of the most relevant energetic components: the potential energy and the kinetic energies of the newly formed species and the ice. AIMD simulations were centered only for the reactions occurring on the amorphous water ice surface model. Simulation results show that about 58%–90% of the energy liberated by the reaction is quickly absorbed (within 1 ps) by the ice, which, as a result, is heated up by about 12–18 K. Since the amount of TKE retained by the newly formed species is small compared to their binding energies, the products are doomed to remain adsorbed on the surface without being released to the gas phase, that is, CD is not expected. Moreover, we estimated that the velocity of the energy dissipation is much faster than the H surface diffusion, indicating that the NH₃ formation through the H-addition adopts a true three-step elementary mechanism, in which, for a given H-addition, the energy is dissipated before the next H-addition. According to these results, we can figure out that the formation of NH₃ takes place through a continuous H-addition until reaching the formation of the final product, NH₃, which remains on the water ice, becoming over time a bulky component of the interstellar ice, in agreement with observational evidence.

Finally, by computing the VDOS of the water molecules of the ice surface along the AIMD trajectories, we elucidated the mechanisms through which the energy transfer from the formed species to the ice surface takes place. We identified a general channel based on the vibrational coupling between the highly excited stretching and bending vibrational modes of the newly formed species and the low-energy libration modes of the water

ice surface. Additionally, a second and more exclusive mechanism (observed in the NH₃ formation trajectories but not in the NH and NH₂ ones) has been found. It is based on the formation of a transient H₃O⁺/NH₂ ion pair due to a proton transfer from NH₃ to a nearby icy H₂O molecule, which permits a faster energy dissipation, thereby relaxing more efficiently the newly born NH₃.

This project has received funding within the European Union's Horizon 2020 research and innovation program from the European Research Council (ERC) for the projects “The Dawn of Organic Chemistry” (DOC), grant agreement No. 741002 and “Quantum Chemistry on Interstellar Grains” (QUANTUMGRAIN), grant agreement No. 865657, and from the Marie Skłodowska-Curie grant for the project “Astro-Chemical Origins” (ACO), grant agreement No. 811312. S.F. wishes to thank Lorenzo Tinacci for the useful and inspiring discussions on the subject. Supplementary material consisting of (i) the energetics of the gas-phase reactions for both the H-additions and H-abstractions; (ii) the evolution with time of the total, potential, and kinetic energies of the studied processes; and (iii) results of the NVE AIMD simulations for the NH₃ formation from Pos2 is available in the AstroChemical Origins Zenodo community doi: 10.5281/zenodo.7115984 (Ferrero 2022).

ORCID iDs

Stefano Ferrero https://orcid.org/0000-0001-7819-7657

Stefano Pantaleone https://orcid.org/0000-0002-2457-1065

Cecilia Ceccarelli https://orcid.org/0000-0001-9664-6292

Piero Ugliengo https://orcid.org/0000-0001-8886-9832

Mariona Sodupe https://orcid.org/0000-0003-0276-0524

Albert Rimola https://orcid.org/0000-0002-9637-4554

References

Ásgeirsson, V., Jónsson, H., & Wikfeldt, K. T. 2017, *J. Phys. Chem. C*, **121**, 1648

Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015, *ARA&A*, **53**, 541

Brehm, M., & Kirchner, B. 2011, *J. Chem. Inf. Model.*, **51**, 2007

Brehm, M., Thomas, M., Gehrke, S., & Kirchner, B. 2020, *JChPh*, **152**, 164105

Brieuc, F., Bronstein, Y., Dammak, H., et al. 2016, *J. Chem. Theory Comput.*, **12**, 5688

Buchholz, M., Fallacara, E., Gottwald, F., et al. 2018, *CP*, **515**, 231

Bussi, G., Donadio, D., & Parrinello, M. 2007, *JChPh*, **126**, 014101

Casassa, S., Ugliengo, P., & Pisani, C. 1997, *JChPh*, **106**, 8030

Ceccarelli, C., Caselli, P., Fontani, F., et al. 2017, *ApJ*, **850**, 176

Ceccarelli, C., Codella, C., Balucani, N., et al. 2022, *arXiv:2206.13270*

Ceriotti, M., Fang, W., Kusalik, P. G., et al. 2016, *Chem. Rev.*, **116**, 7529

De Simone, M., Ceccarelli, C., Codella, C., et al. 2022, *ApJL*, **935**, L14

Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, *A&A*, **512**, A30

Enrique-Romero, J., Rimola, A., Ceccarelli, C., et al. 2019, *ESC*, **3**, 2158

Enrique-Romero, J., Rimola, A., Ceccarelli, C., et al. 2022, *ApJS*, **259**, 39

Ferrero, S. 2022, Supporting Information for “Where does the energy go during the interstellar NH3 formation on water ice? A computational Study”, Version 1.0, Zenodo, doi:10.5281/zenodo.7115984

Fraser, H. J., Collings, M. P., Dever, J. W., & McCoustra, M. R. 2004, *MNRAS*, **353**, 59

Fredon, A., & Cuppen, H. 2018, *PCCP*, **20**, 5569

Fredon, A., Lamberts, T., & Cuppen, H. 2017, *ApJ*, **849**, 125

Fuchs, G., Cuppen, H., Ioppolo, S., et al. 2009, *A&A*, **505**, 629

Garrod, R., Wakelam, V., & Herbst, E. 2007, *A&A*, **467**, 1103

Goedecker, S., Teter, M., & Hutter, J. 1996, *PhRvB*, **54**, 1703

Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, *JChPh*, **132**, 154104

Grimme, S., Ehrlich, S., & Goerigk, L. 2011, *JCoCh*, **32**, 1456

Gutiérrez-Quintanilla, A., Layssac, Y., Butscher, T., et al. 2021, *MNRAS*, **506**, 3734

Habershon, S., & Manolopoulos, D. E. 2009, *JChPh*, **131**, 244518
Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, ApJ, 757, 185

Hama, T., & Watanabe, N. 2013, ChRv, 113, 8783

Hattig, C., Klopper, W., Kohn, A., & Tew, D. P. 2012, ChRv, 112, 4

Henning, T. 2010, ARA&A, 48, 21

Hocuk, S., & Cazaux, S. 2015, A&A, 576, A49

Huppert, S., Plé, T., Bonella, S., Depondt, P., & Finocchi, F. 2022, Appl. Sci., 12, 4756

Jones, A., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Jones, A., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46

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

Kühne, T. D., Iannuzzi, M., Del Ben, M., et al. 2020, JChPh, 152, 194103

Lippert, G., Hutter, J., & Parrinello, M. 1997, MolPh, 92, 477

Mangaud, E., Huppert, S., Plé, T., et al. 2019, J. Chem. Theory Comput., 15, 2863

Markland, T. E., & Ceriotti, M. 2018, Nat. Rev. Chem., 2, 0109

Melani, G., Nagata, Y., & Saalfrank, P. 2021, PCCP, 23, 7714

Meyer, J., & Reuter, K. 2014, AngCh, 53, 4721

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016a, A&A, 585, A24

Minissale, M., Moudens, A., Baouche, S., Chaabouni, H., & Dulieu, F. 2016b, MNRAS, 458, 2953

Mukherjee, S., & Barbatti, M. 2022, J. Chem. Theory Comput., 18, 4109

Najibi, A., & Goerigk, L. 2018, J. Chem. Theory Comput., 14, 5725

Neese, F. 2004, JPCS, 65, 781

Neese, F., Wennmohs, F., Becker, U., & Riplinger, C. 2020, JChPh, 152, 224108

Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, ApJ, 897, 56

Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49

Perdew, J. P., Burke, K., & Ernzerhof, M. 1996, PRL, 77, 3865

Perrero, J., Enrique-Romero, J., Martínez-Bachs, B., et al. 2022, ESC, 6, 496

Peterson, K. A., Adler, T. B., & Werner, H.-J. 2008, JChPh, 128, 084102

Pezzella, M., & Meuwly, M. 2019, PCCP, 21, 6247

Potapov, A., Jäger, C., & Henning, T. 2020, PRL, 124, 221103

Qasim, D., Chuang, K.-J., Fedoseev, G., et al. 2018, A&A, 612, A83

Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, 572, A70

Rittmeyer, S. P., Bukas, V. J., & Reuter, K. 2018, AdPhX, 3, 1381574

Ruaud, M., Loison, J., Hickson, K., et al. 2015, MNRAS, 447, 4004

Shakouri, K., Behler, J., Meyer, J., & Kroes, G.-J. 2018, JPhCh, 122, 23470

Takahashi, J., Masuda, K., & Nagaoka, M. 1999, ApJ, 520, 724

Tinacci, L., Germain, A., Pantaleone, S., et al. 2022, ESC, 6, 1514

Upadhyay, M., & Meuwly, M. 2021, Front. Chem., 9, 2021

Valeev, E. F. 2004, CPL, 395, 190

Van Dishoeck, E. F., Herbst, E., & Neufeld, D. A. 2013, ChRv, 113, 9043

VandeVondele, J., & Hutter, J. 2007, JChPh, 127, 114105

Vasyunin, A., & Herbst, E. 2013, ApJ, 769, 34

Vidali, G. 2013, ChRv, 113, 8762

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

Watanabe, N., Kimura, Y., Kouchi, A., et al. 2010, ApJL, 714, L233

Watanabe, N., & Kouchi, A. 2002, ApJL, 571, L173

Watanabe, N., & Kouchi, A. 2008, PrSS, 83, 439

Watanabe, N., Nagaoka, A., Shiraki, T., & Kouchi, A. 2004, ApJ, 616, 638

Weigend, F., & Ahlrichs, R. 2005, PCCP, 7, 3297

Zamirri, L., Ugliengo, P., Ceccarelli, C., & Rimola, A. 2019, ESC, 3, 1499