Gas phase Elemental abundances in Molecular cloudS (GEMS)

VI. A sulphur journey across star-forming regions: study of thioformaldehyde emission

G. Esplugues
A. Fuente
D. Navarro-Almaida
M. Rodríguez-Baras
L. Majumdar
2, 3,
P. Caselli
V. Wakelam
E. Roueff
R. Bachiller
S. Spezzano
P. Rivière-Marichalar
R. Martín-Doménech
G. M. Muñoz Caro

Observatorio Astronómico Nacional (OAN), Alfonso XII, 3, 28014, Madrid, Spain

School of Earth and Planetary Sciences, National Institute of Science Education and Research, Jatni 752050, Odisha, India
Homi Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai 400094, India
Max-Planck-Institut für extraterrestrische Physik, 85748 Garching, Germany
Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, 92190 Meudon, France
Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA
Centro de Astrobiología, INTA-CSIC, Torrejón de Ardoz, E-28850 Madrid, Spain

In the context of the IRAM 30m Large Program GEMS, we present a study of (deuterated) thioformaldehyde in several starless cores located in star-forming filaments of Taurus, Perseus, and Orion.

We investigate the influence of the environmental conditions on the abundances of these molecules in the cores, and the effect of time evolution.

We have modelled the observed lines of H$_2$CS, HDCS, and D$_2$CS using the radiative transfer code RADEX. We have also used the chemical code Nautilus to model the evolution of these species depending on the characteristics of the starless cores.

We derive column densities and abundances for all the cores. We also derive deuterium fractionation ratios, $D$$_{\mathrm{frac}}$, which allow us to determine and compare the evolutionary stage between different parts of each star-forming region. Our results indicate that the north region of the B 213 filament in Taurus is more evolved than the south, while the north-eastern part of Perseus presents an earlier evolutionary stage than the south-western zone. Model results also show that $D$$_{\mathrm{frac}}$ decreases with the cosmic-ray ionisation rate, while it increases with density and with the degree of sulphur depletion. In particular, we can only reproduce the observations when the initial sulphur abundance in the starless cores is at least one order of magnitude lower than the solar elemental sulphur abundance.

The progressive increase in HDCS/H$_2$CS and D$_2$CS/H$_2$CS with time makes these ratios powerful tools for deriving the chemical evolutionary stage of starless cores. However, they cannot be used to derive the temperature of these regions, since both ratios present a similar evolution at two different temperature ranges ($\sim$7-11 K and $\sim$15-19 K). Regarding chemistry, (deuterated) thioformaldehyde is mainly formed through gas-phase reactions (double-replacement and neutral–neutral displacement reactions), while surface chemistry plays an important role as a destruction mechanism.

A sulphur journey across star-forming regions: study of thioformaldehyde emission
G. Esplugues et al.

§ INTRODUCTION

Molecular clouds are known to be the birthplaces of stars. Observations of different star-forming regions at several wavelengths suggest molecular clouds have complex morphologies, with the dust and gas arranged mostly along elongated and filamentary structures <cit.>. These filaments appear to be key structures that are required to reach the densities necessary for star formation according to Spitzer and Herschel observations, from sub-parsec scales in nearby star-forming regions [André et al., 2010] to tens-of-parsecs scales along spiral arms [Molinari et al., 2010]. In particular, filaments are where the initial conditions for star formation are set since they funnel interstellar gas and dust into increasingly denser concentrations that will contract and fragment, leading to gravitationally bound starless cores that will eventually form stars.

Dust and gas properties are key to determining the evolution of starless cores within filaments. In dense regions, grains are covered by icy mantles, which may affect grain coagulation <cit.>, and the presence of a dust size distribution, which in turn influences the grain emissivity. Moreover, grain coagulation may affect the charge balance in the gas phase, hence modifying the coupling of the gas with the magnetic field. In addition, grain surface chemistry plays an important role in the formation of the molecules that are key to the chemical network (e.g. H$_{{2}}$ and H$_{{2}}$O) and in the freeze-out process of molecules in the cold interiors of starless cores (e.g. CO).

One way to understand the dynamics of starless cores is through the study of molecule deuteration since the deuterium fraction, $D$$_{\mathrm{frac}}$, defined as the abundance ratio between a deuterated molecule and its hydrogenated counterpart, changes with their evolution <cit.>. In particular, $D$$_{\mathrm{frac}}$ is predicted to increase when a core evolves towards the onset of gravitational collapse as the core density profile becomes more and more centrally peaked [Crapsi et al., 2005]. Then, $D$$_{\mathrm{frac}}$ drops when the young stellar object (YSO) formed at the core centre begins to heat its surroundings <cit.>.

In molecular clouds ($T$$\sim$10 K), as neutral-neutral reactions often have activation barriers, the dominant reactions are those involving ions, H$^{+}_{3}$ being the first molecular ion formed as a product of the cosmic-ray (CR) ionisation of H$_2$ and H [Ceccarelli et al., 2014]. This ion is the main one that initiates the deuterium enrichment process through an exothermic reaction with HD <cit.>, and it forms H$_2$D$^+$, D$_2$H$^+$, and D$^{+}_{3}$. The collision of all these multi-deuterated forms of H$^{+}_{3}$ with neutral species produces deuterated molecules, such as N$_2$D$^+$ and DCO$^+$. On the other hand, their dissociative electronic recombination increases the D/H atomic ratio by several orders of magnitude with respect to the D cosmic abundance [Caselli et al., 2019], thus allowing the deuteration of molecules (e.g. methanol) on the surface of dust grains.
In cold environments, HD can also react with CH$^{+}_{3}$ and C$_2$H$^{+}_{2}$, leading to large enhancements of deuterated ions as well, such as CH$_2$D$^+$ and C$_2$HD$^+$, since the back reactions are inhibited [Millar, 2005]. All these newly formed deuterated ions react in turn with other molecules and atoms, transferring the D atoms to all the other species.

In pre-stellar cores <cit.>, CO and O freeze-out onto dust grains, which implies fewer destruction events for H$^{+}_{3}$ and its deuterated forms. This leads to an increase in H$_2$D$^+$, D$_2$H$^+$, and D$^{+}_{3}$ [Dalgarno & Lepp, 1984, Roberts et al., 2003, Walmsley et al., 2004], favouring deuteration processes in the presence of high-density gas <cit.>. Therefore, in dense and cold cores, the deuterated fraction is expected to be much higher than the average [D/H] interstellar abundance ratio <cit.>.

Evolution of starless cores can also be studied through molecule depletion. At temperatures of roughly 10 K and densities above 10$^{4}$ cm$^{-3}$, several molecules, such as CO and CS, condense out onto dust grain surfaces <cit.>, and the amount of depletion increases with time <cit.>. Therefore, depletion can be used as a time marker since evolved cores should be more depleted of certain species than younger cores. Apart from core evolution, the study of depletion is of paramount importance because of its effects on molecular clouds; depletion causes variations in the deuterium fractionation and the degree of ionisation [Dalgarno & Lepp, 1984, Caselli et al., 1998], which is one of the fundamental parameters regulating the star formation rate <cit.>. It also affects the gas-phase chemical composition, leading to chemical variations between similar types of clouds. In addition, the thermal balance of clouds may also be affected by the depletion of major gas coolants, such as CO <cit.>.

Although sulphur is one of the most abundant species in the Universe <cit.> and plays a crucial role in biological systems on Earth <cit.>, S-bearing molecules are not as abundant as expected in the interstellar medium (ISM). In fact, one needs to assume a significant sulphur depletion to reproduce observations not only in cold starless cores, but also in hot corinos and hot cores [Esplugues et al., 2014, Crockett et al., 2014]. In particular, it is thought that sulphur is depleted by a factor of up to 1000 compared to its estimated cosmic abundance [Ruffle et al., 1999, Wakelam et al., 2004]. Several studies have been carried out to shed light on the sulphur reservoir in molecular clouds and the sulphur depletion issue <cit.>. Chemical models predict that, in the dense ISM, atomic sulphur would stick on grains and be mostly hydrogenated to form H$_{{2}}$S [Hatchell et al., 1998, Garrod et al., 2007, Esplugues et al., 2014], which is thought to be the main sulphur reservoir in the ice [Vidal et al., 2017, Navarro-Almaida et al., 2020]. In fact, according to Druard & Wakelam, 2012, the lower the (gas and dust) temperature (<20 K), the greater the H$_{{2}}$S abundance on the grain surfaces. However, H$_{{2}}$S has never been detected in interstellar ices, and its abundance has been estimated to be smaller than 5$\times$10$^{-8}$ (/H) [Jiménez-Escobar & Muñoz Caro, 2011]. Other solid species have been proposed as possible sulphur reservoirs, such as OCS, SO$_{{2}}$, H$_{{2}}$S$_{{2}}$, CS$_{{2}}$, and S$_8$ <cit.>, OCS being the only S-bearing molecule unambiguously detected in ice mantles in the infrared [Geballe et al., 1985, Palumbo et al., 1995] along with, tentatively, SO$_{{2}}$ [Boogert et al., 1997, Zasowski et al., 2009].

Gas phase Elemental abundances in Molecular CloudS (GEMS) is an Institut de Radioastronomie Millimétric (IRAM) 30m Large Program, which aims to estimate S, C, N, and O depletions and the gas ionisation fraction, $X$(e$^{-}$), as a function of visual extinction in a selected set of prototypical star-forming filaments [Fuente et al., 2019]. To achieve this goal, it is first necessary to determine the abundances of the main reservoirs of the elements in the gas phase by selecting a sample of filaments located in several clouds, covering different types of star formation activity. Here we focus on a sample of starless cores of the nearby star-forming regions Taurus, Perseus, and Orion. These molecular cloud complexes were previously observed with $\it{Herschel}$ and SCUBA as part of the Gould Belt Survey [André et al., 2010], and accurate visual extinction $(A_{\mathrm{V}})$ and dust temperature $(T_{\mathrm{d}})$ maps are available [Malinen et al., 2012, Hatchell et al., 2005, Lombardi et al., 2014, Zari et al., 2016]. These regions each have different types of star formation activity and therefore different levels of external illumination, allowing us to investigate the influence of UV radiation on the gas composition. In particular, observing several starless cores within each filament will allow us to investigate the effect of time evolution on the chemistry of these cores, while comparing cores in different regions will let us explore the effect of the environment on the chemistry therein.

To carry out the present study, we selected the species H$_{{2}}$CS and its deuterated counterparts (HDCS and D$_{{2}}$CS). Several studies <cit.> confirm that H$_{{2}}$CS is expected to be involved in many grain-surface reactions and that it plays a key role in the synthesis of larger sulphur-bearing species, such as CH$_3$SH (in an analogous process to the sequential hydrogenation of CO that leads to CH$_3$OH through H$_2$CO). H$_2$CS can be directly formed from the hydrogenation of HCS, which, in turn, is formed through a neutral-neutral reaction between atomic carbon and H$_2$S
Recent results also demonstrate that H$_2$CS can be formed via the reaction between CH and H$_2$S [Doddipatla et al., 2020], as well as via the irradiation of CO:H$_2$S ice [Jiménez-Escobar & Muñoz Caro, 2011]. In any case, these reactions are limited by the available amount of H$_2$S within the ices. Nevertheless, as previously mentioned, H$_2$S is considered the main sulphur reservoir in ices, which would facilitate the formation of H$_2$CS. In this way, theoretical studies <cit.> also conclude that, at high densities, H$_2$CS, in addition to being present at a significant level in the gas phase, might also be an abundant C-bearing sulphur species in the ice. All this makes H$_{{2}}$CS and its deuterated counterparts good candidates for characterising starless cores and studying their evolution.

In this paper the observations of H$_2$CS, HDCS, and D$_2$CS are described in Sect. <ref>, and the considered source sample in Sect. <ref>. In Sect. <ref> we present the data and use the non-local thermodynamic equilibrium radiative transfer code RADEX to derive column densities and abundances. A discussion about fractional abundance differences between the sources of the sample is presented in Sect. <ref>. In that section, we also use the Nautilus time-dependent chemical code to analyse the deuterium fraction evolution and to study the CR impact on this fraction. Section <ref> also provides a comparison between theoretical and observational results, as well as an analysis of the main chemical formation and destruction routes of (deuterated) thioformaldehyde. Finally, we summarise our conclusions in Sect. <ref>.

Cores included in the GEMS sample and the observation cuts associated with them and shown in Figs. <ref>-<ref>.

Region     Cloud     Core     2cCoordinates      Cut

ID         RA (J2000)          Dec (J2000)

B 213$^1$       #1       04:17:41.8           $+$28:08:47.0             C1

#2        04:17:50.6           $+$27:56:01.0           C2

#5       04:18:03.8           $+$28:23:06.0             C5

Taurus          #6       04:18:08.4           $+$28:05:12.0            C6

#7       04:18:11.5           $+$27:35:15.0           C7

#10      04:19:37.6           $+$27:15:31.0           C10

#12      04:19:51.7           $+$27:11:33.0            C12

#16      04:21:21.0           $+$27:00:09.0           C16

L1448$^2$        #32      03:25:49.0           $+$30:42:24.6                C1

NGC 1333$^2$                #46      03:29:11.0            $+$31:18:27.4             C3-1

#60     03:28:39.4             $+$31:18:27.4             C3-14

#51      03:29:08.8            $+$31:15:18.1                C4

Perseus            #53     03:29:04.5           $+$31:20:59.1              C5

#57      03:29:18.2              $+$31:25:10.8             C6

#64      03:29:25.5              $+$31:28:18.1             C7

Barnard 5$^2$        #79      03:47:38.9             $+$32:52:15.0              C1

IC 348$^2$          #1         03:44:01.0          $+$32:01:54.8          C1

#10      03:44:05.7             $+$32:01:53.5              C10

Orion A                     05:35:19.5            $-$05:00:41.5               C1

Orion                               05:35:08.1             $-$05:35:41.5               C2

05:35:23.6          $-$05:12:31.8              C3

$^1$Taurus core IDs are from [Hacar et al., 2013].

$^2$Perseus core IDs are from [Hatchell et al., 2007].

§ OBSERVATIONS

The  data  used  in  this  work  are  taken  from  the  GEMS project. A detailed description of the observations is given by Fuente et al., 2019 and Rodríguez-Baras et al., 2021. For clarity, the main observational parameters are summarised below. The 3 mm and 2 mm observations (covering the frequency range 85-172 GHz) were carried out using the IRAM 30m telescope at Pico Veleta (Spain) during three observing periods in July 2017, August 2017, and February 2018. The observational parameters are listed in Table <ref> (see the appendix) with the beam size varying with the frequency as HPBW($\arcsec$)=2460/$\nu$ where $\nu$ is in GHz. The observing mode was frequency switching with a frequency throw of 6 MHz well adapted to removing standing waves between the secondary mirror and the receivers. The Eight MIxer Receivers (EMIR) and the fast Fourier transform spectrometers with a spectral resolution of 49 kHz (equivalent to 0.165 km s$^{-1}$ and 0.087 km s$^{-1}$ for frequencies at 89 GHz and 168 GHz, respectively) were used for these observations. The achieved rms was $\sim$10-20 mK for $\nu$<150 GHz and $\sim$20-30 mK for $\nu$>150 GHz.

The intensity scale is $T$$_{\mathrm{MB}}$, which is related to $T$$^{\star}_{\mathrm{A}}$ by

\begin{equation}
\end{equation}

