Theoretical Modeling of the Adsorption of Neutral and Charged Sulphur-Bearing Species onto Olivine Nanoclusters

Sulphur depletion in the interstellar medium (ISM) is a long-standing issue, as only 1% of its cosmic abundance is detected in dense molecular clouds (MCs), while it does not appear to be depleted in other environments. In addition to gas phase species, MCs also contain interstellar dust grains, which are irregular, micron-sized, solid aggregates of carbonaceous materials and/or silicates. Grains provide a surface where species can meet, accrete, and react. Although freeze-out of sulphur onto dust grains could explain its depletion, only OCS and, tentatively, SO 2 were observed on their surfaces. Therefore, it is our aim to investigate the interaction between sulphur-containing species and the exposed mineral core of the grains at a stage prior to when sulphur depletion is observed. Here, the grain core is represented by olivine nanoclusters, one of the most abundant minerals in the ISM, with composition Mg 4 Si 2 O 8 and Mg 3 FeSi 2 O 8 . We performed a series of quantum mechanical calculations to characterize the adsorption of 9 S-bearing species, both neutral and charged, onto the nanoclusters. Our calculations reveal that the Fe S interaction is preferred to Mg S, causing sometimes the chemisorption of the adsorbate. These species are more strongly adsorbed on the bare dust grain silicate cores than on water ice mantles, and hence therefore likely sticking on the surface of grains forming part of the grain core. This demonstrates that the interaction of bare grains with sulphur species in cloud envelopes can determine the S-depletion observed in dense molecular clouds.


