The Astrophysical Journal, 938:158 (17pp), 2022 October 20 https://doi.org/10.3847/1538-4357/ac9278 © 2022. The Author(s). Published by the American Astronomical Society. Binding Energies of Interstellar Relevant S-bearing Species on Water Ice Mantles: A Quantum Mechanical Investigation Jessica Perrero1,2 , Joan Enrique-Romero1 , Stefano Ferrero1 , Cecilia Ceccarelli3 , Linda Podio4 , Claudio Codella3,4 , 1 Albert Rimola1 , and Piero Ugliengo2 Departament de Química, Universitat Autònoma de Barcelona, Bellaterra, E-08193, Catalonia, Spain; albert.rimola@uab.cat 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 3 Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), F-38000 Grenoble, France 4 INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy 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 fields are highly diminished, and (ii) the formation of thick icy mantles coating the surfaces of dust grains, due to the Stars like our Sun begin their journey in the interstellar medium adsorption and in situ synthesis of species and their subsequent (ISM), where star formation takes place in dense (103–104 cm−3) hydrogenation. The resulting icy mantles are dominated by and cold (10 K) regions, the so-called molecular clouds. In the water, but they also contain other volatile species like CO, past, it was thought that the ISM was too harsh an environment for CO2, CH3OH, and NH3, and hence are referred to as “dirty molecules to survive and thrive. However, almost 80 yr ago, the ices” (Boogert et al. 2015). first diatomic molecules (the methylidyne radical CH, its cation The phenomena occurring on the icy mantle surfaces are of CH+, and the cyano radical CN) were detected by means of optical paramount importance from a chemical point of view. Indeed, and ultraviolet spectroscopy (Swings & Rosenfeld 1937; McKellar ice mantles are thought to facilitate the occurrence of reactions 1940; Douglas & Herzberg 1942), which indeed indicated that taking place at their surfaces, primarily exhibiting three chemistry is present and has an important role in the ISM. At present, we count more than 270 detected gas-phase functions: (i) as pre-concentrators of chemical species, species (McGuire 2021), a number that constantly increases especially relevant in the low-density environments of given the enhanced performances of the observational facilities molecular clouds, where the rate of collision between gas- (e.g., Yebes 40 m, IRAM 30 m, NOEMA, the Atacama Large phase species is extremely small, (ii) as chemical catalysts, by Millimeter/submillimeter Array (ALMA), and Green Bank decreasing the activation energies of the reactions, and hence Telescope (GBT)), and the expectation of the new results overcoming the energy barriers under the interstellar conditions provided by the James Webb Space Telescope. As a (Zamirri et al. 2019b), and (iii) as third bodies, by absorbing the consequence of these thrilling discoveries, curiosity concerning large energy excess released by exothermic reactions, without the chemistry of the ISM has been growing more and undermining the stability of the newly formed products more keen. (Pantaleone et al. 2020, 2021). The first and the third points The characteristic densities and temperatures of molecular are the keys to successful recombination reactions. clouds allow (i) the molecules to survive, since the external UV 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 Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further binding energy (BE) of the adsorbed species. The latter distribution of this work must maintain attribution to the author(s) and the title quantity defines how strong the interaction between the species of the work, journal citation and DOI. and the surface is, and dictates the chance that a given adsorbed 1 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. species can be ejected into the gas phase: the higher the BE, the missing components of the BEs through a clever fitting higher the temperature required to thermally desorb. In the procedure against selected experimental BE values. In Das current astrochemical models, BE values are important input et al. (2018), the adoption of clusters from one to six water parameters, not only to consider likely desorption events (e.g., molecules allowed them to study the BE of 16 species, with the Penteado et al. 2017) but to simulate diffusion processes on the puzzling result that water tetramer was found to be the best surfaces, as the diffusion barriers are usually assumed to be a model for ice. fraction of the BE (Mispelaer et al. 2013; Karssemeijer & One way to generate computer models for amorphous ices is Cuppen 2014; Ghesquiere et al. 2015; Lauck et al. 2015; by using water clusters cut out from a crystalline ice, heating Cuppen et al. 2017; He et al. 2017, 2018; Cooke et al. 2018; them up by molecular dynamics runs and by quenching at low Kouchi et al. 2020; Maté et al. 2020). Therefore, BEs are temperatures (Shimonishi et al. 2018; Rimola et al. 2021b). crucial parameters in determining the ISM chemistry and its Another approach is to refer to a model repeated in space by the resulting composition (Penteado et al. 2017; Wakelam et al. periodic boundary conditions. This allows for the simulation of 2017; Ferrero et al. 2020; Enrique-Romero et al. 2021). 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 1.1. Previous Studies of the BEs starting from a large unit cell of a crystalline ice and Determining BEs is a challenging task (Minissale et al. performing a number of heating/freezing cycles to arrive to 2022). From an experimental point of view, effective and an amorphous ice model. Another strategy, recently proposed coverage-dependent BEs are obtained via temperature pro- by the ACO-FROST scheme (Germain et al. 2022), consists of grammed desorption (TPD) experiments. TPD consists of two growing icy clusters by a random step-by-step addition of steps: first, the substrate, maintained at a constant temperature, molecular water molecules. is exposed to the adsorption of the species; then, the More recently, Ferrero et al. (2020) computed a large set of temperature is increased until the desorption of the species, BEs for 21 astrochemically relevant species, on both crystalline which is collected and analyzed by a mass spectrometer. The and amorphous surfaces of water ice. Both models ensured that BEs are then extracted by applying the direct inversion method hydrogen-bond cooperativity (of fundamental relevance to through the Polanyi–Wigner equation (e.g., Dohnálek et al. obtain accurate BE values) is fully taken into account, at 2001). Despite its great usefulness, this technique shows two variance with the small ice clusters adopted by Wakelam et al. limitations: (i) the obtained BEs depend not only on the (2017) and Das et al. (2018). For the crystalline icy model, only morphology and the composition of the substrate but also on very few adsorption sites were available, limiting the BEs to the regimes in which the experiments are performed, i.e., one or two distinct values per species. For the amorphous ice submonolayer, monolayer, or multilayer (e.g., Noble et al. model, however, the complex surface morphology allowed the 2012; He et al. 2016; Chaabouni et al. 2018) and will return calculation of multiple BE values for each species, resulting in coverage-dependent BEs, and (ii) it measures a desorption a BE range better describing the variety of binding sites enthalpy, a quantity that is equal to the BE only in the absence expected on real amorphous ice mantles. of other activated processes, like ice restructuring, usually In the present work, we apply the methodology of Ferrero assumed to be negligible (He et al. 2016). et al. (2020) in order to enlarge the set of BEs. The focus here is Despite the numerous literature studies that have investi- on the adsorption of 17 interstellar relevant S-bearing species, gated the desorption processes by means of TPD (e.g., Collings namely: CS, C2S••, C3S, C4S••, CH3SH, H2CS, HS•, H2S, HS •2 , et al. 2004; Noble et al. 2012; Dulieu et al. 2013; Fayolle et al. H2S2, NS•, OCS, S••, S•• •• 2 , SO , and SO2 (dots indicate the 2016; He et al. 2016; Smith et al. 2016), only a limited variety unpaired electrons). The chosen species satisfy two criteria: (i) of molecules has been considered, in particular, stable closed- they are neutral species, since positive ions become neutral shell species, in contrast to the large variety of interstellar when landing on negatively charged dust grains (Walmsley species. Moreover, the employed substrates are almost et al. 2004; Ceccarelli & Dominik 2005; Rimola et al. 2021a) invariably water ices. as long recognized (Draine & Sutin 1987; Rae et al. 2004; A possible alternative to TPD is to adopt state-of-the-art Mason et al. 2014), and (ii) they envisage at most six atoms to computational chemistry to simulate the adsorption process. ensure relatively compact structures needed to properly probe The simulation can handle both closed- and open-shell species, the largest number of adsorption sites on ice models with allotting some extra care for the latter. The main limitation is relatively small unit cells. the trade-off between the system size (the icy grain model) and the accuracy of the calculation. A number of computational 1.2. The Sulfur Depletion studies have reported the BEs of important astrochemical species like H, H2, N, O, HF, CO, and CO2 on icy surfaces The choice of studying the BEs of S-bearing species is due to modeled by periodic/cluster crystalline/amorphous systems a long standing issue in the field of astrochemistry: the sulfur (e.g., Al-Halabi & Van Dishoeck 2007; Karssemeijer & depletion problem. In dense clouds, sulfur is severely depleted Cuppen 2014; Zamirri et al. 2017; Shimonishi et al. 2018; from the gas phase (e.g., Tieftrunk et al. 1994; Ruffle et al. 1999; Zamirri et al. 2019a; Bovolenta et al. 2020; Ferrero et al. 2020). Wakelam et al. 2004; Phuong et al. 2018; Vastel et al. 2018; Other works addressed the BEs of a large number of species, van’t Hoff et al. 2020) by more than 2 orders of magnitude with adopting very approximated models for the substrate. Two respect to the solar system S abundance ([S]/[H] = 1.8 × 10−5; examples are the works of Wakelam et al. (2017) and Das et al. Anders & Grevesse 1989; Woods et al. 2015). Therefore, sulfur (2018), in which the ice surface simulated by water clusters of may freeze out on dust grains after its hydrogenation, or react minimal nuclearity allowed the prediction of the BEs for 100 with other species in a similar way as atomic oxygen. H2S was species. In Wakelam et al. (2017), a single water molecule was expected to be the main sink of sulfur on dust grains (Garrod adopted to represent the ice surface, while recovering the et al. 2007; Jiménez-Escobar & Muñoz Caro 2011), but it has 2 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. never been directly detected on interstellar ices yet. The only (molecules) to three (solid bulks) dimensions, it allows us to species detected in the ice mantles are OCS and, tentatively, SO2 rigorously define true slabs to modeling surfaces without the (Boogert et al. 2015). Accordingly, the most important S-carriers need of artificial replica along the direction perpendicular to the are still unknown. Nowadays, the most likely species that act as slab plane, as usually implemented in codes adopting plane reservoirs of S are thought to be organosulfur compounds (Laas waves basis sets. & Caselli 2019) and polysulphanes (H2Sn; Druard & Wakelam On the crystalline model, the hybrid B3LYP 2012). These are refractory species that, once trapped on the (Becke 1988, 1993; Lee et al. 1988) and the hybrid meta- mantles or in the core of the grains, cannot desorb and, therefore, GGA M06-2X (Zhao & Truhlar 2008) density functionals were become undetectable (Woods et al. 2015). In the last three years, used in the geometry optimizations and in computing the BEs seven new gas-phase S-bearing species were identified, most of of closed-shell and open-shell species, respectively. To account them in TMC-1 (Cernicharo et al. 2021; Rodríguez-Almeida for dispersion interactions, the hybrid B3LYP functional was et al. 2021), and characterized by a C–S chemical bond. combined with the Grimme’s D3 empirical correction with the A recent review of the laboratory experiments performed on Becke–Johnson (BJ) damping scheme (i.e., B3LYP-D3(BJ) the chemistry of sulfur in the condensed phase has highlighted (Grimme et al. 2010, 2011; Sure & Grimme 2013). The hybrid that the formation of SO2, SO3, hydrates of H2SO4, and related meta-GGA M06-2X functional was used together with a spin- species is very common when performing photolysis, proton- unrestricted formalism (Pople et al. 1995) to model the irradiation, and radiolysis of mixed ices. However, neutral– adsorption of open-shell species (with one or two unpaired neutral reactions have yet to be deeply explored. Moreover, electrons) because it better describes their adsorption properties gas-phase chemistry is likely to be a major contributor in sulfur (Ferrero et al. 2020) due to its higher percentage of exact astrochemistry (Mifsud et al. 2021). exchange (54% compared to the 20% of B3LYP). Note that for The formation of organosulphur molecules was observed M06-2X, the D3(BJ) correction for dispersion was not applied experimentally for the first time in the work of Ruf et al. because, for this functional, the description of the dispersion (2019), where a 2:1:1 mixture of H2O:CH3OH:NH3 ice was component is already included in its definition. For this reason, bombarded with S7+ ions, although showing a possible source the contribution of the dispersion component to the BEs can be of compounds relative to more evolved environments, like icy worked out only for the closed-shell species (computed at moons, Kuiper Belt objects, comets, and their building blocks. B3LYP-D3(BJ)). The two density functional theory (DFT) In diffuse clouds, sulfur is mostly present in its ionized methods have been used in combination with the Ahlrichs atomic form S+, and close to the cosmic abundance triple zeta valence quality basis set supplemented with a double (Jenkins 2009), it progressively depletes as the cloud evolves set of polarization functions (i.e., A-VTZ*; Schäfer et al. 1992). (Hily-Blant et al. 2022). Sulfur depletion in protoplanetary Thus, for the crystalline models, the BEs were computed from disks and later stages can be explained by the presence of the DFT energies that resulted from the DFT optimized refractory sulfide minerals, such as FeS, in the grains (Keller structures. This scheme is hereafter referred to as DFT//DFT. et al. 2002; Kama et al. 2019). All of the BEs were corrected for the basis set superposition With the present work, thus, we aim to contribute to this puzzling subject by providing a set of accurate BE values of error (BSSE) through the counterpoise correction method by several S-bearing species, which can be used in numerical Boys & Bernardi (1970). modeling studies aimed at rationalizing observations and Finally, in order to check for the accuracy of the DFT//DFT laboratory experiments. 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 2. Methodology a number of species representative of the set. Given the Since the calculations follow the same procedure as Ferrero dependency of CCSD(T) method on the basis set, the Dunning et al. (2020), we refer to this work for the equations used here (1989) correlation consistent family of cc-pVnZ basis sets to compute the BEs of the S-bearing species set. Moreover, the (where n is 2, 3, 4, 5,... and corresponds to a ζ of D, T, Q, 5,K) reader can find details on these expressions in Appendix A. A was employed. This allows us to perform an extrapolation of thorough guide to the computation of BEs is available in the the desired property, here the BE, by plotting it against 1/n3. Supplementary Material. 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 2.1. Computational Details the exception of SO2, for which the aug-cc-pVnZ basis set We adopted a periodic approach to model the water ice (where one more diffuse function is added to all of the atoms, surfaces and the adsorption of the 17 S-bearing species on including H) was adopted to recover the whole interaction them. The calculations were performed with the periodic between SO2 and the two water molecules of the model system. ab initio code CRYSTAL17 developed by Dovesi et al. (2018), See Appendix A for more details on the BSSE correction and which adopts localized Gaussian functions as basis sets. This on the ONIOM2 scheme. code exploits the local combination of atomic orbitals (LCAO) On the amorphous ice model, the DFT//DFT scheme can approximation, expanding the Bloch functions of periodic hardly be carried out due to the larger unit cell size compared systems as a linear combination of atom-centered Gaussian with the crystalline one. Therefore, we adopted the computa- functions. This software implements both the Hartree–Fock tionally cheaper semiempirical HF-3c method for the adsorp- (HF) and Kohn–Sham self-consistent field methods to solve the tion on the amorphous ice surface model. HF-3c is an HF-based electronic Schrödinger equation, taking advantage of the method adopting a minimal basis set (MINI-1; Tatewaki & symmetry of the system, when present. Thanks to its capability Huzinaga 1980), to which three empirical corrections (3c) are to simulate systems with periodicity ranging from zero added (Sure & Grimme 2013): (i) the dispersion energy 3 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. (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 per- formed 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 diagonaliz- ing 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 Å Figure 1. Top panel: side view of the crystalline (010) water P-ice slab model from the minimum along each Cartesian coordinate. For the (along the b lattice vector). Bottom panel: side view of the amorphous water ice crystalline ice, computed at the DFT level, only a portion of the slab model (along the b lattice vector). Numbers identify the different system constituted by the adsorbed species plus the two closest adsorption sites: those in blue represent dangling hydrogen atoms (dH), while those in red represent dangling oxygen atoms (dO). 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 avoid lateral interactions between the adsorbed molecules in computed at the HF-3c level. For this later case, however, the different replicas of the ice cell, thereby simulating the adsorption frequencies were only used to check that the complexes are of isolated species. The supercell shows two dangling hydrogen actual minima of the potential energy surfaces (all of the (dH) atoms and two dangling oxygen (dO) atoms as binding sites frequencies are real), since they are expected to not be accurate (top panel of Figure 1). The structure of the ice is such that it enough because of the approximated HF-3c methodology. In ensures a null electric dipole along the nonperiodic z-axis. This is contrast, from the DFT frequencies computed for the a direct consequence of the symmetry of the system, which fragments, we computed the zero-point vibrational energy shows two identical faces both at the top and at the bottom of the (ZPE) corrections to: (i) correct the electronic BEs for the model. Therefore, the adsorption was modeled only on the top crystalline ice by including the corresponding ZPE contribu- face of the system. tions, and (ii) derive a scaling factor to apply to the BEs The more realistic amorphous model, already used in Ferrero computed on the amorphous ice model to arrive at ZPE- et al. (2020; see bottom panel of Figure 1) consists of 60 water corrected values without the need to explicitly run a full DFT molecules per unit cell. The cell parameters are ∣ a ∣ = 20.355 Å, frequency calculation (see Appendix B for more details on ∣ b ∣ = 10.028 Å, and ∣ g ∣ = 102°. 983 (at HF-3c). The different these two points). 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 2.2. Ice Models variety of different binding sites. Therefore, to model the Interstellar ices are thought to be mostly formed by amorphous adsorption of the S-bearing species, we selected eight solid water (Watanabe & Kouchi 2008; Boogert et al. 2015). characteristic adsorption sites (highlighted in the bottom panel However, as a starting point of our study, we have chosen a of Figure 1). The sites are chosen as follows: cases 1, 2, and crystalline model (see the top panel of Figure 1). This choice is 5 are located on the top surface of the ice, case 4 and 6 are due to two main reasons: (i) crystalline structures are well defined inside the cavity of the model, and case 3, 7, and 8 are on the due to symmetry constraints and are computationally cheap, and bottom-surface of the model. Binding sites from 1–6 are dH, (ii) no definite structure for amorphous ice models is available. while 7 and 8 are dO. Moreover, regions rich in crystalline ices have been observed in protoplanetary disks and stellar outflows (Molinari et al. 1999; 3. Results Terada & Tokunaga 2012). Our ice surface model is derived from the bulk of the proton-ordered P-ice, which was cut along the In this work, the adsorption of 17 S-bearing species (shown (010) surface defining a 2D periodic slab model (Casassa et al. in Figure 2) has been simulated, of which eight are closed-shell 1997; Zamirri et al. 2018). The thickness of the ice (10.9 Å) was species and nine open-shell species (with unpaired electrons as enough to converge the corresponding surface energy. We indicated by the dots in Figure 2). For the open-shell species, adopted a 2 × 1 supercell slab model, consisting of 12 atomic the electronic ground state of each one was checked and found layers and cell parameters of ∣ a ∣ = 8.980 Å and ∣ b ∣ = 7.082 Å to be in agreement with those reported by Woon & Herbst (at B3LYP-D3(BJ)/A-VTZ* level), which are large enough to (2009) and Das et al. (2018). Therefore, NS•, HS•, HS •2 , and 4 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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: (CH3SH, H2CS, H2S, H2S2, HCS•, Figure 2. Set of 17 S-bearing species considered in this work. Dots represent HS•, and HS •2 ): these molecules contain groups that can act unpaired electrons. as H-bond donors and H-bond acceptors. It is worth noting that in the cases of H2CS and HCS•, the H-bond donor HCS• are doublets (with one unpaired electron) and S••, SO••, groups are weak. •• •• 2 , C2S , and C4S are triplets (with two unpaired electrons). S•• Molecules of Group II: (OCS, SO2, NS•, and SO••): these Each species was placed manually on the ice surface models molecules feature only H-donor acceptor groups because according to the molecule/surface electrostatic potential they contain atoms that are more electronegative than S. complementarity, and from these manually built systems, Molecules of Group III: (CS, C2S••, C3S, and C4S••): geometry optimizations were performed. molecules that can establish one H-bond interaction with their C-end atom and a surface dH. Molecules of Group IV: (S•• and S•• 2 ): species consisting of 3.1. BEs on Crystalline Ice only S atoms so that dispersive interactions are expected to be the main binding driving forces. 3.1.1. DFT//DFT BEs The molecule with the smallest dispersion contribution to its Table 1 reports the computed BE values as described in the BE is H2S (of Group I), which, like water, is capable of acting computational section. as both an H-bond donor and acceptor, hence involving two The role of the dispersion interactions to the BEs was water molecules of the ice surface in its adsorption. H2S, H2S2, assessed by considering the D3(BJ) contribution only to the H2CS, and CH3SH (of the same Group I) show this same BE. In other words, for each full B3LYP-D3(BJ) optimized structure, we split the B3LYP-D3(BJ) total energy in a pure adsorption feature. However, the presence of an S–S or a C–S electronic DFT term coming from the B3LYP and the pure bond enhances the contribution of dispersion forces, which can dispersive D3(BJ) one. This approach maximizes the weight of reach up to circa 70% of the BE, as in H2S2. The open-shell the dispersion because the attractive nature of the D3(BJ) term species HS• and HS •2 (also of Group I) have smaller BEs than brings the adsorbates as close as possible to the icy surface, their hydrogenated counterparts due to the absence of one which, in turn, boosts the D3(BJ) value. A different approach H atom. would be to optimize the geometries at B3LYP only and then For molecules of Group II, the sulfur atom interacts through evaluate the D3(BJ) term through a single-point energy dispersive forces with dO of the nearest surface water evaluation (referred to as the B3LYP-D3(BJ)//B3LYP molecules. SO2, thanks to its two oxygen atoms that act as scheme). The missing D3(BJ) term during the geometry H-bond acceptors, presents the highest BE of the entire set of optimization causes both an elongation of the average distance tested molecules (6880 K). OCS, which shows an almost null and important changes in the spatial arrangement of the dipole moment, presents the largest dispersion contribution to adsorbate compared to the full B3LYP-D3(BJ) cases. There- the BE of the entire set (>90%), hence being one of the species fore, the dispersion contribution for the B3LYP-D3(BJ)// less bound to the ice. B3LYP is expected to be smaller (underestimated) than the one Molecules of Group III indeed form an H-bond between from the full B3LYP-D3(BJ) optimization calculation. We their C-ends and a surface dH. This interaction, moreover, is checked both approaches (data not reported) for two sets of enhanced due to the dipole moment of the C chains and to the molecules, O-bearing (CO, H2O, CO2, CH3OH, and H2CO) carbene-like character of the C-end atom. Furthermore, these and S-bearing (CS, H2S, OCS, CH3SH, and H2CS) adsorbed on species can bind strongly to the surface thanks to the onset of the P-ice crystalline surface. We found the B3LYP structures dispersive forces that become stronger with the size of the involving CS and, especially OCS, to be rather different from molecule. Despite these premises, the BEs of these species are the corresponding B3LYP-D3(BJ) ones. The same was found similar (4000–4500 K), with the exception of C3S, which has a for the CO and CH3OH cases. Nonetheless, the percentage higher BE (6723 K at the most) because of its adsorption average difference in the dispersion between the two geometry. That is, once C3S adsorbs on the surface, it is approaches is a modest 11% in favor of the B3LYP-D3(BJ), arranged in rows with a distance of only 3.1 Å between the as expected. In this work, we were interested in establishing the S-end and the C-end of two neighboring molecules, therefore relative weight of dispersion (rwd) between S- and maximizing favorable head-to-tail interactions between dipoles. O-containing species, rather than its absolute value. For the This results in an energetic gain due to the lateral interactions 5 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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 BEs on the Crystalline Water Ice DFT//DFT DFT//HF-3c Species BE Disp BE No Disp Disp (%) BE Disp BE No Disp Disp (%) CS 4582 1600 2982(65%) 3861 1239 2622(68%) C 3S 6038/6723 1564/2129 4474(74%)/4594(68%) 5713/6146 421/1443 5292(93%)/4703(77%) CH3SH 6182 2959 3223(52%) 5340 2081 3259(61%) H2CS 4847/6242 2225/2850 2622(54%)/3392(54%) 4258/5641 1792/2417 2466(58%)/3223(57%) H2 S 5268 3332 1936(37%) 4679 2862 1816(39%) H2 S 2 5436/6158 1624/1936 3813(70%)/4222(69%) 5027 1527 3500(70%) OCS 3007/3440 -229/325 3235(107%)/3115(90%) 2790 132 2658 (95%) SO2 6880 3608 3271(48%) 5929 2502 3428(58%) HCS• 4113 3849 HS• 4270 3247 HS •2 4270/4763/4871 3572/3139 NS• 4330 3765 C2S•• 4089/4691 3993/4498 C4S•• 4366/4414 4246 S•• 3921 3247 S •• 2 2393 2009 SO•• 3861 1924 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 3.1.3. DFT//HF-3c BEs that enhance the BE, a phenomenon that does not happen for To cope with the larger size of the amorphous ice model, we C2S•• and C4S••. In fact, we would expect C4S•• to show larger tested the performance of the cheaper DFT//HF-3c scheme on BEs, but its elongated structure makes it challenging to the crystalline ice model to assess its applicability also for the maximize its contact area with the ice surface, hence interacting amorphous ice. mostly through the C-end and pointing the S-end upward. Starting from the previously optimized DFT structures, we Finally, for molecules of Group IV, even though S•• 2 is re-optimized the geometries at HF-3c level. In general, the HF- heavier and larger than S••, its interaction with the ice is the 3c optimized geometries are comparable with the DFT ones, weakest of the entire set of the tested molecules. S•• 2 has an with the exception of very few cases, such as the SO2 electronic triplet ground state (like O2) and is the only adsorption complex (see below). However, by comparing the homonuclear molecule of the set. Its dipole moment is zero, HF-3c//HF-3c BEs with the DFT//DFT ones, the resulting but not its quadrupole moment, which allows for the onset of a correlation is rather coarse (see Figure 3, panel (b), blue quadrupolar interaction between the S–S bond region and a dH symbols and line). of the ice surface. However, its BE is lower (2393 K) than that However, a much better agreement is shown by the DFT// for S•• (3921 K). 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 3.1.2. ONIOM2-corrected BEs accuracy of the DFT//HF-3c scheme was also exhibited in molecular crystals (Cutini et al. 2016), polypeptides (Cutini For some species, BEs computed at DFT//DFT have been et al. 2017), and pure-silica zeolites (Cutini et al. 2019). These refined at a CCSD(T) level by applying the ONIOM2 results, thus, extend the application of the DFT//HF-3c methodology by Dapprich et al. (1999), considering that most scheme to derive accurate BEs for S-bearing molecules of the BE contribution arises from the local interactions adsorbed on ice surfaces. between the adsorbed species and a small number of nearby HF-3c//HF-3c data (blue symbols) of Figure 3 reveal an water molecules. These calculations were performed on the outlier (for SO2, around 7000 K), the HF-3c//HF-3c BE being following seven S-bearing species: CH3SH, CS, H2CS, OCS, 25% higher than the DFT//DFT one, while the DFT//HF-3c SO2, NS•, and S•• 2 . These molecules have been chosen because value is about 15% smaller than the DFT//DFT value. they well represent both closed- and open-shell species and This effect is because HF-3c optimized geometries show sample the entire range of BEs of the 17 species of this work. shorter distances between the adsorbate and the ice model, Figure 3 (panel (a)), plots the ONIOM2-corrected BEs which, when computing the single-point energy at the DFT against the DFT//DFT ones, showing a very good fit that level, will increase both the dispersive interactions and the confirms the accuracy and reliability of the DFT//DFT scheme exchange repulsion contributions. For instance, HF-3c reduces adopted so far, with the latter slightly overestimated compared the distance between the S atom of SO2 and an O atom of the to the ONIOM2-corrected one. ice surface to 2.25 Å, compared to the 2.48 Å of the B3LYP-D3 6 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. which, when adsorbed, shrinks to 1.633 Å, resulting in a gain in energy of about 10 kJ mol−1 (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), Figure 3. (a) Best linear fit between the DFT//DFT and the ONIOM2- this is not the case for the ASW model. In ASW, the random corrected BEs for the crystalline ice systems. (b) Best linear fit between the organization of the water molecules breaks the H-bonding DFT//DFT BEs and the HF-3c//HF-3c (in blue) and the DFT//HF-3c (in chain extensions, reducing the H-bonding cooperativity and, orange) schemes, for the crystalline ice systems. BEs are in kelvin. therefore, its strength. This, in turn, also affects the H-bonding strength exhibited by the chain terminal ice OH groups (BJ) optimization. Moreover, by looking at the BE disp values involved in the interactions with the adsorbates, decreasing the corresponding BE values. This is shown clearly by the of Table 1, we notice that the percentage of dispersion in higher BE values computed for adsorption on the crystalline DFT//HF-3c BE values is almost always larger than that of P-ice with respect to the ASW model, in agreement with DFT//DFT BEs. previous literature data (Ferrero et al. 2020). However, other Finally, HF-3c can be critical for radical species like NS•. effects, like the easier structure deformation for the amorphous The NS• experimental bond length value is 1.494 Å (Peebles & ice, accommodates the adsorbates getting closer to the ice Marshall 2002), while the calculated ones are 1.480 Å with surface, maximizing both the H-bond strength and the M06-2X and 1.651 Å with HF-3c. The charge and spin density dispersive contribution, overcoming the loss of H-bonding distributions are also different, with 76% and 59% spin density cooperativity. on nitrogen at HF-3c and at M06-2X, respectively. The higher The above interpretation can explain the presence of some spin localization of HF-3c causes the lengthening of the bond, very high BE values, which depart from the rest of the BEs of 7 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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 Amorphous Ice BE Values Species dH1 dH2 dH3 dH4 dH5 dH6 dO7 dO8 〈BE〉 σ CS 4270 1852 3596 3283 3380 3584 1155 1648 2846 1054 C3S 5472 6014 6122 5412 5112 5941 1275 2959 4788 1630 CH3SH 3608 4991 3584 3271 3271 4630 2622 2850 3603 769 H2CS 4161 4967 3416 3488 4378 4053 1960 3115 3692 859 H2 S 3452 5208 3283 4210 4630 2622 2285 2309 3500 1026 H2 S 2 3007 4234 4558 4462 3801 3801 3175 4883 3990 623 OCS 2875 1491 2574 2177 1660 3320 2057 3175 2416 637 SO2 3344 3175 3644 6759 3680 4426 2730 2442 3775 1263 HCS• 2117 2225 2454 2117 1491 3560 2478 2201 2330 545 HS• 3163 2454 2057 2646 5003 2177 1624 1251 2547 1080 HS •2 2562 1527 2586 3031 2670 2550 2778 2694 2550 413 NS• 3283 1431 1828 1768 3416 3428 1768 1876 2350 805 C2S•• 4426 3163 5088 3572 2009 4991 2009 2393 3456 1193 C4S•• 5075 2165 5797 3716 4943 5569 4306 4077 4456 1097 S•• 3247 2081 2550 2321 2526 2321 1852 1852 2344 425 S •• 2 1924 1455 2850 1552 2009 2069 1648 1203 1839 472 SO•• 4041 1082 2634 1419 1696 2345 1672 1696 2073 874 Note. For each species, eight BE values, alongside their mean (〈BE〉) and standard deviation (σ), are provided. the corresponding ensembles. This happens for CH3SH, C3S, H2CS, H2S, H2S2, C2S••, and C4S••, 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); how- ever, 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 Figure 4. Analysis of ZPE corrections. Correlation between DFT//DFT BE(0) adopted the new ice structure as a reference for the calculation and BE for the adsorption at the crystalline ice. The intercept is set to zero. of the BE, which are now coherent with the general BE ensembles. We think this is an important issue, not common in good linear correlation between BE and BE(0), showing the the adsorption phenomena on metallic, ionic, or covalent average ZPE correction is about 14% of the BE. surfaces, which are dominated by stronger forces than those The final BE(0) values for all of the adsorption complexes at occurring with soft matter like water ices. This problem would the amorphous ice have, therefore, been computed by multi- be mitigated only by enlarging the ice model at an increasing plying the DFT//HF-3c BEs by a scale factor of 0.862 (see computational cost. Figure 4, panel (a)). 3.3. Zero-point Energy (ZPE) Corrections 4. Discussion As anticipated in Section 2.1, the ZPE corrections were 4.1. Comparison with BEs for Non-S-bearing Adsorbates computed following two different schemes, depending on the adopted surface models. For the adsorption on the crystalline In the recent work by Ferrero et al. (2020), the adsorption of ice, we computed the harmonic frequencies of a fragment 17 closed-shell and four open-shell species was simulated. The consisting of the S-bearing species and the two closest water list included: (i) nondipolar molecules with significant quadru- molecules of the ice with the assumption that deformation pole moments (i.e., H2, N2, and O2) or higher-order multiple occurs only in the proximity of the adsorption site. The ΔZPE moments (i.e., CH4), hence interacting weakly with the ice; (ii) is computed following the advice described in Appendix B to molecules that, besides exhibiting a quadrupole moment, arrive at a ZPE-corrected BE, i.e., BE(0). Figure 4 shows a established H-bonds with surface dH atoms (i.e., CO, OCS, and 8 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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 ensem- bles, one can argue S-containing species show higher desorption temperatures compared to weakly interacting molecules, like O2, CH4, 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•• 2 (the weakest bounded S-bearing species) also implies its highest mobility. In the hypothesis that polysul- phanes could be a reservoir of sulfur, the high mobility of S•• 2 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 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) In addition to these general considerations, a punctual Distribution of the dispersive (left panel) and nondispersive (right panel) comparison between S-bearing and O-bearing species can be contributions to the crystalline BEs for closed-shell species. drawn for some cases in order to highlight a few exceptions to what we generally understood about our set of molecules. CO2); (iii) amphiprotic molecules, which are both H-bond In Table 3 we compare the contributions to the BE(0) values acceptors and donors (i.e., NH3, H2O, HCl, HCN, and H2S); of some S-bearing species with their oxygen analogs. As a rule (iv) larger molecules that are capable of establishing a greater of thumb, we would expect molecules containing sulfur to number of interactions with the surface due to the presence of always show a BE(0) lower than the corresponding O-bearing several functional groups, hence giving rise to the highest BE ones, with the dispersive interactions being larger for the values of the set (i.e., CH3OH, CH3CN, H2CO, HCONH2, and former group. However, by considering case by case, we notice HCOOH); and (v) the open-shell species OH•, NH •2 , CH •3, and that there are some exceptions to this behavior (see Figures 7–9 HCO•, the two former ones forming strong H-bonds with dH in Appendix C) For the pair CS/CO, the trend is the opposite, and dO sites, at variance with the two latter ones. while the OCS/CO2 pair shows many similarities. Figure 5 compares the BEs of the S-bearing species with Indeed, CS exhibits a larger BE(0) and a lower percentage of those of Ferrero et al. (2020). The distribution of BEs in Ferrero dispersive energy than its O-analog, CO. This can be explained et al. (2020) has a bimodal character (see panel (a)), indicating by looking at the dipole moments of these species. CS shows a two different adsorption behaviors in the set of the considered large dipole moment (1.9 D), at variance with the negligible species: a group of low BEs (characterized by dispersive one of CO. The H-bonding interaction of the two species with a 9 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. Table 3 Contribution to Selected BE(0) (in kelvin) for S-bearing Species vs. O-bearing Species Adsorption at Crystalline Ice Sulfur Bearing Species Oxygen Bearing Species Species BE(0) No Disp Disp (%) μ Species BE(0) No Disp Disp(%) μ CS 2453 785 1668(68%) 1.9 CO 1663 50 1613(97%) 0.1 CH3SH 4603 1795 2808(61%) 1.5 CH3OH 7385 5391 1994(27%) 1.7 H2CS 4267 1792 2475(58%) 1.7 H2CO 5187 3268 1919(37%) 4.6 H2 S 4033 2460 1573(39%) 1.1 H2 O 7200 5832 1368(19%) 2.0 OCS 2405 120 2285(95%) 0.8 CO2 2568 796 1772(69%) 0.0 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 giving rise to an ensemble of BEs. As reported in detail by interaction of the S atom with the water molecule of the surface Ferrero et al. (2020) and not repeated here, extra care should be (see Figure 7 in Appendix C). exerted when comparing the computed BEs with those from Both OCS and CO2 have a linear geometry, an almost null experiments. dipole, and interact through an H-bond with a dangling OH of For instance, while BEs are straightforwardly calculated by the surface. While OCS exhibits a larger dispersion interaction, computer simulations (see Appendix B), experimental BEs are for CO2, the electrostatic component makes up about 30% of never directly observed in the thermal desorption, as TPD the BE, yielding a comparable BE(0) for the two molecules. peaks are usually interpreted through the Redhead (1962) On the other hand, the other pairs behave as we expected, as method or more sophisticated techniques, through which described in Section 4.1. The adsorption geometries show some desorption adsorption energies (DAEs) are worked out. As similarities between the two members of each pair, with the DAEs may be a function of the surface coverage θ, the most important differences due to the presence and strength of comparison with the theoretical BEs is not straightforward (for the H-bonds. Specifically, CH3SH/CH3OH and H2S/H2O a recent review, see Minissale et al. 2022). For instance, the ice exhibit the same number of H-bonds, while for H2CS/H2CO, structure is sensitive to both temperature and adsorption/ only the latter establishes an H-bond with the surface through desorption processes (at variance with surfaces from covalent/ the C=O group. ionic solids) as its layers are held together by weaker interactions (e.g., H-bonds in water ices or dispersion in CO 4.3. Comparison with the Literature ices), which are similar in strength to the BE values. Obviously, ice restructuring processes may affect the final value of Before comparing the present results with literature exper- the DAEs. imental data, some general considerations are at hand. In the In Table 4, our BEs are compared with those from Penteado ISM, the water ice mantle is the result of the in situ water et al. (2017; see Figure 6). These authors showed a list of formation through surface reactions. Therefore, the final ice recommended BEs, collecting data from previous works and structure is the result of the energetic locally hot spots injected providing an uncertainty range for each value. For H2S, OCS, in the mantle where intermediates (e.g., OH• and HOO•) and and SO2, the BEs are derived a posteriori from the experimental final products are formed. As the injected energy is large, it TPD measurements of Collings et al. (2004), using the locally heats up the ice, shaping the final structure of the grain empirical relationship of in a complex way, difficult to attain through laboratory experiments. In the terrestrial laboratories, ice is prepared BE (X ) : Tdes (X ) = BE (H2 O) : Tdes (H2 O) (1 ) directly by deposition from vapor water molecules, with a certain kinetics, on cold metal fingers (Penteado et al. 2017). where X is the considered chemical species, and Tdes is the The final ice can be amorphous (either compact or porous) or desorption temperature of the species deposited on a water ice crystalline, depending on the kinetics of adsorption and extrapolated from the TPD measurements of Collings et al. processing temperature. Unfortunately, the detailed structure (2004). The reference BE(H2O) is set to 4800 K from previous of amorphous ice can only be partly inferred indirectly by works (see Penteado et al. 2017 for more details). The spectroscopic measurements, therefore limiting our knowledge uncertainty assigned to the BE values of H2S, OCS, and SO2 of the atomistic details. The same difficulty occurs for was 3.5%. For the other species, the provided data are based on computer modeling: indeed, we are not aware of computer the work of Hasegawa & Herbst (1993), which in turn were simulations reacting hydrogen and oxygen atoms step by step based on previous works, including Allen & Robinson (1977). up to form a sizable icy grain, even if a simpler attempt to form In the latter, the BE of a molecule was treated as the sum of the a single water molecule on the forsterite surface has been BEs of the constituent atoms to infer the interaction energy of recently reported (Molpeceres et al. 2019). As a function of the adopted ice models in computer several species on SiO2 surfaces, which was subsequently simulation, the external surface may show a variety of different scaled by an empirical factor to obtain BEs on water ice. The kinds of adsorption sites, giving rise to different BEs (Song & uncertainty of these values is set to half the BE for species Kastner 2017). When very few sites are available, like in the whose BE is less than 1000 K, and set to 500 K for all other proton-ordered ice model adopted here, no more than two cases. different BE values resulted. For the amorphous model, The comparison between our computed BEs for the crystal- instead, a significant variability of binding sites resulted, line ice and the literature data shows, as expected, 10 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. Table 4 Summary of the Computed BE(0) Values (in kelvin) in Comparison with Data from the Literature Binding Energies This Work Literature Crystalline Amorphous Computed Database Various Species BE(0) Min Max 〈BE〉 σ Das Wakelam KIDA UMIST Penteado CS 2453 995 3680 2453 909 2217 3200 3200 1900 1800 ± 500 C 3S 4925/5298 1099 5277 4128 1405 3500 3500 3000 ± 500 CH3SH 4603 2260 4302 3106 663 4200 H2CS 3670/4863 1690 4282 3183 740 3110 4400 4400 2700 2025 ± 500 H2 S 4033 1970 4489 3017 884 2556 2900 2700 2743 2290 ± 90 H2 S 2 4333 2592 4209 3439 537 4368 3100 3100 2600 ± 500 OCS 2405 1286 2861 2083 549 1571 2100 2400 2888 2325 ± 95 SO2 5111 2105 5827 3254 1089 3745 5000 3400 5330 3010 ± 110 HCS• 3318 1286 3069 2009 469 2713 2900 2900 2350 2000 ± 500 HS• 2799 1078 4313 2195 931 2221 2700 2700 1500 1350 ± 500 HS •2 2706/3079 1317 2613 2198 356 4014 2650 2650 2300 ± 500 NS• 3245 1234 2955 2026 694 2774 1900 1900 1800 ± 500 C2S•• 3442/3877 1731 4385 2979 1028 2447 2700 2500 ± 500 C4S•• 3660 1866 4997 3841 945 4300 4300 3500 ± 500 S•• 2799 1579 2799 2020 366 1428 2600 2600 1100 985 ± 495 S •• 2 1732 1037 2457 1644 1585 407 2200 2200 2000 ± 500 SO•• 1658 933 3483 2128 1787 754 2900 2800 2600 1800 ± 500 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). 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. 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 Das et al. (2018) computed the BEs for a molecular set model, while blue stars correspond to the ZPE-corrected BEs on the crystalline including 100 species, adopting a water tetramer and, in few ice model. cases, a water hexamer as the ice model, at the MP2/aug-cc- pVDZ theory level, without BSSE and ZPE corrections. discrepancies. However, comparison of the range provided in Among the S-bearing molecules, only H2S and OCS were Penteado et al. (2017) with our BE ensemble reveals that adsorbed on a water hexamer. The ice models employed in Das disagreement is only for S••. The additive principle applied by et al. (2018) still lack the effect of H-bond cooperativity Allen & Robinson (1977) implies that S••2 BE must be twice the (present in our periodic models) and could be the reason of the value of the isolated S•• atom. However, we have seen that S•• small values of BEs for the S-containing species, even lower and S•• 2 showed a very similar distribution of BEs, due to the than those estimated by Wakelam et al. (2017). Moreover, H2S same nature of their interaction with the surface. For the other and OCS BEs show an increment when changing the model species, the ranges provided by Penteado et al. (2017) lay from the water tetramer to the water hexamer, and it is not clear within our ensembles, especially for H2S, OCS, and SO2, for how long this trend would continue before reaching a which have a very narrow distribution of BEs. However, converged value. The Das et al. (2018) values are within the computed BEs for crystalline ice compare poorly with the BE ensembles computed for the amorphous ice. Altogether, 11 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. this indicates that the choices of the ice model and of the the two estimated BEs: 995–3680 K (CS; this work) versus computational method have a strong relevance on the predicted 1109–1869 K (CO; Ferrero et al. 2020). On the other hand, the BEs. Thus, more studies are needed to fully characterize BE observations suggest similar snow lines for CS, H2CO, and distributions, by using larger models and exploring more H2CS, while the estimated BEs of H2CO and H2CS are larger binding sites. than that of CS: 3071–6194 K (H2CO; Ferrero et al. 2020) and 1690–4282 K (H2CS; this work) versus 995–3680 K (CS; this work). We emphasize, however, that the ALMA-DOT 4.4. Binding Energies and Snow Lines in Protoplanetary Disks observations are unable to resolve the inner 25 au of the disk. As mentioned in Section 1, the BEs of the species on ice Therefore, higher spatial resolution observations are needed to mantles have a profound impact on the interstellar chemistry constrain the snow lines of species with large BEs, such as and, more specifically, on the chemical composition of the H2CO and H2CS. regions where planetary systems, like our own, are formed. A In conclusion, the admittedly scarce observations so far very important case is represented by the gaseous versus solid available seem to be in agreement with the estimated BEs chemical composition in protoplanetary disks, which are the obtained through theoretical quantum chemical calculations as sites where planets, asteroids, and comets eventually form. If a provided in Ferrero et al. (2020) and this work, which species is in the solid form, then it will likely be incorporated in encourage us to pursue the work to a larger sample of species, the (aforementioned) formation of rocky objects, whereas, if it both S- and N-bearing and the so-called interstellar complex is in the gaseous form, it will only enrich the giant gaseous organic molecules (Ceccarelli et al. 2017). forming planets of the elements contained in the species (e.g., Öberg & Bergin 2021). Sulfur, on the other hand, seems to be a 5. Conclusions very important element for terrestrial life and, perhaps, even a key element for its emergence (e.g., Shalayel et al. 2020). In the present work, the BEs of 17 astrochemically relevant As discussed in Section 1, the major reservoirs of sulfur in S-bearing species adsorbed on periodic water ice models have cold molecular clouds (Section 1.2, from where the evolution been computed using DFT based on the hybrid B3LYP-D3(BJ) toward a planetary system starts) are unknown, though they are (for closed-shell species) and hybrid meta-GGA M06-2X (for suspected to be in solid organosulfur compounds and open-shell species) with an Ahlrichs-VTZ* double-zeta polar- polysulphates (Druard & Wakelam 2012; Laas & Caselli 2019). ized Gaussian basis sets. We adopted either a crystalline During the evolution, the quantity of sulfur in the gaseous proton-ordered ice model or a more realistic amorphous ice phase changes by orders of magnitude, so that it is very model, both treated within the periodic boundary conditions. important to know it in the protoplanetary disks. So far, there The adopted DFT results have also been refined toward a are no measurements for solid, and only very few for gaseous, CCSD(T) level through the ONIOM2 approach carried out on S-bearing species in protoplanetary disks (see below). The the crystalline model, showing a very good convergence locus where the transition of a species from the solid to the gas between the DFT BEs and those refined with the ONIOM2 phase takes place is usually called “snow-line” in astrophysical scheme. As full DFT is too costly for the simulation of jargon, particularly so in the protoplanetary disks, which are adsorption on the amorphous ice, we checked the reliability of objects that are relatively cold and dense (e.g., Öberg & the cheaper HF-3c method by contrasting the BEs with that at Bergin 2021). Then, one efficient way to measure the solid the DFT level for the crystalline ice. While HF-3c BEs poorly versus gaseous content of S-bearing species in protoplanetary correlate with the DFT ones, the BEs at the single-point DFT// disks is to observe where the species disappear from the gas HF-3c level give an excellent correlation with full DFT values. phase, because they freeze onto the dust grain water-rich The validated DFT//HF-3c method was used to study the mantles, namely, to measure their snow lines. Two major adsorption process on the amorphous ice model, in which the parameters determine where the snow line lies: the temperature richness of adsorption sites (up to eight) implies a BE ensemble of the dust, and the BE of the species. In general, the larger the for each species. We also showed the ZPE correction to the distance from the disk center (heated by the central forming electronic BE to be around 14% of the BE itself irrespective of star), the colder the dust. Therefore, species with lower BEs the adsorbed species, and provided a scale factor to obtain the remain gaseous at larger disk radii. BE(0) for adsorption complexes at the amorphous ice model. The ALMA-DOT survey of young disks provides the first This is important, as it allows for skipping the expensive preliminary observational constraints on the position of the frequency calculation by simple scaling of the electronic BEs. snow lines of some of the S- and O-bearing molecules The analysis of the obtained data showed that, in order to considered in this study (Podio et al. 2020a). In particular, properly describe the behavior of S-bearing species, special ALMA-DOT observed a protoplanetary disk, IRAS 04302 care must be taken in the description of the dispersive +2247, seen edge-on, which makes it ideal to constrain the interactions, which gain more importance as the number of snow lines of different species. These observations show that large and heavy atoms, like sulfur, increases. When considering the molecular emission is bright in an intermediate disk layer the amorphous ice model, which provides enough richness of (the so-called molecular layer; Aikawa & Herbst 1999), while it sites, the resulting BE ensembles match both the experimental decreases in the disk midplane due to molecular freeze-out onto and the literature computed BEs. We posit BE distributions the dust grains. More relevant to this work, CS, H2CO, and rather than single BE values to be essential for providing robust H2CS show a similar spatial distribution, and a sharp decrease parameters to be adopted in numerical models of chemical of their emission is observed at smaller radii than CO (Codella evolution in the universe. et al. 2020; Podio et al. 2020b). The estimated snow line is Regarding the sulfur depletion problem, our calculations do at 25 au for H2CO (and consequently H2CS and CS) and 100 au not solve the open questions. Nonetheless, we provide strong for CO (van’t Hoff et al. 2020). The larger snow-line radius of and relatively accurate BE values, which can elucidate the CS with respect to that of CO is generally in agreement with features of S-bearing species in comparison with S-free 12 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. molecules when adsorbed on icy mantles. An important clue is The Broyden–Fletcher–Goldfarb–Shanno algorithm (Broyden provided by the S•• 2 case in which the lowest BE implies higher 1970; Fletcher 1970; Goldfarb 1970; Shanno 1970) was used to diffusivity and increasing reactivity toward the formation of carry out geometric optimizations. higher-weight S-containing rings strongly bound to the ice through dispersion interactions. The CS molecule also shows Appendix B some counterintuitive results when compared to the analog CO Calculation of the Binding Energies molecule, due to a higher dipole moment and dispersion When using a finite basis set of localized Gaussian functions contribution. The present work also highlights the utility of to describe our systems, BSSEs arise. This is the case in our accurately computing BE values to be correlated with snow BEs computed at both DFT//DFT and DFT//HF-3c. Thus, lines in protoplanetary disks. the a posteriori counterpoise (CP) correction of Boys & Bernardi (1970) was applied to compensate for this error, and This project has received funding within the European our BEs are therefore defined as the opposite of the CP- Union’s Horizon 2020 research and innovation program from corrected interaction energies (ΔECP), that is, the European Research Council (ERC) for the projects “Quantum Chemistry on Interstellar Grains” (Quantumgrain), BE = -DE CP (B1) grant agreement No. 865657 and “The Dawn of Organic DE CP = DE - BSSE (B2) Chemistry” (DOC), grant agreement No. 741002, and from the Marie Sklodowska-Curie for the project “Astro-Chemical where the non-CP-corrected interaction energy ΔE is given by * Origins” (ACO), grant agreement No. 811312. The Italian the sum of the deformation-free interaction energy (ΔE ), the Space Agency for co-funding the Life in Space Project (ASI N. deformation energy of the slab (δES), the deformation energy of 2019-3-U.O), the Italian MUR (PRIN 2020, Astrochemistry the molecule (δEμ), and the lateral interaction (EL) between beyond the second period elements, Prot. 2020AFB3FX) are adsorbate molecules in different replicas of the cell. This also acknowledged for financial support. A.R. is indebted to the quantity corresponds to the common definition of interaction “Ramón y Cajal” program. The authors thank Gretobape for energy, which in this case is the energy of the complex to fruitful and stimulating discussions. The authors wish to thank which energies of the isolated species and isolated ice model the anonymous reviewers for the valuable suggestions. are subtracted. Supplementary Material consisting of (i) the fractional coordinates of HF-3c adsorption structure optimized for both DE = DE * + dES + dEm + EL the crystalline and amorphous ice, (ii) images of the adsorption = Ecomplex - Eice - Especies. (B3) features at crystalline periodic ice models, in which electro- static potential maps, spin density maps (for open-shell For the crystalline systems, inclusion of the ZPE corrections systems), and vibrational features are displayed, (iii) a pdf file to the BE are performed according to the equation BE(0) = BE with a thorough guide to the computation of BEs and the basis − ΔZPE, where sets employed for the calculations, is available on DZPE = ZPE complex - ZPE ice - ZPE species. (B4) Zenodo: https://doi.org/10.5281/zenodo.6798922. For the amorphous case, inclusion of ZPE corrections was done by applying a scaling factor derived from a linear Appendix A correlation between BE and BE(0) in the crystalline systems, Computational Details that is, BE(0) = 0.862*BE. Finally, to check for accuracy and refine DFT BEs for our For optimization and frequency calculations, we refer to the crystalline ice model, the single- and double-electronic computational parameters adopted in the work of Ferrero et al. excitation coupled-cluster method with an added perturbative (2020): the threshold parameters for the evaluation of the description of triple excitations (CCSD(T)) was used, in Coulomb and exchange bielectronic integrals (TOLINTEG) combination with a correlation consistent basis set extrapolated were set to 7, 7, 7, 7, and 14 when working with crystalline ice, to the largest possible basis set. These calculations have been while for the amorphous ice to 7, 7, 7, 7, and 25, which performed with the Gaussian16 software package (Frisch et al. facilitates the convergence of the self-consistent field (SCF), 2016). These refinements were performed through the especially when adsorbing open-shell species. ONIOM2 approach (Dapprich et al. 1999), dividing the Some species required a careful treatment to describe systems into two parts (model and real systems), which are properly their electronic state. The use of the BROYDEN described by two different levels of theory (high and low). The accelerator (Johnson 1988) and the broken-(spin)-symmetry model system (represented by the adsorbate and the two closest ansatz (Neese 2004, 2009) was necessary to obtain an optimal water molecules) was described by the high level of theory, convergence of the SCF procedure. CCSD(T). The real system (that is, the whole system) was For the numerical evaluation of bielectronic integrals described by the corresponding DFT (low) level of theory. (Becke 1988), the standard pruned grid, composed of 75 radial The ONIOM2-corrected BE is: points and a maximum of 974 angular points, was used. To BE(ONIOM2) = BE (low, real) sample the reciprocal space, we adopted the Pack–Monkhorst + BE (high, model) - BE(low, model). mesh (Pack & Monkhorst 1977), with a shrinking factor of 2, which generates 4 k points in the first Brillouin zone. The only (B5) exception is for the adsorption of atomic sulfur on the crystalline In this way, the ΔE = BE(high, model) − BE(low, model) ice model, which required a shrinking factor of 3, generating 5 k represents the correction term to the energy of the real system points, in order to avoid numerical noise in the results. due to the improved description at the high level. This 13 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. procedure ensures that if the low level of theory is improved to where BE(DFT; μ − ice) is the BE computed at DFT//DFT, arrive at the high level of theory or if the model system is while the BEs of the model system (μ − 2H2O) are enlarged to become the real system, the energy of the whole computed through single-point energy calculations at system is that at the higher level of theory. CCSD(T). In this work, for the calculation of the ONIOM2-corrected A thorough guide to the computation of BEs is available in BEs, BE(ONIOM2), Equation (B5) can be rewritten as the Supplementary Material. BE(ONIOM2) = BE (DFT; m - ice) + BE (CCSD (T ) ; m - 2H2 O) Appendix C - BE (DFT; m - 2H2 O) (B6) S-bearing versus O-bearing Species 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) CH3SH vs. CH3OH; (ii) CS vs. CO. 14 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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 H2CS vs. H2CO. 15 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. 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) H2S vs. H2O; (ii) adsorption geometry of OCS vs. CO2. ORCID iDs Becke, A. D. 1993, JChPh, 98, 5648 Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015, ARA&A, 53, 541 Jessica Perrero https://orcid.org/0000-0003-2161-9120 Bovolenta, G., Bovino, S., Vöhringer-Martinez, E., et al. 2020, MolAs, 21, Joan Enrique-Romero https://orcid.org/0000-0002- 100095 2147-7735 Boys, S., & Bernardi, F. 1970, MolPh, 19, 553 Brandenburg, J. G., Alessio, M., Civalleri, B., et al. 2013, JPCA, 117, 9282 Stefano Ferrero https://orcid.org/0000-0001-7819-7657 Broyden, C. G. 1970, JApMa, 6, 76 Cecilia Ceccarelli https://orcid.org/0000-0001-9664-6292 Casassa, S., Ugliengo, P., & Pisani, C. 1997, JChPh, 106, 8030 Linda Podio https://orcid.org/0000-0003-2733-5372 Ceccarelli, C., Caselli, P., Fontani, F., et al. 2017, ApJ, 850, 176 Claudio Codella https://orcid.org/0000-0003-1514-3074 Ceccarelli, C., & Dominik, C. 2005, A&A, 440, 583 Albert Rimola https://orcid.org/0000-0002-9637-4554 Cernicharo, J., Cabezas, C., Agúndez, M., et al. 2021, A&A, 648, L3 Piero Ugliengo https://orcid.org/0000-0001-8886-9832 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, References 354, 1133 Cooke, I. R., Öberg, K. I., Fayolle, E. C., Peeler, Z., & Bergner, J. B. 2018, Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233 ApJ, 852, 75 Al-Halabi, A., & Van Dishoeck, E. F. 2007, MNRAS, 382, 1648 Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, SSRv, 212, 1 Allen, M., & Robinson, G. W. 1977, ApJ, 212, 396 Cutini, M., Civalleri, B., Corno, M., et al. 2016, J. Chem. Theory Comput., Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197 12, 3340 Becke, A. D. 1988, PhRvA, 38, 3098 Cutini, M., Civalleri, B., & Ugliengo, P. 2019, ACS Omega, 4, 1838 16 The Astrophysical Journal, 938:158 (17pp), 2022 October 20 Perrero et al. Cutini, M., Corno, M., & Ugliengo, P. 2017, J. Chem. Theory Comput., Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1 13, 370 Pack, J. D., & Monkhorst, H. J. 1977, PhRvB, 16, 1748 Dapprich, S., Komáromi, I., Byun, K. S., Morokuma, K., & Frisch, M. J. 1999, Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, ApJ, Comput. Theor. Chem., 461, 1 897, 56 Das, A., Sil, M., Gorai, P., Chakrabarti, S. K., & Loison, J. C. 2018, ApJS, Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49 237, 9 Pascale, F., Zicovich-Wilson, C. M., Lopez Gejo, F., et al. 2004, JCoCh, Dohnálek, Z., Kimmel, G. A., Joyce, S. A., et al. 2001, JPCB, 105, 3747 25, 888 Douglas, A. E., & Herzberg, G. 1942, CJRes, 20a, 71 Peebles, L. R., & Marshall, P. 2002, CPL, 366, 520 Dovesi, R., Erba, A., Orlando, R., et al. 2018, WIREs Comput. Mol. Sci., 8, Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71 e1360 Phuong, N. T., Chapillon, E., Majumdar, L., et al. 2018, A&A, 616, L5 Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 Podio, L., Garufi, A., Codella, C., et al. 2020a, A&A, 644, A119 Druard, C., & Wakelam, V. 2012, MNRAS, 426, 354 Podio, L., Garufi, A., Codella, C., et al. 2020b, A&A, 642, L7 Dulieu, F., Congiu, E., Noble, J., et al. 2013, NatSR, 3, 1338 Pople, J. A., Gill, P. M. W., & Handy, N. C. 1995, IJQC, 56, 303 Dunning, T. H. 1989, JChPh, 90, 1007 Rae, J. G. L., Hartquist, T. W., Lepp, S. H., O’Neill, P. T., & Williams, D. A. Enrique-Romero, J., Ceccarelli, C., Rimola, A., et al. 2021, A&A, 655, A9 2004, A&A, 413, 1 Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJL, 816, L28 Redhead, P. 1962, Vacuu, 12, 203 Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11 Rimola, A., Ceccarelli, C., Balucani, N., & Ugliengo, P. 2021a, FrASS, 8, Fletcher, R. 1970, CompJ, 13, 317 655405 Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian 16 Revision Rimola, A., Ferrero, S., Germain, A., Corno, M., & Ugliengo, P. 2021b, Mine, C.01 (Wallingford, CT: Gaussian Inc.), https://gaussian.com 11, 26 Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 Rodríguez-Almeida, L. F., Jiménez-Serra, I., Rivilla, V. M., et al. 2021, ApJL, Germain, A., Tinacci, L., Pantaleone, S., Ceccarelli, C., & Ugliengo, P. 2022, 912, L11 ESC, 6, 1286 Ruf, A., Bouquet, A., Boduch, P., et al. 2019, ApJL, 885, L40 Ghesquiere, P., Mineva, T., Talbi, D., et al. 2015, PCCP, 17, 11455 Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, Goldfarb, D. 1970, MaCom, 24, 23 306, 691 Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, JChPh, 132, 154104 Schäfer, A., Horn, H., & Ahlrichs, R. 1992, JChPh, 97, 2571 Grimme, S., Ehrlich, S., & Goerigk, L. 2011, JCoCh, 32, 1456 Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 Shalayel, I., Youssef-Saliba, S., Vazart, F., et al. 2020, EurJOC, 2020, 3019 He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 825, 89 Shanno, D. F. 1970, MaCom, 24, 647 He, J., Emtiaz, S. M., & Vidali, G. 2017, ApJ, 837, 65 Shimonishi, T., Nakatani, N., Furuya, K., & Hama, T. 2018, ApJ, 855, 27 He, J., Emtiaz, S. M., & Vidali, G. 2018, ApJ, 863, 156 Smith, R. S., May, R. A., & Kay, B. D. 2016, JPCB, 120, 1979 Hily-Blant, P., Pineau des Forêts, G., Faure, A., & Lique, F. 2022, A&A, Song, L., & Kästner, J. 2017, ApJ, 850, 118 658, A168 Sure, R., & Grimme, S. 2013, JCoCh, 34, 1672 Jenkins, E. B. 2009, ApJ, 700, 1299 Swings, P., & Rosenfeld, L. 1937, ApJ, 86, 483 Jiménez-Escobar, A., & Muñoz Caro, G. M. 2011, A&A, 536, A91 Tatewaki, H., & Huzinaga, S. 1980, JCoCh, 1, 205 Johnson, D. D. 1988, PhRvB, 38, 12807 Terada, H., & Tokunaga, A. T. 2012, ApJ, 753, 19 Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114 Tieftrunk, A., Pineau des Forets, G., Schilke, P., & Walmsley, C. M. 1994, Karssemeijer, L. J., & Cuppen, H. M. 2014, A&A, 569, A107 A&A, 289, 579 Keller, L. P., Hony, S., Bradley, J. P., et al. 2002, Natur, 417, 148 Tinacci, L., Germain, A., Pantaleone, S., et al. 2022, ESC, 6, 1514 Kouchi, A., Furuya, K., Hama, T., et al. 2020, ApJL, 891, L22 van’t Hoff, M. L. R., van Dishoeck, E. F., Jørgensen, J. K., & Calcutt, H. 2020, Kruse, H., & Grimme, S. 2012, JChPh, 136, 154101 A&A, 633, A7 Laas, J. C., & Caselli, P. 2019, A&A, 624, A108 Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514 Lauck, T., Karssemeijer, L., Shulenberger, K., et al. 2015, ApJ, 801, 118 Wakelam, V., Caselli, P., Ceccarelli, C., Herbst, E., & Castets, A. 2004, A&A, Lee, C., Yang, W., & Parr, R. G. 1988, PhRvB, 37, 785 422, 159 Mason, N. J., Nair, B., Jheeta, S., & Szymanska, E. 2014, FaDi, 168, 235 Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, ApJS, 217, 20 Maté, B., Cazaux, S., Satorre, M. Á., et al. 2020, A&A, 643, A163 Wakelam, V., Loison, J.-C., Mereau, R., & Ruaud, M. 2017, MolAs, 6, 22 McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, McGuire, B. A. 2022, ApJS, 259, 30 418, 1035 McKellar, A. 1940, PASP, 52, 187 Watanabe, N., & Kouchi, A. 2008, PrSS, 83, 439 Mifsud, D. V., Kaňuchová, Z., Herczku, P., et al. 2021, SSRv, 217, 14 Woods, P. M., Occhiogrosso, A., Viti, S., et al. 2015, MNRAS, 450, 1256 Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ESC, 6, 597 Woon, D. E., & Herbst, E. 2009, ApJS, 185, 273 Mispelaer, F., Theulé, P., Aouididi, H., et al. 2013, A&A, 555, A13 Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427 Molinari, S., Ceccarelli, C., White, G. J., et al. 1999, ApJL, 521, L71 Zamirri, L., Corno, M., Rimola, A., & Ugliengo, P. 2017, ESC, 1, 384 Molpeceres, G., Rimola, A., Ceccarelli, C., et al. 2019, MNRAS, 482, 5389 Zamirri, L., Pantaleone, S., & Ugliengo, P. 2019a, JChPh, 150, 064702 Neese, F. 2004, JPCS, 65, 781 Zamirri, L., Ugliengo, P., Ceccarelli, C., & Rimola, A. 2019b, ESC, 3, 1499 Neese, F. 2009, Coord. Chem. Rev., 253, 526 Zhao, Y., & Truhlar, D. G. 2008, Theor. Chem. Acc., 120, 215 Noble, J. A., Theule, P., Mispelaer, F., et al. 2012, A&A, 543, A5 Zicovich-Wilson, C. M., Pascale, F., Roetti, C., et al. 2004, JCoCh, 25, 1873 17