# The binding energy distribution of H2S: why it is not the major sulphur reservoir of the interstellar ices

Vittorio Bariosco, Stefano Pantaleone, Cecilia Ceccarelli, Albert Rimola, Nadia Balucani, Marta Corno, Piero Ugliengo

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

## To cite this version:

Vittorio Bariosco, Stefano Pantaleone, Cecilia Ceccarelli, Albert Rimola, Nadia Balucani, et al.. The binding energy distribution of H2S: why it is not the major sulphur reservoir of the interstellar ices. Monthly Notices of the Royal Astronomical Society, 2024, 531, pp.1371-1384. 10.1093/mn-ras/stae1210. insu-04836854

HAL Id: insu-04836854

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

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The binding energy distribution of H₂S: why it is not the major sulphur reservoir of the interstellar ices

Vittorio Bariosco<sup>1,2</sup> Stefano Pantaleone<sup>1</sup> Cecilia Ceccarelli<sup>1,3★</sup> Albert Rimola<sup>1,2</sup>
Nadia Balucani<sup>1,4</sup> Marta Corno<sup>1</sup> and Piero Ugliengo<sup>1★</sup>

Accepted 2024 May 2. Received 2024 May 2; in original form 2024 February 28

## ABSTRACT

Despite hydrogen sulphide (H₂S) has been predicted to be the major reservoir of S-bearing species on the icy mantles of interstellar grains, no solid H₂S has been detected so far. A crucial parameter that governs whether or not a species remains frozen on to the grain mantles is its binding energy (BE). We present a new computational study of the H₂S BE on a large amorphous water ice surface, constituted by 200 water molecules. The resulting H₂S BE distribution ranges from 57 K (0.5 kJ mol⁻<fnref n="1" />) to 2406 K (20.0 kJ mol⁻<fnref n="1" />), with the average µ = 984 K (8.2 kJ mol⁻<fnref n="1" />). We discuss the reasons why the low bound of the newly computed BE distribution, which testifies to the very weak interaction of H₂S with the ice surface, has never been found by previous theoretical or experimental works before. In addition, the low H₂S BEs may also explain why frozen H₂S is not detected in interstellar ices. Following previous molecular dynamics studies that show that the energy of reactions occurring on ice surfaces is quickly absorbed by the water molecules of the ice and conservatively assuming that 10 per cent of the HS + H → H₂S formation energy (−369.5 kJ mol⁻<fnref n="1" />) is left to the newly formed H₂S, its energy is more than twice the largest BE and five times the average BE and, hence, H₂S will most likely leave the water surface.

**Key words:** astrochemistry – molecular data – molecular processes – ISM: abundances – ISM: molecules.

## 1 INTRODUCTION

Sulphur is the seventh most abundant heavy (atomic number ≥3) element in the Sun photosphere and the fourth (after O, C, and N) when excluding those largely locked in interstellar refractory dust (Mg and Si) and noble gas (Ne) (Asplund et al. 2009). Remarkably, in terrestrial living beings, sulphur is also the most abundant heavy element after O, C, and N and, consequently, crucial for terrestrial life. Therefore, regardless whether interstellar molecules have been involved in triggering life on Earth or not, understanding the sulphur chemistry in space and its evolution is of paramount importance.

In diffuse clouds, the sulphur abundance is very similar to that measured in the Sun photosphere (S/H ~1.4 × 10⁻⁵; Jenkins 2009). In dense (~10⁴ cm⁻<fnref n="3" />) and cold (~10 K) clouds, it is more than a factor of 100 less abundant (e.g. Tieftrunk et al. 1994; Ruffle et al. 1999; Vastel et al. 2018), as it is the case also for other elements, such as oxygen which is largely locked in frozen water, the major component of the so-called dust icy mantles (Boogert, Gerakines & Whittet 2015). However, while gaseous oxygen abundance goes back to its initial elemental one in warm (≥100 K) objects (e.g. hot cores and hot corinos), where the dust icy mantles sublimate, sulphur

does not. In warm objects, adding up the abundances of all detected gaseous S-bearing species gives a total of ~1 per cent the elemental S abundance (e.g. Blake et al. 1994; Wakelam et al. 2004; Codella et al. 2021). This constitutes a decades-long riddle: where is sulphur in molecular interstellar medium (ISM)?

One could think that sulphur is locked in the refractory dust grain cores (the binding energies of S-bearing species on silicates are indeed larger than on water ices: Perrero et al. 2024), but the abundance in the diffuse clouds (see above) excludes this possibility. Therefore, sulphur has to be in semirefractory or volatile species that are undetectable when they are in the gas phase. A systematic study towards the cold object TMC1 by Fuente et al. (2023) strongly supports this hypothesis (see also Vidal & Wakelam 2018; Laas & Caselli 2019) but still does not answer the question about sulphur major species.

In order to answer the question of where sulphur is in molecular ISM, we need to understand its chemical evolution from diffuse to molecular clouds. In diffuse clouds, sulphur is in the atomic form and, therefore, when the cloud becomes molecular, sulphur should approximately follow the fate of the oxygen atoms (the depletion time being only a factor 1.4 larger, due to the heavier mass of S with respect to O), i.e. being hydrogenated once S atoms freeze-out on to the grain surfaces and mostly become frozen H₂S. Yet, while astrochemical models tend to predict frozen H₂S to be the major reservoir of sulphur in the dust mantles, astronomical observations, even the most recent

<footnote n="1">¹Dipartimento di Chimica, Università degli Studi di Torino, via P. Giuria 7, I-10125 Torino, Italy</footnote>

<footnote n="2">²Departament de Química, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Catalonia, Spain</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4">⁴Dipartimento di Chimica, Biologia e Biotecnologie, Università degli Studi di Perugia, Perugia, I-06123, Italy</footnote>

<footnote>* E-mail: cecilia.ceccarelli@univ-grenoble-alpes.fr (CC); piero.ugliengo@unito.it (PU)</footnote>
ones by JWST, have miserably failed to detect frozen H<sub>2</sub>S at a level at least 10 times lower than the predicted one (e.g. Boogert et al. 2015; McClure et al. 2023). Only OCS and, tentatively, SO<sub>2</sub> have been so far detected in the solid phase and their abundance is too low to account for the total sulphur elemental abundance (e.g. Palumbo, Tielens & Tokunaga 1995; Palumbo, Geballe & Tielens 1997; Boogert et al. 2015; Sturm et al. 2023). Backing this lack of major frozen H<sub>2</sub>S, its abundance in warm objects is always lower than the sulphur elemental abundance (e.g. Kushwahaa et al. 2023). Why sulphur is not mostly frozen H<sub>2</sub>S is, thus, in itself a riddle.

Adding to the sulphur mystery, cometary ices, which are believed to be mostly made up of interstellar ices, do not seem to be depleted of sulphur. Summing up all the S-bearing species detected in the comet 67P/Churyumov–Gerasimenko gives the solar sulphur abundance (Calmonte et al. 2016). What is more, in the same comet, H<sub>2</sub>S is the most abundant S-bearing species (~57 per cent) followed by atomic sulphur (~27 per cent), all other S-species being present in much smaller quantities (Calmonte et al. 2016). Other comets, where remote observations are available in the literature, show a similar behaviour, even though there is a relatively large variability on the H<sub>2</sub>S abundance with respect to other S-bearing species (Biver et al. 2021). Finally, very recent observations by Altwegg et al. (2022) have found very abundant (sometimes as much as water!) ammonium hydro-sulphide (NH<sub>4</sub><sup>+</sup>SH<sup>-</sup>) salts, a semi-refractory species, in the 67P/Churyumov–Gerasimenko comet.

In summary, it seems that, in order to solve the missing sulphur mystery in molecular ISM, we need to understand what happens to H₂S, why it has a small abundance in the interstellar ices and a large one in cometary ices.

The first and most obvious question to answer is in what conditions H₂S remains frozen on the interstellar ices rather than gaseous. The crucial parameter to answer this question is the H₂S binding energy (BE) in interstellar water-rich ices. We remind that the BE enters in two different ways: (1) it governs at what dust temperature the H₂S solid-to-gas transition occurs (e.g. Minissale et al. 2022; Ceccarelli et al. 2023), and (2) it determines whether H₂S remains frozen or goes to the gas phase when it is formed on the grain surfaces by the hydrogenation of S (as for the case of CO or N hydrogenation; Pantaleone et al. 2020; Ferrero et al. 2023).

Two methods can be used to evaluate the BE of a species: laboratory experiments and theoretical quantum mechanical (QM) calculations (e.g. Brown & Bolina 2007; Chaabouni et al. 2018; Nguyen et al. 2020; Ferrero et al. 2022; Minissale et al. 2022; Molpeceres et al. 2022). One of the very recent progresses on this subject concerns the realization that the adsorption of a species on an amorphous water surface cannot be described by a single BE value but rather a distribution of BEs (e.g. Ferrero et al. 2020; Bovolenta et al. 2022; Heyl, Viti & Vermariën 2023; Tinacci et al. 2023), as BE depends on the adsorption site and molecule orientation. In this work, we use state of the art QM calculations to derive the BE distribution of H₂S, where we adopt a large (200 water molecules) model for the ice (Germain et al. 2022) and computed the BE in 141 different adsorbing sites.

The article is organized as follows. In Section 2, we provide an overview of the methods used in the literature to derive the species BE; in Section 3, we described the adopted methodology in detail; in Section 4, we show the results of the calculations and discuss them in the context of the H₂S BE previously reported in the literature; in Section 5, we discuss the astrophysical implications. Section 6 concludes the article summarizing the major results.

## 2 DETERMINATION OF BINDING ENERGIES IN PREVIOUS LITERATURE

### 2.1 Laboratory experiments