INTRODUCTION
On Earth, sulphur is one of the key elements that contributes to the engine of life.It is an essential component of our bodies, as it is found in amino acids such as cysteine and methionine, in coenzyme A, in vitamin B1 and biotin, in antioxidants as glutathione and in many more molecules (Voet et al. 2008).However, sulphur played an important role long before life appeared on Earth, and its chemistry in the interstellar medium (ISM), the environment where stars and planets are formed, is nowadays a challenge for the astronomical and the astrochemical community.Sulphur is the tenth most abundant element in the Universe.Its cosmic abundance, [S]/[H] = 1.8*10 −5 (Anders & Grevesse 1989), equal to an S/O ratio of 1/37, is used as a reference for measuring sulphur abundance during the several stages of the evolution of a planetary system.While in diffuse clouds the majority of sulphur resides in the gas phase as S + (Jenkins 2009), it is not clear which forms it assumes in translucent and dense clouds (also called molecular clouds, MCs) (Woods et al. 2015;Laas & Caselli 2019).Hily-Blant et al. (2022) suggests that in chemically young cores, atomic sulphur is still the main carrier, which however could not be directly observed.Moreover, it would gradually deplete as the cloud evolves, as also inferred in Fuente et al. (2023).The phenomenon of freeze-out that occurs in MCs could be responsible for the adsorption of gas-phase sulphur onto the surface of icy grains (Caselli et al. 1994), a case in which the presence of H 2 S on their surface would be expected.However, only OCS and tentatively SO 2 were observed in MCs and young stellar object (YSOs) ices (Boogert et al. 2015(Boogert et al. , 2022;;McClure et al. 2023) and only upper limits have been estimated for H 2 S (Smith 1991;McClure et al. 2023).More recent studies suggest that sulphur is mostly contained in organic species (Laas & Caselli 2019) or in a refractory residue composed of species characterized by the S S bond (like H 2 S  with 3 ≤ n ≤ 8, and chains of S  with 2 ≤ n ≤ 3), along with complex species such as hexathiepan (S 6 CH 2 ), which are products of ice processing, as experiments indicate (Ferrante et al. 2008;Jiménez-Escobar & Caro 2011;Cazaux et al. 2022).Nevertheless, numerous sulphur-bearing species (S 2 , S 3 and S 4 , CH 3 SH, C 2 H 6 SH, together with the most common H 2 S, OCS, SO, S 2 , SO 2 and CS 2 ) are detected in comets as 67P/Churyumov-Gerasimenko (Calmonte et al. 2016), whose composition seems to be inherited from the prestellar and protostellar evolutionary phases (Drozdovskaya et al. 2019).Comets contain a record of the chemical composition of the primitive Solar Nebula at the place where they formed, 4.6 Gyr ago (Altwegg et al. 2019).The same applies to a particular type of meteorites, in which sulphur is found in the form of minerals containing iron and nickel (Trigo-Rodríguez 2012).Among the variety of existing meteorites, chondrites represent fairly well the material that surrounded the young Sun and from which the planets formed.Ordinary chondrites contain S mostly in the form of troilite FeS (a variety of pyrrhotite, Fe (1− ) S) (Kallemeyn et al. 1989), while carbonaceous chondrites contain also pentlandite (Fe,Ni) 9 S 8 , and some inorganic sulphates and aromatic organic compounds.(Gao & Thiemens 1993;Trigo-Rodríguez 2012;Yabuta et al. 2007) On the other hand, enstatite chondrites are rich in sulphides like niningerite, alabandite and other alkaline sulphides (Sears et al. 1982).
Thus, there seems to be a missing piece of the puzzle between the gas-phase S-bearing species present in diffuse clouds and the solid sulphides found in the remnants of the protoplanetary disks.To understand the fate of sulphur, we aim to investigate what can happen in diffuse clouds, when some simple sulphur species interact with bare dust grains, a system that already contains the atomic ingredients necessary to form the sulphides observed in more evolved environments.Depending on the C/O ratio of the AGB stars where they are formed, dust grains are formed by carbonaceous materials or silicates and represent 1% of the ISM mass (Potapov & McCoustra 2021).The main families of interstellar silicates are pyroxenes and olivines with chemical composition Mg  Fe (1− ) SiO 3 and Mg 2 Fe (2−2 ) SiO 4 (with x = 0-1) (Henning 2010).These nanometre to micrometre-sized nonspherical, porous dust particles lock up nearly 30% of oxygen and 100% of silicon, magnesium and iron available in the interstellar medium (Van Dishoeck et al. 2013).Commonly, interstellar grain silicates are structurally amorphous, although Mg-rich crystalline silicates have also been observed around young stars (Molster & Kemper 2005;Spoon et al. 2022).
In this work, we modelled dust grains as olivine nanoclusters, as in a recent work by Serra-Peralta et al. (2022).In the latter, both cluster and periodic approaches were considered to study the formation of water, and were previously adopted by Navarro-Ruiz et al. (2014, 2016) to study the formation of H 2 .To our knowledge, the only study of S-bearing species interacting with an olivine surface is a paper by some of us (Rimola et al. 2017) in which butan-1-sulphonic acid was adsorbed on the (010) surface of forsterite (Mg 2 SiO 4 ) to seek the correlation between binding energies and measured abundances of a class of soluble organic compounds found in carbonaceous meteorites.Here, we decided to focus on the cluster approach, firstly developing a reliable computational methodology to simulate olivine clusters, using here the adsorption of H 2 S as a test case.We then characterized the adsorption on the Mg and Fe sites of olivine clusters of other sulphur bearing species, both neutral and charged: S, S + , CS, SH, SH + , SO, SO + , H 2 S, and H 2 S + .For each species, we computed the reaction energy, which can be recast as interaction energy, ΔE, which is a negative quantity for an exothermic reaction.Usually, in the context of astrochemical numerical models, it is the binding energy BE (BE=-ΔE) the quantity that defines the strength of the interaction between the species and the surface.The BE, as defined above, is computed by using the electronic energies of the involved species, assuming that the nuclei are immobile.Quantum mechanics imposes that even at 0 K the nuclei undergo a vibrational motion, which contributes to the electronic energy through the zero point energy (ZPE).It is, therefore the BE(0) = BE − ΔZPE, the key parameter to be used in astrochemical models, determining whether desorption phenomena can take place (Penteado et al. 2017), and also governing the diffusion process, as diffusion barriers are usually assumed to be a fraction of the BE(0) (e.g., Karssemeĳer & Cuppen 2014;Cuppen et al. 2017;Kouchi et al. 2020;Maté, Belén et al. 2020).
The work is organized as follows: in Section 2 we describe the methodology used in the study, Section 3 contains the results, and in Section 4 we discuss our finding and their implications in the sulphur depletion problem.Finally, Section 5 summarizes the conclusions.

METHODS
The olivine clusters presented in Serra-Peralta et al. (2022) and their forsterite parent cluster developed by Escatllar et al. (2019) were adopted to model the core of a dust grain.Each model consists of a cluster of two silicate units with composition Mg (4−  ) Fe  Si 2 O 8 (with x = 0-1).Clusters A, B, and C are characterized by the presence of a Fe atom at positions X, Y, and Z, respectively (see Figure 1).Please, note that W and Z positions are equivalent by symmetry.The free Fe 2+ electronic valence configuration is 4s 0 3d 6 , allowing the system to exist in three different spin states: singlet, triplet, and quintet.As reported in the literature (Navarro-Ruiz et al. 2014, 2016), the most stable configuration is that showing the maximum spin multiplicity, therefore the quintet, followed by the triplet and the singlet states at higher energies.
We first performed a benchmark study with the aim of finding a suitable methodology to describe the adsorption of S-bearing species onto olivine clusters.All the calculations were run with the ORCA program v.5.0.3 (Neese 2022).We optimized (using default settings) clusters A, B, and C in the quintet, triplet and singlet spin states with: i) the methodology adopted in Serra-Peralta et al. ( 2022), namely B3LYP-D3(BJ) (Lee et al. 1988;Becke 1988Becke , 1993b) ) combined with the 6-311++G(d,p) basis set (Hehre et al. 1972;Hariharan & Pople 1973); ii) the meta-generalized-gradient approximation r 2 SCAN-3c, which goes with its own tailored triple zeta quality basis set (Grimme et al. 2021), iii) the range-separated hybrid meta-GGA density functional B97M-V (Mardirossian & Head-Gordon 2017), and iv) BHLYP-D3, a hybrid density functional which features 50% of exact Hartree-Fock exchange (Becke 1993a;Lee et al. 1988).The last two functionals were combined with the def2-TZVP basis set (Weigend & Ahlrichs 2005).The calculations have been performed with the spin-unrestricted formalism in case of open-shell systems (Pople et al. 1995).We compared the values obtained with these DFT functionals against the DLPNO-CCSD(T)/aug-cc-pVTZ (Guo et al. 2018) single energy calculations performed on each optimized geometry.The calculations were carried out with a tight-PNO setup and the default settings for the SCF.
The second part of the benchmark involved the characterization of the adsorption of H 2 S and H 2 O on the forsterite grain nanoclusters, starting from four complexes manually built to allow the two species to interact with sites X, W, Y and with neighbouring oxygen atoms of site W (which as mentioned above is equal to site Z).Adsorbing H 2 S with r 2 SCAN-3c and B3LYP-D3(BJ) results in some cases in its spontaneous dissociation.Thus, we used the Nudged Elastic Band (NEB) algorithm to find the minimum energy path connecting the reactant (physisorbed complex) and the product (chemisorbed complex).The algorithm generates a number of configurations of the atoms, referred to as images of the system (in our case we chose 8 images) that link the two minima with the aim of finding the image with the highest energy.Subsequently, the latter structure is optimized as a transition state following the eigenvector associated with the eigenvalue of the Hessian matrix.
Once we selected an appropriate level of theory, we initially characterized the binding sites of the four olivine clusters by adsorbing H 2 S, H 2 O, CS and CO.Because of symmetry, we identified three adsorption sites (X, W, Y) for forsterite cluster, cluster A and cluster B, while in cluster C the Fe substitution at position Z causes a loss of symmetry and the emergence of an additional binding site.We then characterized the adsorption of the nine S-bearing species on cluster C, which is the most stable structure and presents the largest number of binding sites.Each species was described in its ground state: CS, H 2 S, and SH + as singlets, S + , SH, SO + , and H 2 S + as doublets, and finally S and SO as triplets.Combining the spin states of the S-bearing species with the ground state of cluster C, i.e., the quintet, we obtained adsorption complexes of CS, H 2 S, and SH + in the quintet state, S + , SH, SO + , and H 2 S + in the sextet state and S and SO in the septet state.For each complex, the full set of harmonic frequencies was computed to establish each structure to be a true minimum on the potential energy surface (PES).From them, the zero point vibrational energy (ZPE) correction was computed to correct the binding energy, BE, a quantity defined as the opposite of the interaction energy, obtained with the equation: to be corrected for the ZPE as: where ΔZPE = ZPE   -ZPE  -ZPE   .

