Diastereoselective Formation of Trans-HC(O)SH through Hydrogenation of OCS on Interstellar Dust Grains

Sulfur-bearing molecules constitute a fundamental branch of prebiotic chemistry due to the recently gained importance of cysteine C₃H₇NO₂S as a catalyst in the assembly of complex peptides (Foden et al. 2020). In cold, interstellar environments, cysteine is far from being detected, but other complex organic molecules (COMs) of prebiotic importance bearing S have recently been detected, highlighting ethyl mercaptan C₂H₅SH or thioformic acid (HC(O)SH; Kolesniková et al. 2014; Rodríguez-Almeida et al. 2021). In turn, the family of compounds to which HC(O)SH belongs (thioacids) has also been pointed out as possible catalysts in prebiotic processes in a primordial Earth (Chandru et al. 2016). Particularly puzzling is the fact that the presence of these molecules on Earth could have a (at least partial) panspermic origin, owing to both a confirmed presence of sulfur-bearing species in comets (Korth et al. 1986; Krishna Swamy & Wallis 1987; Calmonte et al. 2016; Rubin et al. 2019) and the fact that the simplest thioacid (thioformic acid) has been recently identified in the interstellar medium (ISM; i.e., the giant molecular cloud G+0.693; see Rodríguez-Almeida et al. 2021). The chemical mechanism behind the formation of thioformic acid is unknown, but such a mechanism needs to explain the formation of HC(O)SH from simple, sulfur-bearing precursors. In Rodríguez-Almeida et al. (2021) several possible routes were postulated for the formation of HC(O)SH, namely,

In this work, we have focused on reaction (3) on interstellar dust grains, having two main objectives in mind. First, we wish to check whether or not the reaction proceeds effectively, and, second, we wish to determine the stereoisomerism of the reaction. The study of reaction (3) sparked our interest for two main reasons. First, OCS is a moderately abundant molecule in the astronomical source where HC(O)SH was detected (Armijos-Abendaño et al. 2014; Rodríguez-Almeida et al. 2021); second OCS is the only sulfur-bearing species positively detected on interstellar ices (Palumbo et al. 1997).

Under the harsh conditions operating in molecular clouds, it is expected that a large fraction of the molecular material efficiently depletes on ice-covered dust grains. It is on the surface of these grains that a significant portion of the reaction networks operate, with a prevalence of hydrogenation reactions (Tielens & Hagen 1982; Garrod et al. 2008; Hama & Watanabe 2013) and hydrogen saturated species (with noteworthy exceptions; see, for example, Noble et al. 2015). In interstellar molecular clouds, atomic hydrogen is efficiently formed from the dissociation of H₂, the most abundant interstellar molecule (Wakelam et al. 2017a), by cosmic rays (Padovani et al. 2018). These newly formed H atoms accrete on interstellar dust grains at a rate of ∼1 atom per day (Wakelam et al. 2017a). In addition, note that G+0.693, i.e., the molecular cloud in the Galactic Center where HC(O)SH has been detected, has been proposed to be located in an environment with a large amount of H atoms available (Requena-Torres et al. 2006, 2008). This is due to the fact that a high cosmic ray ionization rate is measured in the Galactic Center (factors of 100–1000 higher than the standard cosmic ray ionization rate in the Galactic disk; Goto 2014) and this induces a strong radiation field of cosmic ray induced secondary UV photons. The prevalence of hydrogenation is due to atomic hydrogen being sufficiently mobile to scan the surface of dust grains in search of possible reaction partners (Hama & Watanabe 2013), but also is due to tunneling, a quantum effect in nature that permits light particles to proceed through kinetic barriers.

