# Binding Energies of Interstellar Relevant S-bearing Species on Water Ice Mantles: A Quantum Mechanical Investigation

Jessica Perrero, Joan Enrique-Romero, Stefano Ferrero, Cecilia Ceccarelli, Linda Podio, Claudio Codella, Albert Rimola, Piero Ugliengo

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

## ► To cite this version:

Jessica Perrero, Joan Enrique-Romero, Stefano Ferrero, Cecilia Ceccarelli, Linda Podio, et al., Binding Energies of Interstellar Relevant S-bearing Species on Water Ice Mantles: A Quantum Mechanical Investigation. The Astrophysical Journal, 2022, 938, 10.3847/1538-4357/ac9278. insu-03860279

HAL Id: insu-03860279

https://insu.hal.science/insu-03860279

Submitted on 19 Nov 2022

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_2.png" /></figure>
OPEN ACCESS

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

CrossMark

# Binding Energies of Interstellar Relevant S-bearing Species on Water Ice Mantles: A Quantum Mechanical Investigation

Jessica Ferrero<fnref n="1" />,<fnref n="2" /> ID, Joan Enrique-Romero<fnref n="1" /> ID, Stefano Ferrero<fnref n="1" /> ID, Cecilia Ceccarelli<fnref n="3" /> ID, Linda Podio<fnref n="4" /> ID, Claudio Codella<fnref n="3" />,<fnref n="4" /> ID, Albert Rimola<fnref n="1" /> ID, and Piero Ugliengo<fnref n="2" /> ID

Received 2022 July 6; revised 2022 September 12; accepted 2022 September 14; published 2022 October 24

## Abstract

Binding energies (BEs) are one of the most important parameters for astrochemical modeling determining, because they govern whether a species stays in the gas phase or is frozen on the grain surfaces. It is currently known that, in the denser and colder regions of the interstellar medium, sulfur is severely depleted in the gas phase. It has been suggested that it may be locked into the grain icy mantles. However, which are the main sulfur carriers is still a matter of debate. This work aims to establish accurate BEs of 17 sulfur-containing species on two validated water ice structural models, the proton-ordered crystalline (010) surface and an amorphous water ice surface. We adopted density functional theory-based methods (the hybrid B3LYP-D3(BJ) and the hybrid meta-GGA M06-2X functionals) to predict structures and energetics of the adsorption complexes. London's dispersion interactions are shown to be crucial for an accurate estimate of the BEs due to the presence of the high polarizable sulfur element. On the crystalline model, the adsorption is restricted to a very limited number of binding sites with single valued BEs, while on the amorphous model, several adsorption structures are predicted, giving a BE distribution for each species. With the exception of a few cases, both experimental and other computational data are in agreement with our calculated BE values. A final discussion on how useful the computed BEs are with respect to the snow lines of the same species in protoplanetary disks is provided.

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

## 1. Introduction

Stars like our Sun begin their journey in the interstellar medium (ISM), where star formation takes place in dense (10³–10⁴ cm⁻<fnref n="3" />) and cold (10 K) regions, the so-called molecular clouds. In the past, it was thought that the ISM was too harsh an environment for molecules to survive and thrive. However, almost 80 yr ago, the first diatomic molecules (the methylidyne radical CH, its cation CH⁺, and the cyano radical CN) were detected by means of optical and ultraviolet spectroscopy (Swings & Rosenfeld 1937; McKellar 1940; Douglas & Herzberg 1942), which indeed indicated that chemistry is present and has an important role in the ISM.

At present, we count more than 270 detected gas-phase species (McGuire 2021), a number that constantly increases given the enhanced performances of the observational facilities (e.g., Yebes 40 m, IRAM 30 m, NOEMA, the Atacama Large Millimeter/submillimeter Array (ALMA), and Green Bank Telescope (GBT)), and the expectation of the new results provided by the James Webb Space Telescope. As a consequence of these thrilling discoveries, curiosity concerning the chemistry of the ISM has been growing more and more keen.

The characteristic densities and temperatures of molecular clouds allow (i) the molecules to survive, since the external UV

fields are highly diminished, and (ii) the formation of thick icy mantles coating the surfaces of dust grains, due to the adsorption and in situ synthesis of species and their subsequent hydrogenation. The resulting icy mantles are dominated by water, but they also contain other volatile species like CO, CO₂, CH₃OH, and NH₃, and hence are referred to as “dirty ices” (Boogert et al. 2015).

The phenomena occurring on the icy mantle surfaces are of paramount importance from a chemical point of view. Indeed, ice mantles are thought to facilitate the occurrence of reactions taking place at their surfaces, primarily exhibiting three functions: (i) as pre-concentrators of chemical species, especially relevant in the low-density environments of molecular clouds, where the rate of collision between gas-phase species is extremely small, (ii) as chemical catalysts, by decreasing the activation energies of the reactions, and hence overcoming the energy barriers under the interstellar conditions (Zamirri et al. 2019b), and (iii) as third bodies, by absorbing the large energy excess released by exothermic reactions, without undermining the stability of the newly formed products (Pantaleone et al. 2020, 2021). The first and the third points are the keys to successful recombination reactions.

For a grain surface reaction to take place, at least one of the two reactants needs to be adsorbed on the grain. The adsorption is regulated by two variables: the dust temperature, and the binding energy (BE) of the adsorbed species. The latter quantity defines how strong the interaction between the species and the surface is, and dictates the chance that a given adsorbed

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

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

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

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

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

<footnote n="4">⁴INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>
species can be ejected into the gas phase: the higher the BE, the higher the temperature required to thermally desorb. In the current astrochemical models, BE values are important input parameters, not only to consider likely desorption events (e.g., Penteado et al. 2017) but to simulate diffusion processes on the surfaces, as the diffusion barriers are usually assumed to be a fraction of the BE (Mispelaer et al. 2013; Karssemeijer & Cuppen 2014; Ghesquiere et al. 2015; Lauck et al. 2015; Cuppen et al. 2017; He et al. 2017, 2018; Cooke et al. 2018; Kouchi et al. 2020; Maté et al. 2020). Therefore, BEs are crucial parameters in determining the ISM chemistry and its resulting composition (Penteado et al. 2017; Wakelam et al. 2017; Ferrero et al. 2020; Enrique-Romero et al. 2021).

## 1.1. Previous Studies of the BEs

Determining BEs is a challenging task (Minissale et al. 2022). From an experimental point of view, effective and coverage-dependent BEs are obtained via temperature pro- grammed desorption (TPD) experiments. TPD consists of two steps: first, the substrate, maintained at a constant temperature, is exposed to the adsorption of the species; then, the temperature is increased until the desorption of the species, which is collected and analyzed by a mass spectrometer. The BEs are then extracted by applying the direct inversion method through the Polanyi–Wigner equation (e.g., Dohnálek et al. 2001). Despite its great usefulness, this technique shows two limitations: (i) the obtained BEs depend not only on the morphology and the composition of the substrate but also on the regimes in which the experiments are performed, i.e., submonolayer, monolayer, or multilayer (e.g., Noble et al. 2012; He et al. 2016; Chaabouni et al. 2018) and will return coverage-dependent BEs, and (ii) it measures a desorption enthalpy, a quantity that is equal to the BE only in the absence of other activated processes, like ice restructuring, usually assumed to be negligible (He et al. 2016).

Despite the numerous literature studies that have investi- gated the desorption processes by means of TPD (e.g., Collings et al. 2004; Noble et al. 2012; Dulieu et al. 2013; Fayolle et al. 2016; He et al. 2016; Smith et al. 2016), only a limited variety of molecules has been considered, in particular, stable closed- shell species, in contrast to the large variety of interstellar species. Moreover, the employed substrates are almost invariably water ices.

A possible alternative to TPD is to adopt state-of-the-art computational chemistry to simulate the adsorption process. The simulation can handle both closed- and open-shell species, allotting some extra care for the latter. The main limitation is the trade-off between the system size (the icy grain model) and the accuracy of the calculation. A number of computational studies have reported the BEs of important astrochemical species like H, H₂, N, O, HF, CO, and CO₂ on icy surfaces modeled by periodic/cluster crystalline/amorphous systems (e.g., Al-Halabi & Van Dishoeck 2007; Karssemeijer & Cuppen 2014; Zamirri et al. 2017; Shimonishi et al. 2018; Zamirri et al. 2019a; Bovolenta et al. 2020; Ferrero et al. 2020). Other works addressed the BEs of a large number of species, adopting very approximated models for the substrate. Two examples are the works of Wakelam et al. (2017) and Das et al. (2018), in which the ice surface simulated by water clusters of minimal nuclearity allowed the prediction of the BEs for 100 species. In Wakelam et al. (2017), a single water molecule was adopted to represent the ice surface, while recovering the

missing components of the BEs through a clever fitting procedure against selected experimental BE values. In Das et al. (2018), the adoption of clusters from one to six water molecules allowed them to study the BE of 16 species, with the puzzling result that water tetramer was found to be the best model for ice.

One way to generate computer models for amorphous ices is by using water clusters cut out from a crystalline ice, heating them up by molecular dynamics runs and by quenching at low temperatures (Shimonishi et al. 2018; Rimola et al. 2021b). Another approach is to refer to a model repeated in space by the periodic boundary conditions. This allows for the simulation of both the crystalline ice (here adopting an already studied proton-ordered model by Casassa et al. 1997) and an amorphous one. In the latter case, the usual procedure is starting from a large unit cell of a crystalline ice and performing a number of heating/freezing cycles to arrive to an amorphous ice model. Another strategy, recently proposed by the ACO-FROST scheme (Germain et al. 2022), consists of growing icy clusters by a random step-by-step addition of molecular water molecules.

More recently, Ferrero et al. (2020) computed a large set of BEs for 21 astrochemically relevant species, on both crystalline and amorphous surfaces of water ice. Both models ensured that hydrogen-bond cooperativity (of fundamental relevance to obtain accurate BE values) is fully taken into account, at variance with the small ice clusters adopted by Wakelam et al. (2017) and Das et al. (2018). For the crystalline icy model, only very few adsorption sites were available, limiting the BEs to one or two distinct values per species. For the amorphous ice model, however, the complex surface morphology allowed the calculation of multiple BE values for each species, resulting in a BE range better describing the variety of binding sites expected on real amorphous ice mantles.