Benchmark
Clusters A, B, and C were optimized to characterize their stability from both a structural and an electronic point of view.Data on the quintet and the singlet states are reported in Table 1, while those concerning the triplet state were excluded due to some convergence problems that arose during the optimization of the structures.We did not aim to compare the energy difference between clusters of different spin states, as we know that DFT is not reliable for treating systems in their excited states, but we compared the relative energies between different clusters with the same spin multiplicity (ΔE).Cluster C was the most stable of the group and it was taken as a reference to calculate the ΔE of cluster A and B. Among the four functionals tested in the benchmark, B97M-V showed the smallest error when compared to DLPNO-CCSD(T), with an average error of 11.1% for the quintet and of 21.9% for the singlet spin states.Although the percentage error appeared to be large, the absolute errors corresponded to a few kJ mol −1 , which are always < 4 kJ mol −1 for B97M-V, therefore, in the limit of chemical accuracy.BHLYP-D3(BJ) performed similarly to B97M-V, followed by B3LYP-D3(BJ), while r 2 SCAN-3c, lacking of Hartree-Fock exchange, was the less suitable functional.However, it is worth pointing out the surprising performance of r 2 SCAN-3c in the description of the quintet spin state, despite the presence of unpaired electrons, with an error of only 9.3%.Nevertheless, the representation of the singlet spin state was not satisfying, showing an error of 72.9%.In general, the energies of the clusters in the singlet spin state were affected by larger errors than those of the structures in the quintet ground state, confirming that DFT does not represent a suitable choice when describing excited states.Nevertheless, the stability order documented in the literature was well reproduced with every functional adopted in the benchmark.In terms of the description of the band gap of the system, the outcome largely depends on the percentage of Hartree-Fock (HF) exchange included in the definition of the functional.Indeed, for r 2 SCAN-3c, in which HF exchange is not accounted for, the band gap was almost zero, sometimes responsible for convergence problems, as in the case of the triplet spin state.
To expand our benchmark, we considered the adsorption of H 2 S and H 2 O on the forsterite (Fo) cluster.The two species were manually placed on the Fo and the obtained complexes were optimized with B97M-V, B3LYP-D3 and r 2 SCAN-3c functionals.When comparing the interaction energies computed with the three different functionals against their respective DLPNO-CCSD(T) single energy points, B97M-V/def2-TZVP resulted again in the best performing level of theory (see Table 2).In this case, the percentage errors were much lower than in the previous benchmark due to the larger adsorption energies involved compared to ΔE.However, at least for B97M-V, the absolute error remained below 4 kJ mol −1 .We also point out the great improvement in the performance of r 2 SCAN-3c, justified by the absence of unpaired electrons in the system.
When modelling the adsorption of H 2 S on Fo, we noticed an interesting phenomenon.A spontaneous dissociation of H 2 S was observed in complexes III and IV at r 2 SCAN-3c level of theory, while only complex IV caused the dissociation of the adsorbate at B3LYP-D3(BJ) level of theory, and no dissociation occurred with the B97M-V functional.In these two complexes (III and IV), H 2 S interacts with the Mg cation located in sites W and X, respectively.The dissociated complexes obtained with r 2 SCAN-3c are stable enough, so that a subsequent optimization with B97M-V did not change their structure.The chemisorbed structures were more stable than the physisorbed ones by about 100 kJ mol −1 .Thus, we computed the energy barrier for the dissociation process at B97M-V (that in which spontaneous dissociation is not observed), obtaining a value of 0.2 kJ mol −1 for complex III and 7.2 kJ mol −1 for complex IV.When considering the ZPE correction, the dissociation barriers at 0 K became -3.4 kJ mol −1 and 3.9 kJ mol −1 for complexes III and IV, respectively.The frequency associated with the imaginary mode of the transition state is 231i cm −1 for complex III and 192i cm −1 for complex IV.In Figure 2, the reactant, the transition state structure and the product are shown in the case of complex III, together with a plot of the PES along the reaction coordinate.The dissociation barriers are indeed very small, in one case even submerged below zero, indicating that the cleavage of H 2 S into HS − and H + is likely to occur onto Fo cluster, once the molecule is adsorbed on the Mg sites.When adsorbing H 2 O on Mg, we observed the spontaneous dissociation of the O H bond in complex IV, in which the oxygen atom interacts with the Mg cation of site X at both r 2 SCAN-3c and B3LYP-D3(BJ) levels of theory.We did not computed the PES of the dissociation at B97M-V level of theory, however, we assume that a small barrier for the process is present, as for H 2 S. Thus, we can conclude that the dissociation barrier depends on the methodology and that can either be submerged below zero or be a few kJ mol −1 .The two benchmark studies indicate B97M-V/def2-TZVP to be the best performing methodology with an average error of 16.5% on the relative stability of olivine clusters and of 2.6% on H 2 S and H 2 O adsorption energies, becoming this functional our choice for the following calculations.