where $F$$_{\mathrm{eff}}$ is the telescope forward efficiency and $B$$_{\mathrm{eff}}$ is the main beam efficiency[http://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies]. The difference between $T$$^{\star}_{\mathrm{A}}$ and $T$$_{\mathrm{MB}}$ is $\sim$17$\%$ at 86 GHz and $\sim$27$\%$ at 145 GHz.

Calibration errors are estimated to be $\sim$10$\%$. The data were reduced and processed using the CLASS and GREG packages from IRAM GILDAS software[http://www.iram.fr/IRAMFR/GILDAS], developed by IRAM.

§ SOURCE SAMPLE

GEMS considers the molecular cloud complexes Taurus, Perseus, and Orion. In particular, observations include several cuts roughly perpendicular to selected filaments, where the separation between one position and another in a given cut is selected
to sample the visual extinction range in regular intervals of $A_{\mathrm{V}}$. In this project, we only focus on the starless cores within filaments of these regions (see Figs. <ref>, <ref>, and <ref>). In the following, we describe the observed positions
in more detail.

The Taurus (L 1495/B 213) molecular cloud, with a total mass of about 1.5$\times$10$^4$ M$_{\sun}$ derived from CO data analysis [Pineda et al., 2010], is one of the closest regions of star formation, at a distance of about 145 pc [Qian et al., 2015, Yan et al., 2019], and is known to contain more than 250 YSOs. Taurus is considered an archetypal low-mass star-forming region. Multiple studies, using different ground-based and space telescopes, have been carried out to analyse its evolution and structure. For instance, Palla & Stahler, 2002 found a young population inside the filaments and a more dispersed and older population outside them, concluding the presence of an age spread in the region. Goldsmith et al., 2008 showed a very complex and highly structured cloud morphology including filaments, cavities, and rings, while Goodman et al., 1992 studied the presence and orientation of magnetic fields with respect to several filaments within the Taurus cloud. We also highlight the $\it{Herschel}$ observations of Taurus obtained in the context of the Gould Belt survey [Kirk et al., 2013, Palmeirim et al., 2013], where a large-scale continuum map of this complex was obtained [Schmalzl et al., 2010], and the recent astrometric studies carried out with Gaia <cit.>.

B 213 molecular hydrogen column density maps as derived by Palmeirim et al., 2013, reconstructed at an angular resolution of 18.2$\arcsec$.
A general view of the region is shown in the top-right panel, and main regions of interest are enlarged. Contours are (3, 6, 9,12, 15, 20, and 25)$\times$10$^{21}$ cm$^{-2}$. Positions observed by GEMS with the 30m telescope are indicated with white triangles. Red circles represent the position of the starless cores. Labels in red indicate the cut IDs. See Table <ref> for further details.

Perseus filament (from right to left and top to bottom: NGC 1333, IC 348, L 1448, and B5) dust opacity maps at 850 µm by Zari et al., 2016, convolved at an angular resolution of 36$\arcsec$. Contours are (0.056, 0.13, 0.24, 0.56, 1.01, and 1.6) $\times$10$^{21}$ cm$^{-2}$, which, according to expression (7) from Zari et al., 2016, corresponds to visual extinctions of $\sim$5, 7.5, 10, 15, 20, and 25 mag, respectively. Positions observed with the 30m telescope are indicated with white triangles. Red circles represent the positions of the starless cores.

Orion dust opacity map at 850 µm by Lombardi et al. (2014), convolved at an angular resolution of 36$\arcsec$. Contours are (0.056, 0.24, 0.56, 1.36, and 1.61)$\times$10$^{21}$ cm$^{-2}$, which, according to Lombardi et al., 2014, correspond to visual extinctions of $\sim$1.3, 5.6, 13.2, 23.8, and 38 mag. Positions observed with the 30m telescope are indicated with white triangles. Red circles represent the positions of the starless cores. Labels in red indicate the cut IDs (see Table <ref>).

One of the main filaments in Taurus is known as the Lynds Dark Nebula 1495 (L1495). L1495 contains several Barnard Dark Nebulae, which are dust-filled regions. Dark nebulae are extremely dense regions of dust that obscure visible light. The central region is known as B10, with B211 and B213 stretching out from the centre. The L 1495/B 213 nebula is a clear example of a star-forming region where the magnetic field lines are perpendicular to the main filament [Soler, 2019]. It has been extensively studied by Hacar et al., 2013 in high-density tracers, such as C$^{18}$O, N$_2$H$^+$, and SO, and by Hacar et al., 2016 in the three main isotopologues of $^{12}$CO, $^{13}$CO, and C$^{18}$O, deducing the presence of several dense cores embedded in B213 [Benson & Myers, 1989, Onishi et al., 2002, Tatematsu et al., 2004, Punanova et al., 2018]. Some of these dense cores are starless, while others are associated with YSOs of different ages, with a density of stars decreasing from north to south [Davis et al., 2010]. This suggests a different dynamical and chemical age along the filament.
Here, we consider the eight (Table <ref>, Fig. <ref>) starless cores #1, #2, #5, #6, #7, #10, #12, and #16 (core numbers from the catalogue of Hacar et al. 2013).

The Perseus molecular cloud (L 1448, NGC 1333, Barnard 5, IC 348) is one of the most active nearby star-forming regions, extending $\sim$10 pc on the sky ([Bachiller & Cernicharo, 1986]). According to recent Gaia parallaxes and photometric data, it is located at a distance ranging from 234 to 331 pc ([Zucker et al., 2020]). The infrared survey of [Ladd et al., 1993] suggests that Perseus is intermediate in its star-forming properties between the Taurus complex, with stars predominantly forming in relative isolation (YSO surface density $\sim$10 pc$^{-2}$), and the Orion complex, with large clusters with densities of 100 pc$^{-2}$ as the main stage for the star formation.

The Perseus complex contains six regions with star-formation activity -L1448, L1455, NGC 1333, Barnard 1 (B1), IC 348, and Barnard 5 (B5)- [Bachiller & Cernicharo, 1986, Ladd et al., 1994] as well as more than 400 YSOs, and $\sim$100 dense cores [Hatchell et al., 2005, Dunham et al., 2013]. Most of the protostars in Perseus are associated with IC 348 and NGC 1333 [Yang et al., 2021], with the latter having many active outflows that may regulate the ongoing star formation [Bachiller & Cernicharo, 1990, Knee & Sandell, 2000, Davis et al., 2008, Curtis et al., 2010, Plunkett et al., 2013]. B1, with a high deuterium enrichment [Marcelino et al., 2005], has been studied in different wavelengths <cit.>, revealing the presence of 10 protostars, multiple molecular outflows, and Herbig-Haro objects. All this suggests that the surrounding region of B1 is very actively forming stars. An energetic outflow can also be found in B5 [Langer et al., 1996] driving into the surrounding cloud material. B5 presents a morphology, with multiple filaments within the velocity coherent region of the core [Pineda et al., 2011]. The starless cores from Perseus considered in this paper are listed in Table <ref> (see also Fig. <ref>). All the cores were classified as starless cores by Hatchell et al., 2007, except #60 (1333-C3-14) that is a young Class 0 object.

The Orion molecular cloud (Orion A), located at $\sim$428 pc [Zucker et al., 2019], is the nearest region with a presence of recent OB star formation <cit.>. From CO maps [Maddalena et al., 1986], two big molecular regions are distinguished, the north molecular complex (Orion B) and the south molecular complex (Orion A). Orion B is an extensive CO emission region associated with the dark cloud L 1630 [Lynds, 1962]. This emission is spread over a region of 4 degrees in the north-south direction from the Horsehead nebula, NGC 2023, and NGC 2024 to the reflection nebulae NGC 2068 and NGC 2071. On the other hand, Orion A is associated with the dark clouds L 1640, L 1641, and L 1647, with emission extending 6 degrees south from the Orion Nebula. The molecular emission of this southern complex is subdivided into three clouds of neutral material: OMC 1, OMC 2, and OMC 3. OMC 1, which is located behind the HII M 42 region, is identified as a dense gas directly associated with Orion KL [Wilson et al., 1970]. OMC 2 is related to the HII M 43 region [Gatley et al., 1974], and OMC 3 is located approximately 16$\arcmin$ north of OMC 2 [Kutner et al., 1976]. These clouds structured as an integral-shaped filament of molecular gas [Bally et al., 1987] present concentrations of sub-millimetre continuum emission in the southern part of the filament, which are referred to as OMC 4 [Johnstone & Bally, 1999] and OMC 5 [Johnstone & Bally, 2006].

In GEMS, we have considered three cuts, along OMC 2 (ORI-C3), OMC 3 (ORI-C1), and
OMC 4 (ORI-C2) (see Table <ref> and Fig. <ref>). The cuts avoid protostars and stars in these active star-forming regions and represent distinct environments because of their different distances from the Orion nebula. We selected the visual extinction peaks along each cut for this study.

§ DATA ANALYSIS AND RESULTS

Interstellar thioformaldehyde (H$_2$CS) is a slightly asymmetric rotor with two interchangeable hydrogen nuclei; therefore, its rotational levels are grouped into ortho ($K_{\mathrm{a}}$ odd) and para ($K_{\mathrm{a}}$ even), with statistical weights 3:1. The ortho ground state lies 14.9 K above the para ground state, and its dipole moment is µ$_{\mathrm{a}}$=1.647 D [Fabricant et al., 1977]. An advantage of the spectrum of a slightly asymmetric rotor molecule is that lines arising between levels of different energies lie nearby in the spectrum.

H$_2$CS was first detected in the observation of the $K$-doublet 2$_{1}$$_{,}$$_{1}$-2$_{1}$$_{,}$$_{2}$ transition at 3 GHz in absorption towards Sgr B2 [Sinclair et al., 1973]. Subsequently, this molecule has been observed towards several interstellar sources, such as interstellar clouds <cit.>, including cold dark clouds (e.g. [Irvine et al., 1989, Vastel et al., 2018]), in circumstellar envelopes <cit.>, in presence of shocked gas <cit.>, and in regions dominated by UV photons <cit.>.

Regarding the deuterated versions of thioformaldehyde, HDCS was first detected in TMC-1 through a spectral survey [Ohishi & Kaifu, 1998, Kaifu et al., 2004], using the frequencies calculated in the laboratory by Minowa et al., 1997.
Twice deuterated thioformaldehyde (D$_2$CS) was observed for the first time in a radio astronomical source (the dark cloud Barnard 1) by Marcelino et al., 2005.

§.§ Line profiles

Figures <ref>-<ref> (see the appendix) show the lines of H$_2$CS, HDCS, and D$_2$CS observed with the IRAM 30m telescope. We have detected two transitions (3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$ and 4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$) of o-H$_2$CS, two transitions (4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$ and 5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$) of p-H$_2$CS, two transitions (3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$ and 3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$) of HDCS, one transition (3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$) of o-D$_2$CS, and one transition (5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$) of p-D$_2$CS (Tables <ref>-<ref>). These transitions span an energy range of $E_{\mathrm{up}}$=8.1-24.7 K.

We first fitted the observed lines with Gaussian profiles using the CLASS software to derive the radial velocity ($V_{\mathrm{LSR}}$), the line-width, and the intensity for each line. Results are shown in Tables <ref>-<ref>. The contribution to the intensity arises from one narrow velocity component in all the cases. Line widths for H$_{2}$CS, HDCS, and D$_{2}$CS vary between 0.18 and $\sim$1.20 km s$^{-1}$, with the widest line profiles found in the starless core 1333-C4-1 (Perseus).
All the detected lines present $T_{\mathrm{MB}}$$\geq$35 mK, with H$_2$CS having $T_{\mathrm{MB}}$$<$0.8 K, $T_{\mathrm{MB}}$$<$0.5 K for HDCS, and D$_2$CS with $T_{\mathrm{MB}}$$<$0.3 K.

§.§ Radiative transfer code

In order to derive column densities and abundances, we carried out a more advanced analysis of the emission of H$_2$CS, HDCS, and D$_2$CS using the molecular excitation and radiative transfer code RADEX (van der Tak et al. 2007). RADEX is a one-dimensional non-local thermodynamic equilibrium radiative transfer code, that uses the escape probability formulation assuming an isothermal and homogeneous medium without large-scale velocity fields.

The collisional rates have been scaled from the (ortho- and para-) H$_2$CO rates derived by Wiesenfeld & Faure, 2013, which were calculated including energy levels up to about 180 cm$^{-1}$ for collisions with H$_2$. For deuterated thioformaldehyde, we also used the same collision rates as for H$_2$CS, with a correction for the reduced mass of
HDCS. To obtain the fit that better reproduces the observed line profiles, we let temperature ($T_{\mathrm{k}}$) and gas density ($n_{\mathrm{H_2}}$) vary as free parameters. However, several good fits can be obtained for different combinations of $n_{\mathrm{H_2}}$ and $T_{\mathrm{k}}$. In order to avoid this degeneracy, we let the gas temperature to vary in a small range around the dust temperature[Obtained by
Palmeirim] et al., 2013, Lombardi et al., 2014, and Zari et al., 2016 on the basis of the Herschel Gould Belt Survey [André et al., 2010] and Planck data [Bernard et al., 2010]., $T$$_{\mathrm{K}}$ = $T$$_{\mathrm{d}}$$\pm$5 K, and let the density and the H$_2$CS column density to vary to constrain accurately the physical conditions and reproduce the observed intensities. The best fit model is obtained by finding the minimum root mean square (rms) value of log$_{10}$($I_{\mathrm{obs}}$/$I_{\mathrm{mod}}$), following Neufeld et al., 2014). This is defined as

\begin{equation}
{\mathrm{rms}}= \displaystyle{ \sqrt { {\frac{1}{n}} \displaystyle\sum_{i=1}^{n} (log_{10} {\frac{ I^{i}_{\mathrm{obs}}}{ I^{i}_{\mathrm{mod}} }}}})^2\,,
\end{equation}

where $n$ is the number of observed lines, $I$$^{i}_{\mathrm{obs}}$ is the observed line intensity calculated from Gaussian fits (Sect. <ref> and Tables <ref>-<ref>), and $I$$^{i}_{\mathrm{mod}}$ is the model line intensity using RADEX. The values obtained for the H$_2$ volume density of the cores in the best fit models are in the same order of magnitude that the ones obtained by Rodríguez-Baras et al., 2021 using the species CS to derive the densities. In particular, the differences between both cases are within a factor of $\sim$3. Column densities for HDCS and D$_2$CS were derived assuming physical conditions obtained for H$_2$CS in each core.

Apart from calibration, given that lines are optically thin, one of the main sources of uncertainty to consider is the low angular resolution of the telescope (between $\sim$14$\arcsec$ and $\sim$29$\arcsec$ depending on frequency), which implies that the emission from the inner region of the pre-stellar core is blended with the outer cold envelope. This mainly affects Orion and Perseus cores at low frequencies. Also, the possible volume density gradients along the line of sight may influence the results as well. Taking this into account, the uncertainties in the obtained column densities are estimated to be 20$\%$ for H$_2$CS, and 25$\%$ for HDCS and D$_2$CS (higher uncertainty for the deuterated versions of H$_2$CS because their lines are weaker). The derived column densities for H$_2$CS, HDCS, and D$_2$CS are shown in Table <ref>.

We also calculated the fractional abundances, $X$, of H$_2$CS, HDCS, and D$_2$CS with respect to total H nuclei using visual extinction data for the cores [Rodríguez-Baras et al., 2021] and the relation between extinction and hydrogen column density in the Galaxy [Güver & Özel, 2009][$N$$_{\mathrm{H}}$ (cm$^{-2}$) = (2.21$\pm$0.09)$\times$10$^{21}$ $A$$_{\mathrm{V}}$ (mag)].
The visual extinction data were calculated for B 213 from the $N$$_{\mathrm{H_2}}$ values provided by Palmeirim et al., 2013. In the Orion case, we used dust opacity maps at 850 $\mu$m [Lombardi et al., 2014] and the expression $A$$_{\mathrm{V}}$=$A$$_{\mathrm{K}}$/0.112, where the $K$-band extinction $A$$_{\mathrm{K}}$=2640$\times$$\tau$$_{850}$+0.012 (see Rodríguez-Baras et al., 2021 for more details). In Perseus, we also used the dust opacity at 850$\mu$m ($\tau$$_{850}$) and expression (7) from Zari et al., 2016.
We have therefore assumed an uncertainty of 30$\%$ for $X$(H$_2$CS) and of 35$\%$ for $X$(HDCS) and $X$(D$_2$CS) due to the uncertainty introduced by the visual extinction values. Fractional abundance results are listed in Table <ref>.

§.§.§ Thioformaldehyde

Figures <ref>-<ref> (in the appendix) show our best fit models for the detected rotational lines of o-H$_2$CS and p-H$_2$CS in B 213, L 1448, NGC 1333, Barnard 5, IC 348, and Orion A, respectively.

For the cores located in the Taurus complex, we find the highest H$_2$CS column density, $N$(H$_{2}$CS)=(7$\pm$1)$\times$10$^{12}$ cm$^{-2}$ (Table <ref>), in the core C5. This core can be found in a region located to the north of the filament L 1495/B 213 (Fig. <ref>), which is characterised by the presence of Class I/Flat objects, but whose stellar population is dominated by more evolved objects. In fact, the region hosting C5 has the largest number of Class III objects <cit.>. The starless core C16, located towards the south-west of B 213, presents  the second highest H$_2$CS column density (Table <ref>). This core is in a region that has recently started to form stars and represents a younger active region compared to C5. Other cores of the sample, such as C1, C2, and C6, which are located north of B 213, although slightly more south from C5 (see maps shown in Figs. 1 and 2 in Hacar et al. 2013), present H$_2$CS column densities up to three times lower than C5.

Regarding Perseus, the lowest H$_2$CS column densities ($<$2$\times$10$^{12}$ cm$^{-2}$) are found in the cores C5, C6, and C7 of NGC 1333, and in IC 348, which is located in the eastern part of the Perseus molecular cloud and associated with a cluster containing a pre-main-sequence star with an age of 0.5-3.5 Myr [Luhman et al., 2003]. By contrast, the highest value of $N$(H$_{2}$CS), (11$\pm$2)$\times$10$^{12}$ cm$^{-2}$, is found in L 1448 situated to the west of the complex.
The core 79-C1-1 (in Barnard 5, north-east of Perseus) presents an intermediate H$_2$CS column density (5$\pm$1$\times$10$^{12}$ cm$^{-2}$).

Fractional abundances for H$_{2}$CS, HDCS, and D$_{2}$CS.

Region      Core     $N$(H$_{2}$CS)/$N$$_{\mathrm{H}}$     $N$(HDCS)/$N$$_{\mathrm{H}}$     $N$(D$_{2}$CS)/$N$$_{\mathrm{H}}$
$\times$10$^{-11}$     $\times$10$^{-11}$     $\times$10$^{-11}$

B 213-C1-1               7$\pm$2                3.0$\pm$1.0             1.8$\pm$0.6

B 213-C2-1               5$\pm$2                1.6$\pm$0.6             1.3$\pm$0.4

B 213-C5-1               14$\pm$4               1.6$\pm$0.6             0.7$\pm$0.3

Taurus     B 213-C6-1      9$\pm$3                2.4$\pm$0.9             2.3$\pm$0.8

B 213-C7-1               10$\pm$3               2.2$\pm$0.8             1.2$\pm$0.4

B 213-C10-1              4$\pm$1                0.5$\pm$0.2             $<$0.6$\pm$0.2

B 213-C12-1              1.3$\pm$0.4            $<$0.4$\pm$0.2          0.6$\pm$0.2

B 213-C16-1              10$\pm$3               1.4$\pm$0.5             1.1$\pm$0.4

L 1448-1                 18$\pm$5               6$\pm$2                 4$\pm$1

1333-C3-1                 4$\pm$1                0.4$\pm$0.1             $<$0.10$\pm$0.04

1333-C4-1                 13$\pm$4               4$\pm$1                 4$\pm$1

1333-C5-1                 1.6$\pm$0.5            $<$0.12$\pm$0.04        $<$0.25$\pm$0.09

Perseus     1333-C6-1       1.4$\pm$0.4     1.1$\pm$0.4                    $<$0.11$\pm$0.04

1333-C3-14                12$\pm$4               1.8$\pm$0.6             1.7$\pm$0.6

1333-C7-1                 5$\pm$1                2.5$\pm$0.9             1.0$\pm$0.4

79-C1-1                   12$\pm$4               3$\pm$1                 1.9$\pm$0.7

IC 348-1                 2.4$\pm$0.7            $<$0.12$\pm$0.04        $<$0.4$\pm$0.1

IC 348-10                2.6$\pm$0.8          $<$0.10$\pm$0.04        $<$0.3$\pm$0.1

ORI-C1-2                  1.3$\pm$0.4            $<$0.03$\pm$0.01        $<$0.02$\pm$0.01

Orion     ORI-C2-3          3.0$\pm$0.9            0.3$\pm$0.1             $<$0.08$\pm$0.03

ORI-C3-1                  0.8$\pm$0.3            $<$0.05$\pm$0.02        $<$0.04$\pm$0.01

In Orion A, we obtain $N$(H$_2$CS)=(2.5-5)$\times$10$^{12}$ cm$^{-2}$ for the three cores (Table <ref>), which are very similar column densities considering their uncertainty values (Table <ref>). Nevertheless, the highest value is found in ORI-C1 <cit.>.

In our analysis, we treated o-H$_2$CS and p-H$_2$CS separately, which allows us to derive the ortho-to-para ratio (OPR) for the different regions. We found that the values of the OPR(H$_2$CS) range from 1.8 to 5.0 with an average value of 2.4$\pm$0.9, without any clear trend with the gas temperature and/or the environment when taking the uncertainties into account.

§.§.§ Deuterated thioformaldehyde

Figures <ref>-<ref> (in the ppendix) show our best fit models for the detected rotational lines of HDCS, o-D$_2$CS, and p-D$_2$CS in the sample of starless cores. As for H$_2$CS, we observe that, although most of the line profiles are very well reproduced by the model, there are a few cases where the best fits underestimate the observed spectra. The presence of residuals may be mainly due to the noise added to the data, which is especially important in the cases where the line intensities are very weak ($<$0.1 K). This is, for instance, the case for the transition 3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$ of HDCS in B213-C2-1, B213-C5-1, and B213-C7-1. The underestimation of the observed spectra by the fits can be also due to the presence of blended species, or simply due to a limited model, since, depending on the physical structure of each source, emission can arise from different possible source components that are not considered in the model.

In the Taurus complex, we detected the deuterated species HDCS in all the starless cores of B 213, except in B 213-C12 (located in the south part of L 1495/B 213, Fig. <ref>) for which we provide an upper limit of $N$(HDCS)$<$(2.2$\pm$0.5)$\times$10$^{11}$ cm$^{-2}$. In fact, the lowest column densities of HDCS, (2.1-8)$\times$10$^{11}$ cm$^{-2}$ (Table <ref>), are also found in this central part of the B 213 filament, in particular in B 213-C10 and B 213-C16. By contrast, the highest values are found in B 213-C1, C6, and C7, located all of them in the north of the filament, suggesting a higher deuterium fraction in the north of this complex than in the south region.

Similar results are found for D$_2$CS, with the lowest column density values (3.0$\times$10$^{11}$ cm$^{-2}$) in the south part of the B 213 filament (core C-12), and the highest ones in cores (C1 and C6) of the north. In particular, we detect o-D$_2$CS in all the cores, except in one (B 213-C10) for which we provide an upper limit. On the other hand, p-D$_2$CS is only detected in B 213-C1 with a column density of (4$\pm$1)$\times$10$^{11}$ cm$^{-2}$. For the rest of the cores, we calculated the column density of p-D$_2$CS from o-D$_2$CS and assuming OPR=2.0$\pm$1.0. This OPR(D$_2$CS) is obtained by averaging the OPRs from the cores B 213-C1-1 and L 1448-1, since only in these two cores we observe o-D$_2$CS and p-D$_2$CS. If there is no detection of o-D$_2$CS, we just provide an upper limit for p-D$_2$CS (as for B 213-C10) in Table <ref>.

In the Perseus complex, we find the highest HDCS column density, (3.5$\pm$0.9)$\times$10$^{12}$ cm$^{-2}$, in L 1448 and NGC 1333-C4 like in the case of H$_2$CS. Non-detection of HDCS nor D$_2$CS is, however, found in IC 348, which is immersed in a more active environment compared to the one where L 1448 is located [Knee & Sandell, 2000, Plunkett et al., 2013].
Values for $N$(HDCS) in 79-C1 and in most of the starless cores of NGC 1333 are found to be lower than in L 1448, with $N$(HDCS)$<$3$\times$10$^{12}$ cm$^{-2}$.
Regarding D$_2$CS, this double deuterated species is only detected in its two versions (o-D$_2$CS and p-D$_2$CS) in L 1448, while in the rest of the cores we provide upper limits or calculate the column density of p-D$_2$CS where o-D$_2$CS has been detected by assuming an OPR=2.0$\pm$1.0 as previously mentioned.

In the three starless cores observed in Orion A, we only detect HDCS, with $N$=(3.7$\pm$0.9)$\times$10$^{11}$ cm$^{-2}$, in ORI-C2, which is located in OMC-4. The OMC-4 region is less luminous and turbulent than OMC-2 and OMC-3 (where the cores ORI-C3 and ORI-C1 are located, respectively), and it does not show evidence of outflow activity [Johnstone & Bally, 1999]. It suggests that OMC-4 might be in a pre-collapse phase of protostar evolution with colder temperatures, explaining thus the presence of deuterated thioformaldehyde. Regarding doubly deuterated thioformaldehyde, we do not detect its presence in any of the sampled starless cores, but we also provide upper limits in Table <ref>.

Fractional abundances of H$_2$CS, HDCS, and D$_2$CS with respect to total H nuclei for the core sample.

§ DISCUSSION

§.§ Fractional abundances

Figure <ref> shows the H$_2$CS, HDCS, and D$_2$CS fractional abundances, $X$, with respect to total H nuclei in the core sample (see also Table <ref>).

In the Taurus cloud, we find the highest H$_2$CS fractional abundance (1.4$\times$10$^{-10}$) in the starless core B 213-C5. The rest of cores present fractional abundances up to ten times lower than C5. As previously mentioned, the region hosting C5 is characterised by the presence of several YSOs, some of them being Class III [Hacar et al., 2013]. This implies that the core B 213-C5 may be affected by nearby star formation activity and, therefore, by the presence of winds and outflows that could excavate molecular material, leading to a more efficient penetration of the interstellar radiation field into the molecular cloud. The presence of low-mass stars nearby can also accelerate energetic particles and increase the local CR flux and ionisation rate <cit.>, leading to an increase in the H$_2$CS abundance as also found in our model results (see Fig. <ref>).
For deuterated thioformaldehyde, we find fractional abundances for HDCS and D$_2$CS in Taurus in the range (0.5-3)$\times$10$^{-11}$ (Table <ref>). In Fig. <ref>, we observe that the highest deuterium fractional abundances are found in the cores located in the north of B 213 (e.g. C1, C6, and C7), while in cores located in the south part of the filament (C10, C12, and C16), $X$(HDCS) and $X$(D$_2$CS) are lower than in the north or simply not detected.

In Perseus complex, we find $X$(H$_2$CS) in the range (1.4-18)$\times$10$^{-11}$, similar to the range found in Taurus, with the highest value in L 1448, which hosts a few young stars associated with extremely long outflows <cit.>. By contrast, some of the lowest $X$(H$_2$CS) values are found in the starless cores of IC 348, a young region where only $<$0.05$\%$ of its mass is contained in outflows [Curtis et al., 2010].
The highest $X$(HDCS) value in the Perseus sample is also found in L 1448, while we do not detect it in IC 348, and the values obtained for most of the cores of NGC 1333 <cit.> are low (0.4-2.5$\times$10$^{-11}$). Regarding D$_2$CS, its obtained fractional abundances are between 1-4$\times$10$^{-11}$, with one of the highest values also found in L 1448.

In Orion A, the highest $X$(H$_2$CS) is (3.0$\pm$0.9)$\times$10$^{-11}$, and it is found in ORI-C2-3. This value is more than double those found in the other cores of the Orion sample. The fractional abundance of HDCS in ORI-C2-3 is 0.3$\pm$0.1)$\times$10$^{-11}$, which is ten times lower than $X$(H$_2$CS).