From the experimental point of view, Temperature Programmed Desorption (TPD) is a commonly used technique to obtain BEs. Although TPD experiments represent a reliable procedure to study desorption phenomena, this technique exhibits some drawbacks, especially when one wants to use its results for astrophysical applications. First, BE values extracted from TPD spectra are significantly influenced by the chemical composition and morphology of the substrate, as well as by the coverage regime of the studied species (submono-, mono-, or multilayer) (e.g. Puletti 2014; Collings et al. 2015; He, Acharyya & Vidali 2016). Furthermore, the observable measured through TPD experiments is the desorption temperature, from which the BE and the pre-exponential factor (also called pre-factor) governing the desorption are extracted by means of simple mathematical relationships (e.g. Luna et al. 2017; Penteado, Walsh & Cuppen 2017; Ferrero et al. 2022; Minissale et al. 2022; Ceccarelli et al. 2023). Besides the dependence of both the pre-factor and the BE on the adopted methods, these values are also affected by processes, like ice restructuring, which are commonly regarded as negligible but they are not (e.g. He et al. 2016). Finally, very likely, the BEs derived by TPD are biased towards the high-energy BE sites, thus providing an upper limit of the BE that are likely correct in the cold interstellar conditions (e.g. Minissale et al. 2022; Tinacci et al. 2023), as we will also discuss later in this work (Section 4).

### 2.2 QM calculations

Due to the aforementioned difficulties in simulating the conditions of dense clouds in terrestrial laboratories (e.g. Potapov & McCoustra 2021), QM calculations assume a crucial role to provide BEs. In the last decade, several works have been devoted to theoretically determine BEs, using different surface models and levels of theory (e.g. Wakelam et al. 2017; Das et al. 2018; Ferrero et al. 2020; Bovolenta et al. 2022; Perrero et al. 2022; Tinacci et al. 2022).

Wakelam et al. (2017) reported the BEs of a large set of species at M06-2X/aug-cc-pVTZ level, in which the ice mantle was approximated with a single water molecule. In order to recover the contribution from a real ice surface, a subset of the computed BEs (16 cases) were fitted against experimental data, and all the BEs computed for the other molecules scaled accordingly. A little improvement on such ‘simple’ model was provided by Das et al. (2018), who expanded the number of water molecules to six, using MP2/aug-cc-pVDZ. They demonstrated that enlarging the size of the ice model can yield BEs more in agreement with the experimental results, eliminating the need for a correction coefficient. Bovolenta et al. (2022) utilized 20 amorphous clusters, each consisting of 22 water molecules, thus introducing the effect of different surface morphologies. For each cluster, an adsorbate samples the surface in a random manner leading to 295 sites, thus obtaining a statistical meaningful distribution of BEs. Ferrero et al. (2020) and Perrero et al. (2022) computed a large set of BEs for astrochemically relevant species, using periodic models of either crystalline and amorphous ice. In this case, periodic boundary conditions ensure the rigidity of a macroscopic grain and the hydrogen bond cooperativity. One drawback was that the modelled amorphous water surface was not exhaustively explored with an unbiased procedure, thus possibly affecting the resulted distribution. Finally, a new way of generating an icy grain mimicking the interstellar ones, called ACO-FROST,
was recently presented by some of us (Germain et al. 2022). In the ACO-FROST scheme, the icy cluster is formed by the random addition of molecular water molecules, reaching up to 1000 units. The cluster was constructed by incrementally adding layers of water and subjecting them to short molecular dynamics at 10K. This allowed for the restructuring of the grain due to the energy released by in situ water formation. Then, on top of this grain, a dense grid of points is distributed, constituting the starting positions of the centre of masses of the adsorbates. A full procedure is defined in the works by Tinacci et al. (2022, 2023), which allows computing the BE distribution at the DLPNO-CCSD(T) high level of theory. The ACO-FROST procedure allows sampling the icy grain surface without any user bias, at variance with the previous works by Ferrero et al. (2020) and Perrero et al. (2022). Moreover, the number of molecules composing the icy grain model takes full into account the hydrogen-bond cooperativity and the complex morphology of the amorphous macroscopic ice. Both of these latter crucial effects were missing in the previous cluster models (vide supra), where only limited cluster sizes have been adopted.

In this work, we apply the ACO-FROST composite procedure by Germain et al. (2022) and Tinacci et al. (2023) to compute an accurate and statistically meaningful distribution of H₂S on the top of a 200 water cluster at DLPNO-CCSD(T) level of theory.

# 3 METHODOLOGY

## 3.1 BE definition

BE is defined as a negative of the interaction energy ($\Delta E$), which in turn is the difference between the complex total energy ($E_c$) and the sum of the isolated adsorbate ($E_{ads}^{iso}$) and the isolated grain total energies ($E_{gra}^{iso}$), hence:

$$BE = -\Delta E = E_{ads}^{iso} + E_{gra}^{iso} - E_c, \qquad (1)$$

$$\Delta E = E_c - E_{ads}^{iso} - E_{gra}^{iso}. \qquad (2)$$

BE can be decomposed in two different terms: (i) the electronic interaction ($BE_e$) corrected for the Basis Set Superposition Error (where $E_{ads}^{iso}$ and $E_{gra}^{iso}$ are calculated at the geometry of the complex), and (ii) the deformation energy contribution ($\delta E_{def}$) arising from the geometrical change of the isolated systems due to their mutual interaction. By subtracting to these terms the zero point energy ($\Delta ZPE$) correction, derived from the harmonic frequency calculation, the BE at 0 K can be defined as

$$BE = \underbrace{BE_e - \delta E_{def}}_{BE*} - \Delta ZPE \qquad (3)$$

For the ease of visualization and interpretation, we refer to the BE corrected for the ZPE as BE, while the same value not corrected is labelled as BE*. The $\Delta ZPE$ is computed as

$$\Delta ZPE = ZPE_c - ZPE_{ads}^{iso} - ZPE_{gra}^{iso}. \qquad (4)$$

## 3.2 Computational details

The grain model and the initial grain-adsorbate structures were first optimized at GFN2-xTB level (Grimme, Bannwarth & Shushkov 2017; Bannwarth, Ehlert & Grimme 2019) to get a preliminary set of BEs. Each preliminary structure was refined using the ONIOM method (Mayhall, Raghavachari & Hratchian 2010) (QM:QM2) implemented in the ORCA program (v.5.0.3)(Neese 2022). Within the ONIOM method, the structure is divided into two fictitious

subsystems: (i) the Real Zone, which is the whole system treated at lower level of theory (QM2); (ii) the Model Zone, which is a smaller part of the structure ensuring that there are sufficient water molecules present to account for the significant local interactions of the adsorbate, treated with a higher level of theory (QM). The final ONIOM energy for each structure is obtained as

$$E(\text{ONIOM}) = E_{\text{Real}}(\text{QM2}) - E_{\text{Model}}(\text{QM2}) + E_{\text{Model}}(\text{QM}). \qquad (5)$$

The B97-3c (Brandenburg et al. 2018) was selected as high level method (QM) for the Model zone in order to compute geometries and harmonic frequencies, while GFN2-xTB was set as low level method (QM2). The atoms outside the Model region are kept fixed and the Model zone is embedded electrostatically (default setting). Tight convergence criteria were set for both the SCF and geometry optimizer, while for the integral evaluation the grid was set to the highest density (keywords: !tightSCF, !tightOPT, and !defgrid3 in the ORCA program). In the few cases in which imaginary frequencies were found, the geometry is restarted using a tighter convergence criterion on the optimization (i.e. !verytightopt), mostly leading to the correct minima. After the optimization procedure, energies were refined applying the DLPNO-CCSD(T) (Guo et al. 2018) as QM method in the ONIOM scheme. Along with this, the aug-cc-pVTZ was selected as primary basis set (Kendall, Dunning & Harrison 1992; Woon & Dunning 1993), while the aug-cc-pVTZ/C (Weigend, Köhn & Hättig 2002) is used as auxiliary basis set for the resolution of the identity approximation. DLPNO-CCSD(T) calculations were performed with a tight-PNO set-up and using the default settings on the SCF convergence. Finally, the energies of the complex ($E_c$(DLPNO-CCSD(T)) were corrected for the Basis Set Superposition Error applying the Counterpoise method (Boys & Bernardi 1970).

A complete benchmark on the H₂S ... H₂O dimer and trimer is available in the Annex (see Table B1).

## 3.3 Water-ice grain model and procedure to compute BEs

The initial sampling and computation of BE is performed following our previous work (Germain et al. 2022). A cluster of 200 water molecules is used as model to mimic the interstellar icy grains (see Fig. 1a). Using the ACO-FROST code, a spherical grid of 162 points evenly distributed around the grain is created. Each of these points is replaced by H₂S molecule with three different random orientations (see Fig. 1b) generating a random distribution of 486 BE sites. After that, the centre of mass of each molecule is projected at a distance of 2.5 Å from the grain surface (see Fig. 1c). The procedure to compute the H₂S BE distribution consists of six subsequent steps:

(i) The geometry of the grain is kept fixed while each individual H₂S molecule is freely left to relax at GFN2 level.

(ii) With the optimized geometry from point (1), all the molecules within a 5 Å radius from the adsorbate centre of mass, adsorbate included, are relaxed again at GFN2 level.

(iii) On the optimized structure from (2), a geometrical check is carried out in order to verify that the number of water molecules has not changed in the 5 Å region, due to the geometry optimization. If this is the case, step (2)-(3) are repeated until convergence on the number of molecules included in the movable region is reached.

(iv) On the previous converged structure, the ONIOM(QM:QM2) model is applied. The Model zone is defined as the 5 Å region from the adsorbate and it also defines the only atoms free to relax (the same as for the GFN2 screening).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Model of the water ice grain and its adsorption sites. a) 200 water molecules grain created with the ACO-FROST procedure (Germain et al. 2022); b) initial grid of 486 points reported as yellow spheres; c) final grid projected onto the grain at 2.5 Å from the surface. Red and white spheres are oxygen and hydrogen atoms, respectively, while yellow spheres are the centres of mass of the H₂S molecule.</figcaption></figure>