Cluster versus periodic approach
The computational approach adopted in this work to study the interaction of S-bearing species with the silicate core of the grains may appear too simplistic due to the limited size of the cluster models.The exposure of the metal centres could be responsible for their higher reactivity, which could influence the outcomes of the calculations.Enlarging the cluster model to increase the coordination of the metal centres is not practical due to the steep increase of the computational cost of the calculations.Therefore, we adopted a periodic approach to model different extended forsterite slab surfaces, with the closest level of theory to that of the cluster.The periodic approach represents the opposite situation compared to the cluster, as the metal center binding sites (albeit at the surfaces) possess the highest coordination due to the extended nature of the slabs.The cluster-based adsorption complexes described in Table 2 were optimized at PBE level of theory, combined with the def2-TZVP basis set.For the periodic calculations, H 2 S was adsorbed onto the (001), ( 021), ( 101), ( 110), (111), and (120) forsterite surfaces, which were previously characterized in Zamirri et al. (2017) and Bancone et al. (2023).The periodic adsorption complexes were also optimized at PBE, and the forsterite atoms were described with the basis set proposed by Bruno et al. (2014), while H 2 S with the same def2-TZVP basis set adopted for the cluster calculations.We did not include dispersion interactions in the description of the system, as it was not possible to reproduce the desired level of theory (D2 with customized dispersion coefficient for the Mg atom) in both the cluster and the periodic calculations due to specific implementation of the D2 coding in the ORCA program.Despite the missing dispersion term, the calculations show that H 2 S still spontaneously dissociates on the Fo cluster, in complexes I, II, III and IV.On the periodic systems, molecular physisorption is found on the most stable forsterite surfaces, ( 120) and ( 001), while on the less stable ( 101), ( 111), (021), and ( 110), H 2 S dissociative chemisorption takes place.Therefore, we can thus conclude that H 2 S chemisorption is not simply an artifact of the size of the cluster that exposes low-coordinated metal centres.Depending on the surface and on the specific adsorption site, it is possible to have spontaneous dissociation of H 2 S also on large and extended surfaces in which the ion coordination is more complete than in the cluster.We remark that the level of theory chosen to characterize the system may be responsible for the presence or absence of a small barrier for the dissociation process.