§.§ Deuterium fractionation

Deuterium was formed at the birth of the Universe with an abundance D/H estimated to be 1.6$\times$10$^{-5}$ and is destroyed in the interiors of stars at $T$$\sim$0.5$\times$10$^6$ K [Molinari & Ricci, 1986, Tsujimoto & Bland-Hawthorn, 2010, Moscoso et al., 2021]. The deuterium fractionation ratio, defined as the ratio of the column density of a deuterated molecule to its hydrogen counterpart, is found to increase greatly in some sources. This is, for instance, the case of cold ($T$$_{\mathrm{g}}$$\sim$10 K) dark clouds, where the increase is a few orders of magnitude compared to the cosmic D/H ratio <cit.>.

Deuterium fractionation ratios.

Region      Core     [HDCS]/[H$_2$CS]     [D$_2$CS]/[H$_2$CS]

B 213-C1-1               0.5$\pm$0.2            0.3$\pm$0.1

B 213-C2-1               0.3$\pm$0.1            0.2$\pm$0.1

B 213-C5-1               0.12$\pm$0.05          0.05$\pm$0.02

Taurus     B 213-C6-1      0.3$\pm$0.1            0.3$\pm$0.1

B 213-C7-1               0.2$\pm$0.1            0.11$\pm$0.05

B 213-C10-1              0.13$\pm$0.06          $<$0.16$\pm$0.07

B 213-C12-1              $<$0.3$\pm$0.2         0.5$\pm$0.2

B 213-C16-1              0.14$\pm$0.06          0.11$\pm$0.05

L 1448-1                 0.3$\pm$0.1            0.2$\pm$0.1

1333-C3-1                 0.10$\pm$0.05          $<$0.03$\pm$0.01

1333-C4-1                 0.3$\pm$0.1            0.3$\pm$0.1

1333-C5-1                 $<$0.08$\pm$0.04       $<$0.15$\pm$0.07

Perseus     1333-C6-1       0.8$\pm$0.4            $<$0.08$\pm$0.04

1333-C3-14                0.15$\pm$0.07          0.14$\pm$0.06

1333-C7-1                 0.5$\pm$0.2            0.2$\pm$0.1

79-C1-1                   0.2$\pm$0.1            0.16$\pm$0.07