(v) Harmonic frequencies are computed only for the atoms inside the QM region, while the other nuclei are kept fixed.

(vi) Finally, the isolated grain surface reference ($E_{gra}^{iso}$) is defined following our previous approach (called ‘TPD’) which somehow mimic the physics of the experimental TPD process (see Tinacci et al. 2022), relaxing the QM region after the desorption of the hydrogen sulphide.

## 3.4 Procedure to obtain unique BE sites

In order to avoid redundancies in the automatic BE site sampling, we developed a pruning procedure. With respect to the previous one reported by our group (Tinacci et al. 2022), we apply the following pruning criteria on the initial GFN2-xTB geometries, instead of the ONIOM optimized ones (see Fig. A2). This new expedient relies on the experience acquired with the previous works from Tinacci et al. (2022) and Tinacci et al. (2023), thus reducing the computational effort. The comparison between the two pruning procedures is reported in the Annex for the case of the H₂O molecule. Finally, the redundant structures are picked out in the following manner: (i) all the combinations between Model zone with the same number of atoms are taken into account; (ii) for each combination, we align the two structures and compute the RMSD and the absolute value of the energy difference $|\Delta E_{Cij}| = |E_{Ci} - E_{Cj}|$ between the i-th and j-th complex structure; (iii) in the end, if the two structures have an RMSD below 0.3 Å and $|\Delta E_{Cij}| < 1.5$ kJ mol⁻¹ they are considered identical (see the Annex for a detailed discussion on the selection of this criteria).

# 4 RESULTS AND DISCUSSION

## 4.1 BE distribution

After the pruning procedure, 338 redundant BE sites were found, returning a total number of 148 unique BE sites. Among these structures, only two BE sites showed one imaginary frequency, each in the range of [-32, -16] cm⁻¹. Since the imaginary frequencies span a low wave-numbers range, we decided to keep these sites in order to improve our statistical sample. In addition, after DLPNO-CCSD(T) refinement, we obtained seven negative BEs that were excluded from the sample, reducing the total number of BE sites to 141. The excluded BEs are minima of the PES according to frequency calculations and all of them are positive values at B97-3c

<figure><img src="image_4.png" /><figcaption>**Figure 2.** BE distribution obtained at DLPNO-CCSD(T) level for the final 141 BE sites (all positive). Geometries and ZPEs are evaluated at the ONIOM(B97-3c:xTB-GFN2) level. The relative population for each bin is reported on the top of the bars. The dotted blue line is the unnormalized Gaussian best-fitting function (see equation 6).</figcaption></figure>

level. However, all the corrections to the electronic energy applied to improve the quality of the results (i.e. ZPE, DLPNO-CCSD(T)/aug-cc-pVTZ and its corresponding BSSE correction) tend overall to shift the B97-3c distribution towards lower values. We decided to delete the seven cases with negative BEs because they do not have a physical meaning and their omission does not alter the distribution. The final BE distribution (equation 3) is reported in Fig. 2 and Table 1. The number of bins and their width are obtained following the Freedman Diaconis estimator (Freedman & Diaconis 1981). With respect to the previous work on H₂S by Perrero et al. (2022), the distribution spans a lower BE range, and shows an asymmetric shape as highlighted by the truncated Gaussian fit reported in Fig. 2:

$$G = A \exp \left( -\frac{(x - \mu)^2}{2\sigma^2} \right), \qquad (6)$$

where A is the pre-factor equal to 29.3, the mean value μ is 7.5 kJ mol⁻¹ (equivalent to 900 K), and the standard deviation σ is 5 kJ mol⁻¹ (equivalent to 600 K).
<table><caption>**Table 1.** Relative fraction of ice with a given BE value and the associated pre-factor computed using the Tait et al. (2005) recipe. The fractions were identified with the number of bins obtained in the distribution using the Freedman Diaconis estimator (Freedman & Diaconis 1981), as in Fig. 2. The given BE value is the average value of each bin, while the pre-factor was computed at *T*<sub>peak</sub> (see the main body).</caption><thead><tr><th>BE<br>kJ/mol (K)</th><th>Pre-factor<br>(s<sup>-1</sup>)</th><th>Fraction<br>of the ice</th></tr></thead><tbody><tr><td>1.7 (204)</td><td>5.25 × 10<sup>11</sup></td><td>0.13</td></tr><tr><td>4.1 (497)</td><td>8.48 × 10<sup>12</sup></td><td>0.15</td></tr><tr><td>6.6 (790)</td><td>3.64 × 10<sup>13</sup></td><td>0.18</td></tr><tr><td>9.2 (1083)</td><td>1.02 × 10<sup>14</sup></td><td>0.23</td></tr><tr><td>11.5 (1376)</td><td>2.14 × 10<sup>14</sup></td><td>0.16</td></tr><tr><td>13.9 (1669)</td><td>3.82 × 10<sup>14</sup></td><td>0.07</td></tr><tr><td>16.4 (1961)</td><td>6.48 × 10<sup>14</sup></td><td>0.05</td></tr><tr><td>18.8 (2254)</td><td>9.97 × 10<sup>14</sup></td><td>0.02</td></tr></tbody></table>

## 4.2 Energetic trends

As discussed previously, the BE is a multicomponent term that can be decomposed in order to better understand the origin of this value, and hence the shape of the distribution. In Fig. 3, the correlation plot between BE with several important constituents such as $BE_e$ (BSSE) (electronic BE without BSSE correction), $\delta E_{def}$, BSSE, and $\Delta$ZPE is shown. Along with, each point is associated with the number of water molecules in the QM zone and sketched in the coloured bar. The main outcomes for each plot can be summarized as follows:

(i) The $BE_e$ (BSSE) is in good correlation with the BE. As expected, $BE_e$ (BSSE) values are all higher than BE. Moreover, $BE_e$ (BSSE) does not depend on the number of water molecules, corroborating the fact that the 5 Å criteria identified by the previous work of our group (Tinacci et al. 2022) represents a converged value, which does not bias the BEs.

(ii) The $\delta E_{def}$ contribution is rather sparse, the majority at low values (<7 kJ mol$^{-1}$), and the BE slightly increases with the $\delta E_{def}$. In the same manner, $\delta E_{def}$ roughly correlates with the number of water molecules, indicating that the morphology of the BE site, such as small cavities, is more important with respect to the number of included molecules.

(iii) BSSE contribution meagrely correlates with BE showing an increasing trend. This is because BSSE is a geometrical factor, meaning that the closer the adsorbate to the surface, the higher the BE and the BSSE correction. Finally, as expected, BSSE correlates with the number of H₂O molecules in the Model Zone.

(iv) $\Delta$ZPE does not show any particular correlation with the site energy or other quantities, exhibiting the smallest values in the lower part of the distribution. In the same way, H₂O molecules do not affect the $\Delta$ZPE contribution, which is in agreement with previous observations in previous works by some of us (e.g. Zamirri et al. 2017). Moreover, internal test on the H₂O ... H₂S interaction was conducted to assess the role of the anharmonicity in the calculation of the ZPE. The results revealed no significant variations with respect to the harmonic $\Delta$ZPE values.

## 4.3 H₂O ... H₂S interaction

Despite the hydrogen bond (H-bond) is commonly defined for second-period elements (N, O, F), also the atoms of the corresponding groups of the third period (i.e. P, S, Cl) can either behave as H-bond donors and acceptors. In our specific case, the electronegativity of sulphur is lower than oxygen, and accordingly

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Correlation plot between BE and several contributions: electronic BE not corrected for BSSE ($BE_e$ (BSSE)), deformation energy ($\delta E_{def}$), Basis Set Superposition Error correction (BSSE), and zero point energy ($\Delta$ZPE). As in equation (3), $BE = BE_e(BSSE) - \delta E_{def} - BSSE - \Delta$ZPE. All values are in kJ mol$^{-1}$ and presented at DLPNO-CCSD(T)/aug-cc-pVTZ//B97-3c level. The colour map shows the number of H₂O molecules included in the Model region.</figcaption></figure>

also its tendency to engage in H-bonding. In Table B1, the interaction between H₂S and H₂O with golden standards methods is reported. As expected, the H₂S and H₂O interaction calculated with the golden standard CCSD(T)/aug-cc-pvtz method is 10.5–11.5 kJ mol$^{-1}$, a factor of 2 lower than H₂O dimer (ca 21 kJ mol$^{-1}$; Halkier et al. 1997). Lamberts & Kästner (2017) carried out a benchmark on all the possible combinations of H₂O and H₂S, finding the following hierarchy of binding energies: H₂O ... H₂O (21 kJ mol$^{-1}$) > H₂O ... H₂S (11.3–12.3 kJ mol$^{-1}$) > H₂S ... H₂S (6.9 kJ mol$^{-1}$). H₂S does not show any preference as H-bond donor or acceptor, indicating a less specific interaction, with respect to second-row elements (i.e. O, N, F) (McQuarrie 2008; Atkins, Atkins & de Paula 2014; Chand & Biswal 2019). This energetic change has a one-to-one correspondence with the specificity of the H-bond geometry. Therefore, a deep geometrical analysis was performed to better understand the subtle chemistry of S-bearing molecules.
To this end, in the BE distribution we identified three clusters based on the number of H-bonds established between H₂S and H₂O (i.e. HSH ... OH₂ and/or H₂S ... HOH), reported in Fig. 4. The H-bond cut-off criteria that we used, selected after an extensive visual check on the structures, are the following: ranges of 1.7–2.9 Å for the H-bond length, and 140°–180° for the H-bond angle. In Fig. 4, the total distribution is divided into three different subclasses based on the number of H-bonds involving the H₂S, namely three H-bonds, two H-bonds, and one H-bond. At each distribution, we associated the correlation plot between distance and angle for all the H-bonds in the specific subclass. The first panel of Fig. 4 shows that only 14 BE sites are involved in 3 H-bonds, with a BE average value of 12.8 kJ mol⁻¹, spanning from 1.4 to 20 kJ mol⁻¹. The most populated class is the one in which H₂S form two H-bonds (the second panel of Fig. 4) with the water surface, with a total number of 64 BE sites. This class goes from 0.7 to 17.5 kJ mol⁻¹, with a mean value of 9.5 kJ mol⁻¹. Finally, the last subclass is reported in the third panel of Fig. 4, illustrating the energy contribution arising from the H₂S adsorbed only with one H-bond. This class spans the lowest part of the distribution with a range of 13 kJ mol⁻¹ and a mean value of 5.8 kJ mol⁻¹. All these outcomes are in agreement with the expectation that the higher the number of H-bonds, the stronger the resulting BE.