In the present work, we apply the methodology of Ferrero et al. (2020) in order to enlarge the set of BEs. The focus here is on the adsorption of 17 interstellar relevant S-bearing species, namely: CS, C₂S`, C₃S, C₄S`, CH₃SH, H₂CS, HS`, H₂S, HS₂, H₂S₂, NS`, OCS, S`, S`, SO`, and SO₂ (dots indicate the unpaired electrons). The chosen species satisfy two criteria: (i) they are neutral species, since positive ions become neutral when landing on negatively charged dust grains (Walmsley et al. 2004; Ceccarelli & Dominik 2005; Rimola et al. 2021a) as long recognized (Draine & Sutin 1987; Rae et al. 2004; Mason et al. 2014), and (ii) they envisage at most six atoms to ensure relatively compact structures needed to properly probe the largest number of adsorption sites on ice models with relatively small unit cells.

## 1.2. The Sulfur Depletion

The choice of studying the BEs of S-bearing species is due to a long standing issue in the field of astrochemistry: the sulfur depletion problem. In dense clouds, sulfur is severely depleted from the gas phase (e.g., Tieftrunk et al. 1994; Ruffle et al. 1999; Wakelam et al. 2004; Phuong et al. 2018; Vastel et al. 2018; van't Hoff et al. 2020) by more than 2 orders of magnitude with respect to the solar system S abundance ([S]/[H] = 1.8 × 10⁻⁵; Anders & Grevesse 1989; Woods et al. 2015). Therefore, sulfur may freeze out on dust grains after its hydrogenation, or react with other species in a similar way as atomic oxygen. H₂S was expected to be the main sink of sulfur on dust grains (Garrod et al. 2007; Jiménez-Escobar & Muñoz Caro 2011), but it has
never been directly detected on interstellar ices yet. The only species detected in the ice mantles are OCS and, tentatively, SO₂ (Boogert et al. 2015). Accordingly, the most important S-carriers are still unknown. Nowadays, the most likely species that act as reservoirs of S are thought to be organosulfur compounds (Laas & Caselli 2019) and polysulphanes (H₂S*n*, Druard & Wakelam 2012). These are refractory species that, once trapped on the mantles or in the core of the grains, cannot desorb and, therefore, become undetectable (Woods et al. 2015). In the last three years, seven new gas-phase S-bearing species were identified, most of them in TMC-1 (Cernicharo et al. 2021; Rodríguez-Almeida et al. 2021), and characterized by a C–S chemical bond.

A recent review of the laboratory experiments performed on the chemistry of sulfur in the condensed phase has highlighted that the formation of SO₂, SO₃, hydrates of H₂SO₄, and related species is very common when performing photolysis, proton-irradiation, and radiolysis of mixed ices. However, neutral-neutral reactions have yet to be deeply explored. Moreover, gas-phase chemistry is likely to be a major contributor in sulfur astrochemistry (Mifsud et al. 2021).

The formation of organosulphur molecules was observed experimentally for the first time in the work of Ruf et al. (2019), where a 2:1:1 mixture of H₂O:CH₃OH:NH₃ ice was bombarded with S⁷⁺ ions, although showing a possible source of compounds relative to more evolved environments, like icy moons, Kuiper Belt objects, comets, and their building blocks.

In diffuse clouds, sulfur is mostly present in its ionized atomic form S⁺, and close to the cosmic abundance (Jenkins 2009), it progressively depletes as the cloud evolves (Hily-Blant et al. 2022). Sulfur depletion in protoplanetary disks and later stages can be explained by the presence of refractory sulfide minerals, such as FeS, in the grains (Keller et al. 2002; Kama et al. 2019).

With the present work, thus, we aim to contribute to this puzzling subject by providing a set of accurate BE values of several S-bearing species, which can be used in numerical modeling studies aimed at rationalizing observations and laboratory experiments.

## 2. Methodology

Since the calculations follow the same procedure as Ferrero et al. (2020), we refer to this work for the equations used here to compute the BEs of the S-bearing species set. Moreover, the reader can find details on these expressions in Appendix A. A thorough guide to the computation of BEs is available in the Supplementary Material.

### 2.1. Computational Details

We adopted a periodic approach to model the water ice surfaces and the adsorption of the 17 S-bearing species on them. The calculations were performed with the periodic ab initio code CRYSTAL17 developed by Dovesi et al. (2018), which adopts localized Gaussian functions as basis sets. This code exploits the local combination of atomic orbitals (LCAO) approximation, expanding the Bloch functions of periodic systems as a linear combination of atom-centered Gaussian functions. This software implements both the Hartree–Fock (HF) and Kohn–Sham self-consistent field methods to solve the electronic Schrödinger equation, taking advantage of the symmetry of the system, when present. Thanks to its capability to simulate systems with periodicity ranging from zero

(molecules) to three (solid bulks) dimensions, it allows us to rigorously define true slabs to modeling surfaces without the need of artificial replica along the direction perpendicular to the slab plane, as usually implemented in codes adopting plane waves basis sets.

On the crystalline model, the hybrid B3LYP (Becke 1988, 1993; Lee et al. 1988) and the hybrid meta-GGA M06-2X (Zhao & Truhlar 2008) density functionals were used in the geometry optimizations and in computing the BEs of closed-shell and open-shell species, respectively. To account for dispersion interactions, the hybrid B3LYP functional was combined with the Grimme's D3 empirical correction with the Becke–Johnson (BJ) damping scheme (i.e., B3LYP-D3(BJ) (Grimme et al. 2010, 2011; Sure & Grimme 2013). The hybrid meta-GGA M06-2X functional was used together with a spin-unrestricted formalism (Pople et al. 1995) to model the adsorption of open-shell species (with one or two unpaired electrons) because it better describes their adsorption properties (Ferrero et al. 2020) due to its higher percentage of exact exchange (54% compared to the 20% of B3LYP). Note that for M06-2X, the D3(BJ) correction for dispersion was not applied because, for this functional, the description of the dispersion component is already included in its definition. For this reason, the contribution of the dispersion component to the BEs can be worked out only for the closed-shell species (computed at B3LYP-D3(BJ)). The two density functional theory (DFT) methods have been used in combination with the Ahlrichs triple zeta valence quality basis set supplemented with a double set of polarization functions (i.e., A-VTZ*; Schäfer et al. 1992). Thus, for the crystalline models, the BEs were computed from the DFT energies that resulted from the DFT optimized structures. This scheme is hereafter referred to as DFT//DFT. All of the BEs were corrected for the basis set superposition error (BSSE) through the counterpoise correction method by Boys & Bernardi (1970).

Finally, in order to check for the accuracy of the DFT//DFT BEs and without the aim of performing a benchmarking study, a refinement at the CCSD(T) level was performed by applying the ONIOM2 scheme as proposed by Dapprich et al. (1999) to a number of species representative of the set. Given the dependency of CCSD(T) method on the basis set, the Dunning (1989) correlation consistent family of cc-pVnZ basis sets (where *n* is 2, 3, 4, 5, ... and corresponds to a ζ of D, T, Q, 5, ...) was employed. This allows us to perform an extrapolation of the desired property, here the BE, by plotting it against 1/*n*³. We checked that the numbers thus obtained approached the values given by DFT//DFT. The Jun-cc-pVnZ basis set (with *n* = 2, 3, and 4) was chosen to perform the calculations, with the exception of SO₂, for which the aug-cc-pVnZ basis set (where one more diffuse function is added to all of the atoms, including H) was adopted to recover the whole interaction between SO₂ and the two water molecules of the *model* system. See Appendix A for more details on the BSSE correction and on the ONIOM2 scheme.

On the amorphous ice model, the DFT//DFT scheme can hardly be carried out due to the larger unit cell size compared with the crystalline one. Therefore, we adopted the computationally cheaper semiempirical HF-3c method for the adsorption on the amorphous ice surface model. HF-3c is an HF-based method adopting a minimal basis set (MINI-1; Tatewaki & Huzinaga 1980), to which three empirical corrections (3c) are added (Sure & Grimme 2013): (i) the dispersion energy
(D3(BJ)) for noncovalent interactions (Grimme et al. 2010); (ii) a short range bond correction to recover the systematically overestimated covalent bond lengths for electronegative elements, due to the MINI-1 basis set (Grimme et al. 2011; Brandenburg et al. 2013); and (iii) the geometrical counterpoise (gCP) method developed by Grimme (Kruse & Grimme 2012) to a priori remove the BSSE. On the optimized HF-3c geometries, DFT single-point energy calculations were performed to compute the BEs (this procedure is hereafter referred to as DFT//HF-3c), which were also corrected for BSSE. Before adopting DFT//HF-3c for the amorphous system, we proved its accuracy compared to DFT//DFT for the crystalline ice, resulting in an excellent match.

For both the crystalline and amorphous systems, geometry optimizations were carried out by relaxing both the internal atomic positions inside the unit cell and the cell parameters.

Each stationary point was characterized by means of harmonic frequency calculations at the Γ point by diagonalizing the mass-weighted Hessian matrix of the second-order energy derivatives with respect to atomic displacements (Pascale et al. 2004; Zicovich-Wilson et al. 2004). Each Hessian matrix element was computed numerically by means of a six-point formula based on two displacements of ±0.003 Å from the minimum along each Cartesian coordinate. For the crystalline ice, computed at the DFT level, only a portion of the system constituted by the adsorbed species plus the two closest water molecules to the adsorbate was considered to build the Hessian matrix. In the case of the amorphous ice, the entire system was considered, and the full Hessian matrix was computed at the HF-3c level. For this later case, however, the frequencies were only used to check that the complexes are actual minima of the potential energy surfaces (all of the frequencies are real), since they are expected to not be accurate enough because of the approximated HF-3c methodology. In contrast, from the DFT frequencies computed for the fragments, we computed the zero-point vibrational energy (ZPE) corrections to: (i) correct the electronic BEs for the crystalline ice by including the corresponding ZPE contributions, and (ii) derive a scaling factor to apply to the BEs computed on the amorphous ice model to arrive at ZPE-corrected values without the need to explicitly run a full DFT frequency calculation (see Appendix B for more details on these two points).

## 2.2. Ice Models

Interstellar ices are thought to be mostly formed by amorphous solid water (Watanabe & Kouchi 2008; Boogert et al. 2015). However, as a starting point of our study, we have chosen a crystalline model (see the top panel of Figure 1). This choice is due to two main reasons: (i) crystalline structures are well defined due to symmetry constraints and are computationally cheap, and (ii) no definite structure for amorphous ice models is available. Moreover, regions rich in crystalline ices have been observed in protoplanetary disks and stellar outflows (Molinari et al. 1999; Terada & Tokunaga 2012). Our ice surface model is derived from the bulk of the proton-ordered P-ice, which was cut along the (010) surface defining a 2D periodic slab model (Casassa et al. 1997; Zamirri et al. 2018). The thickness of the ice (10.9 Å) was enough to converge the corresponding surface energy. We adopted a 2 × 1 supercell slab model, consisting of 12 atomic layers and cell parameters of $|a|=8.980$ Å and $|b|=7.082$ Å (at B3LYP-D3(BJ)/A-VTZ\* level), which are large enough to

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Top panel: side view of the crystalline (010) water P-ice slab model (along the *b* lattice vector). Bottom panel: side view of the amorphous water ice slab model (along the *b* lattice vector). Numbers identify the different adsorption sites: those in blue represent dangling hydrogen atoms (dH), while those in red represent dangling oxygen atoms (dO).</figcaption></figure>

avoid lateral interactions between the adsorbed molecules in different replicas of the ice cell, thereby simulating the adsorption of isolated species. The supercell shows two dangling hydrogen (dH) atoms and two dangling oxygen (dO) atoms as binding sites (top panel of Figure 1). The structure of the ice is such that it ensures a null electric dipole along the nonperiodic *z*-axis. This is a direct consequence of the symmetry of the system, which shows two identical faces both at the top and at the bottom of the model. Therefore, the adsorption was modeled only on the top face of the system.

The more realistic amorphous model, already used in Ferrero et al. (2020; see bottom panel of Figure 1) consists of 60 water molecules per unit cell. The cell parameters are $|a|=20.355$ Å, $|b|=10.028$ Å, and $|\gamma|=102.983$ (at HF-3c). The different structural features between the upper and lower surfaces are responsible for the presence of a small electric dipole moment across the nonperiodic direction, and also for the presence of a variety of different binding sites. Therefore, to model the adsorption of the S-bearing species, we selected eight characteristic adsorption sites (highlighted in the bottom panel of Figure 1). The sites are chosen as follows: cases 1, 2, and 5 are located on the top surface of the ice, case 4 and 6 are inside the cavity of the model, and case 3, 7, and 8 are on the bottom-surface of the model. Binding sites from 1–6 are dH, while 7 and 8 are dO.

## 3. Results

In this work, the adsorption of 17 S-bearing species (shown in Figure 2) has been simulated, of which eight are closed-shell species and nine open-shell species (with unpaired electrons as indicated by the dots in Figure 2). For the open-shell species, the electronic ground state of each one was checked and found to be in agreement with those reported by Woon & Herbst (2009) and Das et al. (2018). Therefore, NS\*, HS\*, HS\*₂, and
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Set of 17 S-bearing species considered in this work. Dots represent unpaired electrons.</figcaption></figure>

HCS\* are doublets (with one unpaired electron) and S\*\*, SO\*\*, S\*\*\*, C\_2S\*\*, and C\_4S\*\*\* are triplets (with two unpaired electrons). Each species was placed manually on the ice surface models according to the molecule/surface electrostatic potential complementarity, and from these manually built systems, geometry optimizations were performed.

## 3.1. BEs on Crystalline Ice

### 3.1.1. DFT/DFT BEs

Table 1 reports the computed BE values as described in the computational section.

The role of the dispersion interactions to the BEs was assessed by considering the D3(BJ) contribution only to the BE. In other words, for each full B3LYP-D3(BJ) optimized structure, we split the B3LYP-D3(BJ) total energy in a pure electronic DFT term coming from the B3LYP and the pure dispersive D3(BJ) one. This approach maximizes the weight of the dispersion because the attractive nature of the D3(BJ) term brings the adsorbates as close as possible to the icy surface, which, in turn, boosts the D3(BJ) value. A different approach would be to optimize the geometries at B3LYP only and then evaluate the D3(BJ) term through a single-point energy evaluation (referred to as the B3LYP-D3(BJ)//B3LYP scheme). The missing D3(BJ) term during the geometry optimization causes both an elongation of the average distance and important changes in the spatial arrangement of the adsorbate compared to the full B3LYP-D3(BJ) cases. Therefore, the dispersion contribution for the B3LYP-D3(BJ)//B3LYP is expected to be smaller (underestimated) than the one from the full B3LYP-D3(BJ) optimization calculation. We checked both approaches (data not reported) for two sets of molecules, O-bearing (CO, H₂O, CO₂, CH₃OH, and H₂CO) and S-bearing (CS, H₂S, OCS, CH₃SH, and H₂CS) adsorbed on the P-ice crystalline surface. We found the B3LYP structures involving CS and, especially OCS, to be rather different from the corresponding B3LYP-D3(BJ) ones. The same was found for the CO and CH₃OH cases. Nonetheless, the percentage average difference in the dispersion between the two approaches is a modest 11% in favor of the B3LYP-D3(BJ), as expected. In this work, we were interested in establishing the relative weight of dispersion (rwd) between S- and O-containing species, rather than its absolute value. For the

cases described above, we computed an rwd percentage average difference between B3LYP-D3(BJ)//B3LYP and a fully optimized B3LYP-D3(BJ) of only 4%, without any order reverting in the dispersion weight. Therefore, in agreement with this test case, the data shown in this work are relative to the full B3LYP-D3(BJ) approach, aware that the role of dispersion may be overestimated.

At the DFT//DFT level, computed BE disp values span a range from 2400–6900 K, depending on the kind of interaction that the adsorbate experiences with the ice surface, i.e., hydrogen bonds, interactions between both permanent and instantaneous dipole moments, and dispersive interactions. According to that, we can distribute the molecules in four different groups:

*Molecules of Group I:* (CH₃SH, H₂CS, H₂S, H₂S₂, HCS\*, HS\*, and HS\*₂): these molecules contain groups that can act as H-bond donors and H-bond acceptors. It is worth noting that in the cases of H₂CS and HCS\*, the H-bond donor groups are weak.

*Molecules of Group II:* (OCS, SO₂, NS\*, and SO\*\*): these molecules feature only H-donor acceptor groups because they contain atoms that are more electronegative than S.

*Molecules of Group III:* (CS, C₂S\*\*, C₃S, and C₄S\*\*): molecules that can establish one H-bond interaction with their C-end atom and a surface dH.

*Molecules of Group IV:* (S\*\*\* and S\*\*₂): species consisting of only S atoms so that dispersive interactions are expected to be the main binding driving forces.

The molecule with the smallest dispersion contribution to its BE is H₂S (of Group I), which, like water, is capable of acting as both an H-bond donor and acceptor, hence involving two water molecules of the ice surface in its adsorption. H₂S, H₂S₂, H₂CS, and CH₃SH (of the same Group I) show this same adsorption feature. However, the presence of an S–S or a C–S bond enhances the contribution of dispersion forces, which can reach up to circa 70% of the BE, as in H₂S₂. The open-shell species HS\* and HS\*₂ (also of Group I) have smaller BEs than their hydrogenated counterparts due to the absence of one H atom.

For molecules of Group II, the sulfur atom interacts through dispersive forces with dO of the nearest surface water molecules. SO₂, thanks to its two oxygen atoms that act as H-bond acceptors, presents the highest BE of the entire set of tested molecules (6880 K). OCS, which shows an almost null dipole moment, presents the largest dispersion contribution to the BE of the entire set (>90%), hence being one of the species less bound to the ice.

Molecules of Group III indeed form an H-bond between their C-ends and a surface dH. This interaction, moreover, is enhanced due to the dipole moment of the C chains and to the carbene-like character of the C-end atom. Furthermore, these species can bind strongly to the surface thanks to the onset of dispersive forces that become stronger with the size of the molecule. Despite these premises, the BEs of these species are similar (4000–4500 K), with the exception of C₃S, which has a higher BE (6723 K at the most) because of its adsorption geometry. That is, once C₃S adsorbs on the surface, it is arranged in rows with a distance of only 3.1 Å between the S-end and the C-end of two neighboring molecules, therefore maximizing favorable head-to-tail interactions between dipoles. This results in an energetic gain due to the lateral interactions
<table><caption>**Table 1**
Computed Binding Energies (BE disp, in kelvin, and without ZPE Corrections) of the S-bearing Species Set at the P-ice (010) (2 × 1) Supercell</caption><caption>BEs on the Crystalline Water Ice</caption><thead><tr><th rowspan="2">Species</th><th colspan="3">DFT//DFT</th><th colspan="3">DFT//HF-3c</th></tr><tr><th>BE Disp</th><th>BE No Disp</th><th>Disp (%)</th><th>BE Disp</th><th>BE No Disp</th><th>Disp (%)</th></tr></thead><tbody><tr><td>CS</td><td>4582</td><td>1600</td><td>2982(65%)</td><td>3861</td><td>1239</td><td>2622(68%)</td></tr><tr><td>C<sub>3</sub>S</td><td>6038/6723</td><td>1564/2129</td><td>4474(74%)/4594(68%)</td><td>5713/6146</td><td>421/1443</td><td>5292(93%)/4703(77%)</td></tr><tr><td>CH<sub>3</sub>SH</td><td>6182</td><td>2959</td><td>3223(52%)</td><td>5340</td><td>2081</td><td>3259(61%)</td></tr><tr><td>H<sub>2</sub>CS</td><td>4847/6242</td><td>2225/2850</td><td>2622(54%)/3392(54%)</td><td>4258/5641</td><td>1792/2417</td><td>2466(58%)/3223(57%)</td></tr><tr><td>H<sub>2</sub>S</td><td>5268</td><td>3332</td><td>1936(37%)</td><td>4679</td><td>2862</td><td>1816(39%)</td></tr><tr><td>H<sub>2</sub>S<sub>2</sub></td><td>5436/6158</td><td>1624/1936</td><td>3813(70%)/4222(69%)</td><td>5027</td><td>1527</td><td>3500(70%)</td></tr><tr><td>OCS</td><td>3007/3440</td><td>-229/325</td><td>3235(107%)/3115(90%)</td><td>2790</td><td>132</td><td>2658 (95%)</td></tr><tr><td>SO<sub>2</sub></td><td>6880</td><td>3608</td><td>3271(48%)</td><td>5929</td><td>2502</td><td>3428(58%)</td></tr><tr><td>HCS*</td><td>4113</td><td></td><td></td><td>3849</td><td></td><td></td></tr><tr><td>HS*</td><td>4270</td><td></td><td></td><td>3247</td><td></td><td></td></tr><tr><td>HS<sub>2</sub></td><td>4270/4763/4871</td><td></td><td></td><td>3572/3139</td><td></td><td></td></tr><tr><td>NS*</td><td>4330</td><td></td><td></td><td>3765</td><td></td><td></td></tr><tr><td>C<sub>2</sub>S**</td><td>4089/4691</td><td></td><td></td><td>3993/4498</td><td></td><td></td></tr><tr><td>C<sub>4</sub>S**</td><td>4366/4414</td><td></td><td></td><td>4246</td><td></td><td></td></tr><tr><td>S**</td><td>3921</td><td></td><td></td><td>3247</td><td></td><td></td></tr><tr><td>S<sub>2</sub>**</td><td>2393</td><td></td><td></td><td>2009</td><td></td><td></td></tr><tr><td>SO**</td><td>3861</td><td></td><td></td><td>1924</td><td></td><td></td></tr></tbody></table>

Note. The first column refers to the adsorbed species. Columns (2)–(4) contain the values computed at the DFT//DFT level of theory. The last three columns report the same values computed at the DFT//HF-3c level of theory. For the closed-shell species (computed at the B3LYP-D3(BJ) level), the BE contributions arising from the nondispersive (no disp) and the dispersive (disp) forces, together with their percentage (%), are listed. For the open-shell species (computed at the M06-2X level), the BE no disp and disp contributions cannot be separated. Whenever more than one adsorption geometry was found, all of the BEs are listed and separated by a slash.

between adsorbate molecules in different replicas of the cell that enhance the BE, a phenomenon that does not happen for C₂S** and C₄S**. In fact, we would expect C₄S** to show larger BEs, but its elongated structure makes it challenging to maximize its contact area with the ice surface, hence interacting mostly through the C-end and pointing the S-end upward.

Finally, for molecules of Group IV, even though S₂ is heavier and larger than S**, its interaction with the ice is the weakest of the entire set of the tested molecules. S₂ has an electronic triplet ground state (like O₂) and is the only homonuclear molecule of the set. Its dipole moment is zero, but not its quadrupole moment, which allows for the onset of a quadrupolar interaction between the S–S bond region and a dH of the ice surface. However, its BE is lower (2393 K) than that for S** (3921 K).

### 3.1.2. ONIOM2-corrected BEs

For some species, BEs computed at DFT//DFT have been refined at a CCSD(T) level by applying the ONIOM2 methodology by Dapprich et al. (1999), considering that most of the BE contribution arises from the local interactions between the adsorbed species and a small number of nearby water molecules. These calculations were performed on the following seven S-bearing species: CH₃SH, CS, H₂CS, OCS, SO₂, NS*, and S₂*. These molecules have been chosen because they well represent both closed- and open-shell species and sample the entire range of BEs of the 17 species of this work.

Figure 3 (panel (a)), plots the ONIOM2-corrected BEs against the DFT//DFT ones, showing a very good fit that confirms the accuracy and reliability of the DFT//DFT scheme adopted so far, with the latter slightly overestimated compared to the ONIOM2-corrected one.

### 3.1.3. DFT//HF-3c BEs

To cope with the larger size of the amorphous ice model, we tested the performance of the cheaper DFT//HF-3c scheme on the crystalline ice model to assess its applicability also for the amorphous ice.

Starting from the previously optimized DFT structures, we re-optimized the geometries at HF-3c level. In general, the HF-3c optimized geometries are comparable with the DFT ones, with the exception of very few cases, such as the SO₂ adsorption complex (see below). However, by comparing the HF-3c//HF-3c BEs with the DFT//DFT ones, the resulting correlation is rather coarse (see Figure 3, panel (b), blue symbols and line).

However, a much better agreement is shown by the DFT//HF-3c BEs (the linear correlation improves notably, see Figure 3, panel (b), orange symbols and line). Such an agreement was already shown in Ferrero et al. (2020), and the accuracy of the DFT//HF-3c scheme was also exhibited in molecular crystals (Cutini et al. 2016), polypeptides (Cutini et al. 2017), and pure-silica zeolites (Cutini et al. 2019). These results, thus, extend the application of the DFT//HF-3c scheme to derive accurate BEs for S-bearing molecules adsorbed on ice surfaces.

HF-3c//HF-3c data (blue symbols) of Figure 3 reveal an outlier (for SO₂, around 7000 K), the HF-3c//HF-3c BE being 25% higher than the DFT//DFT one, while the DFT//HF-3c value is about 15% smaller than the DFT//DFT value.

This effect is because HF-3c optimized geometries show shorter distances between the adsorbate and the ice model, which, when computing the single-point energy at the DFT level, will increase both the dispersive interactions and the exchange repulsion contributions. For instance, HF-3c reduces the distance between the S atom of SO₂ and an O atom of the ice surface to 2.25 Å, compared to the 2.48 Å of the B3LYP-D3
<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a) Best linear fit between the DFT//DFT and the ONIOM2-corrected BEs for the crystalline ice systems. (b) Best linear fit between the DFT//DFT BEs and the HF-3c//HF-3c (in blue) and the DFT//HF-3c (in orange) schemes, for the crystalline ice systems. BEs are in kelvin.</figcaption></figure>

(BJ) optimization. Moreover, by looking at the *BE* *disp* values of Table 1, we notice that the percentage of dispersion in DFT//HF-3c BE values is almost always larger than that of DFT//DFT BEs.

Finally, HF-3c can be critical for radical species like NS*. The NS* experimental bond length value is 1.494 Å (Peebles & Marshall 2002), while the calculated ones are 1.480 Å with M06-2X and 1.651 Å with HF-3c. The charge and spin density distributions are also different, with 76% and 59% spin density on nitrogen at HF-3c and at M06-2X, respectively. The higher spin localization of HF-3c causes the lengthening of the bond,

which, when adsorbed, shrinks to 1.633 Å, resulting in a gain in energy of about 10 kJ mol⁻¹ (at DFT//HF-3c). As this behavior is systematic (for both free and adsorbed molecules), it cancels out when computing the BE. In conclusion, HF-3c can be used as a cheap level of theory provided that a careful check is exerted for species with uncommon structures and electronic configurations (especially open-shell ones).

## 3.2. BEs on Amorphous Ice

The amorphous ice shows a richer variety of adsorption sites causing a distribution of BE values for the considered molecules. We characterized eight surface binding sites (dH1-6 as H-bond donors and dO7-8 as H-bond acceptors), which were chosen as a function of their H-bonding ability. We manually built up the starting geometry, exposing the amorphous solid water (ASW) electrophilic regions (due to dangling OH bonds, while their counterparts correspond to exposed O atoms) to the nucleophilic ones of the adsorbates and vice versa (following the principle of electrostatic complementarity). The distance between the adsorbate and the surface was set up as the sum of the van der Waals radii of the corresponding closest atoms, by means of computer graphics manipulation. By proceeding this way, we obtained a total of 136 adsorption complexes. As explained in Section 2.1, the complexes were first optimized at the HF-3c level followed by a harmonic frequency calculation to confirm all structures as minima. Table 2 shows the BEs value at the final DFT//HF-3c level. The mean (⟨BE⟩) and standard deviation (σ) for the BEs are also provided, although we recommend that the reader not consider them as representative of the ensemble of BEs here listed, given the limited sampling of binding sites performed in this study, far from being suitable for a statistical treatment.

From the analysis of our data, it stands out that dH4 and dH6 sites, residing in the cavity domain of the ice model, exhibit the highest BE values due to the close proximity of the surrounding water molecules. The reverse happens for dO7 and dO8, placed on the bottom and as the outermost water molecules of the model, which are the least favorable binding sites.

On average, the dispersion contributions to the BEs are higher when the adsorption takes place at the amorphous ice than at the crystalline ice due to closer contact of the adsorbate with the ice surface. While in the crystalline proton-ordered P-ice model, the H-bond chains extend to infinity ensuring a reinforcement of the H-bonding strength (cooperativity effect), this is not the case for the ASW model. In ASW, the random organization of the water molecules breaks the H-bonding chain extensions, reducing the H-bonding cooperativity and, therefore, its strength. This, in turn, also affects the H-bonding strength exhibited by the chain terminal ice OH groups involved in the interactions with the adsorbates, decreasing the corresponding BE values. This is shown clearly by the higher BE values computed for adsorption on the crystalline P-ice with respect to the ASW model, in agreement with previous literature data (Ferrero et al. 2020). However, other effects, like the easier structure deformation for the amorphous ice, accommodates the adsorbates getting closer to the ice surface, maximizing both the H-bond strength and the dispersive contribution, overcoming the loss of H-bonding cooperativity.

The above interpretation can explain the presence of some very high BE values, which depart from the rest of the BEs of
<table><caption>**Table 2**
Summary of the DFT//HF-3c BEs (in kelvin, and without ZPE Corrections) of the 17 S-bearing Species at the Amorphous Ice Periodic Model</caption><thead><tr><th rowspan="2">Species</th><th colspan="9">Amorphous Ice BE Values</th></tr><tr><th>dH1</th><th>dH2</th><th>dH3</th><th>dH4</th><th>dH5</th><th>dH6</th><th>dO7</th><th>dO8</th><th>⟨BE⟩</th><th>σ</th></tr></thead><tbody><tr><td>CS</td><td>4270</td><td>1852</td><td>3596</td><td>3283</td><td>3380</td><td>3584</td><td>1155</td><td>1648</td><td>2846</td><td>1054</td></tr><tr><td>C<sub>3</sub>S</td><td>5472</td><td>6014</td><td>6122</td><td>5412</td><td>5112</td><td>5941</td><td>1275</td><td>2959</td><td>4788</td><td>1630</td></tr><tr><td>CH<sub>3</sub>SH</td><td>3608</td><td>4991</td><td>3584</td><td>3271</td><td>3271</td><td>4630</td><td>2622</td><td>2850</td><td>3603</td><td>769</td></tr><tr><td>H<sub>2</sub>CS</td><td>4161</td><td>4967</td><td>3416</td><td>3488</td><td>4378</td><td>4053</td><td>1960</td><td>3115</td><td>3692</td><td>859</td></tr><tr><td>H<sub>2</sub>S</td><td>3452</td><td>5208</td><td>3283</td><td>4210</td><td>4630</td><td>2622</td><td>2285</td><td>2309</td><td>3500</td><td>1026</td></tr><tr><td>H<sub>2</sub>S<sub>2</sub></td><td>3007</td><td>4234</td><td>4558</td><td>4462</td><td>3801</td><td>3801</td><td>3175</td><td>4883</td><td>3990</td><td>623</td></tr><tr><td>OCS</td><td>2875</td><td>1491</td><td>2574</td><td>2177</td><td>1660</td><td>3320</td><td>2057</td><td>3175</td><td>2416</td><td>637</td></tr><tr><td>SO<sub>2</sub></td><td>3344</td><td>3175</td><td>3644</td><td>6759</td><td>3680</td><td>4426</td><td>2730</td><td>2442</td><td>3775</td><td>1263</td></tr><tr><td>HCS*</td><td>2117</td><td>2225</td><td>2454</td><td>2117</td><td>1491</td><td>3560</td><td>2478</td><td>2201</td><td>2330</td><td>545</td></tr><tr><td>HS*</td><td>3163</td><td>2454</td><td>2057</td><td>2646</td><td>5003</td><td>2177</td><td>1624</td><td>1251</td><td>2547</td><td>1080</td></tr><tr><td>HS<sub>2</sub>*</td><td>2562</td><td>1527</td><td>2586</td><td>3031</td><td>2670</td><td>2550</td><td>2778</td><td>2694</td><td>2550</td><td>413</td></tr><tr><td>NS*</td><td>3283</td><td>1431</td><td>1828</td><td>1768</td><td>3416</td><td>3428</td><td>1768</td><td>1876</td><td>2350</td><td>805</td></tr><tr><td>C<sub>2</sub>S**</td><td>4426</td><td>3163</td><td>5088</td><td>3572</td><td>2009</td><td>4991</td><td>2009</td><td>2393</td><td>3456</td><td>1193</td></tr><tr><td>C<sub>4</sub>S**</td><td>5075</td><td>2165</td><td>5797</td><td>3716</td><td>4943</td><td>5569</td><td>4306</td><td>4077</td><td>4456</td><td>1097</td></tr><tr><td>S**</td><td>3247</td><td>2081</td><td>2550</td><td>2321</td><td>2526</td><td>2321</td><td>1852</td><td>1852</td><td>2344</td><td>425</td></tr><tr><td>S<sub>2</sub>**</td><td>1924</td><td>1455</td><td>2850</td><td>1552</td><td>2009</td><td>2069</td><td>1648</td><td>1203</td><td>1839</td><td>472</td></tr><tr><td>SO**</td><td>4041</td><td>1082</td><td>2634</td><td>1419</td><td>1696</td><td>2345</td><td>1672</td><td>1696</td><td>2073</td><td>874</td></tr></tbody></table>

Note. For each species, eight BE values, alongside their mean (⟨BE⟩) and standard deviation (σ), are provided.

the corresponding ensembles. This happens for CH<sub>3</sub>SH, C<sub>3</sub>S, H<sub>2</sub>CS, H<sub>2</sub>S, H<sub>2</sub>S<sub>2</sub>, C<sub>2</sub>S**, and C<sub>4</sub>S**, adsorbed in the dH2 and dH5 sites (located at the edges of the top surface). In these cases, a dramatic structural rearrangement of the amorphous ice takes place, with a global energy gain that can reach 2400 K, contributing to the “anomalous” high BEs. In practice, the species/ice interactions cause a large restructuring of the amorphous ice, whose new structure has a lower energy than the initial ice taken as a reference.

This phenomenon is not usual when dealing with rigid structures (e.g., metal surfaces, oxides, and graphene); however, soft matter like water ice is held together by the same kind of interactions keeping the adsorbate attached at the surface. Therefore, the geometrical relaxation induced by the adsorption can rearrange the structure in such a way that, once optimized, it falls in a new minimum with a lower energy with respect to the reference system. A full analysis of this phenomenon has been recently described in Tinacci et al. (2022).

To cope with this “artificial” increment of the BE, we adopted the new ice structure as a reference for the calculation of the BE, which are now coherent with the general BE ensembles. We think this is an important issue, not common in the adsorption phenomena on metallic, ionic, or covalent surfaces, which are dominated by stronger forces than those occurring with soft matter like water ices. This problem would be mitigated only by enlarging the ice model at an increasing computational cost.

### 3.3. Zero-point Energy (ZPE) Corrections

As anticipated in Section 2.1, the ZPE corrections were computed following two different schemes, depending on the adopted surface models. For the adsorption on the crystalline ice, we computed the harmonic frequencies of a fragment consisting of the S-bearing species and the two closest water molecules of the ice with the assumption that deformation occurs only in the proximity of the adsorption site. The ΔZPE is computed following the advice described in Appendix B to arrive at a ZPE-corrected BE, i.e., BE(0). Figure 4 shows a

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Analysis of ZPE corrections. Correlation between DFT//DFT BE(0) and BE for the adsorption at the crystalline ice. The intercept is set to zero.</figcaption></figure>

good linear correlation between BE and BE(0), showing the average ZPE correction is about 14% of the BE.

The final BE(0) values for all of the adsorption complexes at the amorphous ice have, therefore, been computed by multiplying the DFT//HF-3c BEs by a scale factor of 0.862 (see Figure 4, panel (a)).

## 4. Discussion

### 4.1. Comparison with BEs for Non-S-bearing Adsorbates

In the recent work by Ferrero et al. (2020), the adsorption of 17 closed-shell and four open-shell species was simulated. The list included: (i) nondipolar molecules with significant quadrupole moments (i.e., H<sub>2</sub>, N<sub>2</sub>, and O<sub>2</sub>) or higher-order multiple moments (i.e., CH<sub>4</sub>), hence interacting weakly with the ice; (ii) molecules that, besides exhibiting a quadrupole moment, established H-bonds with surface dH atoms (i.e., CO, OCS, and
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** (a) Distribution of BEs computed on the amorphous ice model for S-bearing species (orange) and the species of Ferrero et al. (2020; blue). (b) Distribution of the dispersive (left panel) and nondispersive (right panel) contributions to the crystalline BEs for closed-shell species.</figcaption></figure>

CO₂); (iii) amphiprotic molecules, which are both H-bond acceptors and donors (i.e., NH₃, H₂O, HCl, HCN, and H₂S); (iv) larger molecules that are capable of establishing a greater number of interactions with the surface due to the presence of several functional groups, hence giving rise to the highest BE values of the set (i.e., CH₃OH, CH₃CN, H₂CO, HCONH₂, and HCOOH); and (v) the open-shell species OH*, NH₂*, CH₃, and HCO*, the two former ones forming strong H-bonds with dH and dO sites, at variance with the two latter ones.

Figure 5 compares the BEs of the S-bearing species with those of Ferrero et al. (2020). The distribution of BEs in Ferrero et al. (2020) has a bimodal character (see panel (a)), indicating two different adsorption behaviors in the set of the considered species: a group of low BEs (characterized by dispersive

interactions), and a group of high BEs (amphiprotic and larger species involved in the H-bond). In contrast, S-containing species span a much narrow range of BEs, suggesting a common adsorption feature, despite the diverse chemical structures.

Deeper insights into the contributions to the BE of these two sets of species arise from the analysis of the closed-shell species adsorbed at the crystalline ice model. Each BE can be decomposed in dispersive and nondispersive contributions (see panel (b) of Figure 5) for all considered species.

The distribution of the dispersive contributions highlights two main facts: (i) this term is more important when heavy atoms (like sulfur) belong to the adsorbed species; (ii) the mean value of the distribution for S-bearing species (ca. 3500 K) is higher than that of the nondispersive terms (ca. 2000 K), as dispersion accounts for more than 50% of the BEs provided in this work. Therefore, weak interactions are the most relevant and crucial forces in the description of the adsorption of S-containing species on ices and, thus, they need to be carefully accounted for to properly describe the sulfur chemistry on icy mantles. The above considerations are not surprising due to higher atomic polarizability of second-row atoms compared to C-, N-, and O-containing molecules.

The different behaviors of the two sets of species will have consequences on their chemistry. Considering the BE ensembles, one can argue S-containing species show higher desorption temperatures compared to weakly interacting molecules, like O₂, CH₄, or CO. On the contrary, S-bearing species interact more weakly than amphiprotic molecules with ice. The same inference applies to their surface mobility. As usually, the barriers for diffusion are taken as an arbitrary fraction of the BEs (Cuppen et al. 2017; Kouchi et al. 2020), so we can assume that species with lower BEs will diffuse faster than those with higher BEs, setting up constrains on the reactions that can take place on an icy surface. The particularly low BEs of S₂⁺ (the weakest bounded S-bearing species) also implies its highest mobility. In the hypothesis that polysulphanes could be a reservoir of sulfur, the high mobility of S₂⁺ could favor their encounter and reaction, forming chains and cyclic sulfur species (Druard & Wakelam 2012; Mifsud et al. 2021) that hardly desorb due to the very large dispersive contribution to their BE.

## 4.2. S-bearing versus O-bearing Species

In addition to these general considerations, a punctual comparison between S-bearing and O-bearing species can be drawn for some cases in order to highlight a few exceptions to what we generally understood about our set of molecules.

In Table 3 we compare the contributions to the BE(0) values of some S-bearing species with their oxygen analogs. As a rule of thumb, we would expect molecules containing sulfur to always show a BE(0) lower than the corresponding O-bearing ones, with the dispersive interactions being larger for the former group. However, by considering case by case, we notice that there are some exceptions to this behavior (see Figures 7–9 in Appendix C) For the pair CS/CO, the trend is the opposite, while the OCS/CO₂ pair shows many similarities.

Indeed, CS exhibits a larger BE(0) and a lower percentage of dispersive energy than its O-analog, CO. This can be explained by looking at the dipole moments of these species. CS shows a large dipole moment (1.9 D), at variance with the negligible one of CO. The H-bonding interaction of the two species with a
<table><caption>**Table 3**
Contribution to Selected BE(0) (in kelvin) for S-bearing Species vs. O-bearing Species Adsorption at Crystalline Ice</caption><thead><tr><th colspan="5">Sulfur Bearing Species</th><th colspan="5">Oxygen Bearing Species</th></tr><tr><th>Species</th><th>BE(0)</th><th>No Disp</th><th>Disp (%)</th><th>μ</th><th>Species</th><th>BE(0)</th><th>No Disp</th><th>Disp(%)</th><th>μ</th></tr></thead><tbody><tr><td>CS</td><td>2453</td><td>785</td><td>1668(68%)</td><td>1.9</td><td>CO</td><td>1663</td><td>50</td><td>1613(97%)</td><td>0.1</td></tr><tr><td>CH<sub>3</sub>SH</td><td>4603</td><td>1795</td><td>2808(61%)</td><td>1.5</td><td>CH<sub>3</sub>OH</td><td>7385</td><td>5391</td><td>1994(27%)</td><td>1.7</td></tr><tr><td>H<sub>2</sub>CS</td><td>4267</td><td>1792</td><td>2475(58%)</td><td>1.7</td><td>H<sub>2</sub>CO</td><td>5187</td><td>3268</td><td>1919(37%)</td><td>4.6</td></tr><tr><td>H<sub>2</sub>S</td><td>4033</td><td>2460</td><td>1573(39%)</td><td>1.1</td><td>H<sub>2</sub>O</td><td>7200</td><td>5832</td><td>1368(19%)</td><td>2.0</td></tr><tr><td>OCS</td><td>2405</td><td>120</td><td>2285(95%)</td><td>0.8</td><td>CO<sub>2</sub></td><td>2568</td><td>796</td><td>1772(69%)</td><td>0.0</td></tr></tbody></table>

**Note.** The latter are taken from Ferrero et al. (2020). The BE(0) are decomposed into *no disp* and *disp* contributions, applying the percentage of dispersion (%) obtained for each species. The dipole moment *μ* (in Debye) of the free molecules is computed at the B3LYP-D3(BJ)/A-VTZ level of theory. When more than one BE was present, the average value is reported.

dangling OH is therefore stronger for CS, as it is the dispersive interaction of the S atom with the water molecule of the surface (see Figure 7 in Appendix C).

Both OCS and CO₂ have a linear geometry, an almost null dipole, and interact through an H-bond with a dangling OH of the surface. While OCS exhibits a larger dispersion interaction, for CO₂, the electrostatic component makes up about 30% of the BE, yielding a comparable BE(0) for the two molecules.

On the other hand, the other pairs behave as we expected, as described in Section 4.1. The adsorption geometries show some similarities between the two members of each pair, with the most important differences due to the presence and strength of the H-bonds. Specifically, CH₃SH/CH₃OH and H₂S/H₂O exhibit the same number of H-bonds, while for H₂CS/H₂CO, only the latter establishes an H-bond with the surface through the C=O group.

### 4.3. Comparison with the Literature

Before comparing the present results with literature experimental data, some general considerations are at hand. In the ISM, the water ice mantle is the result of the in situ water formation through surface reactions. Therefore, the final ice structure is the result of the energetic locally hot spots injected in the mantle where intermediates (e.g., OH⁺ and HOO⁻) and final products are formed. As the injected energy is large, it locally heats up the ice, shaping the final structure of the grain in a complex way, difficult to attain through laboratory experiments. In the terrestrial laboratories, ice is prepared directly by deposition from vapor water molecules, with a certain kinetics, on cold metal fingers (Penteado et al. 2017). The final ice can be amorphous (either compact or porous) or crystalline, depending on the kinetics of adsorption and processing temperature. Unfortunately, the detailed structure of amorphous ice can only be partly inferred indirectly by spectroscopic measurements, therefore limiting our knowledge of the atomistic details. The same difficulty occurs for computer modeling: indeed, we are not aware of computer simulations reacting hydrogen and oxygen atoms step by step up to form a sizable icy grain, even if a simpler attempt to form a single water molecule on the forsterite surface has been recently reported (Molpeceres et al. 2019).

As a function of the adopted ice models in computer simulation, the external surface may show a variety of different kinds of adsorption sites, giving rise to different BEs (Song & Kastner 2017). When very few sites are available, like in the proton-ordered ice model adopted here, no more than two different BE values resulted. For the amorphous model, instead, a significant variability of binding sites resulted,

giving rise to an ensemble of BEs. As reported in detail by Ferrero et al. (2020) and not repeated here, extra care should be exerted when comparing the computed BEs with those from experiments.

For instance, while BEs are straightforwardly calculated by computer simulations (see Appendix B), experimental BEs are never directly observed in the thermal desorption, as TPD peaks are usually interpreted through the Redhead (1962) method or more sophisticated techniques, through which desorption adsorption energies (DAEs) are worked out. As DAEs may be a function of the surface coverage *θ*, the comparison with the theoretical BEs is not straightforward (for a recent review, see Minissale et al. 2022). For instance, the ice structure is sensitive to both temperature and adsorption/desorption processes (at variance with surfaces from covalent/ionic solids) as its layers are held together by weaker interactions (e.g., H-bonds in water ices or dispersion in CO ices), which are similar in strength to the BE values. Obviously, ice restructuring processes may affect the final value of the DAEs.

In Table 4, our BEs are compared with those from Penteado et al. (2017; see Figure 6). These authors showed a list of recommended BEs, collecting data from previous works and providing an uncertainty range for each value. For H₂S, OCS, and SO₂, the BEs are derived a posteriori from the experimental TPD measurements of Collings et al. (2004), using the empirical relationship of

$$BE(X): T_{\text{des}}(X) = BE(H_2O): T_{\text{des}}(H_2O) \quad (1)$$

where X is the considered chemical species, and Tdes is the desorption temperature of the species deposited on a water ice extrapolated from the TPD measurements of Collings et al. (2004). The reference BE(H₂O) is set to 4800 K from previous works (see Penteado et al. 2017 for more details). The uncertainty assigned to the BE values of H₂S, OCS, and SO₂ was 3.5%. For the other species, the provided data are based on the work of Hasegawa & Herbst (1993), which in turn were based on previous works, including Allen & Robinson (1977). In the latter, the BE of a molecule was treated as the sum of the BEs of the constituent atoms to infer the interaction energy of several species on SiO₂ surfaces, which was subsequently scaled by an empirical factor to obtain BEs on water ice. The uncertainty of these values is set to half the BE for species whose BE is less than 1000 K, and set to 500 K for all other cases.

The comparison between our computed BEs for the crystalline ice and the literature data shows, as expected,
<table><caption>**Table 4**
Summary of the Computed BE(0) Values (in kelvin) in Comparison with Data from the Literature</caption><thead><tr><th rowspan="3">Species</th><th colspan="5">This Work</th><th colspan="4">Literature</th></tr><tr><th rowspan="2">Crystalline<br>BE(0)</th><th colspan="4">Amorphous</th><th colspan="2">Computed</th><th colspan="2">Database</th><th rowspan="2">Various<br>Penteado</th></tr><tr><th>Min</th><th>Max</th><th>⟨BE⟩</th><th>σ</th><th>Das</th><th>Wakelam</th><th>KIDA</th><th>UMIST</th></tr></thead><tbody><tr><td>CS</td><td>2453</td><td>995</td><td>3680</td><td>2453</td><td>909</td><td>2217</td><td>3200</td><td>3200</td><td>1900</td><td>1800 ± 500</td></tr><tr><td>C<sub>3</sub>S</td><td>4925/5298</td><td>1099</td><td>5277</td><td>4128</td><td>1405</td><td></td><td></td><td>3500</td><td>3500</td><td>3000 ± 500</td></tr><tr><td>CH<sub>3</sub>SH</td><td>4603</td><td>2260</td><td>4302</td><td>3106</td><td>663</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>H<sub>2</sub>CS</td><td>3670/4863</td><td>1690</td><td>4282</td><td>3183</td><td>740</td><td>3110</td><td>4400</td><td>4400</td><td>2700</td><td>2025 ± 500</td></tr><tr><td>H<sub>2</sub>S</td><td>4033</td><td>1970</td><td>4489</td><td>3017</td><td>884</td><td>2556</td><td>2900</td><td>2700</td><td>2743</td><td>2290 ± 90</td></tr><tr><td>H<sub>2</sub>S<sub>2</sub></td><td>4333</td><td>2592</td><td>4209</td><td>3439</td><td>537</td><td>4368</td><td></td><td>3100</td><td>3100</td><td>2600 ± 500</td></tr><tr><td>OCS</td><td>2405</td><td>1286</td><td>2861</td><td>2083</td><td>549</td><td>1571</td><td>2100</td><td>2400</td><td>2888</td><td>2325 ± 95</td></tr><tr><td>SO<sub>2</sub></td><td>5111</td><td>2105</td><td>5827</td><td>3254</td><td>1089</td><td>3745</td><td>5000</td><td>3400</td><td>5330</td><td>3010 ± 110</td></tr><tr><td>HCS*</td><td>3318</td><td>1286</td><td>3069</td><td>2009</td><td>469</td><td>2713</td><td>2900</td><td>2900</td><td>2350</td><td>2000 ± 500</td></tr><tr><td>HS*</td><td>2799</td><td>1078</td><td>4313</td><td>2195</td><td>931</td><td>2221</td><td>2700</td><td>2700</td><td>1500</td><td>1350 ± 500</td></tr><tr><td>HS<sub>2</sub></td><td>2706/3079</td><td>1317</td><td>2613</td><td>2198</td><td>356</td><td>4014</td><td></td><td>2650</td><td>2650</td><td>2300 ± 500</td></tr><tr><td>NS*</td><td>3245</td><td>1234</td><td>2955</td><td>2026</td><td>694</td><td>2774</td><td></td><td>1900</td><td>1900</td><td>1800 ± 500</td></tr><tr><td>C<sub>2</sub>S**</td><td>3442/3877</td><td>1731</td><td>4385</td><td>2979</td><td>1028</td><td>2447</td><td></td><td>2700</td><td></td><td>2500 ± 500</td></tr><tr><td>C<sub>4</sub>S**</td><td>3660</td><td>1866</td><td>4997</td><td>3841</td><td>945</td><td></td><td></td><td>4300</td><td>4300</td><td>3500 ± 500</td></tr><tr><td>S**</td><td>2799</td><td>1579</td><td>2799</td><td>2020</td><td>366</td><td>1428</td><td>2600</td><td>2600</td><td>1100</td><td>985 ± 495</td></tr><tr><td>S<sub>2</sub>**</td><td>1732</td><td>1037</td><td>2457</td><td>1644</td><td>1585</td><td>407</td><td></td><td>2200</td><td>2200</td><td>2000 ± 500</td></tr><tr><td>SO**</td><td>1658</td><td>933</td><td>3483</td><td>2128</td><td>1787</td><td>754</td><td>2900</td><td>2800</td><td>2600</td><td>1800 ± 500</td></tr></tbody></table>

**Note.** The second column contains the BEs on the crystalline P-ice (010) model, the third and the fourth columns contain the minimum and maximum BE values, respectively, found for the amorphous ice, the fifth and the sixth columns list the mean (⟨BE⟩) and the standard deviation (σ), respectively, relative to the amorphous BEs, while columns (7)–(11) report BEs from the literature, respectively, from computed data (Das et al. 2018 and Wakelam et al. 2017), from databases (KIDA; Wakelam et al. 2015; and UMIST; McElroy et al. 2013) and from experiments (Penteado et al. 2017). The BE(0) values were obtained by applying the Equation (0).862*BE = BE(0) (see the text for more details).

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Computed and literature BEs (in kelvin). For each species, the solid black line covers the min/max ZPE-corrected BE values for the amorphous ice model, while blue stars correspond to the ZPE-corrected BEs on the crystalline ice model.</figcaption></figure>

discrepancies. However, comparison of the range provided in Penteado et al. (2017) with our BE ensemble reveals that disagreement is only for S**. The additive principle applied by Allen & Robinson (1977) implies that S<sub>2</sub>** BE must be twice the value of the isolated S** atom. However, we have seen that S** and S<sub>2</sub>** showed a very similar distribution of BEs, due to the same nature of their interaction with the surface. For the other species, the ranges provided by Penteado et al. (2017) lay within our ensembles, especially for H<sub>2</sub>S, OCS, and SO<sub>2</sub>, which have a very narrow distribution of BEs. However, computed BEs for crystalline ice compare poorly with the

average values provided by Penteado et al. (2017), showing that the additive principle is not suitable for the description of the interactions between polyatomic species and extended ice models.

In addition to the experimental data, comparison is carried out with recent BEs from computational works for a large number of species, some of which contain sulfur (see Figure 6).

In Wakelam et al. (2017), a single water molecule was adopted as an ice model (the M06-2X/aug-cc-pVTZ theory level; with neither ZPE nor BSSE corrections). The lost contributions from a real ice surface were somehow recovered by a clever fit between computed BEs and experimental ones for a subset of species. The resultant fitted parameters were used to derive all other BE values. Comparison between our crystalline BEs with those of Wakelam et al. (2017) shows a difference of about 30% of the value. A better agreement with Wakelam et al. (2017) values is shown by the BE range spanned by our amorphous BE ensembles.

Das et al. (2018) computed the BEs for a molecular set including 100 species, adopting a water tetramer and, in few cases, a water hexamer as the ice model, at the MP2/aug-cc-pVDZ theory level, without BSSE and ZPE corrections. Among the S-bearing molecules, only H<sub>2</sub>S and OCS were adsorbed on a water hexamer. The ice models employed in Das et al. (2018) still lack the effect of H-bond cooperativity (present in our periodic models) and could be the reason of the small values of BEs for the S-containing species, even lower than those estimated by Wakelam et al. (2017). Moreover, H<sub>2</sub>S and OCS BEs show an increment when changing the model from the water tetramer to the water hexamer, and it is not clear for how long this trend would continue before reaching a converged value. The Das et al. (2018) values are within the BE ensembles computed for the amorphous ice. Altogether,
this indicates that the choices of the ice model and of the
computational method have a strong relevance on the predicted
BEs. Thus, more studies are needed to fully characterize BE
distributions, by using larger models and exploring more
binding sites.

4.4. *Binding Energies and Snow Lines in Protoplanetary Disks*

As mentioned in Section 1, the BEs of the species on ice
mantles have a profound impact on the interstellar chemistry
and, more specifically, on the chemical composition of the
regions where planetary systems, like our own, are formed. A
very important case is represented by the gaseous versus solid
chemical composition in protoplanetary disks, which are the
sites where planets, asteroids, and comets eventually form. If a
species is in the solid form, then it will likely be incorporated in
the (aforementioned) formation of rocky objects, whereas, if it
is in the gaseous form, it will only enrich the giant gaseous
forming planets of the elements contained in the species (e.g.,
Öberg & Bergin 2021). Sulfur, on the other hand, seems to be a
very important element for terrestrial life and, perhaps, even a
key element for its emergence (e.g., Shalayel et al. 2020).

As discussed in Section 1, the major reservoirs of sulfur in
cold molecular clouds (Section 1.2, from where the evolution
toward a planetary system starts) are unknown, though they are
suspected to be in solid organosulfur compounds and
polysulphates (Druard & Wakelam 2012; Laas & Caselli 2019).
During the evolution, the quantity of sulfur in the gaseous
phase changes by orders of magnitude, so that it is very
important to know it in the protoplanetary disks. So far, there
are no measurements for solid, and only very few for gaseous,
S-bearing species in protoplanetary disks (see below). The
locus where the transition of a species from the solid to the gas
phase takes place is usually called “snow-line” in astrophysical
jargon, particularly so in the protoplanetary disks, which are
objects that are relatively cold and dense (e.g., Öberg &
Bergin 2021). Then, one efficient way to measure the solid
versus gaseous content of S-bearing species in protoplanetary
disks is to observe where the species disappear from the gas
phase, because they freeze onto the dust grain water-rich
mantles, namely, to measure their snow lines. Two major
parameters determine where the snow line lies: the temperature
of the dust, and the BE of the species. In general, the larger the
distance from the disk center (heated by the central forming
star), the colder the dust. Therefore, species with lower BEs
remain gaseous at larger disk radii.

The ALMA-DOT survey of young disks provides the first
preliminary observational constraints on the position of the
snow lines of some of the S- and O-bearing molecules
considered in this study (Podio et al. 2020a). In particular,
ALMA-DOT observed a protoplanetary disk, IRAS 04302
+2247, seen edge-on, which makes it ideal to constrain the
snow lines of different species. These observations show that
the molecular emission is bright in an intermediate disk layer
(the so-called molecular layer; Aikawa & Herbst 1999), while it
decreases in the disk midplane due to molecular freeze-out onto
the dust grains. More relevant to this work, CS, H₂CO, and
H₂CS show a similar spatial distribution, and a sharp decrease
of their emission is observed at smaller radii than CO (Codella
et al. 2020; Podio et al. 2020b). The estimated snow line is
at 25 au for H₂CO (and consequently H₂CS and CS) and 100 au
for CO (van't Hoff et al. 2020). The larger snow-line radius of
CS with respect to that of CO is generally in agreement with

the two estimated BEs: 995–3680 K (CS; this work) versus
1109–1869 K (CO; Ferrero et al. 2020). On the other hand, the
observations suggest similar snow lines for CS, H₂CO, and
H₂CS, while the estimated BEs of H₂CO and H₂CS are larger
than that of CS: 3071–6194 K (H₂CO; Ferrero et al. 2020) and
1690–4282 K (H₂CS; this work) versus 995–3680 K (CS; this
work). We emphasize, however, that the ALMA-DOT
observations are unable to resolve the inner 25 au of the disk.
Therefore, higher spatial resolution observations are needed to
constrain the snow lines of species with large BEs, such as
H₂CO and H₂CS.

In conclusion, the admittedly scarce observations so far
available seem to be in agreement with the estimated BEs
obtained through theoretical quantum chemical calculations as
provided in Ferrero et al. (2020) and this work, which
encourage us to pursue the work to a larger sample of species,
both S- and N-bearing and the so-called interstellar complex
organic molecules (Ceccarelli et al. 2017).

5. Conclusions

In the present work, the BEs of 17 astrochemically relevant
S-bearing species adsorbed on periodic water ice models have
been computed using DFT based on the hybrid B3LYP-D3(BJ)
(for closed-shell species) and hybrid meta-GGA M06-2X (for
open-shell species) with an Ahlrichs-VTZ<sup>a</sup> double-zeta polar-
ized Gaussian basis sets. We adopted either a crystalline
proton-ordered ice model or a more realistic amorphous ice
model, both treated within the periodic boundary conditions.
The adopted DFT results have also been refined toward a
CCSD(T) level through the ONIOM2 approach carried out on
the crystalline model, showing a very good convergence
between the DFT BEs and those refined with the ONIOM2
scheme. As full DFT is too costly for the simulation of
adsorption on the amorphous ice, we checked the reliability of
the cheaper HF-3c method by contrasting the BEs with that at
the DFT level for the crystalline ice. While HF-3c BEs poorly
correlate with the DFT ones, the BEs at the single-point DFT//
HF-3c level give an excellent correlation with full DFT values.
The validated DFT//HF-3c method was used to study the
adsorption process on the amorphous ice model, in which the
richness of adsorption sites (up to eight) implies a BE ensemble
for each species. We also showed the ZPE correction to the
electronic BE to be around 14% of the BE itself irrespective of
the adsorbed species, and provided a scale factor to obtain the
BE(0) for adsorption complexes at the amorphous ice model.
This is important, as it allows for skipping the expensive
frequency calculation by simple scaling of the electronic BEs.
The analysis of the obtained data showed that, in order to
properly describe the behavior of S-bearing species, special
care must be taken in the description of the dispersive
interactions, which gain more importance as the number of
large and heavy atoms, like sulfur, increases. When considering
the amorphous ice model, which provides enough richness of
sites, the resulting BE ensembles match both the experimental
and the literature computed BEs. We posit BE distributions
rather than single BE values to be essential for providing robust
parameters to be adopted in numerical models of chemical
evolution in the universe.

Regarding the sulfur depletion problem, our calculations do
not solve the open questions. Nonetheless, we provide strong
and relatively accurate BE values, which can elucidate the
features of S-bearing species in comparison with S-free
molecules when adsorbed on icy mantles. An important clue is provided by the $S_2^*$ case in which the lowest BE implies higher diffusivity and increasing reactivity toward the formation of higher-weight S-containing rings strongly bound to the ice through dispersion interactions. The CS molecule also shows some counterintuitive results when compared to the analog CO molecule, due to a higher dipole moment and dispersion contribution. The present work also highlights the utility of accurately computing BE values to be correlated with snow lines in protoplanetary disks.

This project has received funding within the European Union's Horizon 2020 research and innovation program from the European Research Council (ERC) for the projects “Quantum Chemistry on Interstellar Grains” (Quantumgrain), grant agreement No. 865657 and “The Dawn of Organic Chemistry” (DOC), grant agreement No. 741002, and from the Marie Skłodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement No. 811312. The Italian Space Agency for co-funding 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. A.R. is indebted to the “Ramón y Cajal” program. The authors thank Gretobape for fruitful and stimulating discussions. The authors wish to thank the anonymous reviewers for the valuable suggestions. Supplementary Material consisting of (i) the fractional coordinates of HF-3c adsorption structure optimized for both the crystalline and amorphous ice, (ii) images of the adsorption features at crystalline periodic ice models, in which electrostatic potential maps, spin density maps (for open-shell systems), and vibrational features are displayed, (iii) a pdf file with a thorough guide to the computation of BEs and the basis sets employed for the calculations, is available on Zenodo: https://doi.org/10.5281/zenodo.6798922.

## Appendix A
### Computational Details

For optimization and frequency calculations, we refer to the computational parameters adopted in the work of Ferrero et al. (2020): the threshold parameters for the evaluation of the Coulomb and exchange bielectronic integrals (TOLINTEG) were set to 7, 7, 7, 7, and 14 when working with crystalline ice, while for the amorphous ice to 7, 7, 7, 7, and 25, which facilitates the convergence of the self-consistent field (SCF), especially when adsorbing open-shell species.

Some species required a careful treatment to describe properly their electronic state. The use of the BROYDEN accelerator (Johnson 1988) and the broken-(spin)-symmetry ansatz (Neese 2004, 2009) was necessary to obtain an optimal convergence of the SCF procedure.

For the numerical evaluation of bielectronic integrals (Becke 1988), the standard pruned grid, composed of 75 radial points and a maximum of 974 angular points, was used. To sample the reciprocal space, we adopted the Pack–Monkhorst mesh (Pack & Monkhorst 1977), with a shrinking factor of 2, which generates 4 k points in the first Brillouin zone. The only exception is for the adsorption of atomic sulfur on the crystalline ice model, which required a shrinking factor of 3, generating 5 k points, in order to avoid numerical noise in the results.

The Broyden–Fletcher–Goldfarb–Shanno algorithm (Broyden 1970; Fletcher 1970; Goldfarb 1970; Shanno 1970) was used to carry out geometric optimizations.

## Appendix B
### Calculation of the Binding Energies

When using a finite basis set of localized Gaussian functions to describe our systems, BSSEs arise. This is the case in our BEs computed at both DFT//DFT and DFT//HF-3c. Thus, the a posteriori counterpoise (CP) correction of Boys & Bernardi (1970) was applied to compensate for this error, and our BEs are therefore defined as the opposite of the CP-corrected interaction energies ($\Delta E^{CP}$), that is,

$$BE = -\Delta E^{CP} \qquad (B1)$$

$$\Delta E^{CP} = \Delta E - \text{BSSE} \qquad (B2)$$

where the non-CP-corrected interaction energy $\Delta E$ is given by the sum of the deformation-free interaction energy ($\Delta E^*$), the deformation energy of the slab ($\delta E_S$), the deformation energy of the molecule ($\delta E_\mu$), and the lateral interaction ($E_L$) between adsorbate molecules in different replicas of the cell. This quantity corresponds to the common definition of interaction energy, which in this case is the energy of the complex to which energies of the isolated species and isolated ice model are subtracted.

$$\begin{aligned} \Delta E &= \Delta E^* + \delta E_S + \delta E_\mu + E_L \\ &= E_{\text{complex}} - E_{\text{ice}} - E_{\text{species}}. \end{aligned} \qquad (B3)$$

For the crystalline systems, inclusion of the ZPE corrections to the BE are performed according to the equation $BE(0) = BE - \Delta ZPE$, where

$$\Delta ZPE = ZPE_{\text{complex}} - ZPE_{\text{ice}} - ZPE_{\text{species}} \qquad (B4)$$

For the amorphous case, inclusion of ZPE corrections was done by applying a scaling factor derived from a linear correlation between BE and $BE(0)$ in the crystalline systems, that is, $BE(0) = 0.862^*BE$.

Finally, to check for accuracy and refine DFT BEs for our crystalline ice model, the single- and double-electronic excitation coupled-cluster method with an added perturbative description of triple excitations (CCSD(T)) was used, in combination with a correlation consistent basis set extrapolated to the largest possible basis set. These calculations have been performed with the Gaussian16 software package (Frisch et al. 2016). These refinements were performed through the ONIOM2 approach (Dapprich et al. 1999), dividing the systems into two parts (model and real systems), which are described by two different levels of theory (high and low). The model system (represented by the adsorbate and the two closest water molecules) was described by the high level of theory, CCSD(T). The real system (that is, the whole system) was described by the corresponding DFT (low) level of theory.

The ONIOM2-corrected BE is:

$$\begin{aligned} \text{BE}(\text{ONIOM2}) = \text{BE}(\text{low, real}) \\ + \text{BE}(\text{high, model}) - \text{BE}(\text{low, model}). \end{aligned} \qquad (B5)$$

In this way, the $\Delta E = \text{BE(high, model)} - \text{BE(low, model)}$ represents the correction term to the energy of the real system due to the improved description at the high level. This
procedure ensures that if the low level of theory is improved to arrive at the high level of theory or if the model system is enlarged to become the real system, the energy of the whole system is that at the higher level of theory.

In this work, for the calculation of the ONIOM2-corrected BEs, BE(ONIOM2), Equation (B5) can be rewritten as

$$ BE(ONIOM2) = BE(DFT; \mu - \text{ice}) + BE(CCSD(T); \mu - 2H_2O) - BE(DFT; \mu - 2H_2O) \quad (B6) $$

where BE(DFT; μ - ice) is the BE computed at DFT//DFT, while the BEs of the model system (μ - 2H₂O) are computed through single-point energy calculations at CCSD(T).

A thorough guide to the computation of BEs is available in the Supplementary Material.

## Appendix C
S-bearing versus O-bearing Species

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Adsorption geometry of S-bearing vs. O-bearing species at P-ice (010) surface. Side view is at the top; the top view is at the bottom. Complexes on the right (O-bearing species) are adapted from Ferrero et al. (2020). Distances are given in angstrom units. (i) CH₃SH vs. CH₃OH; (ii) CS vs. CO.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Adsorption geometry of S-bearing vs. O-bearing species at P-ice (010) surface. Side view is at the top; the top view is at the bottom. Adsorption complexes on the right (O-bearing species) are adapted from Ferrero et al. (2020). Distances are given in angstrom units. (i) First and (ii) second adsorption geometry of H₂CS vs. H₂CO.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Adsorption geometry of S-bearing vs. O-bearing species at P-ice (010) surface. Side view is at the top; the top view is at the bottom. Adsorption complexes on the right (O-bearing species) are adapted from Ferrero et al. (2020). Distances are given in angstrom units. (i) H₂S vs. H₂O; (ii) adsorption geometry of OCS vs. CO₂.</figcaption></figure>

ORCID iDs

Jessica Ferrero https://orcid.org/0000-0003-2161-9120

Joan Enrique-Romero https://orcid.org/0000-0002-2147-7735

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

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

Linda Podio https://orcid.org/0000-0003-2733-5372

Claudio Codella https://orcid.org/0000-0003-1514-3074

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

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

References

Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233

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

Allen, M., & Robinson, G. W. 1977, ApJ, 212, 396

Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197

Becke, A. D. 1988, PhRvA, 38, 3098

Becke, A. D. 1993, JChPh, 98, 5648

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

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

Brandenburg, J. G., Alessio, M., Civalleri, B., et al. 2013, JPCA, 117, 9282

Broyden, C. G. 1970, JApMa, 6, 76

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

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

Ceccarelli, C., & Dominik, C. 2005, A&A, 440, 583

Cernicharo, J., Cabezas, C., Agúndez, M., et al. 2021, A&A, 648, L3

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

Codella, C., Podio, L., Garufi, A., et al. 2020, A&A, 644, A120

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

Cooke, I. R., Öberg, K. I., Fayolle, E. C., Peeler, Z., & Bergner, J. B. 2018, ApJ, 852, 75

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

Cutini, M., Civalleri, B., Corno, M., et al. 2016, J. Chem. Theory Comput., 12, 3340

Cutini, M., Civalleri, B., & Ugliengo, P. 2019, ACS Omega, 4, 1838
Cutini, M., Corno, M., & Ugliengo, P. 2017, J. Chem. Theory Comput., 13, 370

Dapprich, S., Komáromi, I., Byun, K. S., Morokuma, K., & Frisch, M. J. 1999, Comput. Theor. Chem., 461, 1

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

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

Douglas, A. E., & Herzberg, G. 1942, CJRes, 20a, 71

Dovesi, R., Erba, A., Orlando, R., et al. 2018, WIREs Comput. Mol. Sci., 8, e1360

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Druard, C., & Wakelam, V. 2012, MNRAS, 426, 354

Dulieu, F., Congiu, E., Noble, J., et al. 2013, NatSR, 3, 1338

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

Enrique-Romero, J., Ceccarelli, C., Rimola, A., et al. 2021, A&A, 655, A9

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

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

Fletcher, R. 1970, CompJ, 13, 317

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

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

Germain, A., Tinacci, L., Pantaleone, S., Ceccarelli, C., & Ugliengo, P. 2022, ESC, 6, 1286

Ghesquiere, P., Mineva, T., Talbi, D., et al. 2015, PCCP, 17, 11455

Goldfarb, D. 1970, MaCom, 24, 23

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

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

Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83

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

He, J., Emtiaz, S. M., & Vidali, G. 2017, ApJ, 837, 65

He, J., Emtiaz, S. M., & Vidali, G. 2018, ApJ, 863, 156

Hily-Blant, P., Pineau des Forêts, G., Faure, A., & Lique, F. 2022, A&A, 658, A168

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

Jiménez-Escobar, A., & Muñoz Caro, G. M. 2011, A&A, 536, A91

Johnson, D. D. 1988, PhRvB, 38, 12807

Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114

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

Keller, L. P., Hony, S., Bradley, J. P., et al. 2002, Natur, 417, 148

Kouchi, A., Furuya, K., Hama, T., et al. 2020, ApJL, 891, L22

Kruse, H., & Grimme, S. 2012, JChPh, 136, 154101

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

Lauck, T., Karssemeijer, L., Shulenberger, K., et al. 2015, ApJ, 801, 118

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

Mason, N. J., Nair, B., Jheeta, S., & Szymanska, E. 2014, FaDi, 168, 235

Maté, B., Cazaux, S., Satorre, M. Á., et al. 2020, A&A, 643, A163

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

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

McKellar, A. 1940, PASP, 52, 187

Mifsud, D. V., Kaňuchová, Z., Herczku, P., et al. 2021, SSRv, 217, 14

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

Mispelaer, F., Theulé, P., Aouididi, H., et al. 2013, A&A, 555, A13

Molinari, S., Ceccarelli, C., White, G. J., et al. 1999, ApJL, 521, L71

Molpeceres, G., Rimola, A., Ceccarelli, C., et al. 2019, MNRAS, 482, 5389

Neese, F. 2004, JPCS, 65, 781

Neese, F. 2009, Coord. Chem. Rev., 253, 526

Noble, J. A., Theule, P., Mispelaer, F., et al. 2012, A&A, 543, A5

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

Pack, J. D., & Monkhorst, H. J. 1977, PhRvB, 16, 1748

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

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

Pascale, F., Zicovich-Wilson, C. M., Lopez Gejo, F., et al. 2004, JCoCh, 25, 888

Peebles, L. R., & Marshall, P. 2002, CPL, 366, 520

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

Phuong, N. T., Chapillon, E., Majumdar, L., et al. 2018, A&A, 616, L5

Podio, L., Garufi, A., Codella, C., et al. 2020a, A&A, 644, A119

Podio, L., Garufi, A., Codella, C., et al. 2020b, A&A, 642, L7

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

Rae, J. G. L., Hartquist, T. W., Lepp, S. H., O'Neill, P. T., & Williams, D. A. 2004, A&A, 413, 1

Redhead, P. 1962, Vacuu, 12, 203

Rimola, A., Ceccarelli, C., Balucani, N., & Ugliengo, P. 2021a, FrASS, 8, 655405

Rimola, A., Ferrero, S., Germain, A., Corno, M., & Ugliengo, P. 2021b, Mine, 11, 26

Rodríguez-Almeida, L. F., Jiménez-Serra, I., Rivilla, V. M., et al. 2021, ApJL, 912, L11

Ruf, A., Bouquet, A., Boduch, P., et al. 2019, ApJL, 885, L40

Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691

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

Shalayel, I., Youssef-Saliba, S., Vazart, F., et al. 2020, EurJOC, 2020, 3019

Shanno, D. F. 1970, MaCom, 24, 647

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

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

Song, L., & Kästner, J. 2017, ApJ, 850, 118

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

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

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

Terada, H., & Tokunaga, A. T. 2012, ApJ, 753, 19

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

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

van't Hoff, M. L. R., van Dishoeck, E. F., Jørgensen, J. K., & Calcutt, H. 2020, A&A, 633, A7

Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514

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

Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, ApJS, 217, 20

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

Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035

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

Woods, P. M., Occhiogrosso, A., Viti, S., et al. 2015, MNRAS, 450, 1256

Woon, D. E., & Herbst, E. 2009, ApJS, 185, 273

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. 2019a, JChPh, 150, 064702

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

Zhao, Y., & Truhlar, D. G. 2008, Theor. Chem. Acc., 120, 215

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