IC 348-1                 $<$0.05$\pm$0.02       $<$0.18$\pm$0.08

IC 348-10                $<$0.04$\pm$0.02       $<$0.12$\pm$0.05

ORI-C1-2                  $<$0.02$\pm$0.01       $<$0.013$\pm$0.006

Orion     ORI-C2-3          0.10$\pm$0.04          $<$0.03$\pm$0.01

ORI-C3-1                  $<$0.07$\pm$0.03       $<$0.05$\pm$0.02

Deuterium fractionation ratios in different interstellar sources. References are
[Minowa et al., 1997] for TMC-1, [Irvine et al., 1989] for L 183, [Marcelino et al., 2005] for Barnard 1, [Agúndez et al., 2019] for L 483, [Drozdovskaya et al., 2018] for IRAS 16293-2422, [Rivière-Marichalar et al., 2019] for the Horeshead PDR, and [Minh et al., 2018] for G33.92+0.11.

Table <ref> shows the deuterium fractionation ratios obtained in our starless core sample. Figure <ref> shows a comparison between these results and the ones obtained for other starless cores from the literature, as well as a comparison with other type of interstellar sources, such as Class 0 objects and photodissociation regions (PDRs). The HDCS/H$_2$CS ratios obtained in our sample of starless cores are in the range $\sim$0.1-0.8, with Orion showing the lowest value (HDCS/H$_2$CS=0.10$\pm$0.04) compared to those obtained in Taurus and Perseus. The obtained HDCS/H$_2$CS ratios are up to $\sim$4 orders of magnitude higher than the cosmic D/H ratio, but they are similar to those found in other starless cores (see Fig. <ref>), such as TMC-1 and L 183 <cit.>, where the low ($\sim$7-13 K) temperatures of these regions, as well as the large fraction of neutral heavy species freeze-out onto dust grains, favour the formation and enhancement of deuterated species [Millar et al., 1989, Ceccarelli et al., 2014]. In particular, the deuterium enrichment in cold regions occurs through the reaction between the molecular ion H$^{+}_{3}$ with HD to form H$_2$D$^+$, which reacts with other molecules and atoms transferring deuterium to other species <cit.>. Nevertheless, in Fig. <ref>, we also observe that similar HDCS/H$_2$CS ratios ($\sim$0.1-0.5) are found as well in interstellar regions characterised by higher temperatures than those found in starless cores, such as Class 0 protostar <cit.>, high-mass star-forming clouds <cit.>, and PDRs <cit.>. A rich deuterium chemistry is observed in PDRs, such as the Orion Bar [Parise et al., 2009] and Monoceros R2 [Treviño-Morales et al., 2014]. In the regions characterised by warm gas temperatures (30$\lesssim$$T$$\lesssim$100 K), the deuterium fractionation is driven through a chemistry that is different to that from cold regions. In particular, the D atoms are transferred to molecules by CH$_2$D$^+$, whose activation barrier when reacting with H$_2$ is larger than that of H$_2$D$^+$ [Roberts & Millar, 2000, Roueff et al., 2007]. In  the  case  of  single-dish  observations of dense cores around Class 0 protostars, the origin of the deuterated compounds emission could be  in  the  cold  and dense region of the core that still remains unaffected by the heating produced by the recently born star [Agúndez et al., 2019]. In hot corinos, with temperatures $T$$>$100 K, the ice mantles evaporate releasing deuterated compounds formed in the gas and/or the ice during the cold phase [Drozdovskaya et al., 2018].

The derived fractionation ratio for the double deuterated thioformaldehyde, D$_2$CS/H$_2$CS, is in the range $\sim$0.05-0.4 (Table <ref> and Fig. <ref>) for the starless cores in Taurus and Perseus, while only upper limits are reported for those in Orion, where D$_2$CS has not been detected. For the cases where D$_2$CS is detected, the double deuterated ratio is on average a factor $\sim$1.9 lower than the single deuterated ratio (HDCS/H$_2$CS). This is a similar range to that found for Class 0 objects, such as Barnard 1 and L 483, where D$_2$CS/H$_2$CS is $\sim$1.3 and $\sim$3 times lower than HDCS/H$_2$CS, respectively [Marcelino et al., 2005, Agúndez et al., 2019]. For the case of Orion where only upper limits are provided (i.e. cores ORI-C1-2 and ORI-C3-1), we obtain HDCS/H$_2$CS values of $<$0.02 and $<$0.07, respectively, which are about $\sim$10 times lower than those found in Taurus and Perseus. In the case of D$_2$CS/H$_2$CS, we also observe this difference between the Orion cores and the cores in Taurus and Perseus. Only in 1333-C3-1 and B 213-C5-1 (affected by the presence of nearby star formation activity as ORI-C3-1; see Sects. <ref> and <ref>), we derive a deuterium fractionation ratio similar to the ones obtained in Orion. In regions characterised by more extreme environments and higher temperatures, such as PDRs, D$_2$CS has not been detected either.

Our results, therefore, make it evident that the deuterium fractionation is higher in the pre-stellar cores of Taurus and Perseus than those of Orion. Nevertheless, we should also take the different distances of the studied regions (Taurus, Perseus, and Orion) into account. In particular, typical pre-stellar core sizes are of about $\sim$10,000 AU <cit.>, which it is equivalent to $\sim$0.05 pc or $\sim$24$\arcsec$ at a distance of $d$$\sim$428 pc (average distance to Orion). In the case of Perseus ($d$$\sim$300 pc) and Taurus ($d$$\sim$145 pc), the angular sizes would be $\sim$30$\arcsec$ and $\sim$71$\arcsec$, respectively. These angular sizes are significantly larger than the telescope beam (up to $\sim$29”). However, in the Orion case, the angular size is smaller than the beam for the lowest frequencies. This implies that in the nearest regions (in particular in Taurus), the observed emission may be limited to the high-density central region of the cores, while, in the Orion case, we are not able to distinguish the emission coming from the densest zone from that coming from the envelope. This could cause the observed trend in the deuterium fractionation, which shows the lowest values in the Orion cores and the highest values in Taurus. In particular, the D/H values obtained in Orion are a factor of $\sim$10 lower than those obtained in Taurus. This suggests that the D/H ratio in the envelope of the Orion cores should be (much) lower than a factor of 10 with respect to the D/H in the high-density central regions. We would need observations at higher frequencies or interferometry to resolve the emission from the pre-stellar cores in Orion.

§.§ Deuterium fractionation evolution

Apart from temperature, deuterium fractionation is also thought to vary along the dynamical evolutionary stages of starless cores. In particular, when starless cores evolve towards star formation, the deuterium fraction of molecules formed in the gas phase increases, reaching their maximum at the onset of star formation [Crapsi et al., 2005, Hirota & Yamamoto, 2006, Emprechtinger et al., 2009, Feng et al., 2019], while after the star birth, the deuterium fraction decreases [Fontani et al., 2011, Sakai et al., 2012, Gerner et al., 2015]. Deuterium fractionation ratios can therefore be used as evolutionary tracers.

Abundances with respect to total hydrogen nuclei (the initial H$_2$ OPR is 3).

Species      Abundance         Species     Abundance

HD          1.6$\times$10$^{-5}$             S$^+$       1.5$\times$10$^{-5}$/$\times$10$^{-6}$/$\times$10$^{-7}$

He          9.0$\times$10$^{-2}$             Na$^+$      2.0$\times$10$^{-9}$

N           6.2$\times$10$^{-5}$             Mg$^+$      7.0$\times$10$^{-9}$

O           2.4$\times$10$^{-4}$             P$^+$       2.0$\times$10$^{-10}$

C$^+$       1.7$\times$10$^{-4}$             Cl$^+$      1.0$\times$10$^{-9}$

Si$^+$      8.0$\times$10$^{-9}$             F$^+$       6.7$\times$10$^{-9}$

Fe$^+$      3.0$\times$10$^{-9}$

In order to study the evolution of the HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios in starless cores, we used the Nautilus gas-grain chemical code [Ruaud et al., 2016]. Nautilus is a three-phase model in which gas, grain surface, grain mantle phases, and their interactions are considered. Nautilus solves the kinetic equations for both the gas phase and the surface of interstellar dust grains and computes the evolution with time of chemical abundances for a given physical structure. We use the chemical network presented by Majumdar et al., 2017, which considers multiple deuterated molecules and includes the spin chemistry in Nautilus. Moreover, we used this state-of-the-art chemical model to explore the influence that poorly known parameters, such as the sulphur elemental abundance and the CR ionisation rate for H$_2$ molecules, would have on the H$_2$CS deuterium fractionation.

Evolution of the HDCS/H$_2$CS ratio (colour bar) as a function of time and temperature for a CR ionisation rate $\xi$=1.3$\times$10$^{-17}$ s$^{-1}$, a hydrogen number density $n_{\mathrm{H}}$=10$^5$, 10$^6$, and 10$^7$ cm$^{-3}$, and an initial sulphur abundance $S$$^+$=1.5$\times$10$^{-5}$, 1.5$\times$10$^{-6}$, and 1.5$\times$10$^{-7}$.

Evolution of the D$_2$CS/H$_2$CS ratio (colour bar) as a function of time and temperature for a CR ionisation rate $\xi$=1.3$\times$10$^{-17}$ s$^{-1}$, a hydrogen number density $n_{\mathrm{H}}$=10$^5$, 10$^6$, and 10$^7$ cm$^{-3}$, and an initial sulphur abundance $S$$^+$=1.5$\times$10$^{-5}$, 1.5$\times$10$^{-6}$, and 1.5$\times$10$^{-7}$.

In particular, we calculated the evolution of HDCS/H$_2$CS and D$_2$CS/H$_2$CS as a function of time (considering a time range of 10$^3$-10$^7$ yr) and the core temperature (ranging from 7 K to 20 K). We also considered a visual extinction $A$$_{\mathrm{V}}$=15 mag, a UV field strength $\chi$$_{\mathrm{UV}}$=1 (in Draine units), different volume densities (10$^5$, 10$^6$, and 10$^7$ cm$^{-3}$), and also different initial sulphur abundances (1.5$\times$10$^{-5}$, 1.5$\times$10$^{-6}$, and 1.5$\times$10$^{-7}$; see Table <ref>). We considered different initial sulphur abundances
since the value of elemental sulphur abundance is uncertain. In particular, while the observed gaseous sulphur accounts for its total cosmic abundance in the primeval diffuse clouds and HII  regions [Neufeld et al., 2015, Goicoechea & Cuadrado, 2021], there is an unexpected paucity of S-bearing molecules within dense molecular clouds <cit.>. In starless cores (i.e. the cold and dense cores where stars are formed), the sulphur atoms locked within the observed molecules constitute only $<$1$\%$ of the expected amount <cit.>). One could think that most of the sulphur is locked on the icy mantles that cover dust grains but, surprisingly, a similar trend is encountered within the solid phase. Nowadays, solid OCS and tentatively solid SO$_2$, are the unique sulphur compounds observed, while only upper limits to the solid H$_2$S abundance have been derived [Palumbo et al., 1995, Palumbo et al., 1997, Boogert et al., 1997, Jiménez-Escobar & Muñoz Caro, 2011]. According to these measurements, the abundances of the icy species could account for $<$5$\%$ of the total sulphur. It has been suggested that this so-called depleted sulphur may be locked as refractory material <cit.>, claiming sulphur depletion of more than two orders of magnitude. Intense observational and theoretical work has been carried out in the last decade to constrain the value of sulphur depletion, and there is increasing evidence for a moderate factor of 1-20 of sulphur depletion in cold and warm clouds [Esplugues et al., 2013, Esplugues et al., 2014, Fuente et al., 2019, Navarro-Almaida et al., 2020, Navarro-Almaida et al., 2021, Bulut et al., 2021, Codella et al., 2021]. For that, we consider in our models the un-depleted case (1.5$\times$10$^{-5}$), as well as sulphur depletions of a factor of 10 and 100 (Table <ref>).

Figures <ref> and <ref> show the evolution of HDCS/H$_2$CS and D$_2$CS/H$_2$CS, respectively. For the case of HDCS/H$_2$CS, when the number density is $n_{\mathrm{H}}$=10$^5$ cm$^{-3}$, its maximum reached value at any time is $\sim$0.2. However, for higher densities, HDCS/H$_2$CS increases up to $\sim$0.5 at $T_{\mathrm{gas}}$$<$12 K. Nevertheless, this maximum value is reached with a difference of about 5 million years depending on if the considered density is $n_{\mathrm{H}}$=10$^6$ cm$^{-3}$ or 10$^7$ cm$^{-3}$ (the smaller the density, the longer the evolution time). Apart from density, we also observe in Fig. <ref> a variation in the deuterated ratio with sulphur fractional abundance in such a way that the lower the initial $X$(S), the higher the HDCS/H$_2$CS ratio, especially at low temperatures. Similar trend is also observed in Fig. <ref> for D$_2$CS/HDCS.
In particular, we find that decreasing the sulphur initial fractional abundance from 1.5$\times$10$^{-5}$ <cit.> to 1.5$\times$10$^{-7}$, leads to an increase by a factor of up to $\sim$1.5 and $\sim$2.5 in the maximum values of the HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios, respectively.

From Figs. <ref> and <ref>, we also obtain a strong dependence of HDCS/H$_2$CS and D$_2$CS/H$_2$CS with time, since their values progressively increase with the evolution of the starless cores. These ratios, therefore, can be used as powerful tools to derive the evolutionary stage of starless cores. However, they cannot be used to derive the temperature of the starless cores since their evolution presents a double peak in two ranges of temperature ($T_{\mathrm{gas}}$$\sim$7-11 K and $\sim$15-19 K).

§.§ Cosmic-ray impact on deuterium fractionation

Evolution of the HDCS/H$_2$CS (left) and D$_2$CS ratio (right) as a function of time and temperature for an initial sulphur fractional abundance $S$$^+$=1.5$\times$10$^{-7}$, a hydrogen number density $n_{\mathrm{H}}$=10$^6$ cm$^{-3}$, and a CR ionisation rate $\zeta$=1$\times$10$^{-16}$ (bottom row) and 1.3$\times$10$^{-17}$ s$^{-1}$ (top row).

Another important parameter regulating deuterium fractionation in the cold phases of the ISM is the ubiquitous presence of CRs, since in the denser regions where UV photons are efficiently absorbed and cannot penetrate, CRs become the main ionising agent of the gas. In fact, its interaction with H$_2$ produces H$^{+}_{2}$, which rapidly reacts with another hydrogen molecule leading to H$^{+}_{3}$. This ion is the starting point for reactions between charged and neutral species, producing several other key molecules, such as HCO$^+$ and H$_2$D$^+$ [Dalgarno & Lepp, 1984]. In spite of its role in deuterium chemistry, the CR ionisation rate is still uncertain and spreads over a range of values. While <cit.> determined a CR ionisation rate of $\sim$3$\times$10$^{-17}$ s$^{-1}$ in dense clouds, Indriolo & McCall, 2012 obtained a range (1.7-10.6)$\times$10$^{-16}$ s$^{-1}$ of values considering a sample of diffuse clouds, and Neufeld & Wolfire, 2017 derived a CR ionisation rate of the order of a few 10$^{-16}$ s$^{-1}$ in the Galactic disc. A very recent study of the pre-stellar core L 1544 shows that observations are better reproduced with a CR ionisation rate $\zeta$=3$\times$10$^{-17}$ s$^{-1}$ than with $\zeta$=1.3$\times$10$^{-17}$ s$^{-1}$ [Redaelli et al., 2021].