Characterization of the nanocluster binding sites
We characterized the binding sites of each nanocluster (Fo and the olivines) by adsorbing H 2 S, CS, and their oxygen analogues.H 2 S and H 2 O are species with an electronegative atom (S and O) that preferentially interacts with the positively charged metal center, and two hydrogen atoms that can interact via H-bond with the oxygen atoms of the silicate unit, when geometrically feasible.On the other hand, CS and CO are species whose dipole shows a concentration of electronic charge on the carbon atom, and accordingly this atom is the one interacting with the metal center.The bond length of CS in the gas phase is 1.523 Å, and when it is adsorbed on Mg it shortens to 1.500 Å, due to the electron donation coming from the  orbitals with moderate antibonding character to Mg.When CS is adsorbed on Fe, the bond length becomes 1.509 Å, due to the electron backdonation effect of Fe on the antibonding orbitals of CS.The same occurs with CO, in which the CO bond lengths in the Mgand Fe-adsorbed complexes are 1.117 Å and 1.120 Å, respectively, shorter than that in the gas-phase molecule (1.123 Å).Remarkably, the largest CS bond shortening compared to CO, together with the differences in their dipole moments (0.1 and 1.9 Debye for CO and CS, respectively), causes the BE of CS to be almost doubled with respect to CO.We previously observed the same phenomenon when characterizing the interaction of CO and CS on water ice (Perrero et al. 2022), where the larger dipole of CS caused an increase in BE of 50% with respect to the CO BE value.This further strikes the difference in the chemical behaviour of CO and CS.
Table 3 shows a summary of the results.We obtained a range of values for the adsorption onto Mg centers and one value for each Fe.It appears that CS and CO prefer to interact with Fe rather than Mg, while H 2 O clearly shows a larger affinity with Mg rather than with Fe.However, the behaviour of H 2 S is not straightforward.The interaction with both metals seems to be comparable in clusters A and B, but cluster C represents an exception.In fact, there is a notable difference between the 74-85 kJ mol −1 BE values computed when H 2 S interacts with Mg and the BE value of 97 kJ mol −1 computed for the adsorption onto Fe.An overall look at the values, particularly the BEs obtained for the Fe site at different positions (X, W, Y, and Z) hints at a possible role of the coordination of the metal center in determining the binding strength.Indeed, sites W and Z (cluster C) are the most exposed and less coordinated, and they tend to give the highest BE, while site X (cluster A) is the most coordinated and is usually the weakest one.In general, cluster C, in addition to being the most stable structure and offering the largest number of binding sites, is also covering the entire set of binding energies obtained with the entire set of nanoclusters.All of these reasons are in favour of choosing this structure (cluster C) for the investigation of the adsorption of the whole set of S-bearing species.

