# Binding Energies of N-bearing Species on Interstellar Water Ice Mantles by Quantum Chemical Calculations

Berta Martínez-Bachs, Stefano Ferrero, Cecilia Ceccarelli, Piero Ugliengo, Albert Rimola

► To cite this version:

Berta Martínez-Bachs, Stefano Ferrero, Cecilia Ceccarelli, Piero Ugliengo, Albert Rimola. Binding Energies of N-bearing Species on Interstellar Water Ice Mantles by Quantum Chemical Calculations. The Astrophysical Journal, 2024, 969, 10.3847/1538-4357/ad47b6. insu-04836847

HAL Id: insu-04836847

https://insu.hal.science/insu-04836847v1

Submitted on 14 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# Binding Energies of N-bearing Species on Interstellar Water Ice Mantles by Quantum Chemical Calculations

Berta Martínez-Bachs<fnref n="1" />, ∗, Stefano Ferrero<fnref n="1" />, †, Cecilia Ceccarelli<fnref n="2" />, ‡, Piero Ugliengo<fnref n="3" />, ‡, and Albert Rimola<fnref n="1" />, ‡

Received 2024 February 14; revised 2024 April 10; accepted 2024 May 3; published 2024 July 1

## Abstract

Of the about 300 gas-phase molecular species so far detected in the interstellar medium (ISM), mostly via observations of their rotational lines, around 40% contain nitrogen (N) atoms. Likewise, of the less than a dozen interstellar molecules, firmly or likely detected in the solid-state water-dominated icy matrix by means of infrared observations, two bear N. A crucial parameter that regulates whether a species is in the gas or adsorbed on the icy phase is their binding energy (BE) toward the icy grain. Therefore, an accurate quantification of the BE is of paramount importance to properly model the ISM chemistry through numerical models. However, very few BEs are available in the literature, either determined experimentally or theoretically. In the present study, we calculate the BEs of 21 among the most abundant interstellar N-bearing species. We adopted two structural water ice models, representing a crystalline and an amorphous surface, using a reliable cost-effective procedure based on the density functional theory. While on the crystalline surface model only one BE per species is obtained due to the high symmetry of the unit cell, on the amorphous model from 5 to 10 BEs are obtained, due to its richer surface morphological variety. Most of our computed BEs agree with available experimental and other computational values. Finally, we discuss how the newly computed BEs can help estimate which N-bearing species can be frozen at the water snow line and, therefore, incorporated in water-rich ice planetesimals.

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

The interest to understand the chemical evolution of the interstellar medium (ISM) started more than 80 yr ago, when the first diatomic species CH, CN, and CH⁺ were detected (Dunham 1937; Swings & Rosenfeld 1937; McKellar 1940; Douglas & Herzberg 1941). Currently, more than 290 gas-phase species have been detected in the ISM through their rotational and vibrational transitions (McGuire 2022; Endres et al. 2016). Of these, about 40% contain N-bearing atoms, making this class of interstellar detected molecules the fourth most numerous, after H-, C-, and O-bearing species. Obviously, this is a consequence of the relative nitrogen abundance in the ISM (the fifth after H, He, O, and C) and its chemical propensity to combine with other atoms. In addition, nitrogen is a key element for life, in which N-bearing compounds played a crucial role in chemical processes that could have eventually led to the emergence of life (e.g., synthesis of amino acids and nucleotides; Saladino et al. 2012a; Saitta & Saija 2014; Santalucia et al. 2022).

In molecular clouds, the cold and dense regions of the ISM, the gas (consisting of atoms and molecules) is mixed with solid-state dust grains. The latter have cores of either silicates or carbonaceous materials, which are covered by thick ice mantles, made predominantly of water (formed in situ), but also of other volatile species such as CO, CO₂, NH₃, CH₃OH, and CH₄. For this reason, these ices are commonly referred to as

“dirty ices” (Boogert et al. 2015). Regarding N-bearing species, in addition to ammonia, the likely detection of solid OCN⁻ and NH₄⁺ has also been reported (e.g., McClure et al. 2023).

Observational measurements based on infrared spectroscopy reveal that interstellar ices are likely to be in an amorphous state (Watanabe & Kouchi 2008; Oba et al. 2009; Boogert et al. 2015; McClure et al. 2023). The composition of the ices is the result of the adsorption of species coming from the external gas environment as well as of chemical reactions occurring on their surfaces. Thus, the binding energy (BE) of each species to the icy dust mantle and the dust temperature are the most important parameters that regulate whether the species is in the gas or solid phase.

Observational measurements are usually interpreted via astrochemical numerical models, which aim to reproduce and understand the chemical behavior of the ISM. In these models, various parameters need to be introduced as input data, those describing the physical status of the object to model and others describing the chemical processes (e.g., Ceccarelli 2023). BEs are crucial parameters, since they are directly associated with the ice-to-gas transition in desorption events. Additionally, other surface processes also depend on the BEs of adsorbed species like the diffusion, in which diffusion barriers are commonly assumed to be a fraction of its BE (Cuppen et al. 2017; Penteado et al. 2017). Therefore, having accurate BE values of interstellar species is of uttermost importance for the reliability of the astrochemical modeling outcomes and, hence, to properly know and understand the chemical evolution and composition of the ISM.

Experimentally, BEs are determined by means of temperature-programmed desorption (TPD) experiments. They

<figure><img src="image_3.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, 08193, Catalonia, Spain; albert.rimola@uab.cat</footnote>

<footnote n="2">²Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), 38000 Grenoble, France</footnote>

<footnote n="3">³Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Set of N-bearing species considered in this work. Green represents carbon atoms, blue nitrogen, red oxygen, and gray hydrogen.</figcaption></figure>

basically consist of three consecutive steps: (i) the substrate is carefully prepared by condensing water in a phase as close as possible to that observed in the ISM icy grains; (ii) the substrate is exposed to the adsorption of the species at a (very low) constant temperature, and (iii) the temperature below the substrate is increased through a defined ramp in time to measure the desorption rate of the species as a function of the temperature, usually detected with mass spectrometry (Noble et al. 2012; Penteado et al. 2017). Then, by a proper mathematical manipulation of the Polanyi–Wigner equation (Kolasinski 2012), which describes the desorption process, the BE and the desorption rate preexponential factor are worked out from the experimental TPD data. There is no unique way to derive the BE and the preexponential factor from the experiments, causing some degree of variability in these data and difficulties in making proper comparisons. To be noticed, TPD experiments do not measure the BEs of the species but its desorption enthalpy, which is usually equivalent to the BE, if no other activated processes take place (He et al. 2016). In general, BE values obtained through TPD experiments are strongly affected not only by the composition and morphology of the icy surface but also by the adsorbate coverage. Despite the great amount of work done investigating desorption processes via TPD experiments, BEs on amorphous water ices have been derived only for less than two dozen molecules (e.g., Dohnálek et al. 2001; Collings et al. 2004; Noble et al. 2012; Dulieu et al. 2013; Fayolle et al. 2016; He et al. 2016; Smith et al. 2016; Chaabouni et al. 2018; Ferrero et al. 2022; Minissale et al. 2022). In addition, TPD experiments cannot deal with radical species due to their high reactivity (Cuppen et al. 2017; Minissale et al. 2022).

BEs can also be obtained theoretically by means of quantum chemical calculations. In this way, many TPD limitations can be overcome, like dealing with open-shell species (e.g., radicals). Several computational studies have derived the BEs of species such as H, H₂, C, N, O, OH, CO, CO₂, HCO, H₂O, NH₃, and HF on periodic/cluster crystalline/amorphous model systems of different sizes (e.g., Al-Halabi & Van Dishoeck 2007;

Karssemeijer & Cuppen 2014; Karssemeijer et al. 2014; Ásgeirsson et al. 2017; Sameera et al. 2017; Senevirathne et al. 2017; Zamirri et al. 2017, 2019; Shimonishi et al. 2018; Bovolenta et al. 2020; Duflot et al. 2021; Germain et al. 2022; Tinacci et al. 2022, 2023; Hendrix et al. 2023). Other computational studies calculated the BEs of a large number of interstellar species, but adopting as a substrate very small ice models. For instance, in Wakelam et al. (2017), it was assumed a single water molecule is representative of the whole grain, despite some clever a posteriori corrections to overcome such a limitation. A different strategy was carried out by Das et al. (2018) and Hendrix et al. (2023), in which water ice surfaces were mimicked by clusters made up by up to four/six water molecules. However, none of these two surface approaches can fairly reproduce the amorphous nature of water icy dust grains. Moreover, these models are not large enough to ensure the occurrence of H-bond cooperativity, which affect the binding capacity of the surface adsorption sites. In contrast, Duflot et al. (2021) performed BE calculations by combining force field molecular dynamics (MD) simulations with quantum chemical calculations on some configurations emerged from the MD sampling for a set of atoms and molecules on hexagonal ice and amorphous cluster models. The use of the ONIOM procedure (by combining density functional theory, hereafter DFT, with semiempirical methods for the description of the binding site and environmental regions, respectively) allowed the authors to use cluster models of 120–150 atoms, this way overcoming the cluster size and the lack of H-bonding cooperativity limitations mentioned above.

In a recent work, Ferrero et al. (2020) computed a large set of BEs for 21 astrochemical-relevant species, both on a crystalline and on amorphous periodic water ice surface models. Both models properly reproduce the behavior of interstellar water ices and ensure the cooperativity of hydrogen-bonded interacting systems. For the crystalline model, since only a few possible adsorption sites exist, just one or two BEs were computed for each species. In contrast, on the amorphous water ice model, a broad distribution of BEs for each species was calculated because multiple adsorption sites were considered. In a more recent work, Ferrero et al. (2022) followed a methodology similar to that of Ferrero et al. (2020) to compute the BE distributions of 17 S-bearing species absorbed on the same crystalline and amorphous periodic water ice models.

The present work adopts the same approach, i.e., ice models and methodology of Ferrero et al. (2020) and Ferrero et al. (2022), to study the adsorption features of a set of 21 N-bearing species (summarized in Figure 1). The 21 chosen N-bearing species include 14 closed-shell species (NH₃, HCN, HNC, N₂, HNO, CH₃NH₂, CH₃CN, NH₂CN, HNCO, N₂O, NH₂CHO, HC₃N, CH₂CHCN, and HC₅N) and seven open-shell species (N, NH, NH₂, CN, NO, OCN, and NO₂). Six species (NH₃, HCN, N₂, CH₃CN, NH₂CHO, and NH₂) were already present in the Ferrero et al. (2020) work, and we repeated the computations of their BEs to ensure internal consistency, as the method applied here may slightly differ in some aspects (see Section 2). Moreover, for each species, the desorption rate preexponential factors were computed, as they are crucial to work out the BEs from the TPD experiments.

All the studied species have been detected in the ISM with one exception, NO₂. In particular, the sample consists of detected neutral N-bearing species not containing atoms
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top view (A) and lateral view (B) of the crystalline water surface model; (C) side view of the amorphous solid water surface. The black lines indicate the unit cell sizes.</figcaption></figure>

that belong to the third period (e.g., NS—computed in Ferrero et al. 2022, PN and SiN, etc.). With this rule, all the di-atomic and tri-atomic N-species, with one exception ($C_2N$), are considered in the present study, whereas we included half of the neutral N-bearing species with four atoms and a smaller fraction with a larger number of atoms. Generally speaking, we privileged the polyatomic species, which are more often detected in star-forming regions, which in turn also implies having the largest abundances.

The goal of the present study is to provide the community with accurate BE values and distributions of these species on water-rich ices so that they can serve as useful input data in astrochemical models, as well as the corresponding desorption rate preexponential factors. As already mentioned, the correct modeling of the species desorption is a basic need for rationalizing the observations to ultimately understand the ISM chemistry.

# 2. Methodology

## 2.1. Water Ice Surface Models

The surfaces mimicking the interstellar water ice mantles were modeled adopting a periodic approach. In particular, two model systems were adopted: a crystalline system and an amorphous system.

As a first approach, a crystalline water ice surface model based on the P-ice structure was used for the adsorption of the N-bearing compounds. P-ice is the proton-ordered analog to the hexagonal (proton-disordered) water ice, which has been demonstrated to reproduce fairly well the physicochemical features of crystalline water ice (Pisani et al. 1996). The surface model was generated by cutting out the P-ice 3D periodic bulk structure perpendicular to the [010] direction, resulting in the (010) slab surface model. This choice is because the (010) slab model is one of the most stable planes of the crystalline P-ice, and it has no dipole component

across the surface (Casassa et al. 1997; Zamirri et al. 2018). The use of P-ice offers an advantage due to its small unit cell, which efficiently treats the adsorption of small species on its surfaces. Furthermore, when dealing with molecules of large size, it is easy and computationally efficient to adopt the supercell strategy, as explained in the following paragraph. A non-proton-ordered model of I<sub>h</sub> ice would be, for the present purposes, too expensive computationally, as the unit cell should be large enough to represent enough local environment of each water molecules. For a crystalline ice, it is important that the long-range H-bonding cooperativity is taken into account while ensuring zero dipole moment across the considered slab models, both features exhibited by the P-ice model. This crystalline water surface (CWS) model consists of 12 atomic layers, with a whole thickness of 10.9 Å, and includes 24 water molecules (72 atoms) per unit cell (see Figures 2(A) and (B)). The cell parameters are |a|=9.066 Å and |b|=7.155 Å, resulting in a surface area of A = 64.863 Å<sup>2</sup>. For this CWS model, to simulate the adsorption of the large HC<sub>3</sub>N and CH<sub>2</sub>CHCN species, we adopted a 1 × 2 unit cell (namely, the b lattice parameter was doubled) giving rise to a surface that contains 48 water molecules (144 atoms) to avoid lateral interactions between unit cell replicas. For this 1 × 2 surface model, the cell parameters are |a|=9.066 Å and |b|=14.310 Å, resulting in a surface area of A = 129.726 Å<sup>2</sup>. Certainly, the CWS model does not represent the structural features of interstellar water ices since they are mostly in an amorphous form (Watanabe & Kouchi 2008; Boogert et al. 2015). However, this model has been very useful to assess the accuracy of the employed methodology due to the cheaper cost of the calculations compared to the amorphous model (see below).