In the same line, [Fuente et al., 2019] derived a CR ionisation rate of $\zeta$=(0.5-1.8)$\times$10$^{-16}$ s$^{-1}$ in the translucent phase of TMC1, probing that high CR ionisation rate is found in the outer regions of molecular clouds. A lower value of $\zeta$=1.3$\times$10$^{-17}$ s$^{-1}$ was derived by Navarro-Almaida et al., 2021 in the dense cores TMC1 (C) and TMC1 (CP), in qualitative agreement with the expected behaviour of $\zeta$ decreasing with visual extinction along the molecular cloud <cit.>.

Classification of the core sample according to the evolutionary stage.

2c$t$<1 My       2c$t$$\sim$1-5 My

Region      Core             HDCS/H$_2$CS        Core           HDCS/H$_2$CS

Taurus      B 213-C7-1      0.2$\pm$0.1        B 213-C1-1     0.5$\pm$0.2

B 213-C16-1     0.14$\pm$0.06      B 213-C2-1     0.3$\pm$0.1

B 213-C10-1     0.13$\pm$0.06      B 213-C6-1     0.3$\pm$0.1

B 213-C5-1      0.12$\pm$0.05

Perseus     79-C1-1          0.2$\pm$0.1        1333-C6-1       0.8$\pm$0.4

1333-C3-14       0.15$\pm$0.07      1333-C7-1       0.5$\pm$0.2

1333-C3-1        0.10$\pm$0.05      1333-C4-1       0.3$\pm$0.1

-                -                  L 1448         0.3$\pm$0.1

Orion       ORI-C2-3         0.10$\pm$0.04          -                 -

In Fig. <ref>, we show the impact of varying  the CR ionisation rate (from $\zeta$=10$^{-16}$ s$^{-1}$ to $\zeta$=1.3$\times$10$^{-17}$ s$^{-1}$, which is considered as the standard CR rate for  dense cores) on the deuterium fractionation for the HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios along 10$^7$ yr. We obtain that, for both models with different $\zeta$, the evolutionary trend is similar for each deuterium fractionation ratio, with them progressively increasing with time and showing two peaks at the range temperatures previously mentioned. Nevertheless, for HDCS/H$_2$CS, we also observe that, although the maximum reached value is the same ($\sim$0.54) independently on the considered $\zeta$, in the model with the higher CR ionisation rate the maximum HDCS/H$_2$CS ratio is reached about 3 million years later than in the model with the low CR ionisation rate. For D$_2$CS/H$_2$CS, we also obtain lower values for the higher $\zeta$ in each evolution time with respect to the model with the low CR ionisation rate. These results indicate that the main effect on the deuterium fractionation when increasing by one order of magnitude the CR ionisation rate is to slow down its evolution by a few million years.

§.§ Comparison with observations

In the Taurus complex, we find the highest HDCS/H$_2$CS ratios ($\sim$0.3-0.5, Table   <ref>) in the starless cores located in the north of the B 213  filament (C1, C2, and C6). Models shown in Fig. <ref> reproduce  these values at $t$$\gtrsim$5$\times$10$^5$ yr, densities $n_{\mathrm{H}}$$>$10$^5$  cm$^{-3}$, and with a high sulphur initial depletion. In fact, the highest deuterium ratio  (HDCS/H$_2$CS=0.5) from the starless core B 213-C1-1 is only reproduced when the sulphur  fractional abundance is as low as 1.5$\times$10$^{-7}$ and for a long evolution time of 1-5 Myr.
Regarding the double deuterated ratio, the starless cores C1, C2, and C6 also present the highest values ($\sim$0.2-0.3) obtained observationally. These D$_2$CS/H$_2$CS ratios are reproduced by a model with high density, $t$$\gtrsim$10$^6$ yr, and also with a low sulphur initial fractional abundance ($\leq$1.5$\times$10$^{-6}$; see Fig. <ref>).
By contrast, the cores C10 and C16 located in the south part of the B 213 filament present low HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios ($\lesssim$0.14). These values are reproduced at earlier times ($t$$<$10$^6$ yr) than those for the cores C1, C2, and C6, indicating that cores located in the south part of B 213 are less chemically evolved than the ones located in the north of the filament. In any case, we should also consider the uncertainties associated with the D/H ratios of these cores (except for the core C1 with the highest D/H even considering its uncertainty; Table <ref>), which makes it difficult to clearly determine the evolutionary stages of the cores since some of these uncertainties represent up to the 50$\%$ of the D/H value. Nevertheless, comparing with results from other authors, we find that our results are in agreement with the ones found by Hacar et al., 2013 through observations of CO and N$_2$H$^+$. In particular, they found that the region harbouring the cores C1, C2, and C6 is weak in C$^{18}$O and has a number of N$_2$H$^+$-bright dense cores, which suggests that some CO depletion has already taken place, and that this region must therefore be more chemically evolved than other regions located more to the south, such as the one hosting C7,
which is more bright in C$^{18}$O. In fact, we obtain a low deuterium ratio in C7 (0.2$\pm$0.1). Table <ref> summarises the evolutionary stages of our core sample obtained by comparing the observations and the model results.

In Perseus, the starless core NGC 1333-C7 presents HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios as high as 0.5 and 0.2, respectively. As we observe in Figs. <ref> and <ref>, we reproduce these values for a model with a density $>$10$^5$ cm$^{-3}$, a $T_{\mathrm{K}}$$<$12 K, a sulphur fractional abundance of at least 1.5$\times$10$^{-6}$, and an evolutionary time between 1 and 5 Myr. This is a slightly older age than that derived (0.5-1.5 Myr) by Lada et al., 1996 and Wilking et al., 2004 for NGC 1333.
For the cores NGC 1333-C4 and L 1448-1, located in the south-western part of the Perseus complex like NGC 1333-C7, we also obtain high ratios of HDCS/H$_2$CS (0.3) and, in particular, of D$_2$CS/H$_2$CS (0.3 and 0.2, respectively). We reproduce them for a chemical age of $\sim$2 Myr for L 1448-1, and of $\sim$5 Myr for NGC 1333-C4 (Fig.  <ref>).
On the other hand, the core 79-C1-1, which is located farthest north-east of the Perseus molecular cloud, presents low deuterium ratios (0.2 and 0.16 for HDCS and D$_2$CS, respectively). For the core IC 348 (also situated in the north-east region), we do not detect emission of HDCS nor D$_2$CS. All this suggests that, in general terms, the north-eastern part of Perseus is in an earlier chemical evolutionary stage ($<$1 Myr) than the south-western part ($\geq$1-5 Myr; see Table <ref>).  In any case, Perseus observations are reproduced considering a sulphur initial fractional abundance $\leq$1.5$\times$10$^{-6}$ as in the case of Taurus.

Main chemical reaction rates forming (top) and destroying (bottom) H$_2$CS.

Regarding Orion, we detect H$_2$CS in the three cores of the sample ($N$=1.9-3.6$\times$10$^{12}$ cm$^{-2}$, Table <ref>), but we  only detect HDCS in ORI-C2-3. This is probably due to the more quiescent properties of OMC-4 (where ORI-C2 is located) compared to the other two regions harbouring the cores ORI-C1 and ORI-C3, as mentioned in Sect. <ref>. The deuterium fractionation in ORI-C2-3 is as low as 0.1, which, for a high-density ($n_{\mathrm{H}}$$\geq$10$^6$ cm$^{-3}$) model characteristic of the Orion region, is reproduced at $t$$\lesssim$0.5 Myr.

§.§ Thioformaldehyde formation and destruction

Out of more than 200 molecules detected in space, about 10% contain sulphur. In spite of the use of astronomical observations coupled with laboratory experiments and astrochemical modelling to study interstellar sulphur chemistry, results have been unsuccessful in predicting fractional abundances of species, such as carbonyl monosulfide (CS), falling short by up to two orders of magnitude compared to astronomical observations <cit.>.
In this way, the analysis of formation mechanisms of (organo)sulphur molecules in different astrochemical environments is key to shed light on the origin and evolution of sulphur in our Galaxy.

Main chemical reaction rates forming (top) and destroying (bottom) HDCS.

Main chemical reaction rates forming (top) and destroying (bottom) D$_2$CS.

Figures <ref>-<ref> show chemical rates for the main chemical reactions forming and destroying H$_2$CS, HDCS, and D$_2$CS calculated with the Nautilus code, considering a sulphur fractional abundance of 1.5$\times$10$^{-7}$, core temperatures of 10 K and 20 K, and densities ranging from 10$^{5}$ to 10$^{7}$ cm$^{-3}$. In the case of H$_2$CS (Fig. <ref>), the main gas-phase routes (i.e. the reactions with the highest reaction rates) forming this species at any evolutionary time are

\begin{equation}
\mathrm{S + CH_3 \rightarrow H + H_2CS}
\end{equation}

\begin{equation}
\mathrm{H_3CS^+ + e^- \rightarrow H + H_2CS}
\end{equation}

in agreement with Prasad & Huntress, 1982 and Laas & Caselli, 2019. Nevertheless, we note that while the ion H$_3$CS$^+$ is mainly formed through the reaction between CH$_4$ and S$^+$ at $t$$\lesssim$10$^4$ yr, for longer times H$_3$CS$^+$ is formed in turn through H$_2$CS.
In the case of a very high-density (10$^{7}$ cm$^{-3}$) starless core (right top panel) and $T$$_{\mathrm{g}}$=20 K (dashed lines), the surface reactions

\begin{equation}
\mathrm{JH + JHCS \rightarrow H_2CS}
\end{equation}

\begin{equation}
\mathrm{JH_2CS \rightarrow H_2CS}
\end{equation}

also become dominant at 5$\times$10$^4$$\lesssim$$t$$\lesssim$5$\times$10$^5$ yr. These chemical and thermal desorption reactions barely influence instead the formation of H$_2$CS when the density of the core is $\leq$10$^6$ cm$^{-3}$.
Regarding destruction, H$_2$CS is mainly destroyed through its reaction with C$^+$ during the early timescales ($t$$\lesssim$10$^4$, $\lesssim$10$^3$, and $\lesssim$10$^2$ yr for $n_{\mathrm{H}}$=10$^5$, 10$^6$, and 10$^7$ cm$^{-3}$, respectively). For longer times, H$_2$CS depletion becomes one of the main destruction routes, especially at very high densities, together with the reaction between H$_2$CS and H$^+$.

Formation and destruction reaction rates for HDCS are shown in Fig. <ref>. HDCS is mainly formed in all the considered cases by the reactions

\begin{equation}
\mathrm{H_2DCS^+ + e^- \rightarrow H + HDCS}
\end{equation}

\begin{equation}
\mathrm{S + CH_2D \rightarrow H + HDCS}.
\end{equation}

Nevertheless, in the low temperature case (10 K), the chemical desorption

\begin{equation}
\mathrm{JH + JDCS \rightarrow HDCS}
\end{equation}

also plays an important role in forming HDCS at $t$$<$10$^3$ yr. Regarding HDCS destruction, its reaction with C$^+$ represents the main destroying mechanism for $t$$\lesssim$10$^4$ yr with $n$$\leq$10$^6$ cm$^{-3}$. At longer times, HDCS is mainly destroyed by H$^+$. In a very high-density core (right bottom panel), depletion of HDCS becomes also important removing this deuterated species from the gas phase.

With respect to the double deuterated thioformaldehyde, for an early time ($t$$\lesssim$10$^3$ yr) and a cold core ($T$$_{\mathrm{g}}$=10 K), D$_2$CS is efficiently formed both on dust grains (through chemical desorption from the reaction between JD and JDCS) and through the gas-phase reaction

\begin{equation}
\mathrm{HD_2CS^+ + e^- \rightarrow H + D_2CS},
\end{equation}

while for a warmer core ($T$$_{\mathrm{g}}$=20 K) the chemical desorption becomes less efficient in forming D$_2$CS, in favour of the neutral–neutral gas-phase reaction

\begin{equation}
\mathrm{S + CHD_2 \rightarrow H + D_2CS}.
\end{equation}

At late timescales ($>$10$^5$ yr), D$_2$CS is mainly formed in the gas phase through the reactions

\begin{equation}
\mathrm{HD_2CS^+ + e^- \rightarrow H + D_2CS},
\end{equation}

\begin{equation}
\mathrm{S + CHD_2 \rightarrow H + D_2CS},
\end{equation}

\begin{equation}
\mathrm{DC_3S^+ + e^- \rightarrow D + D_2CS}.
\end{equation}

By contrast, surface reactions are irrelevant forming D$_2$CS at long times ($>$10$^5$ yr), independently on the temperature and density of the starless core. Regarding D$_2$CS destruction, as in the case of HDCS and H$_2$CS, the ions C$^+$ and H$^+$ are the main destroyers of D$_2$CS. D$_2$CS depletion has also an important role, especially in cores with very high density, where it becomes the main D$_2$CS destroying mechanism during a period of $\sim$10$^4$ years.

All these results show that (deuterated) thioformaldehyde is mainly formed through gas-phase reactions, where the ions H$_3$CS$^+$, H$_2$DCS$^+$, and HD$_2$CS$^+$ (formed in turn from H$_2$D$^+$ and D$_2$H$^+$), as well as the neutral precursors CH$_3$, CH$_2$D, and CHD$_2$, have a key role in the formation of H$_2$CS, HDCS, and D$_2$CS, respectively.

§ SUMMARY AND CONCLUSIONS

Within the IRAM 30m Large Program GEMS, we have carried out a comprehensive observational and theoretical study of thioformaldehyde and its deuterated counterparts in a sample of starless cores in filaments of the nearby star-forming regions Taurus, Perseus, and Orion. These regions have different degrees of star formation activity, and therefore different physical conditions, providing a possibility to explore the effect of environment on gas chemistry.

Using the molecular excitation and radiative transfer code RADEX, we have obtained column densities and fractional abundances for H$_2$CS, HDCS, and D$_2$CS in each core of the sample, also providing upper limits for cases where these species were not detected. We have also derived deuterium fractionation ratios for both single and double deuterated thioformaldehyde. The obtained ratios are up to four orders of magnitude higher than the cosmic D/H ratio. This result is similar to the ones obtained in other starless cores, such as L 183, and it is due to the low temperatures of these regions, as well as the large fraction of neutral heavy species that freeze-out onto dust grains, favouring the formation of deuterated species. We have also found that the HDCS/H$_2$CS ratio is similar to that found in regions characterised by higher temperatures, such as Class 0 objects, high-mass star-forming clouds, and PDRs. In these warm regions, the deuterium fractionation is driven by a chemistry different from that of cold regions, with D atoms being transferred to molecules by CH$_2$D$^+$ instead of H$_2$D$^+$.

Since the deuterium fractionation ratio can be considered as an evolutionary tracer, we have also used it to derive the chemical evolution stage of the starless core sample through the Nautilus gas-grain chemical code, which computes the evolution of chemical abundances as a function of time. The comparison of model results with observations reveals that the north region of the B 213 filament in Taurus, harbouring cores such as C1, C2, and C6, is more chemically evolved ($t$$\gtrsim$1 Myr) than other regions located more to the south of the filament ($t$$<$1 Myr), where we find the cores C7, C10, and C16. In Perseus, the results indicate that the north-eastern part presents an earlier chemical evolutionary stage ($t$$<$1 Myr) than the south-western part ($t$$\geq$1-5 Myr). In Orion, for the only core where deuterated thioformaldehyde is detected, we have obtained a chemical age of $t$$\lesssim$0.5 Myr.