As a general trend for all the three classes, when H₂S acts as H-bond acceptor (circle marker in the right side of Fig. 4), the majority of the cases is localized in the bond length range between 2.30 Å and 2.55 Å and spans an angle from 150° to 180°. On the contrary, when it is a H-bond donor (triangle marker in the right side of Fig. 4), the points are more scattered, and no clear trend is found. In particular when H₂S is engaged in three H-bonds (top panel of Fig. 4), H-bond donor distances spread all over the chart, indicating that some geometries are borderline in the geometric definition of H-bond, and this also justifies the presence of few low-BE values in the distribution, even for this class, where one may expect to found only strong BEs. Moving to the central-right panel of Fig. 4, H-bond donor interactions are still spread, even if the larger amount of data allows us to find a trend between the H-bond distance and the angle: the shorter the H-bond length, the more linear the angle (i.e. the geometry gets closer to the definition of H-bond). Finally, when H₂S establishes only one H-bond (bottom-right panel of Fig. 4), in only one case out of 61 it acts as H-bond acceptor. Moreover, data are clusterized at shorter distances (2.00–2.25 Å), while the angle is always widely spread (150°–180°), indicating that they are ‘true’ H-bonds.

## 4.4 Comparison with previous theoretical results

In this section, we compare our results with the existing ones in the literature, highlighting the differences and the improvements with respect to previous simulations. As a general consideration, there are three parameters that affect the final BE of a certain species: (i) the choice of the model surface, (ii) the utilized sampling procedure, and (iii) the specific employed computational method. Our computed BE values range from 57 K (0.5 kJ mol⁻¹) to 2406 K (20.0 kJ mol⁻¹), the mean is μ = 900 K (7.5 kJ mol⁻¹), and its standard deviation is σ = 600 K (5.0 kJ mol⁻¹). Table 2 lists the BE values obtained by previous theoretical works following the order of increasing complexity for what concern the ice model, to compare with those presented here.

The first group of theoretical works considered a very limited number (≤7) of water molecules to simulate the ice surface. Wakelam et al. (2017) used one single H₂O molecule and reported a computed

value of 1565 K (13 kJ mol⁻¹), not corrected for BSSE and ZPE, and then fitted against the experimental results of 2743 K to report a final value of 2700 K. Similar outcomes are reported by Lamberts & Kästner (2017) for the H₂S ... H₂O interaction, here reported without including the ZPE correction.

Das et al. (2018) reported values with an increasing number of water molecules from 1 to 4. The water monomer BE values are in agreement with those by Wakelam et al. (2017) and Lamberts & Kästner (2017); however, the BE values increase with the number of considered water molecules up to 27 kJ mol⁻¹, indicating that many more molecules are needed to obtain converged BEs.

Although the previous reported values are comparable with our benchmark reported in Table B1 on the H₂S ... H₂O interaction (BE* = not corrected for ZPE), the size of the model used and the lacks of the corrective procedure, make the comparison misleading and for these reasons we do not compare further.

Oba et al. (2018) computed the BEs on a cluster of seven H₂O molecules. Calculations were carried out using MPWB1K/def2-TZVP, which was benchmarked against the same level of accuracy as the one used in this work, and corrected for the ZPE. In this case, the two derived values fall within the high-BE part of our distribution (794–1936 K).

While the previous works considered a relatively limited number of water molecules to model the ice surface, in Bovolenta et al. (2022) several clusters, each consisting of 22 water molecules, were used to obtain a BE distribution.

Their data, μ = 1794 K and σ = 422 K, are shifted towards larger BEs values with respect to our distribution (μ = 984 K and σ = 532 K; Fig. 2). The difference can be attributed to several factors, including variations in the employed icy surface model (small clusters) and the level of theory utilized. The latter is likely more relevant, as the geometries by Bovolenta et al. (2022) were computed at a lower level of accuracy, namely HF-3c, which can largely deviate from geometries obtained at higher level of theory, as shown in the benchmark reported in Table B1.

A different method to model the ice surface was used by some of us, in the sense of increasing the complexity of the substrate while keeping at the same time a reasonable computational cost. To this end, Ferrero et al. (2020) and Perrero et al. (2022) used a periodic model of amorphous ice made up by 60 water molecules, thus mimicking the rigidity and the H-bond cooperativity of a macroscopic ice.

BE values obtained by Ferrero et al. (2020) and Perrero et al. (2022) are shown in Table 2. Despite they share the same structural model and computational methodology the two works provide different BE values. This is due to the fact that the binding sites were selected by the authors chemical intuition, thus biasing the sampling procedure and number of samples used (5 and 8, respectively). Moreover, the ZPE and CCSD(T) corrections were applied only on a subset of crystalline systems and then extrapolated to all the other cases (amorphous). In Appendix B2, we present an accurate description of the error in the extrapolation of these quantities. Here, we summarize the main outcomes: the error on the extrapolated ZPE is 4 kJ mol⁻¹ and on the CCSD(T) is 9 kJ mol⁻¹, a non-negligible quantity considering the mean value of our B97-3c distribution (16 kJ mol⁻¹, see Fig. B2). A general warning rises from this discussion, in particular about using fitting procedures, either from experimental, but also from computed data (even if calculated at the same level of theory of the value to be corrected with); this may strongly bias in particular small quantities as BE of weakly bound systems (in practice the large majority of small molecules interacting with ice surfaces).
<figure><img src="image_6.png" /><figcaption>Figure 4. Correlation plot between different sub-distributions based on the number of H-bond established by H₂S and the interaction peculiarity of each BE site. On the left panel column the total distribution is decomposed in three differently populated classes based on the number of H-bonds established by the H₂S. From the top to the bottom, three H-bonds, two H-bonds, and one H-bond, respectively. On the right hand column, the associated correlation plot between the H-bond distance and the angle formed by the bond is reported. Different markers indicate H₂S acting as H-bond donor (triangle) or acceptor (circle). The colour map indicates the BE energy reported in kJ mol⁻¹.</figcaption></figure>

## 4.5 Comparison with experimental results

Direct comparison between computed BEs and those derived from TDP experiments is a complex task, as already briefly mentioned above. One of the main differences between the two approaches is that, while computer simulations allow directly measuring the heat of desorption of a single molecule from the grain (BE), TPD experiments report desorption temperature curves and peaks ($T_{peak}$), and the corresponding BE and pre-exponential factor are derived numerically by fitting the data with the Polanyi-Wigner equation (Polanyi & Wigner 1925; Minissale et al. 2022). In addition, there have been cases in the literature in which the BE is extracted from $T_{peak}$ through a linear relationship, where the BE of water is assumed to be the reference (Penteado et al. 2017).

With this caveat in mind, we computed the TPD curves from our BE distribution and compared them with the $T_{peak}$ quoted by Oba et al. (2018) and Collings et al. (2004), respectively. Our simulated TPD curves were obtained using the Polanyi-Wigner equation with the surface heating rate $\beta = \frac{dT}{dt}$ (King 1975):

$$-\frac{dN}{dT} = \frac{N}{\beta}v(T)\exp\left(-\frac{BE}{RT}\right), \qquad (7)$$