With the aim to obtain BE values on more realistic interstellar surfaces, an amorphous solid water (ASW) surface model, first proposed by Ferrero et al. (2020), was adopted. It is a slab model based on a periodic approach, which was constructed by joining small water clusters obtained from Shimonishi et al. (2018). In this case, a large unit cell was considered (cell parameters of |a|=20.352 Å and |b|=10.027 Å, and a surface area of A = 204.071 Å<sup>2</sup>), consisting of 60 water molecules (180 atoms) per unit cell, in order to cover the heterogeneity of surface sites present in ASW systems. Additionally, the model presents different structural features for the upper and the lower surfaces, and therefore, different binding sites for adsorbates can be considered (see Figure 2(C)).

## 2.2. Computational Details

The periodic ab initio CRYSTAL17 code (Dovesi et al. 2018) was employed to carry out all the simulations. This software implements both the Hartree-Fock (HF) and the Kohn-Sham self-consistent field (SCF) methods to solve the electronic Schrödinger equation. CRYSTAL17 can simulate periodic (crystals, surfaces, and polymers) and nonperiodic (molecules) systems since it uses localized Gaussian functions as basis sets. As a result, surfaces can be defined as true slab models, hence avoiding the 3D replica usually implemented in plane waves' basis set codes.

The adsorption of the N-bearing species on the water ice models was computed adopting quantum chemical electronic structure methods. A DFT//HF-3c composite procedure was adopted to simulate the interaction between the adsorbates and the water ice slab models. In this cost-effective procedure, geometries were optimized with the HF-3c method while the
energetics were refined by performing single point calculations at the DFT level on the HF-3c optimized structures.

HF-3c is a semiempirical method in which the electronic energy is based on the HF method using a minimal Gaussian atomic orbital basis set, MINIX, in which three corrections are considered to alleviate the deficiencies caused by the approximations, namely, (i) to include dispersive interactions; (ii) to remove the basis set superposition error (BSSE) using the geometrical counterpoise (CP) correction; and (iii) to correct the short-ranged deficiencies due to using minimal basis sets (Jansen & Ros 1969; Liu & McLean 1973; Tatewaki & Huzinaga 1980; Sure & Grimme 2013).

The single point energy calculations at DFT were performed using two different methods: (i) for closed-shell species, the hybrid B3LYP functional, which includes a 20% of exact exchange, (Lee et al. 1988; Becke 1993), complemented by the D3 Grimme's empirical corrections with the Becke-Johnson (BJ) damping scheme, i.e., D3(BJ), to account for the dispersive forces (Grimme et al. 2011); and (ii) for open-shell species (atoms and radicals), the meta-hybrid M06-2X functional, which includes a 45% of exact exchange, in these cases adopting a spin-unrestricted formalism (Pople et al. 1995). For both DFT functionals, the Ahlrichs triple-zeta valence basis set with polarization functions was employed (Schäfer et al. 1992).

In order to assess the performance of the cost-effective DFT//HF-3c methodology, adsorptions were first simulated on the CWS model in a full DFT treatment (namely, geometry optimizations also at DFT, hereafter referred to as DFT//DFT). The same DFT functionals and basis set used in the DFT//HF-3c methodology were adopted in this DFT//DFT procedure.

Additionally, to ensure the accuracy of the DFT//HF-3c scheme (as well as the DFT//DFT one), an energy refinement adopting an ONIOM approach was carried out on the CWS systems, in line with other literature works citesameera2017oniom, duflot2021theoretical. In this ONIOM approach (inspired from the ONIOM2 strategy; Dapprich et al. 1999), the system is divided in two parts: (i) the model system, which here is a molecular system extracted from the crystalline periodic system that includes the adsorbate and the nearest water molecules (two or three) defining the binding site region; (ii) the real system, which is the whole system, treated as a periodic system. The two parts are computed at different theory levels: the model system at CCSD(T) as a high theory level, and the real system at DFT as a low theory level. This ONIOM approach is hereafter referred to as O[CCSD(T)//DFT]. The selection of these two regions is based on the localized nature of the interaction between the adsorbates and the binding site region. Calculations at the high and low theory levels were performed as single point energy calculations on the CWS DFT-optimized adsorption complex structures. The single point CCSD(T) energy calculations were performed with the Gaussian16 program (Frisch et al. 2016). Since the CCSD(T) method is strongly dependent on the basis set (Cramer & Bickelhaupt 2003), the Dunning correlation consistent family of basis sets (cc-pVnZ) was used to perform the calculations (Dunning 1989). Here, the Jun-cc-pVnZ family of basis set was employed, with n standing for double (“D”), triple (“T”), and quadruple (“Q”; Papajak et al. 2011), followed by extrapolation to the complete basis set by plotting the energy against 1/n³.

Periodic geometry optimizations were carried out by relaxing both the internal atomic positions and the cell parameters. All the adsorption complexes were characterized

by calculating their harmonic frequencies, which in CRYSTAL is done numerically at the Γ point by diagonalizing the mass-weighted Hessian matrix of the second-order energy derivatives with respect to atom displacements (central difference formula), in this case of ±0.003 Å from the minimum along each Cartesian coordinate (Pascale et al. 2004; Zicovich-Wilson et al. 2004). For the HF-3c-optimized structures, the whole periodic systems were considered in the harmonic frequency calculations. For geometries optimized at DFT, only the fragment regarding the binding region was considered (the same model system adopted in the ONIOM approach).

## 2.3. Calculation of the BEs

BEs were calculated from the interaction between the slab surface models and the adsorbed species. Interaction energies were calculated by considering the deformation of the surface and the adsorbates, and also the lateral interactions between adsorbates of adjacent unit cells. For localized basis functions as in this case, BSSE arises (Boys & Bernardi 1970), and, therefore, the CP correction was applied to the interaction energies (Davidson & Feller 1986).

BSSE-corrected interaction energies ($\Delta E^{CP}$) were calculated as

$$ \Delta E^{CP} = \Delta E^* + \delta E_S + \delta E_\mu + \delta E_L - \text{BSSE} \quad (1) $$

where $\Delta E^*$ stands for the deformation-free interaction energy, $\delta E_S$ for the deformation energy of the slab model, $\delta E_\mu$ for the deformation energy of the adsorbate, $\delta E_L$ for the adsorbate-adsorbate lateral interaction energy, and BSSE for the CP correction. Details on the calculation of each energetic term of Equation (1) can be found in the Appendix B.

BEs are defined as the opposite of BSSE-corrected interaction energies,

$$ BE = -\Delta E^{CP}. \quad (2) $$