The theoretical study of deuterium fractionation ($D$$_{\mathrm{frac}}$) evolution has not only allowed us to analyse how the ratios HDCS/H$_2$CS and D$_2$CS/H$_2$CS vary with time, but also with other parameters, such as density, temperature, initial sulphur depletion, and the CR ionisation rate. Our results confirm a strong dependence of $D$$_{\mathrm{frac}}$ with time since it progressively increases with the evolution of the starless cores. This makes HDCS/H$_2$CS and D$_2$CS/H$_2$CS powerful tools for deriving the chemical evolution of these regions. However, they cannot be used to derive the temperature of the starless cores since their value presents a double peak in two temperature ranges ($T$$_{\mathrm{gas}}$=$\sim$7-11 K and $T$$_{\mathrm{gas}}$=$\sim$15-19 K).
On the other hand, our results also show that both deuterium fractionation ratios increase with the density of the region and with the degree of sulphur depletion. In particular, we find that decreasing the sulphur initial abundance from 1.5$\times$10$^{-5}$ (the solar elemental sulphur fractional abundance) to 1.5$\times$10$^{-7}$ leads to an increase by a factor of up to $\sim$1.5 and $\sim$2.5 in the maximum values of the HDCS/H$_2$CS and D$_2$CS/H$_2$CS ratios, respectively. In any case, the comparison with the observations of the starless cores of our sample reveals that they are reproduced when the initial sulphur fractional abundance is as low as $\leq$1.5$\times$10$^{-6}$ (i.e. at least one order of magnitude lower than the solar elemental sulphur abundance), in agreement with previous GEMS results <cit.>. Regarding the CR ionisation rate, we find that its increase by one order of magnitude slows down the deuterium fraction evolution by a few million years.

Finally, we have also studied the main formation and destruction mechanisms of H$_2$CS, HDCS, and D$_2$CS using the Nautilus code. We have obtained that they are mainly formed through gas-phase reactions (double-replacement reactions and neutral–neutral displacement reactions). By contrast, grain chemistry (depletion) has an important role in the destruction of these species, especially at $t$$>$10$^4$ yr.

We finally stress that this work has allowed the sulphur chemistry in cold regions to be studied in detail through thioformaldehyde observations. We have also demonstrated the power of using its deuterated counterparts for the analysis of the properties and chemical evolution of starless cores.

We thank the Spanish MICINN for funding support from PID2019-106235GB-I00. L.M. acknowledges the financial support of DAE and DST-SERB research grants (SRG/2021/002116 and MTR/2021/000864) of the Government of India.

§ TABLES AND FIGURES

Observational parameters.

Telescope     Observed frequency range (GHx)     rms (mK)     HPBW($\arcsec$)

IRAM 30m      85.0-86.8          10-20         29

88.3-91.5          10-20         27

92.0-94.8          10-20         26

95.3-97.0          10-20         24

97.6-99.4          10-20         23

101.0-103.0        10-20         21

105.4-107.2        10-20         20

108.7-110.5        10-20         19

133.7-135.5        10-20         18

137.0-138.8        10-20         17

143.2-145.0        10-20         16

146.4-148.3        10-20         16

158.8-160.7        20-30         15

162.1-163.9        20-30         15

167.2-169.0        20-30         14

170.5-172.3        20-30         14

H$_{2}$CS parameters from Gaussian fits.

Core     Transition     3co-H$_{2}$CS     Transition     3cp-H$_{2}$CS

$J_{K,k}-J'_{K',k'}$      $V$$_{\mathrm{LSR}}$ (km s$^{-1}$)     $\Delta$$V$ (km s$^{-1}$)     $T$$_{\mathrm{MB}}$ (K)     $J_{K,k}-J'_{K',k'}$
$V$$_{\mathrm{LSR}}$ (km s$^{-1}$)     $\Delta$$V$ (km s$^{-1}$)     $T$$_{\mathrm{MB}}$ (K)

B 213-C1-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       5.89$\pm$0.01     0.38$\pm$0.08      0.5$\pm$0.1     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   6.02$\pm$0.01     0.27$\pm$0.05                 0.37$\pm$0.07

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 6.23$\pm$0.01       0.27$\pm$0.05      0.37$\pm$0.07     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 5.63$\pm$0.03       0.35$\pm$0.07               0.08$\pm$0.02

B 213-C2-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.98$\pm$0.01     0.42$\pm$0.08      0.15$\pm$0.03     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   7.11$\pm$0.04       0.44$\pm$0.09               0.06$\pm$0.01

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.28$\pm$0.03       0.35$\pm$0.07      0.07$\pm$0.01     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   7.02$\pm$0.05       0.28$\pm$0.06               0.04$\pm$0.01

B 213-C5-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.32$\pm$0.01     0.6$\pm$0.1      0.27$\pm$0.05     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   6.42$\pm$0.03       0.7$\pm$0.1               0.08$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 6.60$\pm$0.02       0.5$\pm$0.1      0.10$\pm$0.02     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   6.94$\pm$0.05       $<$0.25$\pm$0.05              $<$0.07$\pm$0.01

B 213-C6-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.90$\pm$0.01     0.36$\pm$0.07     0.30$\pm$0.06     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   6.96$\pm$0.02       0.31$\pm$0.06               0.09$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 6.72$\pm$0.02       0.34$\pm$0.06      0.13$\pm$0.03     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   6.78$\pm$0.05       $<$0.21$\pm$0.04      $<$0.05$\pm$0.01

B 213-C7-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.80$\pm$0.01     0.39$\pm$0.08      0.45$\pm$0.09     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   6.83$\pm$0.01       0.25$\pm$0.05               0.23$\pm$0.05

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 6.60$\pm$0.01       0.31$\pm$0.06      0.26$\pm$0.06     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   7.14$\pm$0.09       $<$0.35$\pm$0.07      $<$0.04$\pm$0.01

B 213-C10-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      6.74$\pm$0.02       0.45$\pm$0.09      0.13$\pm$0.03     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                 6.72$\pm$0.03       0.23$\pm$0.05      0.05$\pm$0.01

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                  6.01$\pm$0.03       0.21$\pm$0.04      0.07$\pm$0.01      5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 6.02$\pm$0.04       $<$0.18$\pm$0.04          $<$0.06$\pm$0.01

B 213-C12-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      6.64$\pm$0.02       0.42$\pm$0.08      0.09$\pm$0.02     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                 6.65$\pm$0.04       0.24$\pm$0.05      0.06$\pm$0.01

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                  5.82$\pm$0.04        0.20$\pm$0.04      0.06$\pm$0.01       5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 6.58$\pm$0.02       $<$0.12$\pm$0.02      $<$0.07$\pm$0.01

B 213-C16-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      6.64$\pm$0.01     0.41$\pm$0.08      0.6$\pm$0.1     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   6.77$\pm$0.01       0.39$\pm$0.08      0.22$\pm$0.04

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                   6.72$\pm$0.01     0.30$\pm$0.06      0.33$\pm$0.06      5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                6.75$\pm$0.04       $<$0.12$\pm$0.02      $<$0.04$\pm$0.01

L 1448-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$         4.30$\pm$0.01     0.6$\pm$0.1      0.7$\pm$0.1     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   4.44$\pm$0.01     0.7$\pm$0.1      0.31$\pm$0.06

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                   4.42$\pm$0.01     0.6$\pm$0.1      0.35$\pm$0.07     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$               4.79$\pm$0.03       0.41$\pm$0.08      0.25$\pm$0.05

1333-C3-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       8.49$\pm$0.01     0.7$\pm$0.1      0.24$\pm$0.05     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                     8.65$\pm$0.01       0.42$\pm$0.08      0.11$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 8.87$\pm$0.01       0.6$\pm$0.2      0.22$\pm$0.04      5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   8.58$\pm$0.04       0.32$\pm$0.06      0.11$\pm$0.02

1333-C4-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.67$\pm$0.01     1.0$\pm$0.2      0.7$\pm$0.1     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                     7.73$\pm$0.01       0.7$\pm$0.1       0.45$\pm$0.09

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.79$\pm$0.08      0.9$\pm$0.2      0.6$\pm$0.1      5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   7.83$\pm$0.05       0.6$\pm$0.1      0.17$\pm$0.03

1333-C5-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.63$\pm$0.03       0.8$\pm$0.2      0.09$\pm$0.02     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   7.62$\pm$0.02       0.18$\pm$0.04      0.08$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.00$\pm$0.04       0.45$\pm$0.09      0.11$\pm$0.02      5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                  6.49$\pm$0.05       $<$0.10$\pm$0.02      $<$0.11$\pm$0.02

1333-C6-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.42$\pm$0.04       0.36$\pm$0.07      0.08$\pm$0.02     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   7.36$\pm$0.05       $<$0.17$\pm$0.03      $<$0.05$\pm$0.01

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.06$\pm$0.04       0.31$\pm$0.06      0.08$\pm$0.02     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   7.50$\pm$0.09       $<$0.15$\pm$0.03      $<$0.5$\pm$0.01

1333-C3-14     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      8.21$\pm$0.01     0.6$\pm$0.1      0.5$\pm$0.1     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                     8.28$\pm$0.01       0.6$\pm$0.1      0.23$\pm$0.05

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 8.26$\pm$0.01     0.5$\pm$0.1      0.32$\pm$0.06     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                     8.31$\pm$0.03       0.30$\pm$0.06      0.17$\pm$0.03

1333-C7-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.56$\pm$0.02       0.8$\pm$0.2      0.17$\pm$0.03     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   7.66$\pm$0.02       0.25$\pm$0.05      0.13$\pm$0.03

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.04$\pm$0.02       0.33$\pm$0.07      0.16$\pm$0.03     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   7.59$\pm$0.06       $<$0.10$\pm$0.02      $<$0.06$\pm$ 0.01

79-C1-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$         10.33$\pm$0.01     0.6$\pm$0.1      0.36$\pm$0.07     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                    10.06$\pm$0.01       0.46$\pm$0.09       0.19$\pm$0.04

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 10.20$\pm$0.01     0.5$\pm$0.1      0.23$\pm$0.05     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 10.45$\pm$0.03       0.28$\pm$0.06      0.09$\pm$0.02

IC 348-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       8.94$\pm$0.02       0.7$\pm$0.1      0.11$\pm$0.02     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                 8.94$\pm$0.02       0.43$\pm$0.09      0.09$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 8.65$\pm$0.02       0.33$\pm$0.07      0.11$\pm$0.02       5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 7.91$\pm$0.04       $<$0.13$\pm$0.03      $<$0.06$\pm$0.01

IC 348-10     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      8.46$\pm$0.02       0.6$\pm$0.1      0.13$\pm$0.03     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                 8.46$\pm$0.02       0.29$\pm$0.06      0.08$\pm$0.02

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 8.15$\pm$0.02       0.36$\pm$0.07      0.11$\pm$0.02       5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                 8.70$\pm$0.04       $<$0.10$\pm$0.02      $<$0.05$\pm$0.01

ORI-C1-2     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$        10.96$\pm$0.02      0.8$\pm$0.2      0.36$\pm$0.07     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   11.05$\pm$0.07      0.6$\pm$0.1      0.15$\pm$0.03

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 11.10$\pm$0.03      0.7$\pm$0.1      0.27$\pm$0.05     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   11.31$\pm$0.05      0.5$\pm$0.1      0.09$\pm$0.02

ORI-C2-3     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$        8.16$\pm$0.01       0.41$\pm$0.08      0.24$\pm$0.05     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                   8.13$\pm$0.02       0.29$\pm$0.06      0.15$\pm$ 0.03

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 7.30$\pm$0.03       0.40$\pm$0.08      0.13$\pm$0.03     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                   8.54$\pm$0.03       $<$0.22$\pm$0.04      $<$0.12$\pm$0.02

ORI-C3-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$        10.63$\pm$0.02       1.0$\pm$0.2      0.18$\pm$0.04     4$_{0}$$_{,}$$_{4}$-3$_{0}$$_{,}$$_{3}$                  10.53$\pm$0.02      0.33$\pm$0.07      0.14$\pm$0.03

4$_{1}$$_{,}$$_{4}$-3$_{1}$$_{,}$$_{3}$                 10.40$\pm$0.02       0.7$\pm$0.1      0.23$\pm$0.05     5$_{0}$$_{,}$$_{5}$-4$_{0}$$_{,}$$_{4}$                  11.4$\pm$0.05       $<$0.19$\pm$0.04      $<$0.18$\pm$0.04

HDCS parameters from Gaussian fits.

Core     Transition     3cHDCS

$J_{K,k}-J'_{K',k'}$     $V$$_{\mathrm{LSR}}$ (km s$^{-1}$)     $\Delta$$V$ (km s$^{-1}$)     $T$$_{\mathrm{MB}}$ (K)

B 213-C1-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       5.88$\pm$0.01     0.6$\pm$0.2      0.06$\pm$0.02

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 5.94$\pm$0.01     0.37$\pm$0.09      0.29$\pm$0.07

B 213-C2-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.07$\pm$0.04       0.4$\pm$0.1      0.04$\pm$0.01

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 7.07$\pm$0.01     0.4$\pm$0.1      0.10$\pm$0.02

B 213-C5-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.33$\pm$0.05       0.36$\pm$0.09      0.04$\pm$0.01

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 6.40$\pm$0.02       0.5$\pm$0.1      0.11$\pm$0.03

B 213-C6-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       6.98$\pm$0.01     0.36$\pm$0.09      0.19$\pm$0.05

B 213-C7-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       6.90$\pm$0.03       0.5$\pm$0.1      0.07$\pm$0.02

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 6.86$\pm$0.01     0.36$\pm$0.09      0.18$\pm$0.05

B 213-C10-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$      6.68$\pm$0.03       0.30$\pm$0.08      0.04$\pm$0.01

B 213-C16-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$      6.74$\pm$0.03       0.32$\pm$0.08      0.04$\pm$0.01

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 6.72$\pm$0.01       0.4$\pm$0.1      0.15$\pm$0.04

L 1448-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$         4.51$\pm$0.02       0.5$\pm$0.1      0.12$\pm$0.04

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 4.46$\pm$0.01     0.5$\pm$0.1      0.4$\pm$0.1

1333-C3-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       8.57$\pm$0.06       0.4$\pm$0.1      0.05$\pm$0.01

1333-C4-1     3$_{1}$$_{,}$$_{3}$-2$_{1}$$_{,}$$_{2}$       7.79$\pm$0.04       1.2$\pm$0.3      0.09$\pm$0.02

3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$                 7.67$\pm$0.02       1.2$\pm$0.3      0.18$\pm$0.05

1333-C6-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       7.09$\pm$0.04       0.5$\pm$0.1      0.05$\pm$0.01

1333-C3-14     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$      8.23$\pm$0.01       0.6$\pm$0.1      0.16$\pm$0.04

1333-C7-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       7.58$\pm$0.02       0.6$\pm$0.1      0.13$\pm$0.03

79-C1-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$         10.43$\pm$0.01      0.6$\pm$0.1      0.16$\pm$0.04

ORI-C2-3     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$        8.28$\pm$0.05       0.5$\pm$0.1      0.04$\pm$0.01

D$_{2}$CS parameters from Gaussian fits.