where $\beta$ was kept fixed at 5 K min⁻¹, the same condition used by Oba et al. (2018) and very similar to the one used by Collings et al. (2004) (i.e. 4 K min⁻¹). The pre-factor term, $v(T)$, was computed following recipe from Tait et al. (2005). However, we verified that the simulated TPD curve is very similar if the recipe in Hasegawa, Herbst & Leung (1992) is used.
<table><caption>**Table 2.** Comparison of our results with previous theoretical works. First column lists the BE values in K (values in parenthesis are given in kJ mol⁻¹), which have been obtained by the different authors, cited in the last column, using different ice models (second column), methods (third column) and corrective procedure applied (fourth column).</caption><thead><tr><th>BE</th><th>Ice model</th><th>Method</th><th>BSSE/ZPE corrections</th><th>Reference</th></tr></thead><tbody><tr><td>2700 (22.5)</td><td>Water monomer</td><td>M06-2X/aug-cc-pvtz</td><td>No/No</td><td>1</td></tr><tr><td>1479 (12.3)</td><td>Water monomer</td><td>UCCSD(T)-F12/cc-VTZ-F12 // MPWB1K/def2-TZVP</td><td>No/No<sup>a</sup></td><td>2</td></tr><tr><td>1727 (14.4)</td><td>Water monomer</td><td>MP2/A-VDZ</td><td>No/No<sup>a</sup></td><td>3</td></tr><tr><td>2556 (21)</td><td>Water Tetramer</td><td></td><td></td><td></td></tr><tr><td>3232 (27)</td><td>Water Hexamer</td><td></td><td></td><td></td></tr><tr><td>Min: 794 (6.6)</td><td>Water Heptamer</td><td>MPWB1K/def2-TZVP</td><td>No/Yes</td><td>4</td></tr><tr><td>Max: 1936 (16)</td><td></td><td></td><td></td><td></td></tr><tr><td>Max: 2940 (24)</td><td>Amorphous clusters (22 H<sub>2</sub>O)</td><td>ω - PBE/def2-TZVP // HF-3c</td><td>Yes/Yes</td><td>5</td></tr><tr><td>μ: 1794 (15)</td><td></td><td></td><td></td><td></td></tr><tr><td>σ: 422 (3.5)</td><td></td><td></td><td></td><td></td></tr><tr><td>Min: 2291 (19)</td><td>Amorphous periodic (60 H<sub>2</sub>O)</td><td>B3LYP-D3/A-VTZ<sup>b</sup></td><td>Yes/Yes</td><td>6</td></tr><tr><td>Max: 3338 (28)</td><td></td><td></td><td></td><td></td></tr><tr><td>μ: 2753 (23)</td><td></td><td></td><td></td><td></td></tr><tr><td>σ: 437 (3.7)</td><td></td><td></td><td></td><td></td></tr><tr><td>Min: 1970 (16)</td><td>Amorphous periodic (60 H<sub>2</sub>O)</td><td>B3LYP-D3/A-VTZ<sup>b</sup></td><td>Yes/Yes</td><td>7</td></tr><tr><td>Max: 4489 (37)</td><td></td><td></td><td></td><td></td></tr><tr><td>μ: 3017 (25)</td><td></td><td></td><td></td><td></td></tr><tr><td>σ: 884 (7.4)</td><td></td><td></td><td></td><td></td></tr><tr><td>Min: 57 (0.5)</td><td>Iced grain (200 H<sub>2</sub>O)</td><td>ON(DLPNO-CCSD(T)/aug-cc-pvtz:GFN2) // ON(B97-3C:GFN2)<sup>c</sup></td><td>Yes/Yes</td><td>This work</td></tr><tr><td>Max: 2406 (20)</td><td></td><td></td><td></td><td></td></tr><tr><td>μ: 984 (8.2)</td><td></td><td></td><td></td><td></td></tr><tr><td>σ: 532 (4.4)</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. <sup>a</sup>The authors reported also the value corrected for the ZPE but the one not corrected is reported in this table. <sup>b</sup>A-VTZ refers to the Ahlrichs triple zeta basis set augmented with two polarization functions. <sup>c</sup>ON refers to ONIOM calculation method (Mayhall et al. 2010). References: [1] Wakelam et al. (2017); [2] Lamberts & Kästner (2017); [3] Das et al. (2018); [4] Oba et al. (2018); [5] Bovolenta et al. (2022); [6] Ferrero et al. (2020); [7] Ferrero et al. (2022).

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between simulated and experimental TPD spectra. The dashed lines represent the $T_{\text{peak}}$ of the TPD experiments reported by Collings et al. (2004) and Oba et al. (2018), respectively. The continuous lines represent the summed TPD spectra for each value of BE computed in this work, fitted with a Gaussian function.</figcaption></figure>

As evident from Fig. 5, our simulated TPD curve provides a $T_{\text{peak}}$ at about 30 K, which is much lower than the experimental values: ~80 K (Collings et al. 2004) and ~100 K (Oba et al. 2018). The first possible explanation of this large differences is that, as already emphasized, the experimental and computational methods actually produce different quantities (see also Ferrero et al. 2022; Tinacci et al. 2023). Specifically, TPD experiments measure the desorption

rate from an ice continuously heated, which enables the adsorbate to explore all available binding sites before eventually being desorbed (see, for instance, Minissale et al. 2022). Consequently, in submonolayer coverage regimes, the abundance of binding sites exceeds the number of adsorbates, allowing them to diffuse on to the surface. This facilitates a predominant interaction with the most favourable binding sites, thereby establishing an upper limit to the BE value. Furthermore, our methodology entails the interaction of a single H<sub>2</sub>S molecule with a cluster of water molecules, thus operating within an ultra-low coverage regime. Conversely, experimental results may be influenced by the wetting phenomenon, wherein small islands of mutually interacting H<sub>2</sub>S can form even in the submonolayer regime. During TPD experiments, lateral H<sub>2</sub>S···H<sub>2</sub>S interactions, in addition to the H<sub>2</sub>S/ice BE, can significantly elevate the final desorption temperature. Taken together all these effects, unaccounted for in our simulations, could explain the observed shift of the simulated TPD curve towards lower BE values.

Nevertheless, the discrepancy between experiments and simulations prompts the question of which technique is more reliable in simulating the ISM conditions and, hence, which one is better in providing BE for astrochemical numerical models. In the cold ISM where ices are formed, all binding sites are likely populated due to the extremely low temperature, including those mapped by the theoretical calculations with the lowest BE values. Additionally, the ultra-low coverage assumed in our simulations may better reflect the low-density conditions inside dense clouds. In our opinion, therefore, theoretically BE computations may be more realistic in the description of the interstellar ices and, hence, computed BE may provide more reliable values to be used in astrochemical models.
<table><caption>**Table 3.** List of the systems where *ab initio* molecular dynamic simulations have provided the fraction of the reaction energy absorbed by the water molecules of the icy surface. The last line reports the HS + H reaction energy only, for comparison, computed with the same methods adopted for studying the H₂S BE.</caption><thead><tr><th>System</th><th>Reaction energy (kJ mol<sup>-1</sup>)</th><th>Ice absorbed energy (per cent)</th><th>Reference</th></tr></thead><tbody><tr><td>H + H → H<sub>2</sub></td><td>-436 - -440</td><td>45-65</td><td>1</td></tr><tr><td>CO + H → HCO</td><td>-132</td><td>90</td><td>2</td></tr><tr><td>N + H → NH</td><td>-343</td><td>76</td><td>3</td></tr><tr><td>NH + H → NH<sub>2</sub></td><td>-382</td><td>71</td><td>3</td></tr><tr><td>NH<sub>2</sub> + H → NH<sub>3</sub></td><td>-441</td><td>58-80</td><td>3</td></tr><tr><td>HS + H → H<sub>2</sub>S</td><td>-369</td><td></td><td>4</td></tr></tbody></table>

*References:* [1] Pantaleone et al. (2021); [2] Pantaleone et al. (2020); [3] Ferrero et al. (2023); [4] This work.

## 5 ASTROCHEMICAL IMPLICATIONS

The major result of our computations is the relatively low values of the H₂S BE, lower than the previous estimates obtained theoretically and experimentally (see above). This has a major implication for the composition of S-bearing ices. As mentioned in the Introduction, astrochemical models predict H₂S to be the most abundant S-bearing frozen species (e.g. Vidal et al. 2017; Navarro-Almaida et al. 2020). At odd with the predictions, H₂S has never been detected in ices, clearly indicating that it is not a major frozen sulphur carrier (e.g. Boogert et al. 2015; McClure et al. 2023). Various hypothesis have been evoked in the literature to explain this discrepancy, including that H₂S undergoes chemical reactions on the grain surfaces (e.g. Laas & Caselli 2019).

A breakthrough occurred with the detection of a large desorption of H₂S upon hydrogen bombardment of HS and H₂S deposited on amorphous solid water (Oba et al. 2018). However, as discussed in Section 4.5, these authors obtained a TPD curve with a $T_{peak} \sim 100$ K, much larger than those of our BE distribution (Fig. 5), while their theoretical estimates lay in the upper end of our BE distribution (Section 4.4).

Our new H₂S BE distribution actually can explain the Oba et al. (2018) experimental results and also why H₂S is not detected in the interstellar ices. In previous works, *ab initio* molecular dynamics (AIMD) simulations have shown that the energy released by chemical reactions may trigger the diffusion and the desorption of the newly formed species even at very low time-scales (of order of few picoseconds Pantaleone et al. 2020, 2021; Ferrero et al. 2023). We list the results of these computations in Table 3. As it can be seen, the percentage of energy absorbed may depend on the reaction and on the site where the reaction occurs (of order of few picoseconds Ferrero et al. 2023). The remaining energy is conserved by the newly formed molecule, which may or may not leave the surface, depending whether its energy is larger or not than its BE and on the coupling between the internal and translational energies (e.g. Fredon, Lamberts & Cuppen 2017; Fredon, Radchenko & Cuppen 2021). In practice, three parameters enter in the chemical desorption probability: the reaction energy, the percentage which is absorbed by the surface and the BE of the newly formed species. The lower the latter the higher is the chemical desorption probability.

For example, Pantaleone et al. (2021) found that H₂ is liberated into the gas-phase because the energy of the newly formed molecule is about ten times larger than its BE. On the contrary, HCO, formed

by CO + H, does not leave the surface because the molecule energy is lower than its BE (Pantaleone et al. 2020).

In the case of the HS hydrogenation (with a computed reaction energy of 369 kJ mol⁻¹), one can conservatively assume that more than 10 per cent of the reaction energy remains to the formed H₂S, which is about 37 kJ mol⁻¹. This energy is an almost factor 5 larger than the average H₂S BE of our distribution and twice its highest BE. Therefore, the newly formed H₂S will very likely leave the surface after its formation from the HS hydrogenation, in agreement with the experiments and also with the lack of detection of frozen H₂S in interstellar ices. In order to give more precise estimates of the desorption percentage of formed H₂, *ad hoc* AIMD simulations are required at different BE sites, but our crude estimate suggests that the vast majority of newly formed H₂S is doomed to go into the gas phase soon after being formed on the grain. Once in the gas phase, H₂S reacts with the most abundant cations in the gas (e.g. H₃⁺ and H₃O⁺ in cold and warm gas, respectively), becomes protonated and then the electron recombination of H₃S⁺ would lead to SH that reacts with O to produce SO, which is indeed an abundant molecule in molecular clouds (e.g. Wakelam et al. 2004).