Adsorption of S-bearing species
The set of the adsorbed species includes CS, H 2 S, H 2 S + , HS, HS + , S, S + , SO, and SO + .These species have been selected for being the most abundant sulphur carriers in diffuse clouds and the external layers of photon-dominated regions (PDR) where the water ice has not been formed yet, either in the cloud envelopes or the dense PDRs formed in the vicinity of recently formed young massive stars (Fuente et al. 2003;Neufeld et al. 2015;Goicoechea et al. 2021;Fuente et al. 2023).We modelled four adsorption complexes for each species, from which some common characteristics emerged: i) we can draw a colour map (see Figure 3) to highlight the strength of the adsorption sites which is valid for all the species; ii) the Fe S interaction is stronger than the Mg S one (see Table 4).Although we may attribute part of the reason to the coordination of the metal center, it is not the only explanation to such differences; iii) SO and SO + interact with the metal center through oxygen, which in this case is the most electronegative atom.This is responsible for the differences in the adsorption energies listed in Table 4; iv) the presence of a charge dramatically increases the BE, as we observed for the four ions included in this study, (H 2 S + , HS + , S + , and SO + ).
From the calculations, we found that H 2 S + and HS + are physisorbed on Mg, but chemisorbed on Fe, in which the proton transfers to one of the oxygen atoms belonging to the closest silicate unit (similarly to what we observed for H 2 S on Fo), and the charge is distributed onto the entire system.Moreover, when an open-shell species adsorbs on Fe, part of its spin density transfers from the sulphur atom to the metal, increasing the binding energy as in the case of S and HS.On the other hand, this effect is prevented for SO, whose two unpaired electrons remain localized onto the radical species, due to the adsorption through the oxygen atom.The BEs of both SO and SO + can be explained by looking at the behaviour of H 2 O, in which the interaction with the metal center takes place through oxygen and for which the interaction Mg O is favoured over the Fe O one.Finally, the inclusion of the ΔZPE correction in computing the BE(0) does not alter the trend derived from the analysis of the BEs.