Within the O[CCSD(T)//DFT] scheme, the refined BEs ($BE_{O[CCSD(T)//DFT]}$) were calculated as

$$ \begin{aligned} BE_{O[CCSD(T)//DFT]} &= BE_{(low;real)} + BE_{(high;model)} - BE_{(low;model)} \end{aligned} \quad (3) $$

where $BE_{(low;real)}$ stands for the BE of the whole periodic (real) system at the low theory level (DFT), and $BE_{(high;model)}$ and $BE_{(low;model)}$ stand for the BEs of the molecular (model) system at the high (CCSD(T)) and low (DFT) theory levels, respectively. The BSSE correction was applied to these three BE terms computed in the O[CCSD(T)//DFT] framework.

Vibrational zero-point energy (ZPE) corrections of each structure were obtained from their harmonic frequency calculations (see above). The contribution of the ZPE to the BEs ($\Delta ZPE$) was then added to BE to obtain the final BE(0), namely,

$$ BE(0) = BE - ZPE. \quad (4) $$

Within the DFT//HF-3c scheme, harmonic frequencies (computed at HF-3c) are not expected to be accurate enough, and therefore, they were not used to consider $\Delta ZPE$, but only to verify that the optimized structures were really minima stationary points. In contrast, within the DFT//DFT scheme, the harmonic frequencies of the fragment (calculated at DFT) were employed to compute the $\Delta ZPE$ terms. As this procedure was limited to the CWS model, a linear correlation between the
BE and the corresponding BE(0)s calculated at DFT//DFT allowed us to obtain a fitting scaling factor, which was then applied to the DFT//HF-3c BE values to obtain the corresponding BE(0) values for each adsorption complex.

## 2.4. Calculation of the Desorption Rate Preexponential Factor

As mentioned in the Introduction, in TPD experiments, the BE of a species on a given substrate is obtained by applying the Polanyi-Wigner equation (Kolasinski 2012), in which the desorption rate ($k_{des}(T)$) is given by

$$k_{des}(T) \approx \nu(T) e^{-BE(0)/(k_B T)} \quad (5)$$

where $\nu(T)$ is the desorption rate preexponential factor, and $k_B$ is the Boltzmann constant. Thus, accurate determination of BE(0)s relies on an accurate estimate of the $\nu(T)$ value, which considers the influence of entropy in the desorption process.

As in previous articles (Minissale et al. 2022; Tinacci et al. 2022), to calculate $\nu$, we here adopted the formula from Tait et al. (2005):

$$\nu(T) = \frac{k_B T}{h} q_{tr,2D}^{\ddagger} q_{rot,3D}^{\ddagger} \quad (6)$$

where $k_B$ is the Boltzmann constant, $h$ is the Planck constant, $q_{tr,2D}^{\ddagger}$ is the 2D-translational partition function of the adsorbate, excluding the contribution from the translational motion perpendicular to the surface, and $q_{rot,3D}^{\ddagger}$ represents the adsorbate rotational partition function.

The 2D translational partition function is given by

$$q_{tr,2D}^{\ddagger} = \frac{2\pi mk_B T}{h^2} A \quad (7)$$

where $m$ is the mass of the molecule, and $A$ is the surface area per adsorbed molecule as computed using the unit cell parameters for the CWS and ASW cases.

The 3D rotational partition function is given by

$$q_{rot,3D}^{\ddagger} = \frac{\sqrt{\pi}}{\sigma h^3} (8\pi^2 k_B T)^{\frac{3}{2}} \sqrt{I_y I_y I_z} \quad (8)$$

where $I_y$ is the $i$th adsorbate principal moment of inertia, and $\sigma$ is the symmetry adsorbate rotation factor.

For diatomic molecules and linear species, the rotational partition function is reduced to

$$q_{rot,2D}^{\ddagger} = \frac{\sqrt{\pi}}{\sigma h^2} (8\pi^2 k_B T) \sqrt{I_y I_z} \quad (9)$$

For atoms, only the translational partition function survives in the calculation of the desorption rate prefactor.

## 3. Results

This work presents the calculated BEs of the N-containing species shown in Figure 1, both on the CWS and ASW models. To this end, the considered species have been classified in two groups, according to their capability to establish hydrogen bonds (H-bonds): (i) those that can only act as H-bond acceptors (N, CN, N₂, NO, OCN, N₂O, and NO₂), hereafter referred to as Hb-A; and (ii) those that can act as both H-bond donors and acceptors (NH, NH₂, NH₃, HCN, HNC, HNO, CH₃NH₂, CH₃CN, NH₂CN, HNCO, NH₂CHO, HC₃N, CH₂CHCN, and HC₅N), hereafter referred to as Hb-AD. It is important to remark that, although CH₃CN, HC₃N, CH₂CHCN,

<table><caption>**Table 1**
BEs (in Kelvin Units) Computed at DFT//HF-3c, DFT//DFT, and O[CCSD(T)]//DFT Levels on the CWS Model</caption><thead><tr><th>Species</th><th>BE<sub>DFT//HF-3c</sub></th><th>BE<sub>DFT//DFT</sub></th><th>BE<sub>O[CCSD(T)]//DFT</sub></th></tr></thead><tbody><tr><td>N</td><td>1913</td><td>2225</td><td>1586</td></tr><tr><td>NH</td><td>4179</td><td>4643</td><td>4125</td></tr><tr><td>NH<sub>2</sub></td><td>6042</td><td>5939</td><td>5861</td></tr><tr><td>NH<sub>3</sub></td><td>7301</td><td>7373</td><td>7233</td></tr><tr><td>HCN</td><td>5725</td><td>5124</td><td>5642</td></tr><tr><td>HNC</td><td>8482</td><td>8756</td><td>8257</td></tr><tr><td>N<sub>2</sub></td><td>1455</td><td>1564</td><td>1477</td></tr><tr><td>NO</td><td>1369</td><td>1638</td><td>1695</td></tr><tr><td>HNO</td><td>4601</td><td>5150</td><td>4737</td></tr><tr><td>CH<sub>3</sub>NH<sub>2</sub></td><td>8498</td><td>8455</td><td>8167</td></tr><tr><td>CH<sub>3</sub>CN</td><td>7090</td><td>7551</td><td>7532</td></tr><tr><td>NH<sub>2</sub>CN</td><td>10,413</td><td>10,312</td><td>9856</td></tr><tr><td>HNCO</td><td>8232<sup>a</sup></td><td>8340</td><td>8150</td></tr><tr><td>N<sub>2</sub>O</td><td>3092</td><td>3288</td><td>3430</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>9634</td><td>10,079</td><td>9638</td></tr><tr><td>NO<sub>2</sub></td><td>2199</td><td>2681</td><td>2751</td></tr><tr><td>HC<sub>3</sub>N</td><td>3781</td><td>3886<sup>b</sup></td><td>...</td></tr><tr><td>CH<sub>2</sub>CHCN</td><td>5818</td><td>6199</td><td>6000</td></tr></tbody></table>

**Notes.**

<sup>a</sup> BEs have been computed using DFT//PBE procedure; see Section 3.1.1.

<sup>b</sup> BE has been interpolated from the correlation between binding energies at DFT//DFT level vs. DFT//HF-3c method; see Figure 4.

and HC₅N belong to the Hb-AD group, they establish weak donor H-bonds since hydrogen atoms are bonded to carbon atoms.

### 3.1. BEs for the Crystalline Systems

The BEs of all the N-containing species (with the exception of CN, OCN, and HC₅N) have been computed on the CWS model. Table 1 reports the computed BE values with the three different methodologies, namely, DFT//HF-3c, DFT//DFT, and O[CCSD(T)]//DFT (see Section 2.2). The adsorption of the CN, OCN, and HC₅N on the CWS model has not achieved for either self-consistent energy or geometry optimization convergence problems, probably due to the subtle electronic structure (in the cases of CN and OCN radicals) or difficulties in fitting the HC₅N structure onto the relatively small unit cell of the CWS model.

#### 3.1.1. BEs Computed at DFT//HF-3c

The adsorption on the CWS model was first studied at the DFT//HF-3c level. It is worth mentioning that HNCO is artificially described as a linear molecule by the HF-3c method. For this reason, and only for this specific molecule, we used the DFT Perdew-Burke-Ernzerhof (PBE; a gradient generalized approximated functional; Perdew et al. 1996) combined with the D3(BJ) term for the dispersion interactions, as an alternative method to accurately reproduce the structure of HNCO and its adsorption on the water ice surfaces. Accordingly, for the HNCO adsorption, the B3LYP-D3(BJ)//PBE-D3(BJ) procedure was performed. For the sake of simplicity, this alternative scheme is included in the reference group of DFT//HF-3c methodology.

As H-bond and dispersion forces are the main intermolecular interactions dictating the adsorption of the species on the water ice models, we manually guessed the initial structures of the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** HF-3c optimized geometries of NO₂ and NH₂CN (as illustrative examples) on the CWS model (A and B panels) and on the ASW (case 5, C and D panels). Distances are in Angstrom.</figcaption></figure>

adsorbate at the CWS model surface by maximizing the number of H-bond interactions.

As a result, $BE_{DFT//HF-3c}$ values spread in a range from 1369 to 10,413 K (see Table 1).

The general trend is that Hb-AD species present larger $BE_{DFT//HF-3c}$ values than Hb-A species (see Table 1). The $BE_{DFT//HF-3c}$ values for the Hb-AD cases span a higher range [3781–10,413 K] compared to Hb-A cases [1369–3092 K]. Hb-AD species can establish at least two H-bonds (one as H-bond donor and one as H-bond acceptor), and, even in some cases, three H-bonds. Moreover, Hb-AD group can also establish H-bonds reinforced by cooperativity. For Hb-A species, only one or, at the most, two H-bonds, are formed. Figures 3(A) and (B) report, as examples, the adsorption of NO₂ and NH₂CN (as illustrative examples of Hb-A and Hb-AD groups, respectively) on the CWS model optimized at HF-3c.

Dispersive interactions are more significant in species with larger sizes. As, in general, species belonging to Hb-AD group are larger in size than those of the Hb-A group, they give larger dispersion interactions than Hb-A species.

### 3.1.2. BEs Computed at DFT//DFT

To assess the performance of the cost-effective DFT//HF-3c methodology, the adsorption complexes on the CWS model have also been simulated at DFT//DFT, allowing the comparison between the $BE_{DFT//DFT}$ and the $BE_{DFT//HF-3c}$ values.

A general trend is that the $BE_{DFT//DFT}$ values are slightly higher than the $BE_{DFT//HF-3c}$, as shown in the plot of Figure 4(A). The HF-3c-optimized geometries, which are almost free from BSSE by construction, give longer H-bond distances than the DFT-optimized ones. This, in turn, gives lower-BE values for DFT//HF-3c than the DFT//DFT ones, in which the BSSE shrinks the intermolecular distances. Nevertheless, $BE_{DFT//DFT}$ well correlates with $BE_{DFT//HF-3c}$ values, as shown in Figure 4(A).

Remarkably, for HNCO, the DFT//PBE-D3(BJ) procedure is in perfect agreement with the DFT//DFT one, indicating its suitability to be used for this molecule in further calculations.

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Linear fits between (A) the BE values computed at DFT//DFT vs. those at DFT//HF-3c level, (B) the BE values computed at O[CCSD(T)/DFT] level vs. those at DFT//HF-3c level, and (C) the BE values computed at O[CCSD(T)/DFT] level vs. those at DFT//DFT level. The linear fit equations are shown in the graphs.</figcaption></figure>

### 3.1.3. BE Refinement at O[CCSD(T)//DFT]

As explained in Section 2.3, an ONIOM local approximation has been performed (O[CCSD(T)/DFT]) to assess the accuracy of the $BE_{DFT//DFT}$ and $BE_{DFT//HF-3c}$ values and, in the same time, to refine their values. This correction has been performed for all the species, with the exception of the HC₃N
<table><caption>**Table 2**
BE(0) Values (in Kelvin Units) Computed at DFT//HF-3c Level on the ASW Model</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">CWS</th><th colspan="10">ASW</th></tr><tr><th>DFT//HF-3c</th><th>DFT//DFT</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td>N</td><td>1625</td><td>1891</td><td>1744</td><td>1514</td><td>370</td><td>3194</td><td>1643</td><td>905</td><td>1419</td><td>657</td><td>...</td><td>...</td></tr><tr><td>NH</td><td>3551</td><td>3945</td><td>2728</td><td>1167</td><td>2135</td><td>3792</td><td>4692</td><td>1167</td><td>2558</td><td>1469</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub></td><td>5134</td><td>5047</td><td>6722</td><td>3281</td><td>3712</td><td>5567</td><td>5714</td><td>5714</td><td>3792</td><td>3162</td><td>...</td><td>...</td></tr><tr><td>NH<sub>3</sub></td><td>6204</td><td>6265</td><td>7511</td><td>4476</td><td>5795</td><td>5039</td><td>5590</td><td>5181</td><td>4293</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CN</td><td>...</td><td>...</td><td>6690</td><td>3677</td><td>2605</td><td>4171</td><td>7792</td><td>7471</td><td>3014</td><td>3431</td><td>...</td><td>...</td></tr><tr><td>HCN</td><td>4865</td><td>4354</td><td>2484</td><td>4354</td><td>3076</td><td>4364</td><td>2779</td><td>4977</td><td>5222</td><td>6306</td><td>...</td><td>...</td></tr><tr><td>HNC</td><td>7207</td><td>7440</td><td>6603</td><td>5527</td><td>4813</td><td>7932</td><td>10077</td><td>7200</td><td>4663</td><td>3820</td><td>...</td><td>...</td></tr><tr><td>N<sub>2</sub></td><td>1236</td><td>1329</td><td>1145</td><td>1451</td><td>1114</td><td>1012</td><td>756</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NO</td><td>1163</td><td>1392</td><td>965</td><td>629</td><td>512</td><td>291</td><td>931</td><td>945</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HNO</td><td>3909</td><td>4376</td><td>2346</td><td>2680</td><td>2214</td><td>2704</td><td>3144</td><td>3915</td><td>3329</td><td>2715</td><td>...</td><td>...</td></tr><tr><td>CH<sub>3</sub>NH<sub>2</sub></td><td>7221</td><td>7184</td><td>5760</td><td>5649</td><td>6041</td><td>6932</td><td>4373</td><td>8293</td><td>6147</td><td>5039</td><td>...</td><td>...</td></tr><tr><td>CH<sub>3</sub>CN</td><td>6025</td><td>6416</td><td>4920</td><td>6416</td><td>4328</td><td>7577</td><td>4956</td><td>7585</td><td>2623</td><td>3423</td><td>...</td><td>...</td></tr><tr><td>OCN</td><td>...</td><td>...</td><td>2635</td><td>1471</td><td>2216</td><td>1155</td><td>1816</td><td>1759</td><td>2370</td><td>1708</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CN</td><td>8848</td><td>8763</td><td>7076</td><td>6556</td><td>7931</td><td>10,707</td><td>2503</td><td>9955</td><td>10,704</td><td>6554</td><td>5230</td><td>7367</td></tr><tr><td>HNCO<sup>a</sup></td><td>6994</td><td>7087</td><td>5154</td><td>6932</td><td>4404</td><td>4386</td><td>7977</td><td>6485</td><td>4149</td><td>3390</td><td>...</td><td>...</td></tr><tr><td>N<sub>2</sub>O</td><td>2628</td><td>2794</td><td>2162</td><td>1349</td><td>1989</td><td>3042</td><td>3923</td><td>2389</td><td>3043</td><td>1376</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>8186</td><td>8564</td><td>10,904</td><td>6357</td><td>8789</td><td>5795</td><td>7194</td><td>6009</td><td>5764</td><td>6214</td><td>...</td><td>...</td></tr><tr><td>NO<sub>2</sub></td><td>1869</td><td>2278</td><td>1833</td><td>1286</td><td>1467</td><td>1244</td><td>658</td><td>1088</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>3</sub>N</td><td>3213</td><td>3287<sup>b</sup></td><td>4092</td><td>2688</td><td>2954</td><td>6635</td><td>4784</td><td>6466</td><td>3715</td><td>1889</td><td>...</td><td>...</td></tr><tr><td>CH<sub>2</sub>CHCN</td><td>4944</td><td>5268</td><td>4377</td><td>6411</td><td>4956</td><td>3978</td><td>6547</td><td>4956</td><td>2512</td><td>3900</td><td>...</td><td>...</td></tr><tr><td>HC<sub>5</sub>N</td><td>...</td><td>...</td><td>8044</td><td>605</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.**

BE(0) values at DFT//HF-3c and at DFT//DFT level for the CWS model are also included.

<sup>a</sup> The BE(0) values at DFT//HF-3c level have been computed using the DFT//PBE procedure; see Section 3.1.1.

<sup>b</sup> BE(0) value from the correlation between the BE(0) at DFT//DFT vs. DFT//HF-3c ones; see Figure 4.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Linear fit between between BE(0)s and BEs computed at DFT//DFT level for adsorption at the CWS ice model. The linear fit equation is shown in the graph.</figcaption></figure>

molecule, which gave energy convergence problems (vide supra).

Figures 4(B) and (C) show the correlation between the BE<sub>O[CCSD(T)//DFT]</sub> values against those computed at DFT//HF-3c and DFT//DFT, respectively. The linear fits show an overall perfect agreement. In conclusion, the DFT//DFT procedure is a reliable method to compute the BEs. In the same way, these results indicate that the cost-effective DFT//HF-3c procedure is a robust approach to obtain accurate BEs, which is relevant when dealing with the adsorption on ASW model, whose size prevent us from applying the DFT//DFT approach.

### 3.1.4. Zero-point Energy Correction

As mentioned in Section 2.3, for adsorption on the CWS model, we performed harmonic frequency calculations of a fragment consisting of the adsorbate and the binding site regions on the DFT-optimized structures. The resulting ZPE-corrections are used to correct the BE giving the final BE(0)s by means of Equation (4). The resulting BE(0) values are shown in Table 2, while Figure 5 shows the correlation between BE and BE(0) values. The linear fit is very good and shows that, to get BE(0), the raw BE values should be scaled by a factor of 0.85. The same scaling factor has been applied to correct the DFT//HF-3c BE values for the adsorption at the CWS model (see also Table 2) as well as to estimate the final BE(0) values for adsorption on the ASW model (see below).

### 3.2. BEs on the ASW

On the CWS model, the adsorption of the species was performed on only one (the unique) binding site due to the high symmetry of CWS structure. This is at a variance with the ASW model, which exhibits a variety of binding sites, therefore representing an ice model closer to the amorphous interstellar ice.

As already mentioned, for ASW, the BEs have been obtained at DFT//HF-3c level, and therefore, the final BE(0)s were computed by scaling the BE values by 0.85 (see Section 3.1.4), with the exception of the HNCO molecule, computed at the DFT//PBE level (see Section 3.1.1). Figures 3(C) and (D) show, as an example, the optimized geometries of adsorption at site 5 of the ASW for NO<sub>2</sub> (Hb-A group) and NH<sub>2</sub>CN (Hb-AD group), respectively.
<figure><img src="image_11.png" /><figcaption>**Figure 6.** BE(0)s for all considered species on CWS (filled circles) and on ASW (empty circles) models. For ASW, BE(0) values at DFT//HF-3c. Multiple values for each species are for different adsorption sites.</figcaption></figure>

On the ASW model, the adsorption of the species that were critical on the CWS model (i.e., CN, OCN, and HC₅N) converged smoothly.

For each of the considered N-containing species, at least five different BE values have been computed. For most of the cases, eight different BE values have been obtained, and for the particular case of NH₂CN, a total of 10 cases resulted. The exception is the HC₅N molecule (structure already problematic on the CWS model), for which only two BEs have been computed (which can be considered the upper and lower limits in its BE distribution). Table 2 shows the computed BE(0)s on the ASW model, to be compared with those from the CWS model.

Obviously, the BEs on the ASW model are different from those on the CWS model, as a result of a richer local environment for the ASW model (see Figure 6). This results in a wider range of BE values for each species on the ASW, compared to the CWS. Considering that a nonexhaustive unbiased search for all possible adsorbate positions at the ASW surface was not carried out, the present range of the BEs is somehow limited.

As a general trend, BEs computed on the CWS model are among the largest ones obtained for each species. Surface relaxation is more significant on the ASW model than on the CWS one due to the higher rigidity of the latter. The cost of surface relaxation decreases the final BE for ASW compared to CWS, in agreement with previous literature data (Ferrero et al. 2020; Perrero et al. 2022). Nevertheless, in most cases, the average BE value for the ASW is higher than the crystalline value, a trend that was already observed in Duflot et al. (2021).

Despite these differences, the general trend obtained in the BE values on the ASW is similar to that on the CWS: Hb-AD species develop larger BE values (between 605 and 10,904 K) than Hb-A species (between 291 and 4171 K). The reasons for these BE trends do not differ from those already raised for the CWS case (vide supra). The only exception to this trend is the CN molecule: despite belonging to the Hb-A group, three adsorption complexes exhibit higher BEs compared to the rest of the species of the same group (see entries 1, 5, and 6 of

Table 2). For these cases, CN, rather than establishing H-bonds with the surface, forms a hemibonded system with an outermost water of the ice model, establishing a non-classical two-center three-electron bond (hemibond) between the unpaired electron on the carbon atom of CN and the lone pair electrons of the oxygen atom of the water surface (see the optimized geometries for these cases in Appendix E). These hemibonded interactions between CN and H₂O were already found as a first step toward the formation of interstellar formamide on water ice surfaces (Rimola et al. 2018).

### 3.3. Desorption Rate Preexponential Factors

As detailed in Section 2.4, the desorption rate preexponential factors as a function of the temperature for each N-bearing species have been determined using Equation (5). The preexponential factor curves, as a function of T, are shown in Appendix F. The needed parameters entering Equation (5) are shown in Table 3. Moreover, when the monolayer peak desorption temperature ($T_{\text{peak}}$) from experimental measurements is available, we computed the desorption rate preexponential factors for both CWS and ASW models. The principal moments of inertia for the optimized structures (at the very same theory level as for the adsorption complexes) of each species have been computed with the MOLDRAW program (Ugliengo et al. 1993), while using the symmetry point group to assess the symmetry degeneracy. In the CWS, the area per adsorbate species is $6.49 \times 10^{-19} \text{ m}^2$, with the exception of HC₃N and CH₂CHCN species, for which the area is $1.30 \times 10^{-18} \text{ m}^2$ due to the adoption of the 1 × 2 unit cell. In the ASW, the area per adsorbate species was fixed at the value of $2.04 \times 10^{-18} \text{ m}^2$.

For the species N, NH₃, HCN, N₂, CH₃CN, and NH₂CHO, we have calculated the desorption rate preexponential factors for both models at the $T_{\text{peak}}$ values reported by Minissale et al. (2022). The obtained values range from $10^{13}$ to $10^{20}$, with variations depending on the specific species. In both CWS and ASW models, the preexponential factors differ by a factor of 3 due to the variation in the surface area A between the two surface models.

## 4. Discussion

### 4.1. Comparison of the BEs between N-bearing and N-free Species

In a previous work, Ferrero et al. (2020) computed the BEs of 21 species (15 of which not containing N atoms) absorbed on the same CWS and ASW models adopted here and using the same DFT//HF-3c procedure. The considered species span a wide range of chemical compounds, including simple molecules (H₂, N₂, O₂, CO, CO₂, OCS, HCl, HCN, CH₄, NH₃, H₂O, and H₂S), organic compounds (HCONH₂, H₂CO, HCOOH, CH₃OH, and CH₃CN), and four radical species (OH, NH₂, HCO, and CH₃). As a whole, more than half of them were oxygen-bearing (O-bearing) species. Likewise, Perrero et al. (2022) computed a distribution of BEs using the same icy models and the DFT//HF-3c approach for 17 S-bearing species, of which only one contains N-atoms (NS). The set of S-bearing species included eight closed-shell species (H₂S, H₂S₂, CS, CH₃SH, SO₂, OCS, H₂CS, and C₃), and nine radicals (NS, HS, HS₂, HCS, SO, S₂, C₄S, C₂S, and S). As we have adopted the same icy models and computational methods,
<table><caption>**Table 3**

Parameters Used for the Calculation of the Desorption Rate Preexponential Factors ($\nu$, in per second) for Considered Species Adsorbed at CWS and ASW Ice Models: Mass *m* (in Atomic Mas Units), Symmetry Number (*σ*), Principal Moments of Inertia (Ix, Iy, and Iz, in Amu Å²), and Experimental Monolayer Peak Desorption Temperatures (*T*<sub>peak</sub>, in Kelvins) from Minissale et al. (2022) for Submonolayer Regimes on Compact ASW</caption><thead><tr><th>Species</th><th>m</th><th>σ</th><th>I<sub>x</sub></th><th>I<sub>y</sub></th><th>I<sub>z</sub></th><th>T<sub>peak</sub></th><th>ν<sub>CSW</sub></th><th>ν<sub>ASW</sub></th><th>ν<sub>Minissale</sub></th></tr></thead><tbody><tr><td>N</td><td>14.01</td><td>1</td><td>...</td><td>...</td><td>...</td><td>35</td><td>7.61 × 10<sup>13</sup></td><td>2.39 × 10<sup>14</sup></td><td>1.17 × 10<sup>13</sup></td></tr><tr><td>NH</td><td>15.02</td><td>1</td><td>...</td><td>1.63</td><td>1.63</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub></td><td>16.02</td><td>2</td><td>1.09</td><td>2.25</td><td>3.34</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>3</sub></td><td>17.03</td><td>3</td><td>2.86</td><td>2.86</td><td>4.50</td><td>105</td><td>2.69 × 10<sup>16</sup></td><td>8.47 × 10<sup>16</sup></td><td>1.94 × 10<sup>15</sup></td></tr><tr><td>CN</td><td>26.02</td><td>1</td><td>...</td><td>17.34</td><td>17.34</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HCN</td><td>27.03</td><td>1</td><td>...</td><td>17.97</td><td>17.97</td><td>137</td><td>4.05 × 10<sup>17</sup></td><td>1.27 × 10<sup>18</sup></td><td>1.63 × 10<sup>17</sup></td></tr><tr><td>HNC</td><td>27.03</td><td>1</td><td>...</td><td>39.44</td><td>39.44</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>N<sub>2</sub></td><td>28.01</td><td>2</td><td>...</td><td>14.12</td><td>14.12</td><td>35</td><td>2.75 × 10<sup>15</sup></td><td>8.65 × 10<sup>15</sup></td><td>4.51 × 10<sup>14</sup></td></tr><tr><td>NO</td><td>30.01</td><td>1</td><td>...</td><td>12.40</td><td>12.40</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HNO</td><td>31.04</td><td>1</td><td>0.95</td><td>18.87</td><td>19.81</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CH<sub>3</sub>NH<sub>2</sub></td><td>31.06</td><td>1</td><td>8.30</td><td>36.39</td><td>37.90</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CH<sub>3</sub>CN</td><td>41.05</td><td>3</td><td>5.30</td><td>84.06</td><td>84.06</td><td>120</td><td>3.30 × 10<sup>18</sup></td><td>1.04 × 10<sup>19</sup></td><td>2.37 × 10<sup>17</sup></td></tr><tr><td>OCN</td><td>42.02</td><td>1</td><td>...</td><td>71.68</td><td>71.68</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CN</td><td>42.04</td><td>3</td><td>3.06</td><td>77.95</td><td>79.41</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HNCO</td><td>43.03</td><td>1</td><td>1.00</td><td>73.91</td><td>74.91</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>N<sub>2</sub>O</td><td>44.01</td><td>1</td><td>...</td><td>66.32</td><td>66.32</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>45.04</td><td>3</td><td>11.17</td><td>66.22</td><td>77.39</td><td>176</td><td>5.13 × 10<sup>19</sup></td><td>1.61 × 10<sup>20</sup></td><td>3.69 × 10<sup>18</sup></td></tr><tr><td>NO<sub>2</sub></td><td>46.01</td><td>2</td><td>5.49</td><td>54.09</td><td>59.58</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>3</sub>N</td><td>51.05</td><td>1</td><td>...</td><td>161.36</td><td>161.36</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CH<sub>2</sub>CHCN</td><td>53.06</td><td>1</td><td>15.88</td><td>150.35</td><td>166.22</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>5</sub>N</td><td>75.07</td><td>1</td><td>...</td><td>546.48</td><td>546.48</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Note. The desorption rate preexponential factors calculated in this study for both icy models are compared with those by Minissale et al. (2022).

we can consistently compare the BE distributions for oxygen-, sulfur-, and nitrogen-bearing molecules.

Figure 7 shows the distributions of the BE(0)s from Ferrero et al. (2020) and Perrero et al. (2022), which are compared with those from this work. To avoid bias in the comparison, the BEs of the N-bearing species reported by Ferrero et al. were excluded in the Ferrero’s BE(0) plot.

The distribution of BEs provided by Ferrero et al. (2020) exhibits a bimodal character reminiscent of two distinct adsorption modes: (i) the group with low BEs, which includes simple molecules characterized by weak dispersive interactions; and (ii) the group with high BEs, which includes large size species involving both dispersive and strong H-bond interactions. In contrast, the distribution of BEs computed by Perrero et al. (2022) for the S-bearing species exhibits a unimodal-like behavior, with a relatively narrow range of BEs, due to a dominance of the dispersive interactions.

The distribution of BEs for N-bearing species spreads over a broad range of values, with a rather flat distribution, showing H-bonds as well as dispersive interactions as both relevant. Differences in the BEs are due to the different number and strength of H-bond interactions with the ice surface, resulting in a broad BE distribution.

A more detailed comparison has been performed for four N-bearing species (CN, OCN, NH, and NH<sub>2</sub>) with their O-bearing (CO, CO<sub>2</sub>, OH, and H<sub>2</sub>O) and S-bearing (CS, OCS, SH, and H<sub>2</sub>S) analogs. It is worth noting that all these N-bearing species are radicals while, in contrast, the oxygen and sulfur ones are closed-shell species, with the exception of OH and SH (analogs for NH), which are also radicals. Table 4 shows the calculated BE(0)s for these species, while Figure 8 represents the range of BE(0)s calculated for these species.

As a general trend, the BEs of the selected N-bearing species are similar to those for O- and S-containing analogs, as the nature of the interaction is similar. However, despite this

general trend, some discrepancies emerge. H<sub>2</sub>S molecule exhibits lower BEs compared to H<sub>2</sub>O and NH<sub>2</sub> (the O- and N-analogs). This is because, in H<sub>2</sub>S, dispersion interactions are the main driving forces in the adsorption (with weaker H-bonding, Perrero et al. 2022) while, in contrast, H<sub>2</sub>O and NH<sub>2</sub> establish strong H-bonds with the icy surface, in addition to dispersion.

Regarding the CN species and its O- and N-analogs (CO and CS), the former species presents significantly larger BEs than CO and CS. This difference can be attributed to the formation, in three cases, of hemibonded systems between the CN species and icy water molecules, leading to increased BEs (as mentioned in Section 3.2). Beside the hemibonded complexes, the set of CN BEs on the ASW follows a BE distribution similar to that of CO and CS.

The above considerations raise a warning: it is always risky to guess a BE of a species by using the computed BE values of analogous species (CS versus CO as an example), as it is difficult to extrapolate dispersion and H-bonding contributions from the values of analogous species to guess those for the unknown species.

## 4.2. Comparison of the BEs with Literature Values

The computed BE(0)s are compared with available data from the literature, including experimental and computational studies. Table 5 and Figure 9 summarize the available data. Please note that we considered here the compilations where more than one BE is reported, namely, the work by Das et al. (2018), Bovolenta et al. (2022), Minissale et al. (2022, 2022), and Wakelam et al. (2017).

### 4.2.1. Comparison with Experimental Works

Let us first compare our computed BE(0)s with those obtained experimentally. It is worth reminding that, while
<table><caption>**Table 4**
BE(0)s (in Kelvin Units) for Selected Di-atomic and Tri-atomic N-bearing Species and Their Oxygen- and Sulfur-bearing Analogs (Ferrero et al. 2020; Perrero et al. 2022; See Text), Computed at DFT//HF-3c</caption><thead><tr><th colspan="4">N-bearing Species</th><th colspan="4">O-bearing Species</th><th colspan="4">S-bearing Species</th></tr><tr><th>Species</th><th>CWS</th><th>ASW<sub>min</sub></th><th>ASW<sub>max</sub></th><th>Species</th><th>CWS</th><th>ASW<sub>min</sub></th><th>ASW<sub>max</sub></th><th>Species</th><th>CWS</th><th>ASW<sub>min</sub></th><th>ASW<sub>max</sub></th></tr></thead><tbody><tr><td>CN</td><td>...</td><td>2605</td><td>7792</td><td>CO</td><td>1663</td><td>1109</td><td>1869</td><td>CS</td><td>2453</td><td>995</td><td>3680</td></tr><tr><td>OCN</td><td>...</td><td>1155</td><td>2635</td><td>CO<sub>2</sub></td><td>2568</td><td>1489</td><td>2948</td><td>OCS</td><td>2405</td><td>1286</td><td>2861</td></tr><tr><td>NH</td><td>3551</td><td>1167</td><td>4692</td><td>OH</td><td>5588</td><td>1551</td><td>5321</td><td>SH</td><td>2799</td><td>1078</td><td>4313</td></tr><tr><td>NH<sub>2</sub></td><td>5134</td><td>3162</td><td>6722</td><td>H<sub>2</sub>O</td><td>7200</td><td>3605</td><td>6111</td><td>H<sub>2</sub>S</td><td>4033</td><td>1970</td><td>4489</td></tr></tbody></table>

Note. For the ASW model, only the minimum and maximum BE(0)s are reported.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Distribution of the BE(0)s computed on the ASW model using the DFT//HF-3c methodology for N-bearing species (blue), the non-N-bearing species considered in Ferrero et al. (2020; red), and the S-bearing in Perrero et al. (2022; green).</figcaption></figure>

computational calculations directly provide BE(0) values, TPD experiments estimate desorption energies, which are derived from the desorption rate preexponential factor curves. Among

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Comparison of BE(0)s for selected di-atomic and tri-atomic N-bearing species and their oxygen and sulfur analogs (Ferrero et al. 2020; Perrero et al. 2022; see text). Solid lines cover the ranges between the minimum and maximum BE(0)s computed on the ASW at DFT//HF-3c for N- (blue), O- (red), and for S- (green) bearing species. The range of BE(0)s for the CN species is divided between cases in which a hemibond is formed (top range) and those in which a hemibond is not formed (bottom range). For reference, also the BE(0)s on the CWS computed at DFT//HF-3c level are reported (filled black circles).</figcaption></figure>

other problems, the difference arises because desorption energies are surface coverage dependent, and hence, the comparison between computational BE(0)s and desorption energies is not totally straightforward (e.g., Ferrero et al. 2022; Minissale et al. 2022). In our case, we always simulate a zero coverage limit, at least for the ASW model, while higher coverages are implied by the smaller unit cell of the CWS model.

Considering these caveats, we compare our BE(0)s with those provided in the compilation taken from Minissale et al. (2022), where the authors report a list of BE(0)s and preexponential factors extracted by experimental works and carefully reviewed for 21 relevant species. Of these, only N, NH<sub>3</sub>, HCN, N<sub>2</sub>, CH<sub>3</sub>CN, and NH<sub>2</sub>CHO are comparable with the present results. As already pointed out, we focus on the recommended BEs in the submonolayer regimes on compact ASW. As one can see from Table 5, all the Minissale's recommended BE(0)s fall within the range of our computed BE(0)s.

### 4.2.2. Comparison with Previous Theoretical Works

The easiest and cheapest method for computing the BEs was done in Wakelam et al. (2017). They assumed that a real BE value can be worked out assuming a single water molecule as an icy model. Neither ZPE nor BSSE corrections were added to the BEs. To compensate for the limitation of using a single water molecule to mimic the ice, a scaling factor was derived
<table><caption>**Table 5**
Summary of the Binding Energies (in Kelvins) Computed in the Present Work and Those from Literature Compilations</caption><thead><tr><th rowspan="2">Species</th><th colspan="3">BE(0)s from This Work</th><th colspan="5">Binding Energies from Literature</th></tr><tr><th>CWS</th><th>ASW<sub>min</sub></th><th>ASW<sub>max</sub></th><th>BE(0)<sub>Duflot</sub><sup>c</sup></th><th>BE(0)<sub>Bovolenta</sub><sup>d</sup></th><th>BE(0)<sub>Minissale</sub><sup>e</sup></th><th>BE<sub>Das</sub><sup>f</sup></th><th>BE<sub>Wakelam</sub><sup>g</sup></th></tr></thead><tbody><tr><td>N</td><td>1625</td><td>370</td><td>3194</td><td>⟨359⟩<sub>cry</sub> / ⟨580⟩<sub>am</sub></td><td>...</td><td>806</td><td>619</td><td>720</td></tr><tr><td>NH</td><td>3551</td><td>1167</td><td>4692</td><td>⟨2715⟩<sub>cry</sub> / ⟨2437⟩<sub>am</sub></td><td>...</td><td>...</td><td>...</td><td>2600</td></tr><tr><td>NH<sub>2</sub></td><td>5134</td><td>3162</td><td>6722</td><td>...</td><td>3488 ± 265</td><td>...</td><td>3240</td><td>3200</td></tr><tr><td>NH<sub>3</sub></td><td>6204</td><td>4293</td><td>7511</td><td>⟨4468⟩<sub>cry</sub> / ⟨4305⟩<sub>am</sub></td><td>3347 ± 219 // 1104 ± 197</td><td>5362</td><td>5163</td><td>5600</td></tr><tr><td>CN</td><td>...</td><td>2605</td><td>7792</td><td>...</td><td>...</td><td>...</td><td>1736</td><td>2800</td></tr><tr><td>HCN</td><td>4865</td><td>2484</td><td>6306</td><td>...</td><td>2425 ± 483 // 1899 ± 167</td><td>5344</td><td>2352</td><td>3500</td></tr><tr><td>HNC</td><td>7207</td><td>3820</td><td>10077</td><td>...</td><td>4628 ± 706 // 2552 ± 117</td><td>...</td><td>5225</td><td>3800</td></tr><tr><td>N<sub>2</sub></td><td>1236</td><td>756</td><td>1451</td><td>...</td><td>637 ± 114<sup>b</sup></td><td>1074</td><td>1161</td><td>1100</td></tr><tr><td>NO</td><td>1163</td><td>291</td><td>965</td><td>...</td><td>...</td><td>...</td><td>1988</td><td>1600</td></tr><tr><td>HNO</td><td>3909</td><td>2214</td><td>3915</td><td>...</td><td>...</td><td>...</td><td>2988</td><td>3000</td></tr><tr><td>CH<sub>3</sub>NH<sub>2</sub></td><td>7221</td><td>4373</td><td>8293</td><td>...</td><td>...</td><td>...</td><td>4434</td><td>...</td></tr><tr><td>CH<sub>3</sub>CN</td><td>6025</td><td>2623</td><td>7585</td><td>...</td><td>...</td><td>6253</td><td>3786</td><td>4680</td></tr><tr><td>OCN</td><td>...</td><td>1155</td><td>2635</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CN</td><td>8848</td><td>2503</td><td>10707</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HNCO</td><td>6994<sup>a</sup></td><td>3390<sup>a</sup></td><td>7977<sup>a</sup></td><td>...</td><td>...</td><td>...</td><td>5554</td><td>4400</td></tr><tr><td>N<sub>2</sub>O</td><td>2628</td><td>1349</td><td>3923</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>8186</td><td>5764</td><td>10904</td><td>...</td><td>...</td><td>9561</td><td>...</td><td>6300</td></tr><tr><td>NO<sub>2</sub></td><td>1869</td><td>658</td><td>1833</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>3</sub>N</td><td>3213</td><td>1889</td><td>6636</td><td>...</td><td>...</td><td>...</td><td>2925</td><td>...</td></tr><tr><td>CH<sub>2</sub>CHCN</td><td>4944</td><td>2512</td><td>6547</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>5</sub>N</td><td>...</td><td>605</td><td>8044</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** Second column reports the BE(0)s computed in this work for the CWS model at DFT//HF-3c level. Third and fourth columns report the minimum and maximum BE(0)s computed in this work for the ASW model. Fifth to eighth columns report the binding energies from the literature, as indicated in the header.

<sup>a</sup> BE(0)s have been computed using the DFT//PBE method.

<sup>b</sup> It does not include ZPE correction.

<sup>c</sup> Duflot et al. (2021), calculated average values on crystalline (cry, left) and amorphous (am, right) surface models.

<sup>d</sup> Bovolenta et al. (2022), two values indicate bimodal character in the BE(0) distribution.

<sup>e</sup> Minissale et al. (2022).

<sup>f</sup> Das et al. (2018).

<sup>g</sup> Wakelam et al. (2017).

from fitting the experimental BEs versus the computed BE ones for 16 species. In general, the final BEs by Wakelam et al. lay within the range of our computed BEs, although in the lower limit of our range (see Figure 9). The adoption of a single water molecule as an ice surface model cannot account for the double behavior of H-bond donor and H-bond acceptor of some species, neglecting the H-bond cooperativity.

Das et al. (2018) presented BEs by computing interaction energies of species on icy cluster models of four or six water molecules. As for the Wakelam et al. study, neither ZPE nor BSSE corrections were considered. Overall, the Das et al. BEs are within the range of the BEs computed in the present work, with similar considerations for the Wakelam et al. case. We believe these clusters are still too small to allow H-bonding cooperativity to be at its best.

Duflot et al. (2021) computed a set of BEs for different astrochemical species, among them N, NH, and NH<sub>3</sub> ones. These authors combined force field MD simulations with quantum chemical calculations for various snapshots taken along the MD trajectory. With the ONIOM procedure (combining DFT with semiempirical methods), the authors computed BE(0)s on crystalline and amorphous cluster models comprising 120–150 atoms. The averaged set of calculated BE values, are also reported in Table 5. For a proper comparison with the present data, we averaged the BE(0)s for the N, NH, and NH<sub>3</sub> species, getting ⟨1431⟩, ⟨2464⟩, and ⟨5412⟩ K, showing the Duflot’s BEs to differ by some amount with the present ones. Although the cluster models adopted by Duflot et al. (2021) are

more realistic and bigger than the minimal ones used in Wakelam et al. (2017) and Das et al. (2018), the same trend (although less pronounced) is observed: the Duflot’s values lay in the lower part of our range. The reasons of these differences are difficult to assess, as different surface models (cluster versus periodic), quantum chemical methods, and computational strategies (ONIOM approach versus DFT//HF-3c) were adopted. Furthermore, as the comparison is done for only three species, it is difficult to arrive at definite conclusions.

Bovolenta et al. (2022) adopted a different approach by presenting BE(0) distributions consisting of 225–250 values for 21 interstellar molecules on different amorphous cluster models made up by 22 water molecules. A three-step computational procedure was employed, involving sampling, geometry optimization at HF-3c level, and BE(0) calculation at DFT level of theory. Average values and standard deviations were obtained from a Gaussian fit of the set of BE(0)s. If a bimodal character was present in the distribution of BE(0)s, two Gaussian fits were used. Only five of the Bovolenta et al. species (that is, N<sub>2</sub>, HCN, NH<sub>2</sub>, NH<sub>3</sub>, and HNC) coincide with our set of compounds (see Table 5), making it hard to assess a general trend. Additionally, it is noteworthy that, for the purpose of comparison, we are providing the minimum/ maximum of each BE(0) value, while Bovolenta et al. provides average values and standard deviations. For species presenting low BE(0)s (namely, N<sub>2</sub>, HCN, NH<sub>2</sub>), our results align very well with the Bovolenta’s values. However, discrepancies arise for species presenting large BE(0)s (namely, NH<sub>3</sub> and HNC), in
<figure><img src="image_14.png" /></figure>

**Figure 9.** Comparison of the binding energies computed in this work and those from literature compilations. The solid gray lines show the range between the minimum and maximum BE(0)s computed on the ASW model while the filled purple circles are for the BE(0)s computed on the CWS model at DFT//HF-3c level. The light blue crosses correspond to the mean values from the Gaussian fit of BE(0)s distributions reported by Bovolenta et al. (2022), showing also the standard deviations. The pink lines correspond to the peaks of Maxwell-Boltzmann distribution functions of BE(0)s distribution reported by Tinacci et al. (2023), along with the standard deviations. For the other cases, only one value has been given: BE(0)s reported by Minissale et al. (2022; red diamonds), BEs reported by Das et al. (2018; black squares), and BEs reported by Wakelam et al. (2017; green triangles).

which the Bovolenta's BEs fall below the range of our computed BE(0)s. This could be due to the large size and amorphous nature of our ASW model, showing richer morphological features (including cavities), compared to the compact nature of the cluster models adopted by Bovolenta et al., thereby affecting the limit values of the BE distributions. As already pointed out in Bovolenta et al., adsorption on sites presenting a cavity will give larger BE(0)s compared to those for flat adsorption sites. An additional factor that can affect the comparison is the effect of the ZPE correction, which decreases the BE values by about 25% in the Bovolenta et al. work, at a variance with the 15% as a maximum in our case.

Within a context of prebiotic chemistry, Tóth Ugyonka et al. (2024) recently computed the adsorption of HCN on low density amorphous water ice model by means of grand canonical Monte Carlo simulations based on classical force fields at the temperatures of 50, 100, and 200 K. Their heat of adsorption computed at an infinitely low surface coverage were $-49.4 \pm 3.9 \text{ kJ mol}^{-1}$ (BE(0) of 5941 ± 469 K). This range is well within our computed one for HCN, despite the differences between the surface models and methodological aspects.

Finally, Tinacci et al. (2022) reported a study of the NH₃ BE(0)s obtained on an amorphous water ice model (roughly spherical) composed by 200 water molecules, the largest so far reported in the literature (Germain et al. 2022). These authors characterized 77 unique NH₃ structures on the icy grain, and, consequently, a BE(0) distribution was derived covering a wide range of values, from 1500 to 6000 K. Remarkably, the NH₃ BE(0) distribution exhibited a bimodal shape, with a first peak at around 4000 K and the second one, unexpected and never seen before by other studies, at 1800 K. As in previous works (Ferrero et al. 2020; Duflot et al. 2021), the range of NH₃ BEs

(0) computed here, 4293–7511 K, only covers the highest values of the BE distribution by Tinacci et al., likely because the ASW model is too small to have adsorbing sites with such a low BE (see also the discussion in Tinacci et al. 2022).

In a subsequent work, Tinacci et al. (2023) computed the H₂O BE(0) distribution on the same 200 water molecule grain model. Contrarily to NH₃, the water BE(0)s distribution exhibits a single peak at 4230 K and covers the range 1700–7400 K, while Ferrero et al. (2020), with the same ice model employed here, estimated a range, 3600–6111 K, rather shifted toward the high BEs. Also, in this case, therefore, the Ferrero et al. model missed sites with small BE values, even though the disagreement is much less important than that with ammonia, where the whole low-BE tail is missed. Therefore, the case of ammonia should not be generalized, and not all species have a low-BE second peak. The above considerations provide an important take-home message: the size of the icy model matters, especially when providing adsorption sites characterized by small BE values. We believe this is because large water clusters tend to maximize the number of internal H-bonds, leaving less available strong adsorption sites for the adsorbates. Small icy models, on the contrary, show too many “unsatisfied” water molecules at their surface, pushing up (higher-BE(0) values) their affinity toward adsorbates.

### 4.3. Comparison of the Desorption Rate Preexponential Factors

Our computed desorption rate preexponential factors are compared with those reported in the literature by Minissale et al. (2022), in Table 3.

In Minissale et al. (2022), the reported preexponential factors were computed using the Tait et al. (2005) formula, as in this work. However, some parameters are different. These are (i) the area per absorbed site, and (ii) the inertia moments. Regarding the area per adsorbed site, those for the CWS and ASW models were calculated according to their unit cell parameters, while in Minissale et al. a fixed value of $10^{-19} \text{ m}^2$ was used. This area is of the same order of magnitude of that of our CWS model ($6.49 \times 10^{-19} \text{ m}^2$), but 1 order of magnitude different from that of the ASW model ($2.04 \times 10^{-18} \text{ m}^2$). While the inertia moments were computed with the MOL-DRAW program (Ugliengo et al. 1993; vide supra), Minissale et al. used the structures available in the ChemSpider<fnref n="4" /> database, giving rise to some differences with our results.

Table 3 shows that the Minissale et al. derived $\nu$ values differ by up 2 orders of magnitude (our values being always larger than the Minissale's ones) due to the above described differences. Notably, smaller species exhibit lower differences between our and Minissale's preexponential factors, while larger molecules present more pronounced differences.

Given the key role of the desorption rate preexponential factors in determining the BE(0)s by inversion of Equation (5), accuracy in its quantification becomes imperative. Ultimately, even minor differences in the input parameters in computing a preexponential factor can lead to substantial variations in its value, and consequently to the experimentally derived BE(0) value.

<footnote n="4">⁴ www.chemspider.com</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Comparison of the N-bearing species BEs computed in the present work with the distribution of water BEs computed by Tinacci et al. (2023; red curve). The bars represent the BE ranges of species belonging to groups 1 (magenta), 2 (blue), and 3 (green), defined in the text.</figcaption></figure>

## 4.4. Astrophysical Implications

As emphasized in the Introduction, BEs have a direct impact on the interstellar chemistry, i.e., the chemical composition of the molecular ISM and its evolution. Indeed, BEs directly determine the dust temperature at which species, frozen on dust grains, sublimate. Such a process is mathematically expressed as being exponentially proportional to BEs (see, e.g., the reviews by Minissale et al. 2022; Ceccarelli et al. 2023). Therefore, even small variations can result in significant differences in sublimation temperatures, which can lead to large differences in species abundances in the gas-phase and on the dust grains.

Particularly important in the context of planetary formation is the case of protoplanetary disks: if a species is in the solid form, it will be available for planetesimal accretion, contributing to the formation of rocky objects such as planets or asteroids. However, if a species is in the gas-phase, elements contained within the species are more likely to become part of the forming of giant gaseous planets (e.g., Dodson-Robinson et al. 2009; Öberg & Bergin 2021).

In protoplanetary disks, snow lines are the regions where gaseous species freeze out onto dust grains (e.g., Hartmann et al. 2017; Öberg & Bergin 2021). The locations of snow lines depend on both the dust temperature and the BE of the frozen species. In general, moving farther from the center of the disk, the dust temperature decreases, and consequently, species with low BEs remain gaseous at larger distances from the disk center, while those with larger BEs freeze out onto dust grains closer to the center disk. As a result, snow lines (and therefore BEs) have a direct relation on the gaseous versus solid chemical composition in protoplanetary disks, playing a crucial role in the chemical composition of the rocky/gaseous objects being formed in those regions (e.g., Öberg et al. 2011; Öberg & Wordsworth 2019; Mollière et al. 2022).

Within this context, it is interesting to compare the BEs of the N-containing species with those of water, since water is crucial for the emergence of life. In practice, the point is to know which N-bearing species are frozen or gaseous within the water snow line. As mentioned in the previous section, Tinacci et al. (2023) computed the BE distribution on a 200 water molecule amorphous ice model and found that it follows a

Gaussian distribution with a peak at 4230 K. Tinacci et al. (2023) also introduced the water BE distribution in a model of a protoplanetary disk and showed the influence of the BE distribution on the water snow line with respect to the radius assuming a single value. They also estimated the water content of possible planetesimals at a different distance from the disk center. Explicitly modeling the N-chemistry in a protoplanetary disk is beyond the scope of the present work, but we can still give order of magnitude arguments about what N-species will likely be part of planetesimals that are also enriched with water.

Figure 10 shows the ranges of the calculated BEs for the N-bearing species compared to the BE distribution of water (Tinacci et al. 2023). We can identify three groups of N-bearing species:

1. *The N-species that will completely sublimate in regions where 96% of water is frozen, namely, whose maximum BE is lower than ~1700 K: N₂ and NO.* These species will not be incorporated into the planetesimals. Note that N₂ is expected to be among the major nitrogen reservoirs.

2. *The N-species that (at least a fraction) remains frozen also in regions where 96% of water has completely sublimated, namely, whose maximum BE is larger than ~6500 K: CN, NH₂, HNC, NH₃, HNCO, HC₃N, NH₂CN, CH₃CN, HCONH₂ (indicated as NH₂CHO in astrochemical works). CH₃NH₂, and HC₅N.* These species will very likely be incorporated into water-rich planetesimals. It is particularly important for the case of ammonia, whose large fraction remains frozen in water-rich ices, and which is also a major nitrogen reservoir, with the caveat of a small low-BE tail not caught by our computations (see the discussion in Section 4; and Tinacci et al. 2022). Also interesting are the cases of HNCO and NH₂CHO, which contain the four elements CHON important for biological chemistry and are fully conserved as frozen molecules in the eventual planetesimals. In particular, formamide has been evoked as a crucial molecule that can lead to the formation of a multitude of prebiotic species (e.g., Saladino et al. 2012b), and, for this reason, its chemistry has been the focus of numerous works (see, e.g., the review by López-Sepulcre et al. 2019).

3. *All the other species in between, namely, those that have at least a fraction frozen in a region where more than 4% and less than 96% of water is frozen: N, NH, HCN, HNO, OCN, N₂O, NO₂, and CH₂CHCN.* Among these species, it is worth noticing that much of HCN and CH₂CHCN will be preserved in water-rich ices incorporating eventual planetesimals. In particular, HCN has been invoked as the most important brick leading to the formation of amino acids and nucleic acid bases (Oró 1960; Santalucia et al. 2022).

To sum up, the accurate quantification of BEs is crucial to accurately describe snow line positions when modeling planet formation. The distribution of BE directly affects the chemical composition of the forming rocky objects, as it determines which species remain in the gaseous phase and which ones are frozen out onto dust grains.

## 5. Conclusions

In this work, we computed quantum mechanically the BEs of 21 N-bearing species of astrochemical relevance (14 closed-shell
species and seven radicals) absorbed on two water ice models mimicking the mantle of interstellar ice grains. Two icy models are defined within the periodic boundary conditions, simulating a crystalline proton-ordered ice (CWS model) and an amorphous one (ASW model). The BEs have been computed on the larger ASW model with the cost-effective DFT//HF-3c procedure, in which the structure of the adsorbate/(icy grain) is relaxed at HF-3c level, followed by the punctual energy evaluation at full DFT level, i.e., B3LYP-D3(BJ)//HF-3c (closed-shell species) or M06-2X//HF-3c (open-shell species). The performance of this methodological scheme was validated on the CWS model, showing an excellent correlation with the BEs computed at DFT//DFT level and even when comparing BEs resulting from a more refined ONIOM correction at CCSD(T) level.

We proved that molecules exhibiting both H-bond acceptor and donor character (Hb-AD group) bind with higher-BE values than those with only H-bond acceptor character (Hb-A group). The CWS model has only one adsorption site due to symmetry constraints, resulting in a single BE value for each species. This is in contrast with the ASW model, exhibiting a sensitive number of different surface sites, giving rise to a range of BEs for each species. BEs computed for the CWS model are usually in the higher side of the BEs distributions computed for the ASW model. This is due to a better hydrogen bond cooperativity for the crystalline CWS model compared to the ASW one and, due to the rigidity of the CWS model, to a smaller deformation cost suffered by the icy surface due to the adsorption.

The smallest DFT//HF-3c computed BE value is for NO while the higher are for NH₂CHO and NH₂CN, in agreement with their double H-bonding donor/acceptor character. The BE distribution (500–12,500 K) is featureless and broader than the analogous oxygen-bearing molecules computed by Ferrero et al. (2020), in which a bimodal shape resulted. For sulfur-containing molecules, as computed by Ferrero et al. (2022), the BE distribution is single peaked spanning much narrower BE values (1500–6500 K) than the present one. Computed BEs have been compared with experimental and computational data available from the literature.

Previous BEs from theory overlap the range of our BEs, although they correspond to the lower part of our range due to the adoption of different icy models, which span from minimal clusters (Wakelam et al. 2017; Das et al. 2018; Hendrix et al. 2023) to more developed clusters and periodic systems (Duflot et al. 2021; Bovolenta et al. 2022; Germain et al. 2022; Tinacci et al. 2022; Tóth Ugyonka et al. 2024). Experimentally, the BEs obtained through the mathematical manipulation of the TPD spectra (Ferrero et al. 2022; Minissale et al. 2022) are in fairly good agreement with our computed results.

Finally, the computed BEs have been used to predict which N-bearing molecules can remain frozen in the protoplanetary disk region where also water is frozen and can potentially be incorporated into the water-rich planetesimals that will eventually form a rocky planet. We found that two important N-species, HCN and NH₂CHO, which have a particularly important prebiotic potential (Oró 1960; Saladino et al. 2012b; Santalucia et al. 2022), will likely be incorporated in those water-rich planetesimal ices.

## Acknowledgments

This project has received funding from the Marie Skłodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement No. 811312, and 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. MICIN (projects PID2021-126427NB-100 and CNS2023-144902) is also acknowledged. The Italian Space Agency for cofunding the Life in Space Project (ASI N. 2019-3-U.O), the Italian MUR (PRIN 2020, Astrochemistry beyond the second period elements, Prot. 2020AFB3FX) are also acknowledged for financial support. Support from the Project CH4.0 under the MUR program “Dipartimenti di Eccellenza 2023-2027” (CUP: D13C22003520001) is also acknowledged. The authors thankfully acknowledge RES resources provided by BSC in MareNostrum to activities QHS-2022-3-0007 and QHS-2023-2-0011, and the supercomputational facilities provided by CSUC.

## Appendix A
### Computational Details

We adopted computational parameters used in Ferrero et al. (2020) for optimization and frequency calculations. To sample the reciprocal space, the Pack-Monkhorst mesh was adopted. For the crystalline ice, the SHRINK factor was set to 2 resulting in four K-points in the first Brillouin zone. For the amorphous ice, the SHRINK factor was set to 1, resulting in two K-points in the first Brillouin zone. The threshold parameters for the evaluation of the Coulomb and exchange bielectronic integrals (TOLINTEG) were set to 7, 7, 7, 7, and 14 for both models, the crystalline and the amorphous ices. Some species required the use of BROYDEN accelerator for the convergence of the SCF. Integration grid (XLGRID) was added and fixed to default values.

## Appendix B
### BE Detailed Calculation

#### B.1. BBSE-corrected Interaction Energies

BEs are calculated from the interaction between a surface, S, in this case represented by a slab model, and the adsorbed molecular species, $\mu$,

$$ \Delta E = E_P(S \cdot \mu // S \cdot \mu) - E_P(S // S) - E_M(\mu // \mu) \quad (B1) $$

where $E_Z$ (xy//x'y') stands for the energy of the species xy in the optimized geometry x'y', and M and P stand for molecular or periodic calculations. In this case, $E_P(S \cdot \mu // S \cdot \mu)$ is the energy of the species adsorbed on the surface in its optimized geometry, $E_P(S // S)$ is the energy of surface, and $E_M(\mu // \mu)$ stands for the energy of the optimized single species.

This equation does not consider some factors such as the deformation of S and the $\mu$, the BSSE, or lateral interactions between $\mu$ of different unit cells. Some corrections must be done to obtain the corrected interaction energy. First, the deformation of both the slab model and the species must be quantified.

$$ \delta E_S = E_P(S // S \cdot \mu) - E_P(S // S). \quad (B2) $$

$$ \delta E_{\mu} = E_{P}(\mu // S \cdot \mu) - E_{M}(\mu // \mu). \quad (B3) $$
So, the deformation-free interaction energy $\Delta E^*$ is defined as

$$ \Delta E^* = E_P(S \cdot \mu // S \cdot \mu) - E_P(S // S \cdot \mu) - E_P(\mu // S \cdot \mu). \quad (B4) $$

Due to the use of a finite basis set, BSSE arises, which must be corrected, in this case, adopting the Boys-Bernardi CP correction (Boys & Bernardi 1970).

$$ \mathrm{BSSE}(S) = E_P(S \cdot [\mu] // S \cdot \mu) - E_P(S \cdot \mu // S \cdot \mu), \quad (B5) $$

$$ \mathrm{BSSE}(\mu) = E_P([S] \cdot \mu // S \cdot \mu) - E_P(S \cdot \mu // S \cdot \mu), \quad (B6) $$

where $E_Z$ (x[y]//x'y') is the energy of the x species in the presence of the ghost functions of the y species computed in the geometry optimization of the x'y'.

The last correction corresponds to the lateral interaction $\delta E_L$ due to interaction between species of different unit cells.

$$ \delta E_L = E_P(\mu // S \cdot \mu) - E_M(\mu // S \cdot \mu). \quad (B7) $$

Thus, BSSE-corrected interaction energies considering all these corrections can be expressed as

$$ \Delta E^{CP} = \Delta E^* + \delta E_S + \delta E_\mu + \delta E_L - \mathrm{BSSE}(S) - \mathrm{BSSE}(\mu). \quad (B8) $$

BEs are positive values and are defined as the absolute value of $\Delta E^{CP}$ and including the corresponding ZPE correction.

$$ \mathrm{BE} = |\Delta E^{CP}| - \Delta \mathrm{ZPE} \quad (B9) $$

in which the ZPE correction $\Delta ZPE$ is computed as

$$ \Delta ZPE = ZPE(S \cdot \mu\mu) - ZPE(S) - ZPE(E_M). \quad (B10) $$

## Appendix C: Dispersion Contribution

Dispersive contributions for BSSE-corrected interaction energies can be quantified for closed-shell species. For open-shell species and due to the nature of M06-2X species, dispersive contributions cannot be separated from total BEs. Table 6 shows BSSE-corrected interaction energies, dispersion contributions, and their percentages for all closed-shell species. All cases computed on the ASW are included as well as the ones computed on the CWS at DFT//HF-3c and at DFT//DFT methodologies.

<table><caption>**Table 6**
Computed BEs on the ASW at DFT//HF-3c Level and on the CWS at DFT//HF-3c and at DFT//DFT Level for Closed-shell Species</caption><thead><tr><th rowspan="2">Species</th><th rowspan="2"></th><th colspan="2">CWS</th><th colspan="10">ASW</th></tr><tr><th>DFT//HF-3c</th><th>DFT//DFT</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td rowspan="3">NH<sub>3</sub></td><td>BE</td><td>7301</td><td>7373</td><td>8840</td><td>5268</td><td>6820</td><td>5930</td><td>6579</td><td>6098</td><td>5052</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>1816</td><td>1852</td><td>2917</td><td>843</td><td>2864</td><td>1838</td><td>1118</td><td>1768</td><td>758</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>25</td><td>25</td><td>33</td><td>16</td><td>42</td><td>31</td><td>17</td><td>29</td><td>15</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HCN</td><td>BE</td><td>5725</td><td>5124</td><td>2923</td><td>5124</td><td>3620</td><td>5136</td><td>3271</td><td>5857</td><td>6146</td><td>7421</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3043</td><td>2057</td><td>1111</td><td>2357</td><td>1050</td><td>2208</td><td>1308</td><td>3456</td><td>2458</td><td>1633</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>53</td><td>40</td><td>38</td><td>46</td><td>29</td><td>43</td><td>40</td><td>59</td><td>40</td><td>22</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HNC</td><td>BE</td><td>8482</td><td>8756</td><td>7772</td><td>6504</td><td>5664</td><td>9335</td><td>11,860</td><td>8473</td><td>5488</td><td>4496</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3280</td><td>2936</td><td>2019</td><td>3372</td><td>2844</td><td>3042</td><td>2293</td><td>3037</td><td>2488</td><td>766</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>39</td><td>34</td><td>26</td><td>52</td><td>50</td><td>33</td><td>19</td><td>36</td><td>45</td><td>17</td><td>...</td><td>...</td></tr><tr><td rowspan="3">N<sub>2</sub></td><td>BE</td><td>1455</td><td>1564</td><td>1347</td><td>1708</td><td>1311</td><td>1191</td><td>890</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>1636</td><td>1636</td><td>1859</td><td>2255</td><td>1377</td><td>1763</td><td>409</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>112</td><td>105</td><td>138</td><td>132</td><td>105</td><td>148</td><td>46</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HNO</td><td>BE</td><td>4601</td><td>5150</td><td>2761</td><td>3154</td><td>2606</td><td>3182</td><td>3700</td><td>4608</td><td>3917</td><td>3196</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>1556</td><td>1339</td><td>1524</td><td>649</td><td>705</td><td>1135</td><td>1231</td><td>1748</td><td>1321</td><td>1135</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>34</td><td>26</td><td>55</td><td>21</td><td>27</td><td>36</td><td>33</td><td>38</td><td>34</td><td>36</td><td>...</td><td>...</td></tr><tr><td rowspan="3">CH<sub>3</sub>NH<sub>2</sub></td><td>BE</td><td>8498</td><td>8455</td><td>6779</td><td>6648</td><td>7110</td><td>8158</td><td>5146</td><td>9760</td><td>7234</td><td>5930</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3060</td><td>2918</td><td>2545</td><td>1345</td><td>1597</td><td>2882</td><td>992</td><td>4095</td><td>2300</td><td>1764</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>36</td><td>35</td><td>38</td><td>20</td><td>22</td><td>35</td><td>19</td><td>42</td><td>32</td><td>30</td><td>...</td><td>...</td></tr><tr><td rowspan="3">CH<sub>3</sub>CN</td><td>BE</td><td>7090</td><td>7551</td><td>5790</td><td>7551</td><td>5094</td><td>8918</td><td>5465</td><td>8926</td><td>3088</td><td>4029</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3373</td><td>3090</td><td>3381</td><td>1962</td><td>2537</td><td>4612</td><td>2038</td><td>4601</td><td>2818</td><td>1326</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>48</td><td>41</td><td>58</td><td>26</td><td>50</td><td>52</td><td>37</td><td>52</td><td>91</td><td>33</td><td>...</td><td>...</td></tr><tr><td rowspan="3">NH<sub>2</sub>CN</td><td>BE</td><td>10,413</td><td>10,312</td><td>8327</td><td>7716</td><td>9334</td><td>12,601</td><td>2946</td><td>11,716</td><td>12,597</td><td>7713</td><td>6155</td><td>8670</td></tr><tr><td>Disp</td><td>3478</td><td>3164</td><td>1759</td><td>1293</td><td>2034</td><td>5560</td><td>619</td><td>1811</td><td>5571</td><td>3975</td><td>1480</td><td>2991</td></tr><tr><td>Disp (%)</td><td>33</td><td>31</td><td>21</td><td>17</td><td>22</td><td>44</td><td>21</td><td>15</td><td>44</td><td>52</td><td>24</td><td>34</td></tr><tr><td rowspan="3">HNCO<sup>a</sup></td><td>BE</td><td>8232</td><td>8340</td><td>6066</td><td>8158</td><td>5183</td><td>5162</td><td>9388</td><td>7632</td><td>4883</td><td>3989</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3079</td><td>2944</td><td>1692</td><td>2251</td><td>2682</td><td>2727</td><td>3956</td><td>2486</td><td>2187</td><td>1073</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>37</td><td>35</td><td>28</td><td>28</td><td>52</td><td>53</td><td>42</td><td>33</td><td>45</td><td>27</td><td>...</td><td>...</td></tr></tbody></table>
<table><caption>**Table 6**
(Continued)</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">CWS</th><th colspan="10">ASW</th></tr><tr><th>DFT//HF-3c</th><th>DFT//DFT</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td rowspan="3">N<sub>2</sub>O</td><td>BE</td><td>3092</td><td>3288</td><td>2544</td><td>1588</td><td>2340</td><td>3580</td><td>4617</td><td>2811</td><td>3581</td><td>1620</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>2552</td><td>2484</td><td>2690</td><td>1015</td><td>1689</td><td>2802</td><td>2148</td><td>2807</td><td>2802</td><td>1084</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>83</td><td>76</td><td>106</td><td>64</td><td>72</td><td>78</td><td>47</td><td>100</td><td>78</td><td>67</td><td>...</td><td>...</td></tr><tr><td rowspan="3">NH<sub>2</sub>CHO</td><td>BE</td><td>9634</td><td>10,079</td><td>12,833</td><td>7481</td><td>10,344</td><td>6820</td><td>8467</td><td>7072</td><td>6783</td><td>7313</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3560</td><td>3608</td><td>4107</td><td>2244</td><td>4862</td><td>1296</td><td>3979</td><td>2051</td><td>3052</td><td>1389</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>37</td><td>36</td><td>32</td><td>30</td><td>47</td><td>19</td><td>47</td><td>29</td><td>45</td><td>19</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HC<sub>3</sub>N</td><td>BE</td><td>3781</td><td>...</td><td>4816</td><td>3163</td><td>3477</td><td>7808</td><td>5630</td><td>7610</td><td>4373</td><td>2223</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3707</td><td>...</td><td>3941</td><td>2340</td><td>1224</td><td>4566</td><td>1764</td><td>4550</td><td>4327</td><td>710</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>98</td><td>...</td><td>82</td><td>74</td><td>35</td><td>58</td><td>31</td><td>60</td><td>99</td><td>32</td><td>...</td><td>...</td></tr><tr><td rowspan="3">CH<sub>2</sub>CHCN</td><td>BE</td><td>5818</td><td>6199</td><td>5151</td><td>7545</td><td>5832</td><td>4682</td><td>7705</td><td>5832</td><td>2957</td><td>4590</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>3356</td><td>3350</td><td>3448</td><td>3032</td><td>4028</td><td>4718</td><td>1869</td><td>2065</td><td>2197</td><td>3448</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>58</td><td>54</td><td>67</td><td>40</td><td>69</td><td>101</td><td>24</td><td>35</td><td>74</td><td>75</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HC<sub>5</sub>N</td><td>BE</td><td>...</td><td>...</td><td>9467</td><td>712</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp</td><td>...</td><td>...</td><td>4489</td><td>726</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Disp (%)</td><td>...</td><td>...</td><td>47</td><td>102</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.**

## Appendix D
### Deformation Energies

The deformation energies of the surface and of the species are reported in Table 7.

<table><caption>**Table 7**
Deformation Energies of the Surface (δE<sub>S</sub>) and of the Species (δE<sub>μ</sub>) Computed on the ASW at DFT//HF-3c Level and on the CWS at DFT//HF-3c and DFT//DFT Level</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">CWS</th><th colspan="10">ASW</th></tr><tr><th>DFT//HF-3c</th><th>DFT//DFT</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td rowspan="2">N</td><td>δE<sub>S</sub></td><td>-84</td><td>143</td><td>-96</td><td>168</td><td>36</td><td>-373</td><td>-397</td><td>289</td><td>-12</td><td>-24</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">NH</td><td>δE<sub>S</sub></td><td>229</td><td>156</td><td>180</td><td>457</td><td>216</td><td>1034</td><td>-493</td><td>577</td><td>241</td><td>108</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>12</td><td>-84</td><td>-12</td><td>24</td><td>-36</td><td>12</td><td>36</td><td>12</td><td>-12</td><td>12</td><td>...</td><td>...</td></tr><tr><td rowspan="2">NH<sub>2</sub></td><td>δE<sub>S</sub></td><td>409</td><td>481</td><td>-1058</td><td>674</td><td>529</td><td>662</td><td>168</td><td>156</td><td>758</td><td>325</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>108</td><td>84</td><td>36</td><td>60</td><td>24</td><td>96</td><td>-12</td><td>36</td><td>108</td><td>24</td><td>...</td><td>...</td></tr><tr><td rowspan="2">NH<sub>3</sub></td><td>δE<sub>S</sub></td><td>738</td><td>451</td><td>1831</td><td>437</td><td>471</td><td>1916</td><td>750</td><td>988</td><td>591</td><td>...</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>-329</td><td>-41</td><td>-124</td><td>-65</td><td>-84</td><td>-137</td><td>-7</td><td>22</td><td>-43</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">CN</td><td>δE<sub>S</sub></td><td>...</td><td>...</td><td>1600</td><td>361</td><td>649</td><td>1323</td><td>120</td><td>902</td><td>637</td><td>337</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>...</td><td>...</td><td>-469</td><td>-325</td><td>-397</td><td>-337</td><td>-313</td><td>277</td><td>361</td><td>-469</td><td>...</td><td>...</td></tr><tr><td rowspan="2">HCN</td><td>δE<sub>S</sub></td><td>593</td><td>308</td><td>2355</td><td>1131</td><td>703</td><td>281</td><td>2368</td><td>1635</td><td>253</td><td>1629</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>94</td><td>-795</td><td>216</td><td>114</td><td>-19</td><td>-83</td><td>53</td><td>-168</td><td>38</td><td>-168</td><td>...</td><td>...</td></tr><tr><td rowspan="2">HNC</td><td>δE<sub>S</sub></td><td>1127</td><td>853</td><td>2345</td><td>2021</td><td>2189</td><td>1768</td><td>-589</td><td>2430</td><td>2333</td><td>505</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>671</td><td>494</td><td>698</td><td>541</td><td>589</td><td>830</td><td>758</td><td>625</td><td>613</td><td>421</td><td>...</td><td>...</td></tr><tr><td rowspan="2">N<sub>2</sub></td><td>δE<sub>S</sub></td><td>50</td><td>25</td><td>328</td><td>-54</td><td>65</td><td>74</td><td>-21</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>-86</td><td>-65</td><td>-1</td><td>2</td><td>4</td><td>1</td><td>2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">NO</td><td>δE<sub>S</sub></td><td>-36</td><td>60</td><td>-192</td><td>373</td><td>48</td><td>-24</td><td>457</td><td>24</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>-84</td><td>12</td><td>96</td><td>-241</td><td>1</td><td>12</td><td>-132</td><td>-156</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">HNO</td><td>δE<sub>S</sub></td><td>337</td><td>253</td><td>565</td><td>1636</td><td>120</td><td>541</td><td>1612</td><td>1119</td><td>529</td><td>541</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>144</td><td>72</td><td>73</td><td>84</td><td>36</td><td>132</td><td>144</td><td>96</td><td>24</td><td>-12</td><td>...</td><td>...</td></tr><tr><td rowspan="2">CH<sub>3</sub>NH<sub>2</sub></td><td>δE<sub>S</sub></td><td>575</td><td>747</td><td>509</td><td>1082</td><td>649</td><td>361</td><td>373</td><td>1672</td><td>1143</td><td>830</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>-12</td><td>35</td><td>-72</td><td>48</td><td>-24</td><td>60</td><td>-36</td><td>-23</td><td>12</td><td>-47</td><td>...</td><td>...</td></tr><tr><td rowspan="2">CH<sub>3</sub>CN</td><td>δE<sub>S</sub></td><td>293</td><td>382</td><td>216</td><td>-1395</td><td>589</td><td>1539</td><td>433</td><td>-1058</td><td>926</td><td>192</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>52</td><td>95</td><td>36</td><td>37</td><td>47</td><td>-253</td><td>96</td><td>35</td><td>49</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">OCN</td><td>δE<sub>S</sub></td><td>...</td><td>...</td><td>96</td><td>373</td><td>168</td><td>553</td><td>-253</td><td>926</td><td>120</td><td>119</td><td>...</td><td>...</td></tr><tr><td>δE<sub>μ</sub></td><td>...</td><td>...</td><td>-132</td><td>-180</td><td>-72</td><td>-60</td><td>108</td><td>36</td><td>-120</td><td>34</td><td>...</td><td>...</td></tr><tr><td rowspan="2">NH<sub>2</sub>CN</td><td>δE<sub>S</sub></td><td>625</td><td>682</td><td>1203</td><td>1732</td><td>349</td><td>1960</td><td>241</td><td>445</td><td>1852</td><td>1335</td><td>541</td><td>1131</td></tr></tbody></table>

<footnote>ª BEs have been computed using the DFT//PBE method as an alternative to the DFT//HF-3c procedure; see Section 3.1.1.</footnote>
<table><caption>**Table 7**
(Continued)</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">CWS</th><th colspan="10">ASW</th></tr><tr><th>DFT//HF-3c</th><th>DFT//DFT</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td>HNCO</td><td>δE<sub>μ</sub></td><td>36</td><td>250</td><td>-24</td><td>216</td><td>93</td><td>-156</td><td>204</td><td>108</td><td>-156</td><td>-108</td><td>120</td><td>13</td></tr><tr><td></td><td>δE<sub>s</sub></td><td>553</td><td>806</td><td>1720</td><td>1828</td><td>1407</td><td>613</td><td>1335</td><td>2117</td><td>854</td><td>481</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>565</td><td>349</td><td>-216</td><td>120</td><td>541</td><td>144</td><td>120</td><td>142</td><td>565</td><td>203</td><td>...</td><td>...</td></tr><tr><td>N<sub>2</sub>O</td><td>δE<sub>s</sub></td><td>79</td><td>281</td><td>-21</td><td>1287</td><td>192</td><td>1070</td><td>119</td><td>601</td><td>710</td><td>192</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>153</td><td>23</td><td>108</td><td>132</td><td>119</td><td>325</td><td>108</td><td>253</td><td>204</td><td>122</td><td>...</td><td>...</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>δE<sub>s</sub></td><td>1264</td><td>545</td><td>2769</td><td>937</td><td>1472</td><td>-719</td><td>1026</td><td>1044</td><td>3318</td><td>-1045</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>-72</td><td>-478</td><td>91</td><td>99</td><td>65</td><td>98</td><td>-81</td><td>21</td><td>249</td><td>67</td><td>...</td><td>...</td></tr><tr><td>NO<sub>2</sub></td><td>δE<sub>s</sub></td><td>72</td><td>301</td><td>-120</td><td>36</td><td>96</td><td>842</td><td>-36</td><td>96</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>-36</td><td>48</td><td>-49</td><td>-12</td><td>-25</td><td>-13</td><td>-25</td><td>59</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>HC<sub>3</sub>N</td><td>δE<sub>s</sub></td><td>127</td><td>...</td><td>24</td><td>1034</td><td>842</td><td>229</td><td>385</td><td>854</td><td>457</td><td>216</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>-14</td><td>...</td><td>36</td><td>144</td><td>229</td><td>37</td><td>-227</td><td>216</td><td>60</td><td>117</td><td>...</td><td>...</td></tr><tr><td>CH<sub>2</sub>CHCN</td><td>δE<sub>s</sub></td><td>370</td><td>396</td><td>61</td><td>253</td><td>1696</td><td>2261</td><td>-1588</td><td>251</td><td>1828</td><td>1143</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>85</td><td>38</td><td>36</td><td>71</td><td>-253</td><td>23</td><td>47</td><td>73</td><td>97</td><td>-49</td><td>...</td><td>...</td></tr><tr><td>HC<sub>5</sub>N</td><td>δE<sub>s</sub></td><td>...</td><td>...</td><td>385</td><td>204</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td></td><td>δE<sub>μ</sub></td><td>...</td><td>...</td><td>168</td><td>96</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Note. Units are in kelvins.
## Appendix E
### Hemibond Geometries for the CN Species

Figure 11 provides the adsorption geometries of the CN species on the ASW for the cases when a hemibond system is formed.

<figure><img src="image_16.png" /><figcaption>**Figure 11.** CN adsorption geometries for cases when a hemibond is formed. I for case 1; II for case 5; and III for case 6. Distances of H-bonds and hemibonds are included. Units are in Å.</figcaption></figure>

## Appendix F
### Desorption Prefactors

Figure 12 provides the temperature dependence of the desorption rate prefactors computed for all the species in both models, the CWS and the ASW.
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Temperature dependence of the desorption rate prefactors computed using Equation (5) (Tait et al. 2005). Black line corresponds to the ASW model, and orange line to the CWS model. Parameters used computed the desorption rate prefactors temperature dependence are provided in Table 3.</figcaption></figure>
ORCID iDs

Berta Martínez-Bachs https://orcid.org/0000-0002-1290-7019

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

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

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

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

References

Al-Halabi, A., & Van Dishoeck, E. F. 2007, MNRAS, 382, 1648

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

Becke, A. D. 1993, CP, 98, 1372

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

Bovolenta, G., Bovino, S., Vöhringer-Martinez, E., et al. 2020, MolAs, 21, 100095

Bovolenta, G. M., Vogt-Geisse, S., Bovino, S., & Grassi, T. 2022, ApJS, 262, 17

Boys, S., & Bernardi, F. 1970, MolPh, 19, 553

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

Ceccarelli, C., Codella, C., Balucani, N., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 379

Ceccarelli, C. 2023, Faraday Discuss., 245, 11

Chaabouni, H., Diana, S., Nguyen, T., & Dulieu, F. 2018, A&A, 612, A47

Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133

Cramer, C. J., & Bickelhaupt, F. 2003, Angew. Chem. Int. Ed., 42, 381

Cuppen, H., Walsh, C., Lamberts, T., et al. 2017, SSRv, 212, 1

Dapprich, S., Komáromi, I., Byun, K., Morokuma, K., & Frisch, M. J. 1999, JMoSt, 461-462, 1

Das, A., Sil, M., Gorai, P., Chakrabarti, S., & Loison, J. 2018, ApJS, 237, 9

Davidson, E. R., & Feller, D. 1986, ChRv, 86, 681

Dodson-Robinson, S. E., Willacy, K., Bodenheimer, P., Turner, N. J., & Beichman, C. A. 2009, Icar, 200, 672

Dohnálek, Z., Kimmel, G. A., Joyce, S. A., et al. 2001, JPhCB, 105, 3747

Douglas, A., & Herzberg, G. 1941, ApJ, 94, 381

Dovesi, R., Erba, A., Orlando, R., et al. 2018, Wiley Interdiscip. Rev. Comput. Mol. Sci., 8, e1360

Duflot, D., Toubin, C., & Monnerville, M. 2021, FrASS, 8, 645243

Dulieu, F., Congiul, F., Noble, J., et al. 2013, NatSR, 3, 2045

Dunham, T., Jr 1937, PASP, 49, 26

Dunning, T. H. 1989, CP, 90, 1007

Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. 2016, JMoSp, 327, 95

Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJ, 816, L28

Ferrero, S., Grieco, F., Ibrahim Mohamed, A. S., et al. 2022, MNRAS, 516, 2586

Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11

Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian16 Revision C.01 (Wallingford, CT: Gaussian Inc.)

Germain, A., Tinacci, L., Pantaleone, S., Ceccarelli, C., & Ugliengo, P. 2022, ACS Earth Space Chem., 6, 1286

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

Hartmann, L., Ciesla, F., Gressel, O., & Alexander, R. 2017, SSRv, 212, 813

He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 825, 89

Hendrix, J., Bera, P. P., Lee, T. J., & Head-Gordon, M 2023, MolPh, 122, e2252100

Jansen, H., & Ros, P. 1969, CPL, 3, 140

Karssemeijer, L. J., & Cuppen, H. M. 2014, A&A, 569, A107

Karssemeijer, L. J., de Wijs, G. A., & Cuppen, H. M. 2014, Phys. Chem. Chem. Phys., 16, 15630

Kolasinski, K. W. 2012, Surface Science: Foundations of Catalysis and Nanoscience (New York: John Wiley & Sons),

Lee, C., Yang, W., & Parr, R. G. 1988, PhRvB, 37, 785

Liu, B., & McLean, A. 1973, JChPh, 59, 4557

López-Sepulcre, A., Balucani, N., Ceccarelli, C., et al. 2019, ESC, 3, 2122

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, NatAs, 7, 431

McGuire, B. A. 2022, ApJS, 259, 30

McKellar, A. 1940, PASP, 52, 187

Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ESC, 6, 597

Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74

Noble, J. A., Congiu, E., Dulieu, F., & Fraser, H. J. 2012, MNRAS, 421, 768

Oba, Y., Miyauchi, N., Hidaka, H., et al. 2009, ApJ, 701, 464

Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194

Oró, J. 1960, BBRC, 2, 407

Papajak, E., Zheng, J., Xu, X., Leverentz, H. R., & Truhlar, D. G. 2011, J. Chem. Theory Comput., 7, 3027

Pascale, F., Zicovich-Wilson, C. M., López Gejo, F., et al. 2004, J. Comput. Chem., 25, 888

Penteado, E., Walsh, C., & Cuppen, H. 2017, ApJ, 844, 71

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

Ferrero, J., Enrique-Romero, J., Ferrero, S., et al. 2022, ApJ, 938, 158

Pisani, C., Casassa, S., & Ugliengo, P. 1996, CPL, 253, 201

Pople, J. A., Gill, P. M. W., & Handy, N. C. 1995, IJQC, 56, 303

Rimola, A., Skouteris, D., Balucani, N., et al. 2018, ESC, 2, 720

Saitta, A. M., & Saija, F. 2014, PNAS, 111, 13768

Saladino, R., Botta, G., Pino, S., Costanzo, G., & Di Mauro, E. 2012a, Chem. Soc. Rev., 41, 5526

Saladino, R., Crestini, C., Pino, S., Costanzo, G., & Di Mauro, E. 2012b, Phys. Life Rev., 9, 84

Sameera, W., Senevirathne, B., Andersson, S., Maseras, F., & Nyman, G. 2017, J. Phys. Chem. C, 121, 15223

Santalucia, R., Pazzi, M., Bonino, F., et al. 2022, PCCP, 24, 7224

Schäfer, A., Horn, H., & Ahlrichs, R. 1992, CP, 97, 2571

Senevirathne, B., Andersson, S., Dulieu, F., & Nyman, G. 2017, MolAs, 6, 59

Shimonishi, T., Nakatani, N., Furuya, K., & Hama, T. 2018, AJ, 855, 27

Smith, R. S., May, R. A., & Kay, B. D. 2016, JPCB, 120, 1979

Sure, R., & Grimme, S. 2013, JCoCh, 34, 1672

Swings, P., & Rosenfeld, L. 1937, ApJ, 86, 483

Tait, S. L., Dohnálek, Z., Campbell, C. T., & Kay, B. D. 2005, CP, 122, 164708

Tinacci, L., Germain, A., Pantaleone, S., et al. 2023, ApJ, 951, 32

Tatewaki, H., & Huzinaga, S. 1980, JCoCh, 1, 205

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

Tóth Ugyonka, H., Hantal, G., Szőri, M., & Jedlovszky, P. 2024, ESC, 8, 520

Ugliengo, P., Viterbo, D., & Chiari, G. 1993, Z. Kristallogr. Cryst. Mater., 207, 9

Wakelam, V., Loison, J., Mereau, R., & Ruaud, M. 2017, MolAs, 6, 22

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

Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427

Zamirri, L., Corno, M., Rimola, A., & Ugliengo, P. 2017, ESC, 1, 384

Zamirri, L., Pantaleone, S., & Ugliengo, P. 2019, JChPh, 150, 064702

Zicovich-Wilson, C. M., Pascale, F., Roetti, C., et al. 2004, JCoCh, 25, 1873