Of course, the H₂S formation is a two-step process where also the first one, S + H could be subject to strong chemical desorption leaving little HS on the grain surfaces and, consequently, even less H₂S. As reported in Table D1, S + H has a reaction energy of 339 kJ mol⁻¹ which, again assuming the conservative 10 per cent, as for H₂S, would release more than 33 kJ mol⁻¹ on the formed product. Considering that HS is expected to have a smaller BE than H₂S (Wakelam et al. 2017; Das et al. 2018; Ferrero et al. 2022), the amount of retained energy will be very likely enough for the desorption of the radical molecule.

In the end, we notice that simulations obtained via a microscopic kinetic Monte Carlo method (which uses classical dynamics) suggest a probability of ~3 per cent (Furuya, Oba & Shimonishi 2022) and a percentage of reaction energy remaining in the H₂S product of 7 per cent. Nevertheless, these simulations adopted a very high BE, namely 2700 K (Collings et al. 2004), in addition to lacking a precise atomistic treatment, so that their probability may be largely underestimated. As said, it would be very important to carry out AIMD simulations in a sample of sites covered by our newly computed BE distribution and, based on that, estimate the chemical desorption probability.

## 6 CONCLUSIONS

In the present work, we report a detailed computational investigation on the BE of H₂S adsorbed on a large icy model comprising 200 water molecules representative of the ice mantle of the interstellar grain dust. The computer model was built following the ACO-FROST recipe (Germain et al. 2022). A multistep procedure, similar to the one reported by Tinacci et al. (2023) for the adsorption of NH₃ on ice, was adopted here to refine the computed BEs at DLPNO-CCSD(T) level of theory. Starting from the initial generated 486 adsorbing sites, a pruning procedure gives 141 unique binding sites, whose corresponding BE values define a distribution limited from 57.4 K (0.5 kJ mol⁻¹) to 2406.3 K (20.0 kJ mol⁻¹) and with $\mu = 984.4$ K (8.2 kJ mol⁻¹).

The resulting distribution spans a new range of low BEs never reported before, neither by TPD experiments nor theoretical calculations. The new H₂S BE distribution may finally explain its absence in the interstellar ices. In that respect, our recent AIMD results for the formation of H₂, HCO, and NH₃ on icy grain models (Pantaleone et al. 2020, 2021; Ferrero et al. 2023), showed that more than 70
per cent of the formation energy is transferred to the phonons of the ice, while about 30 per cent of it remains as rotational, translational, and vibrational modes of the adsorbate. Even reducing this fraction to 10 per cent for the process of H₂S formation (formation energy of 369 kJ mol⁻¹) on the icy grain, the corresponding surplus energy (ca. 37 kJ mol⁻¹) exceeds the H₂S BE distribution's mean value by an almost factor of 5 and its upper limit by a factor of 2. Consequently, the newly generated H₂S will likely desorb from the ice surface after its formation. This is consistent with the fact that, to date, not even JWST has detected frozen H₂S in interstellar ices.

That said, due to the pivotal importance of H₂S as a sink for S-bearing species, *ad hoc* AIMD simulations should be carried out to establish how the H₂S formation energy is partitioned between the substrate and the molecule itself and, therefore, to more precisely estimate its chemical desorption probability.

## ACKNOWLEDGEMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation programme 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, the Marie Skłodowska-Curie for the project 'Astro-Chemical Origin' (ACO), the Italian MUR (PRIN 2020, Astrochemistry Beyond the Second Period Elements, Prot. 2020AFB3FX), and the Spanish MICIN (project PID2021-126427NB-I00). We acknowledge the CINECA award under the ISCRA initiative (project HP10BGB202: BE distribution of interstellar molecules beyond the second period of the periodic table: a computational approach), for the availability of high performance computing resources and support. The authors acknowledge also support from Project CH4.0 under the MUR programme 'Dipartimenti di Eccellenza 2023-2027' (CUP: D13C22003520001). V.B. expresses deep gratitude to Lorenzo Tinacci and Aurèle Germain for their support and assistance.

## DATA AVAILABILITY

The data underlying this article are freely available at Zenodo.

## REFERENCES

* Altwegg K. et al., 2022, MNRAS, 516, 3900

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Atkins P., de Paula J., 2014, Atkins' Physical Chemistry. Oxford Univ. Press, Oxford

* Bannwarth C., Ehlert S., Grimme S., 2019, J. Chem. Theory Comput., 15, 1652

* Biver N. et al., 2021, A&A, 648, A49

* Blake G. A., van Dishoeck E. F., Jansen D. J., Groesbeck T. D., Mundy L. G., 1994, ApJ, 428, 680

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

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

* Boys S. F., Bernardi F., 1970, Mol. Phys., 19, 553

* Brandenburg J. G., Bannwarth C., Hansen A., Grimme S., 2018, J. Chem. Phys., 148, 064104

* Brown W. A., Bolina A. S., 2007, MNRAS, 374, 1006

* Calmonte U. et al., 2016, MNRAS, 462, S253

* Ceccarelli C. et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 379

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

* Chand A., Biswal H. S., 2019, J. Indian Inst. Sci., 100, 77

* Codella C. et al., 2021, A&A, 654, A52

* Collings M. P., Anderson M. A., Chen R., Dever J. W., Viti S., Williams D. A., McCoustra M. R. S., 2004, MNRAS, 354, 1133

* Collings M. P., Frankland V. L., Lasne J., Marchione D., Rosu-Finsen A., McCoustra M. R. S., 2015, MNRAS, 449, 1826

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

* Ferrero S., Zamirri L., Ceccarelli C., Witzel A., Rimola A., Ugliengo P., 2020, ApJ, 904, 11

* Ferrero S. et al., 2022, MNRAS, 516, 2586

* Ferrero S., Pantaleone S., Ceccarelli C., Ugliengo P., Sodupe M., Rimola A., 2023, ApJ, 944, 142

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

* Fredon A., Radchenko A. K., Cuppen H. M., 2021, Acc. Chem. Res., 54, 745

* Freedman D., Diaconis P., 1981, Z. Wahrscheinlichkeitstheorie verw. Gebiete, 57, 453

* Frisch M. et al., 2016, Gaussian 16, Vol. 121. Gaussian, Inc., Wallingford, CT, p. 150

* Fuente A. et al., 2023, A&A, 670, A114

* Furuya K., Oba Y., Shimonishi T., 2022, ApJ, 926, 171

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

* Grimme S., Bannwarth C., Shushkov P., 2017, J. Chem. Theory Comput., 13, 1989

* Guo Y., Riplinger C., Becker U., Liakos D. G., Minenkov Y., Cavallo L., Neese F., 2018, J. Chem. Phys., 148, 011101

* Halkier A., Koch H., Jørgensen P., Christiansen O., Nielsen I. M. B., Helgaker T., 1997, Theor. Chem. Acc., 97, 150

* Hasegawa T. I., Herbst E., Leung C. M., 1992, ApJS, 82, 167

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

* Heyl J., Viti S., Vermariën G., 2023, Faraday Discussions, 245, 569

* Jenkins E. B., 2009, ApJ, 700, 1299

* Kendall R. A., Dunning Thom H. J., Harrison R. J., 1992, J. Chem. Phys., 96, 6796

* King D. A., 1975, Surf. Sci., 47, 384

* Kushwahaa T., Drozdovskaya M. N., Tychoniec Ł., Tabone B., 2023, A&A, 672, A122

* Laas J. C., Caselli P., 2019, A&A, 624, A108

* Lamberts T., Kästner J., 2017, J. Phys. Chem. A, 121, 9736

* Luna R., Luna-Ferrándiz R., Millán C., Domingo M., Caro G. M. M., Santonja C., Satorre M. A., 2017, The Astrophysical Journal, 842, 51

* Mayhall N. J., Raghavachari K., Hratchian H. P., 2010, J. Chem. Phys., 132, 114107

* McClure M. K. et al., 2023, Nat. Astron., 7, 431

* McQuarrie D., 2008, Quantum Chemistry. University Science Books, New York, available at:https://books.google.es/books?id=zzxLTIljQB4C

* Minissale M. et al., 2022, ACS Earth Space Chem., 6, 597

* Molpeceres G., Kästner J., Herrero V. J., Peláez R. J., Maté B., 2022, A&A, 664, A169

* Navarro-Almaida D. et al., 2020, A&A, 637, A39

* Neese F., 2022, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 12, e1606

* Nguyen T., Oba Y., Shimonishi T., Kouchi A., Watanabe N., 2020, ApJ, 898, L52

* Oba Y., Tomaru T., Lamberts T., Kouchi A., Watanabe N., 2018, Nat. Astron., 2, 228

* Palumbo M. E., Tielens A. G. G. M., Tokunaga A. T., 1995, ApJ, 449, 674

* Palumbo M. E., Geballe T. R., Tielens A. G. G. M., 1997, ApJ, 479, 839

* Pantaleone S., Enrique-Romero J., Ceccarelli C., Ugliengo P., Balucani N., Rimola A., 2020, ApJ, 897, 56

* Pantaleone S., Enrique-Romero J., Ceccarelli C., Ferrero S., Balucani N., Rimola A., Ugliengo P., 2021, ApJ, 917, 49

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

* Ferrero J., Enrique-Romero J., Ferrero S., Ceccarelli C., Podio L., Codella C., Rimola A., Ugliengo P., 2022, ApJ, 938, 158

* Ferrero J., Beitia-Antero L., Fuente A., Ugliengo P., Rimola A., 2024, MNRAS, 527, 10697

* Polanyi M., Wigner E., 1925, Z. Angew. Phys., 33, 429

* Potapov A., McCoustra M., 2021, Int. Rev. Phys. Chem., 40, 299