DISCUSSION AND ASTROPHYSICAL IMPLICATIONS
In a previous work by some of us (Perrero et al. 2022), the BEs of 17 S-bearing species adsorbed onto crystalline and amorphous water ice were computationally studied adopting a periodic approach and the DFT//HF-3c level of theory.The B3LYP-D3 functional was used to compute the energy of closed-shell systems (such as CS and H 2 S), while for open-shell species (such as HS, S, and SO) M06-2X was adopted.These functionals were combined with an Ahlrichs triple zeta valence quality basis set, supplemented with a double  set of polarization functions (Schäfer et al. 1992).Although the computational approach adopted in Perrero et al. (2022) differs from this work, it does not justify the large differences between the BE(0) values of the same set of species on the amorphous ice surface and the olivine nanoclusters (see Table 5 and Figure 4).While in the first case the interactions are weak and controlled mainly by the dispersion forces, BE(0)s computed in this work are large and determined mostly by electrostatic interactions, to which occasionally electron transfer phenomena add up.
The consequences of such a strong interaction between sulphur and the mineral surface is that S-containing species have a high chance to freeze out onto the core of dust grains and not be able to desorb into the gas phase even at the temperatures found in a diffuse cloud (around 100 K).In those environments, the presence of charged gas-phase species should also favour the adsorption process as long as the grains are negatively charged (Cazaux et al. 2022;Fuente et al. 2023).Even if the grain particles are mostly neutrals or positively charged, the accretion of neutral sulphur species on the grain surfaces could enhance the number of sulphur atoms in refractories.On the other hand, it should be considered that dust grains in dense clouds are not completely covered by the ice mantle, therefore still exposing part of their mineral surface into the gas Table 5. Binding energies (BE(0), in kJ mol −1 ) of five S-bearing species and water on amorphous ice (Ferrero et al. 2020;Perrero et al. 2022) (Ferrero et al. 2020;Perrero et al. 2022) and on olivine.(Potapov & McCoustra 2021).As it clearly appears in Figure 4, H 2 O adsorption on olivine is more energetic than H 2 S, meaning that water is more likely (and probable) than hydrogen sulphide to form a stable interaction with the mineral.However, the modest interaction between sulphur-containing species and water ice allows these molecules to diffuse on its surface and likely move towards the portion of exposed core of the grains.Therefore, whether considering the physical conditions of a diffuse or those of a dense cloud, Sbearing species would have a certain chance either to directly stick onto the grain core or diffuse from the mantle towards the exposed fraction of the core grain where they become strongly chemisorbed.Another interesting issue is the tendency of certain species, namely H 2 S, HS + , and H 2 S + , to chemisorb onto olivine nanoclusters, planting the seed for the formation of the Fe S bond present in sulphide minerals such as troilite and pentlandite.In the work of Kama et al. (2019), approximately 90% of the sulphur is predicted to be locked in the refractory residue of protoplanetary disks, the main carriers being precisely sulphide minerals, followed by sulphur chains S  , through to the study of the photospheres of a sample of young stars.However, chemical models (Druard & Wakelam 2012) and laboratory experiments (Woods et al. 2015) predict a scarce abundance of sulphur chains in these environments, reason why sulphide minerals are considered the major reservoir of refractory sulphur in the protoplanetary disk.These two minerals are also found in cometary dust and meteoritic rocks (Kallemeyn et al. 1989;Trigo-Rodríguez 2012).From a thermodynamic point of view, Fe S is also more favoured than Mg S, as the latter is only found in the less common enstatite chondrites, that due to their formation in reducing conditions, are rich in alkaline sulphides like the Mg-containing niningerite and alabandite (Sears et al. 1982).
In this work we did not consider charged dust grains, even though the model of Ibáñez-Mejía et al. (2019) predicted them to be positively charged in a diffuse medium and negatively charged in denser environments.If this was the case, the interaction between gas-phase cations and positively charged grains in diffuse clouds would clearly be hampered by electrostatic repulsion.However, our results concerning positively charged sulphur-bearing species adsorbed on neutral dust grain cores hold true for the case in which neutral gas-phase species interact with charged dust grains.Although the atomistic details of the system should vary due to the presence of a charge on the olivine nanocluster, the interaction energies should be of the same order of magnitude, since the presence of a charge in the system would still be the main driving force of the interaction.On the other hand, Fuente et al. (2023) argue that sulphur would remain ionized in the gaseous medium until grains are expected to be negatively charged.If this were the case, the electrostatic interactions would be responsible for an increased sulphur depletion effect.
According to what Hily-Blant et al. ( 2022) suggested, the fact that sulphur depletion increases progressively as the core evolves could be in agreement with the fact that during the evolution, atomic S freezes out on the surface of grains and, depending on which surface it reaches, it can either become part of the mineral core or as Cazaux et al. (2022) propose, it can fall on the ice mantle and be subjected to photoprocessing, thus reacting and forming a refractory residue.
Our calculations are also relevant to undestand the sulphur chemistry in the dense PDRs formed on the molecular cloud/HII regions interfaces.It has been considered that in these dense warm regions (dust temperature ∼ 100 K), the chemistry is well explained with gas-phase networks since freeze-out on grain surfaces is negligible.Our new calculations suggest that this may not be the case since the binding energies of sulphur species is up to a factor of ∼10 in bare grains.This means that grain-gas interactions need to be considered in these warm environments.
Finally, the assumption that S-bearing species become part of the grain core, rather than interacting with the icy mantles, would be in agreement with the fact that they are usually used as shock tracers (e.g., Sato, M. T. et al. 2022;Zhang et al. 2023).In shock regions, which are caused by the impact of stellar outflows with the quiescent surrounding gas, H 2 S can be sputtered off the grains and enter in the gas phase, where it marks the zone that has been affected by such an impact (Woods et al. 2015).The harsh conditions of the shock would therefore provide enough energy to break the strong interaction between sulphur and the mineral surface, allowing its transformation in H 2 S and successive ejection in the gas phase (Holdship et al. 2017).