Core     Transition     3co-D$_{2}$CS     Transition     3cp-D$_{2}$CS

$J_{K,k}-J'_{K',k'}$     $V$$_{\mathrm{LSR}}$ (km s$^{-1}$)     $\Delta$$V$ (km s$^{-1}$)     $T$$_{\mathrm{MB}}$ (K)     $J_{K,k}-J'_{K',k'}$
$V$$_{\mathrm{LSR}}$ (km s$^{-1}$)     $\Delta$$V$ (km s$^{-1}$)     $T$$_{\mathrm{MB}}$ (K)

B 213-C1-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       5.95$\pm$0.01       0.31$\pm$0.08      0.21$\pm$0.05     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                 5.78$\pm$0.04       0.30$\pm$0.09              0.05$\pm$0.01

B 213-C2-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$      7.09$\pm$0.02        0.32$\pm$0.08      0.08$\pm$0.02     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  7.05$\pm$0.03       $<$0.18$\pm$0.05              $<$0.04$\pm$0.01

B 213-C5-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       6.43$\pm$0.04       0.36$\pm$0.09      0.05$\pm$0.01     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  6.37$\pm$0.09       $<$0.20$\pm$0.06       $<$0.02$\pm$0.01

B 213-C6-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       7.02$\pm$0.02       0.36$\pm$0.08      0.14$\pm$0.04     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  6.99$\pm$0.03       $<$0.14$\pm$0.04              $<$0.05$\pm$0.01

B 213-C7-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$       6.93$\pm$0.02       0.30$\pm$0.08      0.10$\pm$0.03     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  6.84$\pm$0.04       $<$0.17$\pm$0.05              $<$0.05$\pm$0.01

B 213-C12-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$      6.77$\pm$0.06       0.5$\pm$0.1      0.04$\pm$0.01     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  6.60$\pm$0.05       $<$0.15$\pm$0.05       $<$0.03$\pm$0.01

B 213-C16-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$      6.79$\pm$0.04       0.4$\pm$0.1      0.07$\pm$0.02     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  6.97$\pm$0.07       $<$0.18$\pm$0.05     $<$0.02$\pm$0.01

L 1448-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$         4.54$\pm$0.01     0.4$\pm$0.1      0.30$\pm$0.07     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                    4.51$\pm$0.03       0.30$\pm$0.09      0.05$\pm$0.02

1333-C4-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$         7.76$\pm$0.01       1.0$\pm$0.3      0.21$\pm$0.05     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  7.37$\pm$0.06       $<$0.25$\pm$0.08      $<$0.04$\pm$0.01

1333-C3-14     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$        8.29$\pm$0.04       0.5$\pm$0.1      0.12$\pm$0.03     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  8.00$\pm$0.04       $<$0.23$\pm$0.07      $<$0.04$\pm$0.01

1333-C7-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$         7.73$\pm$0.03       0.34$\pm$0.09      0.07$\pm$0.02     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  7.84$\pm$0.04       $<$0.22$\pm$0.07      $<$0.06$\pm$0.02

79-C1-1     3$_{0}$$_{,}$$_{3}$-2$_{0}$$_{,}$$_{2}$           10.46$\pm$0.02      0.4$\pm$0.1      0.12$\pm$0.03     5$_{1}$$_{,}$$_{5}$-4$_{1}$$_{,}$$_{4}$                  10.24$\pm$0.06      $<$0.4$\pm$0.1     $<$0.03$\pm$0.01

Column densities, $N$, for H$_{2}$CS, HDCS, and D$_{2}$CS and the OPR obtained from H$_2$CS.

Region     Core     $N$(o-H$_{2}$CS)     $N$(p-H$_{2}$CS )    OPR(H$_2$CS)     $N$(H$_{2}$CS)     $N$(HDCS)      $N$(o-D$_{2}$CS)     $N$(p-D$_{2}$CS)     OPR(D$_2$CS)    $N$(D$_{2}$CS)

$\times$10$^{12}$ cm$^{-2}$     $\times$10$^{12}$ cm$^{-2}$          $\times$10$^{12}$ cm$^{-2}$     $\times$10$^{12}$ cm$^{-2}$      $\times$10$^{12}$ cm$^{-2}$     $\times$10$^{12}$ cm$^{-2}$         $\times$10$^{12}$ cm$^{-2}$

B 213-C1-1               2.9$\pm$0.6       1.0$\pm$0.2            3$\pm$1            3.9$\pm$0.8            1.8$\pm$0.5           0.7$\pm$0.2                 0.4$\pm$0.1                -               1.1$\pm$0.3

B 213-C2-1               1.6$\pm$0.3       0.9$\pm$0.2            1.8$\pm$0.7        2.5$\pm$0.5            0.8$\pm$0.2           0.4$\pm$0.1                 0.19$\pm$0.05$^*$          2$\pm$1         0.6$\pm$0.2$^*$

B 213-C5-1               5$\pm$1           2.2$\pm$0.4            2.2$\pm$0.9        7$\pm$1                0.8$\pm$0.2           0.25$\pm$0.06               0.12$\pm$0.03$^*$          -               0.37$\pm$0.09$^*$

Taurus     B 213-C6-1      3.0$\pm$0.6       1.3$\pm$0.3            2.3$\pm$0.9        4.3$\pm$0.9            1.2$\pm$0.3           0.8$\pm$0.2                 0.38$\pm$0.09$^*$          -               1.1$\pm$0.3$^*$

B 213-C7-1               2.8$\pm$0.6       1.6$\pm$0.3            1.8$\pm$0.7        4.4$\pm$0.9            1.0$\pm$0.3           0.35$\pm$0.09               0.18$\pm$0.04$^*$          -               0.5$\pm$0.1$^*$

B 213-C10-1              1.0$\pm$0.2       0.6$\pm$0.1            1.6$\pm$0.6        1.6$\pm$0.3            0.21$\pm$0.05         $<$0.18$\pm$0.05            $<$0.09$\pm$0.02$^*$       -               $<$0.27$\pm$0.07$^*$

B 213-C12-1              0.42$\pm$0.08     0.23$\pm$0.05          1.8$\pm$0.7        0.7$\pm$0.1            $<$0.22$\pm$0.05      0.20$\pm$0.05               0.10$\pm$0.03$^*$          -               0.30$\pm$0.08$^*$

B 213-C16-1              3.9$\pm$0.8       1.7$\pm$0.3            2.3$\pm$0.9        6$\pm$1                0.8$\pm$0.2           0.4$\pm$0.1                 0.21$\pm$0.05$^*$          -               0.6$\pm$0.2$^*$

L 1448-1                 7$\pm$1           4.0$\pm$0.8            1.8$\pm$0.7        11$\pm$2               3.5$\pm$0.9            1.8$\pm$0.4                0.7$\pm$0.2                 2$\pm$1        2.6$\pm$0.7

1333-C3-1                 2.0$\pm$0.4       0.5$\pm$0.1            4$\pm$2            2.5$\pm$0.5            0.26$\pm$0.06          $<$0.04$\pm$0.01           $<$0.02$\pm$0.01$^*$       -               $<$0.07$\pm$0.02$^*$

1333-C4-1                 8$\pm$2           3.0$\pm$0.6            3$\pm$1            11$\pm$3               3.4$\pm$0.9            2.4$\pm$0.6                1.2$\pm$0.3$^*$            -               3.6$\pm$0.9$^*$

1333-C5-1                 0.8$\pm$0.2       0.16$\pm$0.03          5$\pm$2            1.0$\pm$0.2            $<$0.07$\pm$0.02       $<$0.10$\pm$0.03           $<$0.05$\pm$0.01$^*$       -               $<$0.15$\pm$0.04$^*$

Perseus     1333-C6-1       0.34$\pm$0.07     0.14$\pm$0.03$^{**}$      -               0.5$\pm$0.1$^{**}$     0.4$\pm$0.1            $<$0.027$\pm$0.007         $<$0.014$\pm$0.003$^*$     -               $<$0.04$\pm$0.01$^*$

1333-C3-14                5$\pm$1           2.1$\pm$0.4            2.4$\pm$0.9        7$\pm$1                1.1$\pm$0.3            0.7$\pm$0.2                0.33$\pm$0.08$^*$          -               1.0$\pm$0.3$^*$

1333-C7-1                 1.3$\pm$0.3       0.6$\pm$0.1            2.1$\pm$0.8        1.9$\pm$0.4            1.0$\pm$0.2            0.26$\pm$0.07              0.13$\pm$0.03$^*$          -               0.4$\pm$0.1$^*$

79-C1-1                   3.4$\pm$0.7       1.9$\pm$0.4            1.8$\pm$0.7        5$\pm$1                1.2$\pm$0.3            0.6$\pm$0.1                0.28$\pm$0.07$^*$          -               0.9$\pm$0.2$^*$

IC 348-1                 0.8$\pm$0.2       0.38$\pm$0.08          2.1$\pm$0.8        1.2$\pm$0.3            $<$0.06$\pm$0.01       $<$0.14$\pm$0.03           $<$0.07$\pm$0.02$^*$       -               $<$0.21$\pm$0.05$^*$

IC 348-10                0.8$\pm$0.2       0.35$\pm$0.07          2.3$\pm$0.9        1.2$\pm$0.3            $<$0.04$\pm$0.01       $<$0.09$\pm$0.02           $<$0.05$\pm$0.01$^*$       -               $<$0.14$\pm$0.03$^*$

ORI-C1-2                  3.6$\pm$0.7       1.3$\pm$0.3            3$\pm$1            5$\pm$1                $<$0.11$\pm$0.03       $<$0.04$\pm$0.01           $<$0.021$\pm$0.005$^*$     -               $<$0.06$\pm$0.02$^*$

Orion     ORI-C2-3          2.0$\pm$0.4       1.8$\pm$0.4            1.2$\pm$0.5        3.8$\pm$0.8            0.37$\pm$0.09          $<$0.06$\pm$0.02           $<$0.032$\pm$0.008$^*$     -               $<$0.10$\pm$0.02$^*$

ORI-C3-1                  1.9$\pm$0.4       0.6$\pm$0.1            3$\pm$1            2.5$\pm$0.5            $<$0.16$\pm$0.04       $<$0.08$\pm$0.02           $<$0.04$\pm$0.01$^*$       -               $<$0.12$\pm$0.03$^*$

$^{**}$Assuming an average OPR(H$_2$CS)=2.4$\pm$0.9 obtained from all the cores of the sample, except 1333-C6-1, where we only detect o-H$_2$CS

$^*$Assuming an average OPR(D$_2$CS)=2.0$\pm$1.0 obtained from the cores B 213-C1-1 and L 1448-1, where we detect both o- and p-D$_2$CS

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in B 213-C1, C2, C5, C6, C7, C10, C12, and C16 (black), and the best fit (red).

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in L 1448 (black), and the best fit (red).

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in the core sample of NGC 1333 (black) and the best fit (red).

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in Barnard 5 (black) and the best fit (red).

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in the core sample of IC 348 (black) and the best fit (red).

Observed lines of o-H$_{2}$CS and p-H$_{2}$CS in the core sample of Orion A (black) and the best fit (red).

Observed lines of HDCS in B 213-C1, C2, C5, C6, C7, C10, and C16 (black) and the best fit (red).

Observed lines of HDCS in L1448 (black) and the best fit (red).

Observed lines of HDCS in the core sample of NGC 1333 (black) and the best fit (red).

Observed line of HDCS in Barnard 5 (black) and the best fit (red).

Observed line of HDCS in Orion A (black) and the best fit (red).

Observed lines of o-D$_{2}$CS and p-D$_{2}$CS in L1448 (black) and the best fit (red).

Observed lines of o-D$_{2}$CS in the core sample of NGC 1333 (black) and the best fit (red).

Observed lines of o-D$_{2}$CS and p-D$_{2}$CS in B 213-C1, C2, C5, C6, C7, C12, and C16 (black) and the best fit (red).

Evolution of the H$_2$CS fractional abundance as a function of time for an initial sulphur abundance $S$$^+$=1.5$\times$10$^{-7}$, a hydrogen number density $n_{\mathrm{H}}$=10$^5$ cm$^{-3}$, and a CR ionisation rate $\zeta$=1.3$\times$10$^{-17}$ s$^{-1}$ (left) and $\zeta$=1.3$\times$10$^{-16}$ s$^{-1}$ (right).

[Agúndez et al., 2008]
Agúndez, M., Fonfría, J. P., Cernicharo, J., Pardo, J. R., &
Guélin, M. 2008, , 479, 493

[Agúndez et al., 2019]
Agúndez, M., Marcelino, N., Cernicharo, J., Roueff, E., &
Tafalla, M. 2019, , 625, A147

[Agúndez & Wakelam, 2013]
Agúndez, M. & Wakelam, V. 2013, Chemical Reviews, 113, 8710

[Aikawa et al., 2003]
Aikawa, Y., Ohashi, N., & Herbst, E. 2003, , 593, 906

[André et al., 2010]
André, P., Men'shchikov, A., Bontemps, S., et al. 2010, , 518,

[Asplund et al., 2009]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, , 47,

[Bachiller & Cernicharo, 1986]
Bachiller, R. & Cernicharo, J. 1986, , 166, 283

[Bachiller & Cernicharo, 1990]
Bachiller, R. & Cernicharo, J. 1990, , 239, 276

[Bachiller & Pérez Gutiérrez, 1997]
Bachiller, R. & Pérez Gutiérrez, M. 1997, , 487, L93

[Bacmann et al., 2003]
Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2003, , 585, L55

[Bally et al., 1987]
Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ,
312, L45

[Benson & Myers, 1989]
Benson, P. J. & Myers, P. C. 1989, , 71, 89

[Bergin & Langer, 1997]
Bergin, E. A. & Langer, W. D. 1997, , 486, 316

[Bernard et al., 2010]
Bernard, J. P., Paradis, D., Marshall, D. J., et al. 2010, , 518,

[Blake et al., 1987]
Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987,
, 315, 621

[Boogert et al., 1997]
Boogert, A. C. A., Schutte, W. A., Helmich, F. P., Tielens,
A. G. G. M., & Wooden, D. H. 1997, , 317, 929

[Brown et al., 1995]
Brown, A. G. A., Hartmann, D., & Burton, W. B. 1995, , 300, 903

[Bulut et al., 2021]
Bulut, N., Roncero, O., Aguado, A., et al. 2021, , 646, A5

[Caselli et al., 2019]
Caselli, P., Sipilä, O., & Harju, J. 2019, arXiv e-prints,

[Caselli et al., 1999]
Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C.
1999, , 523, L165

[Caselli et al., 1998]
Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ,
499, 234

[Caselli et al., 2002]
Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, , 565, 331

[Cazaux et al., 2022]
Cazaux, S., Carrascosa, H., Muñoz Caro, G. M., et al. 2022, ,
657, A100

[Ceccarelli et al., 2014]
Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in
Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P.
Dullemond, & T. Henning, 859

[Ceccarelli et al., 2007]
Ceccarelli, C., Caselli, P., Herbst, E., Tielens, A. G. G. M., &
Caux, E. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt,
& K. Keil, 47

[Codella et al., 2021]
Codella, C., Bianchi, E., Podio, L., et al. 2021, , 654, A52

[Crapsi et al., 2005]
Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, , 619, 379

[Crockett et al., 2014]
Crockett, N. R., Bergin, E. A., Neill, J. L., et al. 2014, , 781,