* Puletti F., 2014, PhD thesis, UCL (University College London)
Ruffle D. P., Hartquist T. W., Caselli P., Williams D. A., 1999, MNRAS, 306, 691

Sturm J. A. et al., 2023, A&A, 679, A138

Tait S. L., Dohnálek Z., Campbell C. T., Kay B. D., 2005, J. Chem. Phys., 122, 164708

Tieftrunk A., Pineau des Forets G., Schilke P., Walmsley C. M., 1994, A&A, 289, 579

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

Tinacci L., Germain A., Pantaleone S., Ceccarelli C., Balucani N., Ugliengo P., 2023, ApJ, 951, 32

Vastel C. et al., 2018, MNRAS, 478, 5514

Vidal T. H. G., Wakelam V., 2018, MNRAS, 474, 5575

Vidal T. H. G., Loison J.-C., Jaziri A. Y., Ruaud M., Gratier P., Wakelam V., 2017, MNRAS, 469, 435

Wakelam V., Caselli P., Ceccarelli C., Herbst E., Castets A., 2004, A&A, 422, 159

Wakelam V., Loison J. C., Mereau R., Ruaud M., 2017, Mol. Astrophys., 6, 22

Weigend F., Köhn A., Hättig C., 2002, J. Chem. Phys., 116, 3175

Woon D. E., Dunning Thom H. J., 1993, J. Chem. Phys., 98, 1358

Zamirri L., Corno M., Rimola A., Ugliengo P., 2017, ACS Earth Space Chem., 1, 384

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### h2s_supp.zip

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: PRUNING PROCEDURE

While performing geometry optimization, various initial H₂S positions may converge to a common minimum on the potential energy surface (PES) because of the complexity of the PES, the relatively mild interaction energy, and the large set of starting geometries. Tinacci et al. (2023) demonstrated the significance of trimming redundant structures and suggested using RMSD and the ΔBE calculated at the DFT level as general criteria for identifying similar structures. Considering the substantial computational demands of

<figure><img src="image_8.png" /><figcaption>**Figure A1.** Correlation plot between RMSD and $|\Delta E_{Cij}|$ for all the H₂O binding sites with two different data sets from Tinacci et al. (2023). Left panel: Structures and energies at xTB-GFN2 level. Right panel: Structures and energies at DFT B97-3c level. For both panels, the green rectangle identifies the redundant structures included in the area $-0.1:1.5 \text{ kJ mol}^{-1}$ and $-0.1:0.8 \text{ Å}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure A2.** Correlation plot between $|\Delta E_{Cij}|$ and RMSD for H₂S binding sites at xTB-GFN2 level. In the left panel, the total scattered plot is reported. The green area is then zoomed and reported in the right panel. The orange rectangle in the right panel highlights the structures that the pruning procedure has identified as redundant. Transparency was set for each point to show data clustering.</figcaption></figure>

calculating 486 BE sites at the DFT level, we opted to explore alternative methods to alleviate the computational cost while upholding the desired level of accuracy. To this purpose, we carried out a correlation analysis between the RMSD and $|\Delta E_{Cij}|$ for the H₂O molecule data that had been previously published by Tinacci et al. (2023) at two consecutive levels of theory within our procedure: xTB-GFN2 (the left panel of Fig. A1) and B97-3c (the right panel of Fig. A1). The first glimpse at Fig. A1 suggests that the two scatter plots look similar. At a deeper analysis, for each plot there is a set points clusterized at very low values of RMSD and $|\Delta E_{Cij}|$, which are identified as redundant structures (see the green rectangle of Fig. A1). For the xTB-GFN2 data, the identified redundant structures were 351, while for B97-3c 349. Accordingly, for the purpose to reduce the computational burden without losing accuracy, this small discrepancy can be neglected.

Moving to the H₂S molecule, the selection criteria for the pruning procedure were identified by correlating the RMSD and $|\Delta E_{Cij}|$ energy for each possible combination among Model zones with the same number of atoms. In the left panel of Fig. A2, the correlation plot between RMSD and $|\Delta E_{Cij}|$ is reported for all the possible combinations. The green rectangle is zoomed and plotted in the right panel of Fig. A2. As it can be seen, in the lower left area of Fig. A2 right panel, there is a clear clustering of points that identifies the redundant BE sites. In particular, the orange rectangle area include 338 redundant structures which, once pruned, generate a distribution of 148 BE sites.

## APPENDIX B: BENCHMARK OF THE METHOD

### B1 GFN2 versus B97-3c versus CCSD(T)

As shown in Table B1, the results at B97-3c level closely matches those at CCSD(T)/aug-cc-pVTZ. Conversely, the GFN2 method (QM2), when applied to the (H₂S ... H₂O) dimer, only identifies the SH ... O donor minimum. However, as indicated by the colour bar of Fig. 4, the minimum number of water molecules found in the QM zone is two. Accordingly, we extended the benchmark to include the interactions of H₂S with at least two water molecules ((H₂O)₂ ... H₂S, called as a whole the trimer case), as this mimics more closely the realistic cases on the large icy grain (see Table B1).
In this case, the GFN2 minimum correctly match the B97-3c one, where H₂S behaves both as H-bond donor and acceptor. Furthermore, the comparison of BE values computed as single point energy at DLPNO-CCSD(T) are almost insensitive to the adopted geometry, giving confidence in the adoption of the GFN2 structures/energies as the low-level method in the ONIOM calculations. In conclusion, although the H₂O ... H₂S dimer benchmark reveals some drawbacks in the GFN2 methodology, on the (H₂O)₂ ... H₂S trimer either DFT and semi-empirical levels show good accordance, as well as the energies calculated at DLPNO-CCSD(T).

## B2 Comparison with previous results by our group

Previous BE values obtained by Ferrero et al. (2020) and Ferrero et al. (2022) are shown in Table 2. Despite they share the same structural model and computational methodology the two works provide different BE values. This is due to the fact that binding sites were selected by chemical intuition by the authors, thus biasing the sampling procedure and number of samples used (5 and 8, respectively). Once understood the differences between the two works we can now try to rationalize the discrepancies with this work. Moreover, both the previous works largely overestimate the BE range reported in this work, as shown in Table 2. In order to unravel the origin of this substantial difference we carefully looked into the adopted computational procedure.

The first source of divergence can be found in the ZPE correction. In their procedure, in order to reduce the computational burden, the authors decide to not compute explicitly the frequencies of the amorphous cases, but only the ones for crystalline structures. As a way to account for this approximation, they decided to correlate the BE not corrected for the ZPE (BE\*) and the one corrected for the ZPE (BE) for all the adsorbates on the crystalline surface (see Ferrero et al. 2022; Fig. 4). After this, a correction factor of 0.862 was extracted from the linear fit and applied to the BE, i.e. BE = 0.862 BE\*. For the purpose of verifying the goodness of this approximation, we reproduced the same plot but with our computed values. As reported in Fig. B1, the correction factor obtained with our data is 0.707. This difference significantly affects the final distribution, and can

<table><caption>**Table B2.** Benchmark of the correction ONIOM method used by Ferrero et al. (2020). Values in parenthesis are given in kJ mol⁻¹. BE\* is the BE not corrected for the ZPE (see equation B1).</caption><thead><tr><th>Structure</th><th>BE (Kelvin)</th></tr></thead><tbody><tr><td>BE\* Real Zone [B97-3c, P-Ice]</td><td>5770 (48.08)</td></tr><tr><td>BE Real Zone [B97-3c, P-Ice]</td><td>4722 (39.35)</td></tr><tr><td>BE Model Zone [B97-3c, 2H₂O: H₂S]</td><td>2828 (23.57)</td></tr><tr><td>BE Model Zone [CCSD(T)/extrapolated, 2H₂O: H₂S]</td><td>1716 (14.30)</td></tr><tr><td>BE(ONIOM)</td><td>3610 (30.08)</td></tr></tbody></table>

be observed in Fig. B2 in which our BE distribution at DFT level is compared with the one obtained using Ferrero's rigid factor. The difference between the two mean value is 472.8 K (3.94 kJ mol⁻¹).

The second crucial factor is the way in which the high level refinement was applied. In their procedure, the DFT BEs were benchmarked using an embedded cluster model with the ONIOM approach (Mayhall et al. 2010). The Real Zone is the periodic proton ordered ice (P-ice) with H₂S adsorbed, while the Model Zone is created by extracting a portion from the optimized geometry of the periodic system. It consists of the adsorbed molecule along with the two nearest water molecules from the ice surface in proximity to the adsorbates (2H₂O: H₂S). For the low level method (QM2), they selected the B3LYP-D3/A-VTZ BE, and for the high level method, the CCSD(T)/Jun-cc-pVDZ extrapolated until Jun-cc-pVQZ or Jun-cc-pV5Z when feasible. Finally, the refined BE was calculated as

$$ \text{BE(ONIOM)} = \text{BE(DFT, P-Ice)} - \text{BE(DFT, 2H}_2\text{O: H}_2\text{S}) + \text{BE(CCSD(T), 2H}_2\text{O: H}_2\text{S}) \quad \text{(B1)} $$

Similarly to the ZPE correction, Ferrero et al. fitted the BE computed at DFT level with the ONIOM BE on a subset of adsorbates on P-ice. As claimed by the author (see their fig. 6 of Ferrero et al. 2020), the good correlation between the two BEs confirms the quality of the DFT results. For this reason the authors decided not to refine with the CCSD(T) method, and all the BEs were published at DFT level (Ferrero et al. 2020; Ferrero et al. 2022). As done for the previous case, we checked the validity of this approximation for the adsorption of H₂S. The periodic structure was reoptimized at B97-