Tunneling is often invoked as the mechanism behind hydrogenation reactions in the ISM, both via an external incoming H atom (Lamberts & Kästner 2017; Meisner et al. 2017; Álvarez-Barcia et al. 2018; Oba et al. 2018; Miksch et al. 2021; Molpeceres & Kästner 2021) and in intramolecular hydrogen/proton migrations reactions (Rani & Vikas 2020a, 2020b, 2021; García de la Concepción et al. 2021). An important trait of the reactions at the cold temperatures of the ISM (and specifically in this case of H atoms) on dust grains is that they yield kinetically controlled reaction products, meaning that lower reaction barriers are correlated with molecular abundances (Loomis et al. 2015; Shingledecker et al. 2019). This also has severe implications in the isomerism of the reactions (Shingledecker et al. 2019, 2020). It is important to note that it has been recently proved that in the gas phase of the ISM, thermodynamic equilibrium can be achieved, even at shallow temperatures (García de la Concepción et al. 2021), but such equilibrium cannot be taken for granted on dust grains, which experience constant H accretion, diffusion, and reaction (Wakelam et al. 2017a). For HC(O)SH, Rodríguez-Almeida et al. (2021) were only able to detect trans-HC(O)SH, the most stable structural isomer of the two possible ones, with the cis one lying 3.1 kJ mol⁻¹ above it in energy (Kaur & Vikas 2014).

As mentioned above, OCS was postulated as a possible precursor of HC(O)SH because of its relatively large abundance in G+0.693 (Armijos-Abendaño et al. 2014; Rodríguez-Almeida et al. 2021). The sequence of reactions that we have studied include a two-step reaction initiated by:

and followed by (further reactions on the HOCS are not considered, because they are not competitive):

From all of them, we found that reactions (4) and (7) determine the most likely reaction pathway, locking the reaction product in trans-HC(O)SH. Both reactions constitute a diasteroespecific pathway for the formation of trans-HC(O)SH coherent with the absence of cis-HC(O)SH in astronomical observations (Rodríguez-Almeida et al. 2021). We provide reaction energies and, when applicable, activation energies and rate constants accounting for tunneling for all reactions. Furthermore, for radical−radical reactions, we give a tentative branching ratio for the reaction.

The OCS + 2H → HC(O)SH reaction on amorphous solid water (ASW) was investigated by alternating two different surface models: the explicit and the implicit model. In the explicit model, we incorporated the effects of a surface by placing the reactants on top of a 20 water molecule cluster (Shimonishi et al. 2018; Molpeceres & Kästner 2021). In the implicit model, we took the effect of a surface in the reaction into account by fixing the value of the rotational partition function of reactant and transition states to unity (Meisner et al. 2017). The latter approximation allows us to compute accurate reaction barriers with computationally expensive theoretical methods. Such an approximation breaks down in the presence of significant binding between an adsorbate and a surface.

The protocol we followed in deciding which model to employ in each step involves computing the different adsorbates’ binding energies on a water ice cluster containing 20 molecules, mimicking a fragment of interstellar ice, and starting with OCS as an adsorbate. We did the water ice cluster and adsorbate placement following the same procedure as in Molpeceres & Kästner (2021). Very briefly, molecular dynamics simulations (MD) were used to generate five amorphous clusters, of 20 water molecules each, by a heating (100 ps) and cooling (10 ps) loop using the recently developed GFNFF potential (Spicher & Grimme 2020). The resulting clusters were optimized using density functional theory (DFT) and the PW6B95-D4 exchange and correlation functional (Zhao & Truhlar 2005; Caldeweyher et al. 2019), using the def2-SVP basis set (Weigend & Ahlrichs 2005) for the optimization of the structures and the def2-TZVP basis set (Weigend & Ahlrichs 2005) for energy refinement. The theoretical method is abbreviated as PW6B95-D4/def2-TZVP//PW6B95-D4/def2-SVP. All open-shell DFT calculations were carried out using an unrestricted wave-function formalism.

Once the water clusters were generated, the adsorbates were distributed on its surface by placing the adsorbates’ center of mass (CM) in random points of a spherical grid distorted to an ellipsoidal grid surrounding the cluster (Molpeceres & Kästner 2021). Ten samples were placed per adsorbate and ice pair, at a minimum distance of ∼3 Å between the CM of the adsorbate and any other atom in the cluster. The initial orientation of the adsorbate was randomized.

Binding energies were calculated according to:

using the PW6B95-D4/def2-TZVP//PW6B95-D4/def2-SVP level of theory. A low average binding energy of the adsorbate indicates physisorption and a weak influence of the surface and justifies using the implicit surface approach. When such was the case (for example, for OCS; see below), the associated magnitudes of the reaction, reaction energies, and activation barriers were calculated employing the high-level (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP method (Adler et al. 2007; Knizia et al. 2009). Moreover, reaction rate constants accounting for tunneling were obtained using reduced instanton theory above the crossover temperature for tunneling (McConnell & Kästner 2017) and instanton theory (Langer 1967; Miller 1975; Coleman 1977; Kästner et al. 2009; Rommel & Kästner 2011; Rommel et al. 2011) below the crossover temperature at this level of theory.

Crossover temperatures, understood as temperatures where the tunneling regime starts to dominate, are defined as:

Here, ω_i corresponds to the absolute value of the frequency of the reaction transition mode. When instanton rate constants were calculated, we constructed and optimized a discretized Feynman path consisting of 48 images, starting at T ∼ 0.7 T_c. A sequential cooling scheme was applied until we reached a minimum temperature of 50 K. The number of images at lower temperatures was increased to ensure convergence of the path with the number of images (up to a maximum number of 186 images at 50 K). Symmetry factors (σ), accounting for the degeneracy in the reaction channels were not included because such a degeneracy should break in a surface (Fernández-Ramos et al. 2007).

When high binding energies are obtained during the sequential study of the different hydrogenation reactions, the implicit model approach is no longer a valid model, and explicit inclusion of water molecules is mandatory. Such is the case for the OCSH and OCHS radicals studied in this work. The hydrogenation of these radicals involves the study of radical−radical reactions. We studied these processes in the presence of the previously described water clusters by placing hydrogen atoms in the vicinity of the radicals and optimizing the resulting structures following a similar protocol as for the adsorbate placement on ice (a minimum distance of 2 Å to the surface and of ∼2.5 Å to the adsorbate’s CM). Although a minimum distance of 2.5 Å is relatively short, it allows us to obtain a good amount of reactive events as a function of the initial configurations of the H atoms. For all of the outcomes of this simulation, we tested that there were no long-range entrance barriers. Finally, it is essential to mention that radical−radical recombinations must be modeled using a broken-symmetry formalism (Enrique-Romero et al. 2020), which requires a careful initialization of the potential energy surface for the system. In this study, the generation of biradical, open-shell singlet solutions was ensured in a two-step way. First, we converged our system in the high-spin state (triplet) at long internuclear distances, and second, from this solution, we safely converged to the open-shell singlet solution.

The codes we have employed during this work are: ASE and GFN-xTB (GFNFF model) for the MD simulations employed during the cluster generation (Hjorth Larsen et al. 2017; Spicher & Grimme 2020), ChemShell (Sherwood et al. 2003; Metz et al. 2014) for the geometry optimizations, transition state search and instanton optimization, Turbomole v7.5 (Furche et al. 2014) for the DFT calculations, and Molpro 2019 (Werner et al. 2012) for the CCSD(T)-F12 calculations.

3.1. Hydrogenation of OCS

To confirm that the implicit surface model is applicable in the hydrogenation of OCS we have computed the distribution of binding energies of OCS in ASW. From a total of 50 samples obtained in five different model clusters, we got an average binding energy = 1613 K with a standard deviation of 639 K. This average is within the values reported in Wakelam et al. (2017b) for a model of OCS interacting with one water molecule. Such a weak binding represents a physisorbed state. It has been discussed previously that in situations of weak binding and at low temperatures, several diffusion mechanisms can be invoked to justify the consideration of the tail of the binding energy distribution as the predominant binding energy (Shimonishi et al. 2018; Molpeceres & Kästner 2021). If we consider this argument, the five highest binding energies for OCS are in the range of ∼2800–3100 K, in evident agreement with the experimental results of Collings et al. (2004), 2888 K. In thermal desorption events, the whole distribution of binding energies must be considered because at the desorption temperature, diffusion and desorption are competitive processes, and an interchange of binding sites can be assumed. In this work, we are interested in reactions at low temperatures, and thus configurations with higher values of the binding energy constitute a better approach to the situation found in quiescent clouds.