CONCLUSIONS
In this work, we studied the interaction of nine S-bearing species, CS, H 2 S, H 2 S + , HS, HS + , S, S + , SO, and SO + , which include both neutral and charged, as well as closed-shell and open-shell systems.We characterized their adsorption complexes and computed the corresponding binding energies on olivine nanoclusters (forsterite and the Fe-containing analogues) by means of quantum mechanical calculations, to simulate the interaction of a gas-phase species with the surface of a silicate bare dust grain in the diffuse interstellar medium.Our main findings are that: i) from a thermodynamic point of view, the Fe S interaction is more stable than Mg S, while when the oxygen end is interacting with the metal center we notice the opposite trend; ii) S-bearing species are adsorbed much more strongly on the core of bare dust grains than on their icy mantles (represented by pure amorphous water ice surfaces; iii) the metal centers of the olivine clusters can be responsible for the dissociation of hydrogenated Sbearing species, especially if the latter are positively charged; iv) the dissociation of H 2 S on forsterite cluster presents a small barrier, that is likely to be overcome when considering the energy liberated by the adsorption process.The implications of such findings are that sulphur species have a high probability of sticking onto the surfaces of bare dust grains and becoming part of their core in the form of refractory materials, as predicted by a number of experimental and observational studies, in addition to chemical models.Indeed, more investigations are needed in order to finally unveil sulphur chemistry in the ISM.However, the trapping of sulphur into the grain cores in the early stages of the evolution of a prestellar core could, in part, account for the S-depletion in dense cores.

Figure 1 .
Figure 1.Forsterite and olivine clusters used in this work as dust grains models.The labels in cluster C are useful to identify adsorption sites.

Figure 2 .
Figure 2. The dissociation of H 2 S on the forsterite cluster (complex III) characterized at the B97M-V/def2-TZVP level of theory.It exhibits a small potential energy barrier and becomes barrierless when taking into account the ZPE corrections.

Figure 3 .
Figure3.Colour map showing the strength of each binding site from weaker (green) to stronger (red), depending on which species has been adsorbed.It is valid for all the species characterized in this work, with the exception of SO and SO + , for which the Fe is the weakest binding site.In the case of H 2 S + and HS + , the cation dissociate when adsorbed onto Fe.

Table 1 .
Relative energies (ΔE in the text) of olivine clusters optimized with different DFT functionals and their DLPNO-CCSD(T) single energy point calculations.Cluster C is the most stable one and is taken as a reference to compute the ΔE of clusters A and B.

Table 2 .
Adsorption energies (ΔE  in the text) of H 2 S and H 2 O on the forsterite cluster optimized with different DFT functionals and their DLPNO-CCSD(T) single energy point calculations.Complex I and II differ for the orientation of the adsorbate with respect to the binding site (#).The asterisk (*) indicates the cases in which spontaneous dissociation takes place.The percentage error is obtained as the average over those computed for each DFT method against the DLPNO//DFT energy value.The average is done over the four complexes of both H 2 O and H 2 S together.

Table 3 .
BEs (in kJ mol −1 ) of H 2 S, H 2 O, CS and CO adsorbed on forsterite and olivine clusters.BE ranges result from the presence of different adsorption sites on the same cluster.

Table 4 .
Calculated binding energies (BE and BE(0), in kJ mol −1 ) for all the considered S-bearing species on the cluster C.
Binding energies (BE(0), in kJ mol −1 ) ranges covered by CS,H 2 S, HS, S, SO and H 2 O adsorbed on amorphous ice