[Cuadrado et al., 2017]
Cuadrado, S., Goicoechea, J. R., Cernicharo, J., et al. 2017, ,
603, A124

[Cummins et al., 1986]
Cummins, S. E., Linke, R. A., & Thaddeus, P. 1986, , 60, 819

[Curtis et al., 2010]
Curtis, E. I., Richer, J. S., Swift, J. J., & Williams, J. P. 2010,
, 408, 1516

[Dalgarno & Lepp, 1984]
Dalgarno, A. & Lepp, S. 1984, , 287, L47

[Davis et al., 2010]
Davis, C. J., Chrysostomou, A., Hatchell, J., et al. 2010, , 405,

[Davis et al., 2008]
Davis, C. J., Scholz, P., Lucas, P., Smith, M. D., & Adamson, A.
2008, , 387, 954

[204]
Deeyamulla, M. & Husain, D. ????

[Doddipatla et al., 2020]
Doddipatla, S., He, C., Kaiser, R. I., et al. 2020, Proceedings of the
National Academy of Science, 117, 22712

[Drozdovskaya et al., 2018]
Drozdovskaya, M. N., van Dishoeck, E. F., Jørgensen, J. K., et al.
2018, , 476, 4949

[Druard & Wakelam, 2012]
Druard, C. & Wakelam, V. 2012, , 426, 354

[Dunham et al., 2013]
Dunham, M. M., Arce, H. G., Allen, L. E., et al. 2013, , 145, 94

[Emprechtinger et al., 2009]
Emprechtinger, M., Caselli, P., Volgenau, N. H., Stutzki, J., &
Wiedner, M. C. 2009, , 493, 89

[Esplugues et al., 2013]
Esplugues, G. B., Tercero, B., Cernicharo, J., et al. 2013, , 556,

[Esplugues et al., 2014]
Esplugues, G. B., Viti, S., Goicoechea, J. R., & Cernicharo, J. 2014,
, 567, A95

[Fabricant et al., 1977]
Fabricant, B., Krieger, D., & Muenter, J. S. 1977, , 67, 1576

[Feng et al., 2019]
Feng, S., Caselli, P., Wang, K., et al. 2019, , 883, 202

[Fontani et al., 2011]
Fontani, F., Palau, A., Caselli, P., et al. 2011, , 529, L7

[Francioso et al., 2020]
Francioso, A., Baseggio, A. C., Mosca, L., & Fontana, M. 2020,
Hindawi, 8294158

[Fuente et al., 2019]
Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, , 624, A105

[Garrod et al., 2007]
Garrod, R. T., Wakelam, V., & Herbst, E. 2007, , 467, 1103

[Gatley et al., 1974]
Gatley, I., Becklin, E. E., Mattews, K., et al. 1974, , 191, L121

[Geballe et al., 1985]
Geballe, T. R., Baas, F., Greenberg, J. M., & Schutte, W. 1985, ,
146, L6

[Gerner et al., 2015]
Gerner, T., Shirley, Y. L., Beuther, H., et al. 2015, , 579, A80

[Goicoechea & Cuadrado, 2021]
Goicoechea, J. R. & Cuadrado, S. 2021, , 647, L7

[Goldsmith et al., 2008]
Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, , 680, 428

[Goldsmith & Langer, 1978]
Goldsmith, P. F. & Langer, W. D. 1978, , 222, 881

[Goodman et al., 1992]
Goodman, A. A., Jones, T. J., Lada, E. A., & Myers, P. C. 1992, ,
399, 108

[Güver & Özel, 2009]
Güver, T. & Özel, F. 2009, , 400, 2050

[Hacar et al., 2016]
Hacar, A., Alves, J., Burkert, A., & Goldsmith, P. 2016, , 591,

[Hacar et al., 2013]
Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, ,
554, A55

[Hartmann, 2002]
Hartmann, L. 2002, , 578, 914

[Hatchell et al., 2007]
Hatchell, J., Fuller, G. A., Richer, J. S., Harries, T. J., & Ladd,
E. F. 2007, , 468, 1009

[Hatchell et al., 2005]
Hatchell, J., Richer, J. S., Fuller, G. A., et al. 2005, , 440, 151

[Hatchell et al., 1998]
Hatchell, J., Thompson, M. A., Millar, T. J., & MacDonald, G. H. 1998,
, 133, 29

[Herbst & van Dishoeck, 2009]
Herbst, E. & van Dishoeck, E. F. 2009, , 47, 427

[Hirota & Yamamoto, 2006]
Hirota, T. & Yamamoto, S. 2006, , 646, 258

[Howe & Millar, 1993]
Howe, D. A. & Millar, T. J. 1993, , 262, 868

[Indriolo & McCall, 2012]
Indriolo, N. & McCall, B. J. 2012, , 745, 91

[Irvine et al., 1989]
Irvine, W. M., Friberg, P., Kaifu, N., et al. 1989, , 342, 871

[Ivlev et al., 2018]
Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O., et al. 2018, , 855,

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

[Johnstone & Bally, 1999]
Johnstone, D. & Bally, J. 1999, , 510, L49

[Johnstone & Bally, 2006]
Johnstone, D. & Bally, J. 2006, , 653, 383

[Kaifu et al., 2004]
Kaifu, N., Ohishi, M., Kawaguchi, K., et al. 2004, , 56, 69

[Kimura et al., 2020]
Kimura, H., Wada, K., Kobayashi, H., et al. 2020, , 498, 1801

[Kirk et al., 2013]
Kirk, J. M., Ward-Thompson, D., Palmeirim, P., et al. 2013, ,
432, 1424

[Knee & Sandell, 2000]
Knee, L. B. G. & Sandell, G. 2000, , 361, 671

[Kutner et al., 1976]
Kutner, M. L., Evans, N. J., I., & Tucker, K. D. 1976, , 209, 452

[Laas & Caselli, 2019]
Laas, J. C. & Caselli, P. 2019, , 624, A108

[Lada et al., 1996]
Lada, C. J., Alves, J., & Lada, E. A. 1996, , 111, 1964

[Ladd et al., 1993]
Ladd, E. F., Lada, E. A., & Myers, P. C. 1993, , 410, 168

[Ladd et al., 1994]
Ladd, E. F., Myers, P. C., & Goodman, A. A. 1994, , 433, 117

[Langer et al., 1996]
Langer, W. D., Velusamy, T., & Xie, T. 1996, , 468, L41

[Leustek, 2002]
Leustek, T. 2002, Arabidopsis Book, 1:e0017

[Linsky et al., 2006]
Linsky, J. L., Draine, B. T., Moos, H. W., et al. 2006, , 647, 1106

[Lombardi et al., 2014]
Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, , 568, C1

[Luhman et al., 2003]
Luhman, K. L., Stauffer, J. R., Muench, A. A., et al. 2003, , 593,

[Lynds, 1962]
Lynds, B. T. 1962, , 7, 1

[Maddalena et al., 1986]
Maddalena, R. J., Morris, M., Moscowitz, J., & Thaddeus, P. 1986,
, 303, 375

[Majumdar et al., 2017]
Majumdar, L., Gratier, P., Ruaud, M., et al. 2017, , 466, 4470

[Malinen et al., 2012]
Malinen, J., Juvela, M., Rawlings, M. G., et al. 2012, , 544, A50

[Marcelino et al., 2005]
Marcelino, N., Cernicharo, J., Roueff, E., Gerin, M., &
Mauersberger, R. 2005, , 620, 308

[Martín-Doménech et al., 2016]
Martín-Doménech, R., Jiménez-Serra, I., Muñoz Caro,
G. M., et al. 2016, , 585, A112

[Millar, 2005]
Millar, T. J. 2005, Astronomy and Geophysics, 46, 2.29

[Millar et al., 1989]
Millar, T. J., Bennett, A., & Herbst, E. 1989, , 340, 906

[Minh et al., 1991]
Minh, Y. C., Brewer, M. K., Irvine, W. M., Friberg, P., & Johansson,
L. E. B. 1991, , 244, 470

[Minh et al., 2018]
Minh, Y. C., Liu, H. B., Galvań-Madrid, R., et al. 2018, , 864,

[Minowa et al., 1997]
Minowa, H., Satake, M., Hirota, T., et al. 1997, , 491, L63

[Molinari & Ricci, 1986]
Molinari, A. & Ricci, R. A. 1986, From nuclei to stars. A meeting in
nuclear physics and astrophysics exploring the path opened by H. A. Bethe.
Proceedings of the International School of Physics, “Enrico Fermi”, Course
XCI, held at Varenna, Italy, 18 - 23 June 1984., Vol. 91

[Molinari et al., 2010]
Molinari, S., Swinyard, B., Bally, J., et al. 2010, , 518, L100

[Moscoso et al., 2021]
Moscoso, J., de Souza, R. S., Coc, A., & Iliadis, C. 2021, , 923,

[Myers, 2009]
Myers, P. C. 2009, , 700, 1609

[Navarro-Almaida et al., 2021]
Navarro-Almaida, D., Fuente, A., Majumdar, L., et al. 2021, , 653,

[Navarro-Almaida et al., 2020]
Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, , 637,

[Neufeld et al., 2015]
Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, , 577, A49

[Neufeld et al., 2014]
Neufeld, D. A., Gusdorf, A., Güsten, R., et al. 2014, , 781,

[Neufeld & Wolfire, 2017]
Neufeld, D. A. & Wolfire, M. G. 2017, , 845, 163

[Ohishi & Kaifu, 1998]
Ohishi, M. & Kaifu, N. 1998, Faraday Discussions, 109, 205

[Oliveira et al., 2003]
Oliveira, C. M., Hébrard, G., Howk, J. C., et al. 2003, , 587,

[Onishi et al., 2002]
Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui, Y.
2002, , 575, 950

[Padovani et al., 2013]
Padovani, M., Hennebelle, P., & Galli, D. 2013, , 560, A114

[Palla & Stahler, 2002]
Palla, F. & Stahler, S. W. 2002, , 581, 1194

[Palmeirim et al., 2013]
Palmeirim, P., André, P., Kirk, J., et al. 2013, , 550, A38

[Palumbo et al., 1997]
Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, , 479,

[Palumbo et al., 1995]
Palumbo, M. E., Tielens, A. G. G. M., & Tokunaga, A. T. 1995, , 449,

[Parise et al., 2009]
Parise, B., Leurini, S., Schilke, P., et al. 2009, , 508, 737

[Pineda et al., 2008]
Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, , 679, 481

[Pineda et al., 2010]
Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, , 712,

[Pineda et al., 2011]
Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2011, , 739, L2

[Plunkett et al., 2013]
Plunkett, A. L., Arce, H. G., Corder, S. A., et al. 2013, , 774, 22

[Prasad & Huntress, 1982]
Prasad, S. S. & Huntress, W. T., J. 1982, , 260, 590

[Punanova et al., 2018]
Punanova, A., Caselli, P., Pineda, J. E., et al. 2018, , 617, A27

[Qian et al., 2015]
Qian, L., Li, D., Offner, S., & Pan, Z. 2015, , 811, 71

[Redaelli et al., 2021]
Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, , 656,

[Rivière-Marichalar et al., 2019]
Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019,
, 628, A16

[Roberts et al., 2003]
Roberts, H., Herbst, E., & Millar, T. J. 2003, , 591, L41

[Roberts & Millar, 2000]
Roberts, H. & Millar, T. J. 2000, , 361, 388

[Roccatagliata et al., 2020]
Roccatagliata, V., Franciosini, E., Sacco, G. G., Randich, S., &
Sicilia-Aguilar, A. 2020, , 638, A85

[Rodríguez-Baras et al., 2021]
Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P.,
et al. 2021, , 648, A120

[Roueff & Gerin, 2003]
Roueff, E. & Gerin, M. 2003, , 106, 61

[Roueff et al., 2007]
Roueff, E., Parise, B., & Herbst, E. 2007, , 464, 245

[Ruaud et al., 2016]
Ruaud, M., Wakelam, V., & Hersant, F. 2016, , 459, 3756

[Ruffle et al., 1999]
Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999,
, 306, 691

[Sakai et al., 2012]
Sakai, T., Sakai, N., Furuya, K., et al. 2012, , 747, 140

[Schmalzl et al., 2010]
Schmalzl, M., Kainulainen, J., Quanz, S. P., et al. 2010, , 725,

[Shimajiri et al., 2015]
Shimajiri, Y., Kitamura, Y., Nakamura, F., et al. 2015, , 217, 7

[Shingledecker et al., 2020]
Shingledecker, C. N., Lamberts, T., Laas, J. C., et al. 2020, ,
888, 52

[Shu et al., 1987]
Shu, F. H., Adams, F. C., & Lizano, S. 1987, , 25, 23

[Sinclair et al., 1973]
Sinclair, M. W., Fourikis, N., Ribes, J. C., et al. 1973, Australian
Journal of Physics, 26, 85

[Soler, 2019]
Soler, J. D. 2019, , 629, A96

[Spezzano et al., 2022]
Spezzano, S., Fuente, A., Caselli, P., et al. 2022, , 657, A10

[Tatematsu et al., 2004]
Tatematsu, K., Umemoto, T., Kandori, R., & Sekimoto, Y. 2004, ,
606, 333

[Treviño-Morales et al., 2014]
Treviño-Morales, S. P., Pilleri, P., Fuente, A., et al. 2014, ,
569, A19

[Tsujimoto & Bland-Hawthorn, 2010]
Tsujimoto, T. & Bland-Hawthorn, J. 2010, in Light Elements in the
Universe, ed. C. Charbonnel, M. Tosi, F. Primas, & C. Chiappini,
Vol. 268, 499–500

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

[Vidal et al., 2017]
Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, ,
469, 435

[Wakelam et al., 2004]
Wakelam, V., Castets, A., Ceccarelli, C., et al. 2004, , 413, 609

[Walawender et al., 2005]
Walawender, J., Bally, J., Kirk, H., & Johnstone, D. 2005, , 130,

[Walmsley et al., 2004]
Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, ,
418, 1035

[Wiesenfeld & Faure, 2013]
Wiesenfeld, L. & Faure, A. 2013, , 432, 2573

[Wilking et al., 2004]
Wilking, B. A., Meyer, M. R., Greene, T. P., Mikhail, A., & Carlson,
G. 2004, , 127, 1131

[Wilson et al., 1970]
Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, , 161, L43

[Wootten & Mangum, 2009]
Wootten, A. & Mangum, J. 2009, in Astronomical Society of the Pacific
Conference Series, Vol. 417, Submillimeter Astrophysics and Technology: a
Symposium Honoring Thomas G. Phillips, ed. D. C. Lis, J. E. Vaillancourt,
P. F. Goldsmith, T. A. Bell, N. Z. Scoville, & J. Zmuidzinas, 219

[Yan et al., 2019]
Yan, Q.-Z., Zhang, B., Xu, Y., et al. 2019, , 624, A6

[Yang et al., 2021]
Yang, Y.-L., Sakai, N., Zhang, Y., et al. 2021, , 913, 155

[Zari et al., 2018]
Zari, E., Hashemi, H., Brown, A. G. A., Jardine, K., & de Zeeuw,
P. T. 2018, , 620, A172

[Zari et al., 2016]
Zari, E., Lombardi, M., Alves, J., Lada, C. J., & Bouy, H. 2016,
, 587, A106

[Zasowski et al., 2009]
Zasowski, G., Kemper, F., Watson, D. M., et al. 2009, , 694, 459

[Zucker et al., 2020]
Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, , 633,

[Zucker et al., 2019]
Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, , 879,