In Figure 1 we present two snapshots obtained of optimized structures on the water cluster surface belonging to the highest binding group. From the visualization of the picture, it appears evident that even in the deepest binding sites, there should not be any preferential orientation of the OCS on the surface. Therefore, the choice of the implicit surface models is valid for the OCS + H → OCS(H) reaction, and the different hydrogenation channels have been modeled at the (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP level of theory. In the previous equation, the H atom is in parentheses because from the first hydrogenation reaction, three different structural isomers, namely OCSH, OCHS, and HOCS, are possible. Furthermore, OCSH presents two different diastereomers, ³ cis and trans (see Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The activation barriers and reaction energies for the three structural isomers are summarized in Table 1. From the table, we observe that all proceed with an activation barrier. On the one hand, hydrogenations at the sulfur atom and the carbon atom are exothermic processes presenting moderate barriers worth further investigation. The geometries of the transition states for the hydrogenations on the sulfur and carbon atom are presented in Figure 3. On the other hand, the hydrogenation at the oxygen atom gives an activation barrier that is too high to be competitive with the other two processes, especially considering that no significant orientation in the surface has been found in our binding energy calculations. Moreover, hydrogenation on the O atoms is slightly endothermic. The hydrogenation on the sulfur atoms produces the cis-OCSH isomer exclusively because it proceeds via a slightly bent transition state (see Figure 4 for an intrinsic reaction coordinate (IRC) of the hydrogenation reaction showing this effect). The IRC energies are not corrected with (U)CCSD(T)-F12a/cc-pVTZ-F12 energies, and that is the reason behind the seemingly lower activation energy. The values of Table 1 remain more accurate, and the IRC results are presented for visual purposes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Activation Barriers (ΔU_A) and Reaction Energies (ΔU_r ; in kJ mol⁻¹, ZPE Corrected) for the Different Hydrogenation Channels of the OCS + H Reaction

Reaction	ΔUA	ΔUr
OCS + H → cis-OCSH	21.1	–46.8
OCS + H → OCHS	38.4	–41.1
OCS + H → HOCS	92.5	6.6

Note. Values are obtained in the gas phase to be used under the implicit surface approach for the rate constants at the (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP level.

Download table as:  ASCIITypeset image

The reaction rate constants for the discussed barriers are presented in Figure 5. From the figure, we can see that the hydrogenation at the sulfur atom is several orders of magnitude faster than the same reaction in the carbon atom: k_S(50 K) = 1.3 × 10⁴ s⁻¹, k_C(50 K) = 3.9 × 10¹ s⁻¹, or three orders of magnitude of difference. In this context, the reaction at the sulfur atom should dominate the different hydrogenation channels. However, we will also study the second hydrogenation on the carbon atom in subsequent sections. It is important to note that during the review process of this article, a similar one (Nguyen et al. 2021) has been made public as a preprint that shares some of these conclusions. The agreement of the comparable calculations (binding energies, activation barriers) is excellent between the two works.

Zoom In Zoom Out Reset image size

3.2. Second Hydrogenation

The second hydrogenation processes have been modeled for the reactions cis-OCSH + H → products and for the OCHS + H → products. Due to the high barrier for the hydrogenation on the oxygen atom, we have not further considered the reaction HOCS + H → products in this section. In addition, to confirm the stereoselectivity of the whole process, we have modeled the reaction cis-OCSH → trans-OCSH both in the gas phase and on the surface of the ice.

3.2.1. Cis-OCSH + H

To model the second hydrogenation reaction, the first step is to obtain an approximation of the binding energy distribution and, most notably, an approximation of the upper binding tail of such distribution, as explained in the previous section. We have studied the binding properties of this radical as for OCS, obtaining an average binding energy = 3284 K at the PW6B95-D4/def2-TZVP//PW6B95-D4/def2-SVP level of theory, with a standard deviation of 1256 K. This binding energy is higher than for OCS, but of particular interest are the situations of stronger binding. A visualization of the optimized geometries for binding energies higher than 4000 K (12 samples) shows that in all cases, the binding of the cis-OCSH is highly oriented, with the hydrogen of the radical interacting with one oxygen of the water in the cluster and one hydrogen of an adjacent water molecule interacting with the oxygen atom of the radical, in an “anchor-like” configuration (Figure 6 shows this binding conformation for one of the highest binding situations of the distribution). Despite the relatively large binding energy, a closer look at the bond distances of the radical in the gas phase and on the surface reveals minor structural changes, pointing to a physisorbed nature of the adduct. This behavior in the surface is in contrast to the structural changes associated with sulfur compounds in the solvated phase (e.g., thiourea; Vikas et al. 2007). The radical orientation on the surface is crucial to explain the subsequent reactivity for two reasons. First, such a fixed position of binding sites hinders the internal rotation of the radical, precluding the isomerization. We have quantified the influence of this effect on the reactivity in Section 3.2.3. Second, and more importantly, such a binding mode leaves only a single molecular face for a reaction, fixing the possible outcomes of the second reaction in either cis-OCSH + H → trans-HC(O)SH or cis-OCSH + H → H2S + CO. We have confirmed this effect by determining the second hydrogenation channels of cis-OCSH on the surface, following the procedure described in Section 2.

Zoom In Zoom Out Reset image size

From one of the binding configurations of higher binding (the one depicted in Figure 6) we placed H atoms around the cis-OCSH radical, satisfying the constraints presented in Section 2. The configurations generated this way (70) were subsequently optimized, visualizing the reaction output to obtain an estimation branching ratio of reaction. Figure 7 represents the initial geometries for this second reaction that were subsequently and sequentially optimized. In the complex potential generated by the ice + adsorbate, we found many optimizations to yield pre-reactive complexes or co-deposition of H on other binding sites of the water ice. This is an artifact of the geometry optimization procedure. In actual ISM conditions, it is expected that both pre-reactive complexes and H bindings are short-lived in the presence of a near radical. The branching ratios provided here are therefore subjected to a high uncertainty and should be regarded as qualitative. The computation of accurate branching ratios requires bigger structural models and sampled configurations and, ideally, molecular dynamics simulations to determine the outcomes. Such extensive sampling is unfeasible due to computational limitations. We have estimated the uncertainty in our branching ratios profiting from the binomial nature of the product channels. Hence, we computed the confidence intervals with a confidence level of 90% for a binomial distribution including all reacting events and using the Jeffreys method, as implemented in the statsmodel library (Seabold & Perktold 2010).

Zoom In Zoom Out Reset image size

We have only found two possible reaction products for the hydrogenation reaction, either the title product (trans-HC(O)SH) or a mixture of H₂S + CO, the former being dominant in our simulations. A detailed summary of our results for this reaction can be found in Table 2. In the table, we make the distinction between the reaction energy for the reaction on the ice (ΔU_r,i) and in the gas phase (ΔU_r,g). The former is affected by the binding of the products with the surface and depends on the binding site, whereas the latter is univocal. Combining the rate constants present in Figure 5 with the results of the second hydrogenation, we arrive at the most important conclusion of the current work: the recently detected trans-HC(O)SH molecule (Rodríguez-Almeida et al. 2021) can be effectively synthesized on the surface of interstellar dust grains. This reaction happens via two subsequent H additions to OCS, a relatively abundant, sulfur-bearing interstellar molecule and the only one positively detected in interstellar ices (Palumbo et al. 1997), as we mention in the Introduction. Furthermore, the structure of the primary intermediate of reaction, cis-OCSH, ensures that the reaction is diastereoselective, further agreeing with the observations.

Table 2. Reaction Outcome, Branching Ratios, and Reaction Energies on the Ice (ΔU_r,i, in kJ mol⁻¹, ZPE Corrected) and in the Gas Phase (ΔU_r,g, in kJ mol⁻¹, ZPE Corrected) for the Cis-OCSH + H Reaction

Outcome	Branching Ratio	ΔUr,i	ΔUr,g
trans-HC(O)SH	0.86 (0.71–0.94)	−360.9	–356.5
CO + H2S	0.14 (0.05–0.29)	−336.3	–347.6

Note. The numbers in the table are obtained using PW6B95-D4/def2-TZVP//PW6B95-D4/def2-SVP.

Download table as:  ASCIITypeset image

3.2.2. OCHS + H

Although the formation of cis-OCSH is expected to be the principal outcome of the first hydrogenation reaction evinced by the different rate constants for OCS hydrogenation portrayed in Figure 5, we have further developed the description of the second hydrogenation after an initial reaction in the C atom of the OCS molecule. A similar procedure for obtaining the binding sites and energies was employed for the OCHS radical. We found the binding energy for this species to be = 4272 K with a standard deviation of 1294 K, which is moderately higher than for cis-OCSH. A look at the situations of higher binding (binding energies higher than 5200 K, 11 samples) reveals that, contrary to the case of cis-OCSH, a preferential conformation is not straightforward to assign.

Following the procedure outlined in the previous section, we evaluated the products of second hydrogenation for OCHS in one of the deepest binding sites (E_bin = 6458 K). The same procedure for the construction of the initial configurations rendered 67 samples. The second hydrogenation in this reaction produces both trans-HC(O)SH and cis-HC(O)SH in nearly a 65%/35% ratio (see Table 3). This ratio further shows that trans-HC(O)SH is the most likely outcome. The formation of cis-HC(O)SH is an unlikely scenario that requires two subsequent unfavorable steps. However, it cannot be discounted, as that some small amount of cis-HC(O)SH is synthesized on grains.

Table 3. Reaction Outcome, Branching Ratios, and Reaction Energies on the Ice (ΔU_r,i, in kJ mol⁻¹, ZPE Corrected) and in the Gas Phase (ΔU_r,g, in kJ mol⁻¹, ZPE Corrected) for the OCHS + H Reaction

Outcome	Branching Ratio	ΔUr,i	ΔUr,g
trans-HC(O)SH	0.67 (0.54–0.78)	−350.7	–364.5
cis-HC(O)SH	0.33 (0.22–0.46)	−354.6	–361.0

Note. The numbers in the table are obtained using PW6B95-D4/def2-TZVP//PW6B95-D4/def2-SVP.

Download table as:  ASCIITypeset image

3.2.3. Cis-OCSH → Trans-OCSH

In Section 3.2.1 we determined that the reaction proceeds with diastereoselectivity owing to two factors. First, the formation of cis-OCSH is the sole hydrogenation product on the S atom for OCS, and, second, the particular binding arrangement between cis-OCSH and H₂O locks only one possible face for reaction (the one leading to trans-HC(O)SH). A legitimate question, however, is whether cis-OCSH may be able to isomerize to trans-OCSH on short timescales. In a positive case, the face available for the second hydrogenation reaction (assuming a similar binding arrangement with H₂O) would be different, and the isomerism of the whole process would surely change. Trans-OCSH is more stable than cis-OCSH by 8.5 kJ mol⁻¹ (ZPE corrected, estimated with (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP), so cis-OCSH should be, in principle, prone to isomerize under the right conditions. To determine if the isomerization of cis-OCSH is a plausible outcome, we have computed instanton rate constants in the gas phase using the implicit surface approach and (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP as the level of theory and instanton rate constants for the isomerization using a small, cis-OCSH–(H₂O)₂ model removing the water molecules not directly involved in the binding shown in Figure 6. The (U)CCSD(T)-F12a/cc-pVTZ-F12//PW6B95-D4/def2-TZVP method is significantly more expensive for the cis-OCSH–(H₂O)₂ model. In this particular case, we have left the explicit correlation treatment (F12) out. The level of theory for the model with two water molecules is thus (U)CCSD(T)/cc-pVTZ//PW6B95-D4/def2-TZVP. The transition state geometries, along with the torsion angle of the two different models, can be found in Figure 8.

Zoom In Zoom Out Reset image size

The reaction using the implicit surface model has an activation energy ΔU_A = 20.5 kJ mol⁻¹. The same reaction now explicitly including the two water clusters has an activation energy ΔU_A,cluster = 20.8 kJ mol⁻¹. Therefore the activation energy for the isomerization hardly depends on the binding situation, most likely due to a similar stabilization in the reactant and transition state on the surface. The slightly different levels of theory can have a small impact on the height of the barrier, too. What changes and deeply impacts the kinetics of the process is the frequency of the imaginary transition mode, which is 343.3i cm⁻¹ for the reaction in the gas phase, different from the value using the cluster model, of 213.9i cm⁻¹. ⁴ Such a subtle change dramatically changes the tunneling rate constants. It renders the reaction implausible when cis-OCSH interacts with water. The (reduced) instanton rate constant at 50 K, in this case, is 1.6 × 10⁻⁹ s⁻¹, in contrast to the instanton rate constant using the implicit surface model of 6.6 × 10⁻⁴ s⁻¹ (see Figure 9 for an Arrhenius plot of the rate constants at different temperatures). Translating these rate constants to half-lifetimes, we arrive at 13.6 yr for the former case and 17 minutes for the latter. While a half-lifetime of 13.6 yr is not extremely long on astronomical timescales, an H atom lands on interstellar dust grains approximately once per day (Wakelam et al. 2017a). The second hydrogenation reaction will occur much faster than the isomerization, finally confirming that the occurrence of trans-OCSH, and thus cis-HC(O)SH, on dust grains should be extremely rare.

Zoom In Zoom Out Reset image size

Finally, we dedicate some words to the possible trans-HC(O)SH → cis-HC(O)SH reaction after the formation of trans-HC(O)SH. Very briefly, such a reaction is not feasible on interstellar dust grains, owing to a high isomerization barrier ΔU_A = 35.9 kJ mol⁻¹, 1.8 times higher than the same isomerization for the cis-OCSH → trans-OCSH reaction. With such a barrier, an approximation of the rate constant using a symmetric Eckart correction yields 1.9 × 10⁻¹⁶ s⁻¹, rendering the process negligible under cold ISM conditions.

In this work, we determined that the sequential reaction OCS + 2H on interstellar dust grains primarily leads to trans-HC(O)SH, the only thioacid found in the ISM (Rodríguez-Almeida et al. 2021). In the absence of other competitive reactions, the cis-HC(O)SH isomer abundance must be small based on the observed specificity of the OCS + 2H reaction. Cis-HC(O)SH is the less stable isomer (cis-trans gap of 3.1 kJ mol⁻¹ (261 K), as mentioned in the Introduction) with the direct cis-trans isomerization barrier of 37.2 kJ mol⁻¹ (Kaur & Vikas 2014). Such obstacles (endothermicity and a moderate activation barrier) make it difficult to predict a significant abundance of cis-HC(O)SH in cold molecular clouds. This is the reason why in previous works isomerization transformations of compounds such as imines and other acids such as formic acid have been investigated by means of chemical processing (Shingledecker et al. 2019, 2020) or photoisomerization (Cuadrado et al. 2016). Alternatively, a multidimensional treatment of the quantum tunneling effects could explain the isomerization transformation in molecular clouds, as obtained for the E/Z isomers of imines (García de la Concepción et al. 2021). However, it remains unknown whether the same conclusions can be applied to the case of cis-/trans-HC(O)SH or other acids. This possibility will be explored in a forthcoming paper (J. García de la Concepción et al. 2021, in preparation).

Our results build on top of the studies unveiling the chemistry of sulfur in the ISM. Past and recent interstellar searches of saturated (Linke et al. 1979; Kolesniková et al. 2014; Rodríguez-Almeida et al. 2021) and unsaturated (Fuente et al. 2017; Martin-Drumel et al. 2019; Cernicharo et al. 2021; Rodríguez-Almeida et al. 2021) S-bearing molecules present a scenario with many mysteries in their formation routes, each one of them contributing to the understanding of the amount of elemental sulfur locked in organosulfur compounds (Laas & Caselli 2019). In a previous study, Lamberts (2018) showed that hydrogenation to saturation of CS is possible on astronomical timescales, indicating that highly saturated species (such as CH₃SH or C₂H₅SH) must be prevalent at later stages of a molecular cloud life. Based on these results, we also expect OCS to undergo significant hydrogenation in a mature molecular cloud on the same timescales when OCS is predicted to show its peak abundance on dust grains (Taquet et al.2020).

The values of our tunneling corrected rate constants, summed to the tentative branching ratios of reaction provided here, can be used to improve astrochemical models. Similarly, the data provided here complements recent efforts (Shingledecker et al. 2019, 2020; Zhang et al. 2020; García de la Concepción et al. 2021) to study stereochemistry in the ISM, fundamental in the context of prebiotic chemistry. The high diastereospecificity of this reaction makes trans-HC(O)SH a good case study for the study of isomerization processes in the ISM.

This work has important implications for a physicochemical end for other studies of chemistry on top of dust grains. We have found that the long-lived intermediate radical for the reaction cis-OCSH is present in a particular binding arrangement with the surface, which may facilitate subsequent reactions on the free carbon atom of such radical. In the future, we will explore the reactions between cis-OCSH with radicals different from H, as well as analyze the destruction channels of HC(O)SH with H or OH on ice surfaces via the Langmuir−Hinshelwood and the Eley–Rideal mechanisms. The particular nature of the cis-OCSH radical binding with water ice opens the gate to additional reactions with radicals more complex than H (e.g., CH₃, HCO), and to the advent of ever more sulfur-bearing COMs.

We want to emphasize the subtle effects of the explicit consideration of the interaction of reactants and surfaces in our calculations. In this study, we have found that the interaction of the different intermediates of reactions with water plays several subtle yet critical roles in the reaction: first, the abovementioned orientation effect of cis-OCSH with water, and, second, the reduction of the transition mode for reaction in the cis-OCSH → trans-OCSH reaction. Regarding the latter effect, it is easy to overlook its importance, since both an implicit surface treatment and the inclusion of a couple of water molecules yield very similar activation energies for the reaction. However, the reduction in the frequency of the transition mode has drastic effects, rendering reaction rate constants separated by almost six orders of magnitude and confirming that the reaction proceeds through a metastable intermediate cis-OCSH. In the absence of a surface, the cis-OCSH → trans-OCSH is a fast process under ISM conditions.

The detection of new sulfur-bearing compounds of prebiotic relevance in the ISM has proved to be an invaluable source of inspiration for experimental and theoretical studies seeking to understand the formation and chemical evolution of prebiotic molecules in space. Thioformic acid, HC(O)SH, is particularly puzzling, being the first thioacid ever detected in space, with a seemingly pure diastereomeric excess of its trans isomer. Both facts (the detection and the particular isomerism) are rationalized in the present work. There are still unknowns in the chemistry of acids in general and thioacids in particular that we will address in subsequent works, e.g., return of the acids to the gas phase, chemical interconversions after formation, or alternative mechanisms of isomerization.

Due to the length of the raw data, input scripts supporting these calculations, as well as the whole list of adsorption geometries, will be provided on reasonable request to the corresponding author.

The authors thank Prof. Dr. Johannes Kästner for useful discussions. G.M. acknowledges the support of the Alexander von Humboldt Foundation through a postdoctoral research grant. We would also like to acknowledge the support by the state of Baden-Württemberg through the bwHPC consortium for providing computer time and the German Research Foundation (DFG) through grant No. INST 40/575-1 FUGG (JUSTUS 2 cluster). J.G.C. and I.J.-S. acknowledge financial support from the Spanish State Research Agency (AEI) project numbers PID2019-105552RB-C41 and MDM-2017-0737 Unidad de Excelencia “María de Maeztu”- Centro de Astrobiología (CSIC-INTA). We also acknowledge support from the Spanish National Research Council (CSIC) through the i-Link project number LINKA20353.

Facility: bwHPC (Justus Cluster). -

Software: Turbomole v7.5 (Furche et al. 2014), GFN-FF (xTB Code; Grimme et al. 2017; Spicher & Grimme 2020), Chemshell (Sherwood et al. 2003; Metz et al. 2014), Molpro2015 (Werner et al. 2012, 2015), ASE (Hjorth Larsen et al. 2017).