<table><caption>**Table B1.** H₂S complexes with one and two H₂O molecules as benchmarks for the applied methods. BE (kJ mol⁻¹) are both BSSE and ZPE corrected from B97-3c frequencies, while BE\* values are free form ZPE correction. GFN2 identifies only the donor interaction for H₂O ... H₂S case (see Appendix B for the complete discussion). For the (H₂O)₂ ... H₂S system, the superscript 1 or 2 identifies each water molecule. B97-3c(GFN2) refers to geometry optimization starting from the minimum structures at GFN2. Distances (Å) and angles (degrees) refer to H-bonds only.</caption><thead><tr><th colspan="5">H₂O ... H₂S</th></tr></thead><tbody><tr><th scope="row">Level of theory</th><th scope="col">BE (kJ mol<sup>-1</sup>)</th><th scope="col">BE\* (kJ mol<sup>-1</sup>)</th><th scope="col">OH ... S (SH ... O)</th><th scope="col">OHS (SHO)</th></tr><tr><th scope="row">CCSD(T)/aug-cc-pVTZ</th><td>5.1 (4.6)</td><td>11.5 (10.5)</td><td>2.54 (2.22)</td><td>165.9 (176.9)</td></tr><tr><th scope="row">DLPNO-CCSD(T)/aug-cc-pVTZ</th><td>4.9 (4.5)</td><td>11.3 (10.4)</td><td>2.54 (2.22)</td><td>165.9 (176.8)</td></tr><tr><th scope="row">B97-D3/aug-cc-pVTZ</th><td>7.35 (4.3)</td><td>13.8 (10.2)</td><td>2.53 (2.23)</td><td>166.7 (174.3)</td></tr><tr><th scope="row">B97-3c</th><td>9.9 (8.5)</td><td>16.3 (14.2)</td><td>2.51 (2.17)</td><td>165.9 (179.2)</td></tr><tr><th scope="row">HF-3c</th><td>1.2 (7.3)</td><td>7.7 (13.3)</td><td>2.86 (2.08)</td><td>167.1 (177.2)</td></tr><tr><th scope="row">GFIN2</th><td>/(10.5)</td><td>/(16.9)</td><td>/(1.98)</td><td>/(175.2)</td></tr><tr><th scope="row">DLPNO-CCSD(T)/aug-cc-pVQZ</th><td>5.3 (4.9)</td><td>11.8 (10.9)</td><td>/</td><td>/</td></tr><tr><th scope="row">DLPNO-CCSD(T)/aug-cc-pVTZ//B97-3c</th><td>4.8 (3.6)</td><td>11.2 (9.5)</td><td>/</td><td>/</td></tr></tbody></table><table><thead><tr><th colspan="5">(H₂O)₂ ... H₂S</th></tr></thead><tbody><tr><th scope="row">Level of theory</th><th scope="col">BE (kJ mol<sup>-1</sup>)</th><th scope="col">BE\* (kJ mol<sup>-1</sup>)</th><th scope="col">OH<sup>[1]</sup> ... S (SH ... O<sup>[2]</sup>)</th><th scope="col">[1]OHS (SHO<sup>[2]</sup>)</th></tr><tr><th scope="row">GFIN2</th><td>18.5</td><td>26.6</td><td>2.50 (1.94)</td><td>143.4 (156.2)</td></tr><tr><th scope="row">B97-3c(GFN2)</th><td>24.2</td><td>35.3</td><td>2.45 (2.07)</td><td>151.3 (152.9)</td></tr><tr><th scope="row">DLPNO-CCSD(T)/aug-cc-pVTZ//B97-3c(GFN2)</th><td>12.0</td><td>23.1</td><td>/</td><td>/</td></tr><tr><th scope="row">DLPNO-CCSD(T)/aug-cc-pVTZ//GFIN2</th><td>12.3</td><td>20.4</td><td>/</td><td>/</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure B1.** Correlation plot between BE and BE\* of all H₂S binding sites reported in kJ mol⁻¹. The intercept is set to zero. BE is BSSE free.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure B2.** Histogram plot for two different BE distributions at B97-3c level of theory. The top panel reports the explicit evaluation of the ZPE. The bottom panel reports the BE distribution computed using the rigid factor reported in the work by Perrero et al. (2022) of 0.862, i.e. BE = 0.862 BE\*.</figcaption></figure>

3c and all the DFT calculations were carried out using the same functional (structures taken from Zenodo repository). The post-HF calculations were obtained with Gaussian16 program (Frisch et al. 2016) using the CCSD(T) method extrapolated until the Jun-cc-pV5Z basis set. Table B2 shows the BE computed following the ONIOM procedure from Ferrero et al. (2020) using the B97-3c functional and the CCSD(T) extrapolated until the Jun-cc-pV5Z method. The correction factor of the refinement procedure can be extracted by subtracting the BE[B97-3c, P-ice] to the BE(ONIOM) as computed in equation (B1), i.e. 39.35–30.08 kJ mol⁻¹. Approximately, the

correction factor is 9 kJ mol⁻¹ (1000 K), which in our opinion cannot be considered negligible. In conclusion, although averaging on all the adsorbates is a clever solution to reduce the computational burden, when dealing with more electronically subtle systems, such as third period-bearing molecules, the approximation made could severely deviate the BE from more accurate value.

## APPENDIX C: PRE-EXPONENTIAL FACTOR FORMALISM

The experimental desorption process can be described with the Polanyi–Wigner equation (Polanyi & Wigner 1925):

$$ -\frac{dN}{dt} = N^i \nu(T) \exp\left(-\frac{E_a}{RT}\right), \qquad (C1) $$

where N is the number of occupied sites, i is the order of the process, ν is the pre-exponential factor or pre-factor, E_a is the activation energy, R is the ideal gas constant, T is the temperature, and t is the time. Usually, the physisorption of small molecules is assumed to be a non-activated process, so that E_a is equal to the BE.

In this work, we adopted the ν<sub>Tait</sub> pre-factor (Tait et al. 2005; Minissale et al. 2022):

$$ \nu_{\text{Tait}} = \frac{k_B T}{h} \left( \frac{2\pi m_A k_B T}{h^2} \right) A \frac{\sqrt{\pi}}{\sigma h^3} (8\pi^2 k_B T)^{\frac{3}{2}} \sqrt{I_x I_y I_z}, \qquad (C2) $$

where I_x, I_y, I_z are the principal moment of inertia of H₂S and σ is the symmetry factor that identifies the indistinguishable rotational configurations of H₂S, here set to 2 (since H₂S belongs to the C₂v group). In this work all the values used to compute the pre-factor

<figure><img src="image_12.png" /><figcaption>**Figure C1.** Temperature dependence on the Tait et al. (2005) pre-factor.</figcaption></figure>

<table><caption>**Table C1.** Data used to compute the pre-exponential factor with the associated units.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>h</td><td>Planck constant</td><td>6.62606957e-34</td><td>m<sup>2</sup> kg s<sup>-1</sup></td></tr><tr><td>k<sub>B</sub></td><td>Boltzmann constant</td><td>1.3806488e-23</td><td>m<sup>2</sup> kg s<sup>-2</sup> K<sup>-1</sup></td></tr><tr><td>A</td><td>Surface per adsorbed molecule</td><td>1e-19</td><td>m<sup>2</sup></td></tr><tr><td>I<sub>x</sub></td><td>Moment of inertia along x</td><td>2.747e-47</td><td>m<sup>2</sup> kg</td></tr><tr><td>I<sub>y</sub></td><td>Moment of inertia along y</td><td>3.155e-47</td><td>m<sup>2</sup> kg</td></tr><tr><td>I<sub>z</sub></td><td>Moment of inertia along z</td><td>5.912e-47</td><td>m<sup>2</sup> kg</td></tr><tr><td>m<sub>A</sub></td><td>H<sub>2</sub>S mass</td><td>5.6591e-26</td><td>kg</td></tr><tr><td>R</td><td>Molar gas constant</td><td>8.3144621</td><td>J K<sup>-1</sup> mol<sup>-1</sup></td></tr></tbody></table>
with the associated units are reported in Table C1. In Fig. C1, the computed temperature dependence on the Tait et al. (2005) pre-factor is reported.

APPENDIX D: THERMOCHEMISTRY DATA

In Table D1, the reaction enthalpies at 0 K for the H⁺ + HS⁻ → H₂S reaction are reported.

<table><caption>**Table D1.** Absolute enthalpies at 0 K are reported for the species involved in the reaction H⁺ + S⁻ → HS⁻ and H⁺ + HS⁻ → H₂S. All the values are in kJ mol⁻¹. In the last line the reaction enthalpy computed as ΣH<sub>prod</sub>(0K) – ΣH<sub>react</sub>(0K) is shown. Geometries and frequencies are calculated at B97-3c level. DLPNO-CCSD(T) method is coupled with the aug-cc-pVTZ basis set. The superscript ‘’ identifies the open-shell species.</caption><thead><tr><th>Species</th><th>B97-3C</th><th>DLPNO-CCSD(T)//B97-3C</th></tr></thead><tbody><tr><td>H<sup>+</sup> (doublet)</td><td>-1330.945004</td><td>-1312.280488</td></tr><tr><td>S<sup>+</sup> (triplet)</td><td>-1045062.305</td><td>-1044046.267</td></tr><tr><td>HS<sup>+</sup> (doublet)</td><td>-1046718.877</td><td>-1045698.361</td></tr><tr><td>H<sub>2</sub>S</td><td>-1048401.701</td><td>-1047380.137</td></tr><tr><td>Δ<sub>γ</sub>H<sup>°</sup> (H<sub>2</sub>S)</td><td>-351.9</td><td>-369.5</td></tr><tr><td>Δ<sub>γ</sub>H<sup>°</sup> (HS<sup>+</sup>)</td><td>-325.6</td><td>-339.8</td></tr></tbody></table>

This paper has been typeset from a TeX/Latex file prepared by the author.