© The Authors 2024

1 Introduction

Sulfur (S) is the tenth most abundant element in the Universe and the sixth most important biogenic element (Fontani et al. 2023). It is the fifth most abundant element on our planet. Even though the chemistry related to sulfur in the interstellar medium (ISM) has been continuously investigated for many years (e.g., Pineau des Forets et al. 1986; Turner 1995; Wakelam et al. 2004; Goicoechea et al. 2006; Vastel et al. 2018; Fontani et al. 2023), it still is far from being well understood.

One of the unsolved main problems regarding the sulfur chemistry in the ISM is that its content in dense molecular regions is highly depleted compared to diffuse clouds (Laas & Caselli 2019; Rivière-Marichalar et al. 2019). The abundance of sulfur in the gas phase measured through molecular emission in dense areas of the ISM is estimated to be a few percent of the cosmic value at most (Woods et al. 2015; Fontani et al. 2023). Hily-Blant et al. (2022) found that toward early-type cores, the sulfur abundance is close to the cosmic abundance. The authors proposed that sulfur is not depleted and is atomic in these environments, which agrees with some chemical models. Then, according to the authors, cores that are chemically more evolved show sulfur depletion by factors up to 100 in their densest parts.

The principal hypothesis for explaining the depletion is that most of the sulfur is locked into dust grains in the densest regions of the ISM (Laas & Caselli 2019). This generally involves species in ices, such as organo-sulfur molecules, sulfur chains, or rings, such as S₈ (Shingledecker et al. 2020). However, only OCS has been detected on ices, and only upper limits were obtained for H₂S, which means that the main carrier of sulfur remains unknown (Hily-Blant et al. 2022, and references therein).

Recently, Fuente et al. (2023) proposed that the sulfur depletion in the gas phase might depend on the particularities of each environment. The mechanisms that cause this dependence are not known in full detail. The authors proposed that the influence of the grain charges on the chemistry might be one of the causes. Additionally, shocks associated with massive star formation could erode the grain cores, which can contribute to enhance the sulfur elemental abundance in the hosting molecular cloud.

Despite the considerable efforts performed to determine the location of the missing interstellar sulfur, the primary carrier of this element is still currently unknown in the gas and solid phase (Artur de la Villarmois et al. 2023). Hence, more studies of sulfur-bearing species toward dense regions of the ISM are required.

Additionally, as Fontani et al. (2023) pointed out, another unresolved problem regarding the S chemistry is that the formation of some S-bearing molecules detected in the gas phase is unclear. For instance, some basic reactions that initiate sulfur chemistry are highly endothermic, requiring temperatures much higher than those typically found in cold, dense molecular cores and clumps. Therefore, further studies of the detection of sulfur-bearing species in molecular cores are necessary.

In conclusion, abundance measurements of S-bearing molecules and comparisons among them toward different interstellar environments are essential for understanding the sulfur chemistry. In particular, cores in the early stages of the star-forming evolution are of notable interest to constrain the early sulfur chemistry in these environments. Based on data from the Atacama Large Millimeter Array (ALMA) data archive (see Sect. 2), we therefore present a study of several sulfur species toward 37 dense cores that are embedded in the most massive infrared-quiet molecular clumps from the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL).

2 Source sample and data

To carry out this work, we first searched the ALMA database for data containing several lines of sulfur-bearing species. The data acquired in the ALMA Band 7 fulfill this requirement. We found that project 2017.1.00914 (PI: Csengeri) consisted of observations in Band 7 toward massive infrared-quiet clumps in the inner Galaxy, that is, ATLASGAL sources that likely contain fragments, such as molecular cores probably in the early stages of star formation. Following a thorough examination of the data cubes and a survey of sulfur-bearing molecular species in the entire dataset, we selected the molecular lines listed in Table 1.

Based on the continuum maps at 0.87 mm obtained toward each ATLASGAL source, the molecular cores embedded in these clumps were identified. When multiple cores were present in the clump, the most intense core was selected. The ATLAS-GAL source names and the coordinates of the selected and analyzed molecular cores are presented in Table A.1. The systemic velocities (v_LSR), which are needed to identify the molecular lines, were obtained from the cataloged velocity of the ATLAS-GAL source, and were slightly corrected when necessary using the most intense and ubiquitous line of the selected molecules, thioformaldehyde (H₂CS). The distances were obtained from two catalogs in Vizier: ‘ATLASGAL cold high-mass clumps with NH₃’ (Wienen et al. 2012) and ‘Complete sample of Galactic clump properties’ (Urquhart et al. 2018).

The data cubes from project 2017.1.00914 (PI: Csengeri, T.; Band 7) were obtained from the ALMA Science Archive¹. The single-pointing observations for the 37 targets were carried out using the following telescope configuration: L5BL/L80BL(m): 8.7/27.5 in the 7 m array. The observed frequency range and the spectral resolution are 333.3–349.1 GHz and 1.1 MHz. The angular resolution extends from 3″.3 to 3″.8 for the whole sample. The continuum and line (each 10 km s⁻¹) sensitivity are about 1.2 and 30 mJy beam⁻¹.

It is important to remark that even though the data passed the QA2 quality level, which ensures a reliable calibration for science-ready data, the automatic pipeline-imaging process may give rise to a clean image with some artifacts. For example, an inappropriate setting of the parameters of the clean task in CASA could generate artificial dips in the spectra. Thus, we reprocessed some raw data using the CASA 4.5.1 and 4.7.2 versions and the calibration pipeline scripts. Particular care was taken with the different parameters of the task clean. The images and spectra from our data reprocessing, following multiple runs of the clean task with varied parameters, matched the archival data closely. Consequently, we chose to rely on the archival data.

In some cases, when it was appreciated that the spectral lines were mounted on some continuum emission, the task imcontsub in CASA was used to subtract it. A first-order polynomial was enough for an appropriate continuum subtraction after we selected frequency ranges without molecular line emission in each spectral window.

3 Results

After we identified the cores within each ATLASGAL source, the spectra from each of them were extracted at their continuum peak position using a beam-size circular region. In all cases, the peaks of the molecular emissions coincide with those of the continuum emission. After the six molecular lines were identified, Gaussian fittings to each line were applied to obtain the intensity peaks, the full width half maximum (FWHM) ∆v, and the integrated line intensity (i.e., the line flux, W). In an appendix (available on Zenodo), all these parameters are presented in Tables F1 and F2, and the spectra and Gaussian fittings of a particular core are shown as an example of this procedure.

None of the analyzed lines showed signs of absorption or flattening. All of them were well fit by Gaussians, indicating that opacity effects, if any, likely had little impact on our analysis.

Additionally, the lines present well-defined Gaussian shapes in general, and no signature of large spectral wings. The SO and SO₂ lines present some signatures of these wings more frequently, and, to a lesser extent, NS and H₂CS show them as well. When the lines had spectral wings, only the main central component was fit (see Figs. G1, G2, and G3 available on Zenodo). These small spectral wings might be due to some dynamical process in the gas or to the presence of other molecular lines, probably complex organic molecules (COMs), not identified, and beyond the scope of this work.

To estimate the column densities of each molecular species, assuming local thermodynamic equilibrium (LTE) and following Artur de la Villarmois et al. (2023), we first derive the column densities of the upper levels, N_u (in cm⁻²), of each transition from $${\text{N}}_{\text{u}}=2375\times {10}^6\left(\frac{\text{\hspace{0.17em}W}}{1{\text{\hspace{0.17em}Jy\hspace{0.17em}kms}}^{-1}}\right)\left(\frac{1{\text{s}}^{-1}}{{\text{A}}_{\text{ul}}}\right)\left(\frac{{\text{arcsec}}^2}{{\mathrm{\Theta }}_{{\text{area}}^2}}\right),$$(1)

where W is the measured line flux in Jy beam⁻¹ km s⁻¹ (see Tables F1 and F2 available on Zenodo), A_ul is the line Einstein coefficient (in s⁻¹) (see Table 2), and Θ_area is the integration area in arcsec², which was taken to be equal to the size of the beam in all cases.

Then, the total column density of each molecule was obtained from $${\text{N}}_{\text{tot}}=\frac{{\text{N}}_{\text{u}}\text{\hspace{0.17em}Q}\left({\text{T}}_{\text{ex}}\right)\mathrm{\hspace{0.17em}exp}\left({\text{E}}_{\text{u}}/{\text{T}}_{\text{ex}}\right)}{{\text{g}}_{\text{u}}},$$(2)

where Q(T_ex) is the rotational partition function, T_ex is the excitation temperature (in K), 𝑔_u the upper-level degeneracy, and E_up (in K) the upper energy of the line transition. The parameters are shown in Table 2. A T_ex of 20 K was uniformly assumed for the excitation of all lines in the whole sample. The column densities we obtained and their respective errors are presented in Table D.1.

We assumed for these calculations that all transitions are optically thin, which is supported by the lack of high-opacity evidence, such as absorption dips or flattened components. However, because we cannot completely rule out this phenomenon, τ values were obtained for all the lines in the whole sample of cores (see Appendix B). The optically thin assumption is supported by the values we obtained.

We measured the integrated flux (S_int) from a beam toward the peak of each core from the 0.8 mm continuum emission (Col. 2 in Table D.2). Using these values, we derived the H₂ column density (N(H₂)) of each core, which was used to estimate the molecular abundances. We used the following equation (Kauffmann et al. 2008): $$\begin{array}{l}\text{N}\left({\text{H}}_2\right)=2.02\times {10}^{20}{\text{cm}}^{-2}\left({\text{e}}^{1.439{(\lambda /\text{mm})}^{-1}{(T/10\text{\hspace{0.17em}K})}^{-1}}-1\right)\\ \times {\left(\frac{\lambda }{\text{mm}}\right)}^3{\left(\frac{{\kappa }_{\nu }}{0.01{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}\left(\frac{{\text{S}}_{\text{int}}}{\text{mJy}}\right)\left(\frac{{\theta }_{\text{HPBW}}}{10\mathrm{arcsec}}\right),\end{array}$$(3)

where T is the dust temperature, and κ_ν is the dust opacity per gram of matter at 870 µm, for which we adopted the value of 0.0185 cm²g⁻¹ (Csengeri et al. 2017, and references therein), and θ_HPBW is the beam size. We assumed thermal coupling between dust and gas (i.e., T_dust=T_kin). To obtain a value of T_kin for each core, we searched for other molecular lines in the data set that might be used to estimate temperatures. We found that all cores have a well-defined emission of the CH₃OH 7(1,7)–6(1,6)++ and 12(1,11)–12(0,12)-+ lines at 335.582 and 336.865 GHz (for an example of the spectra of these lines, see Fig. G4 available on Zenodo). Then, by integrating the emission of these lines (presented in Table C.1), we obtained through a typical rotational diagram procedure (Goldsmith & Langer 1999; see Appendix C) the rotational temperatures that were assumed to be the T_kin of each core (Col. 4, and their errors in Col. 5 of Table D.2). The N(H₂) we obtained for each core with the errors are presented in Cols. 6 and 7 of Table D.2.

Additionally, to roughly characterize the cores, we calculated the mass using the continuum fluxes measured from a beam toward the peaks of the cores. The mass was estimated from (Kauffmann et al. 2008): $$\begin{array}{l}{\text{M}}_{\text{gas}}=0.12{\text{\hspace{0.17em}M}}_{\odot }\left[\exp \left(\frac{1.439}{(\lambda /\text{mm})(\text{T}/10\text{\hspace{0.17em}K})}\right)-1\right]\\ \times {\left(\frac{{\kappa }_{\nu }}{0.01{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}\left(\frac{{\text{S}}_{\text{int}}}{\text{mJy}}\right){\left(\frac{\text{d}}{100\text{\hspace{0.17em}pc}}\right)}^2{\left(\frac{\lambda }{\text{\hspace{0.17em}mm}}\right)}^3.\end{array}$$(4)

The volume densities were derived assuming spherical symmetries with a diameter of the beam (i.e., the density was obtained for a uniform probe volume at the peak of each core). The masses and densities are included in the last columns of Table D.2. In these cases, errors are not quoted because not all cores had errors in their distances, and these parameters were taken as a rough core characterization.

3.1 Molecular abundances

By defining the abundance of a molecular species as X(molecule) = N(molecule)/N(H₂), we calculated these values and compared them. Additionally, they are compared with the physical parameters for each core.

The temperature can have deep implications for the chemistry (e.g., Gorai et al. 2020). In Fig. 1, we therefore present the abundances versus temperature for each core and show their relations with the density on a color scale. Figure 2 displays the X(molecule) versus N(molecule) relations with a color scale representing the temperature. Finally, Fig. E.1 shows the relation among all the sulfur-bearing molecular abundances with the color scale representing the temperature. All figures present the results of a linear fit with a determination coefficient R².

It was proposed that the abundance ratio X(SO₂)/X(SO) might be used as a chemical clock (e.g., Wakelam et al. 2004, 2011), and we therefore calculated this ratio for our sample of cores. We obtained a media of 1.6, with a maximum value of 9.29 and a minimum of 0.06.

3.2 Molecular line widths

The molecular line widths ∆v (FWHM) are usually associated with the kinematics of the molecular gas, which can reflect the consequences induced by multiple events, such as shocks, outflows, rotation, and turbulence. When we assume that the line widths in our sample of cores primarily result from turbulence related to early star formation (e.g., Sanhueza et al. 2012; Yu & Wang 2015), a comparison of this parameter for different molecular species could provide valuable physical and chemical insights. For instance, it could give information about the location or distribution of the molecular species in different gaseous layers within the core. As reported and pointed out by Fontani et al. (2023), these comparisons may allow us to understand the species that are more likely associated with similar gas, and to search for variations among the molecular lines.

In Fig. 3, we show the velocity width distribution (FWHM) of the six molecular species for the full sample of sources. We identify H₂CS as the molecule with the narrower line widths (the distribution peaks at ~5 km s⁻¹ ), so that we can use it as a tracer of the quiescent-gas component. On the other hand, SO₂ displays wider widths, indicating that it traces the most turbulent gas.

In Fig. E.2, we present the comparisons of all the line widths of the analyzed molecular species. The color code represents the temperature. To better discern whether the emission of some molecular species is wider than the other, we also plot the unity line, and the Pearson coefficient is included in each panel to evaluate the correlation among the line width points.

Fig. 1 Molecular abundances (X(x) = N(x)/N(H2)) (in logarithmic scale) vs. kinetic temperature. The lines are the result of linear fits whose results are included in the top left corner of each panel. The colors correspond to the volume density of the cores (in cm−3), and the scale is presented in the bar at the right of each panel.

Fig. 2 Molecular abundances (X(x) = N(x)/N(H2)) vs. column density (N(x)) presented on a logarithmic scale. The lines are the result of linear fits whose results are included in the top left corner of each panel. The colors correspond to the kinetic temperature measured in the cores, and the scale is presented in the bar at the right of each panel.

4 Discussion

We first discuss the results we obtained from the comparisons among the molecular abundances, column densities, and physical parameters such as the kinetic temperature and the volume density. The relation between the abundances and temperature (see Fig. 1) shows that the abundance of all analyzed sulfur-bearing species increases with the temperature, that is, the abundances are positively correlated with T_k. This agrees with the results by Fontani et al. (2023) for the SO, SO₂, SO⁺, and H₂CS in a sample of high-mass starless cores, high-mass proto-stellar objects, and ultracompact HII regions that were observed with a single-dish telescope.

Based on the linear fittings (see the obtained slopes in Fig. 1), the increment of the molecular abundances with T_k is almost the same for all cases. Thus, we suggest that the chemistry involved in each of the analyzed sulfur-bearing species may have a similar dependence on T_k. This result, obtained for a temperature range from about 20 to 100 K, supports the hypothesis that sulfur may be frozen in the dust grains in dense and cold interstellar regions (Laas & Caselli 2019), and hence, as the temperature rises, the abundance of sulfur-containing molecules increases as well. We cannot discard that the sulfur is already in the gaseous phase and that the temperature increase activates chemical reactions that favor the production of these sulfur-bearing species. The relation we obtained for H₂CS is very similar to what was obtained by Chen et al. (2024) from a large sample of massive protoclusters with diverse physical conditions, but in general, they are warmer than the cores investigated here. The core density does not appear to have an effect on the abundance-temperature relation, at least in the density range estimated here, approximately between 10⁶ and 10⁷ cm⁻³.

The relations between the abundance of a molecular species and its column density (Fig. 2) show an interesting result. The abundance increases with the column density in all cases, which indicates that the column density of the sulfur-containing molecule increases more strongly than N(H₂). However, we cannot rule out that N(H₂) might be underestimated as a result of dust opacity effects. On the other hand, the relation we obtained for thioformaldehyde from the linear fitting is almost the same as was obtained by Chen et al. (2024) in their sample of massive protoclusters.

Additionally, the panels displayed in Fig. 2 show that the temperature plays an important role in the increasing abundances. Warm cores (orange and yellow points) lie in general above the fitted lines in all cases. This indicates that warm cores have higher molecular abundances than cold cores at the same column density. Our results show that the results by Chen et al. (2024) for the thioformaldehyde also apply to the other sulfur-bearing species, which could be an important point to be taken into account in chemical models for studying the formation of these molecules.

Fig. 3 Distribution of the measured line widths ∆v (FWHM) of the sample. The vertical dashed lines are the average values of ∆v of each molecule. In the case of SO and SO+, these averages are almost the same, and the dashed lines appear to overlap. The horizontal axis is truncated at 20 km s−1 for a better visualization.

4.1 Comparisons of the molecular abundances

As presented in Fig. E.1, all comparisons of the abundances of the analyzed sulfur-bearing species present a positive correlation. Additionally, the increase in X versus X correlates with the kinetic temperature. This phenomenon shows that as the temperature increases, the abundances also increase, which supports, as mentioned above, the hypothesis that the sulfur is trapped in dust-grain ice mantles in the cold regions. Moreover, based on the good, and in some cases very good, correlation among the abundances (a determination coefficient R² higher than 0.75, and in some cases about 0.9 in the linear regressions), we can conclude that the chemical processes involved in the formation of all analyzed sulfur-bearing species are probably related.

Interestingly, the correlation between X(SO⁺) and X(SO) is lowest for all pairs of compared abundances. Thus, when we take into account that a possible route of formation of the SO⁺ in this type of region is through the SO ionization via cosmic rays (Rivilla et al. 2022), it is likely that as the SO abundance increases, some of these molecules contribute to generating SO⁺, which would be reflected in a loss of correlation between the abundances of both molecules. Moreover, the fact that the average ∆v of both molecules are the same (see Fig. 3) suggests that both species coincide in the same gas layer.

Fontani et al. (2023) presented relations of the abundances of SO⁺, NS, SO₂, and H₂CS with SO and obtained Pearson correlation coefficients (ρ) of 0.62, 0.70, 0.95, and 0.73, respectively. In our case, the corresponding Pearson coefficients ($\rho =\sqrt{{\text{R}}^2}$, where R² is the determination coefficient of the line regressions presented in Fig. E.1) for these abundance relations are 0.83, 0.87, 0.93, and 0.93, respectively. All are higher than those presented by Fontani and collaborators, except for ρ of the X(SO₂)/X(SO) relation. The two values are almost the same in this case. Our sample contains sources in earlier evolutionary stages than the sample studied by Fontani et al. (2023), which might explain the better correlation among the abundances we found. This might indicate that these molecules can have a similar chemical origin in these early stages (e.g., the sulfur freeze-out from the dust grains). These correlations among the abundances would then decrease as the cores evolve because the chemistry becomes more complex.

Finally, the obtained values of the X(SO₂)/X(SO) ratio for our sample of cores agree with what is yielded by chemical models of a core that is chemically evolved about 10⁵ yr for a temperature of 60 K (see Wakelam et al. 2011). The average measured temperature of our sample is 65 K. Even though we cannot know the chemical ages of the cores in the analyzed sample accurately, we can suggest that our observational results may support the use of this abundance ratio as a chemical clock, as was proposed by the authors.

4.2 Line width comparisons

The obtained ∆v values for the sulfur-bearing molecules are quite similar to those obtained, also from interferometric observations of sulfur species, toward others molecular cores (el Akel et al. 2022; Chen et al. 2024). In this large sample of sources, we must be careful with comparisons with other sources in the literature in different stages of evolution. The clumps in which the analyzed cores are embedded are cataloged as early star-forming regions. However, it is likely that this sample of clumps is not fully uniform. It might exhibit variations in evolutionary stage. As mention above, the ∆v (FWHM) can be used to analyze the location or distribution of the molecular species in different gaseous layers with different dynamics within the cores (Sanhueza et al. 2012; Yu & Wang 2015; Fontani et al. 2023; Martinez & Paron 2024).

Figure E.2 shows that the line widths ∆v (FWHM) of all molecular species are positively correlated, and although it is not well defined, this correlation is in general accompanied by a temperature increase. None of the displayed ∆v comparisons are close to the y = x locus, and only the relation between the line widths of ³⁴SO and NS approaches the unity line slightly. This behavior indicates that the different molecules trace gas with different dynamics or turbulence conditions in general.

In particular, the correlation between H₂CS and NS (ρ = 0.93) is excellent, but the NS line widths are systematically larger. When we assume that the ∆v in our sample of cores mainly reflects turbulence, this would indicate that NS traces the same gaseous layer as H₂CS, but also receives some contribution from more turbulent material. This agrees with the fact that the correlation between the abundances of these two molecules is one of the highest in the sample (ρ = 0.97; see Fig. E.1). The correlations of ∆v of H₂CS and NS with the O-bearing species are generally poor (ρ < 0.7), except for the relation between NS and SO, which have ρ = 0.78. We note that the four O-bearing species have larger ∆v than the other molecular species, indicating that they may be associated with warmer and more turbulent gas. Figs. 3 and E.2 show that the oxygen-containing species trace the most turbulent gas. However, the correlations among them are quite poor, indicating that they may arise from different gaseous layers in the cores. In particular, SO₂ is probably located in the innermost region, which is more affected by a possible incipient central compact source, while SO and ³⁴SO could be located farther away. The nature of SO⁺ is elusive. It presents the poorest correlation with the remaining oxygen-containing species. In addition, its ∆v are generally larger than those of ³⁴SO, but this is not evident when comparing with SO and SO₂.

We suggest that H₂CS and NS are associated with quiescent gas, likely from the outer envelope of the cores, and NS probably located somewhat more internally than H₂CS. The comparison between the FWHM ∆v of these species presents the best correlation, which coincides with the high correlation in the comparison of their abundances. This indicates that the physics and chemistry traced by the two molecular species is probably similar.

The oxygen-containing species may trace inner gas and are probably more affected by incipient star-forming processes. In particular, SO₂ traces the warmer and innermost regions. Our results agree well with those obtained by Fontani et al. (2023), who found that the C- and N-bearing species are associated with the quiescent gas in the extended envelope of the sources, while O-bearing species are located in inner regions of the cores. The phenomenon that S- and O-bearing molecules are associated with more turbulent gas was also observed by other authors (see, e.g., Artur de la Villarmois et al. 2019; Garufi et al. 2022; Artur de la Villarmois et al. 2023).

5 Summary and conclusions

We investigated the sulfur-bearing molecules ³⁴SO, SO₂, NS, SO, SO⁺, and H₂CS toward 37 molecular cores that are embedded in a sample of infrared-quiet ATLASGAL sources, presumably in very early stages of star formation. Based on the submillimeter continuum emission and methanol lines, the cores were characterized in density and temperature, and the corresponding H₂ column densities were derived. The high detection rate of the sulfur-bearing species allowed us to perform a statistical study that is useful for investigating the early sulfur chemistry in dense regions of the ISM. The main results and conclusions are summarized below.

(1) We observed that all the estimated abundances of the analyzed sulfur-bearing species increase in a similar way with increasing gas temperature. This supports the hypothesis that sulfur may be frozen in the dust grains in dense and cold interstellar regions, and hence, as the temperature rises, the abundance of sulfur-containing molecules also increases. It cannot be discarded that the sulfur is already present in the gaseous phase and that the temperature increase activates some new chemical routes of formation of these sulfur-bearing species that increase their abundances. Additionally, we suggest that the chemistry involved in the formation of each of the analyzed sulfur-bearing species may have a similar dependence on T_k, at least in the range from 20 to 100 K.

(2) We found that the warmer cores of our sample have a higher molecular abundances than the colder cores at the same column density for all molecular species. This complements a previous study of thioformaldehyde and indicates that this behavior between the abundances and column densities also applies to the other sulfur-bearing species. This might be a point to be taken into account in chemical models studying the formation of these molecules.

(3) The molecular abundance comparisons presented high correlations in general. When we compare our results with previous similar works done toward more evolved sources, we can conclude that the sulfur-bearing species we analyzed probably have the same chemical origin in the very early stages of the cores. This explains the high correlations. As the cores evolve, the correlations among the molecular abundances decrease, indicating that chemistry becomes more complex.

(4) Based on the analysis of the FWHM ∆v of the molecular lines, we suggest that the analyzed molecules with oxygen content (³⁴SO, SO, SO⁺, and SO₂) may be associated with warmer and more turbulent gas than the other molecules. Our results indicate that the H₂CS and NS are associated with more quiescent gas, probably in the external envelopes of the cores, and the O-bearing species are probably located in the inner regions.

(5) In particular, NS and H₂CS are strongly correlated in the FWHM ∆v and abundance comparisons. This may indicate that the physics and chemistry traced by the two molecular species are similar.

Finally, we point out that we presented a work with an important piece of information regarding the sulfur chemistry in the earliest stages of the molecular cores, probably pre-stellar cores. This work not only complements recent previous similar works done toward more evolved sources, but also provides quantitative information about abundances that might be useful in chemical models for explaining the sulfur chemistry in the ISM.

Data availability

Some appendices are available at Zenodo: https://doi.org/10.5281/zenodo.14036845

Acknowledgements

We would like to thank to the anonymous referee for her/his useful suggestions and corrections that improved our manuscript. N.C.M. is a doctoral fellow of CONICET. M.O., A.P., and S.P. are members of the Carrera del Investigador Científico of CONICET, Argentina. A.A., M.B., C.C., and T.H. are part of our citizen science project Ciencia Popular². This work was partially supported by the Argentina grants PIP 2021 11220200100012 and PICT 2021-GRF-TII-00061 awarded by CONICET and ANPCYT. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00914.S ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Appendix A Source sample

Table A.1 presents the ATLASGAL sources sample and the analyzed cores.

Appendix B Calculation of optical depths

To evaluate the optically thin assumption used to calculate the column densities, we computed the optical depth (τ) of the molecular lines in all cores from the following equation (e.g. Mangum & Shirley 2015): $$\tau =-\ln \left(1-\frac{{\text{T}}_{\text{peak\hspace{0.17em}}}}{\varphi \left(\text{J}\left({\text{T}}_{\text{ex}}\right)-\text{J}\left({\text{T}}_{\text{bg}}\right)\right)}\right),$$(B.1)

where T_peak is the peak of the emission (values in the ‘Peak’ columns in Tables F1 and F2 converted to K; available on Zenodo³), ϕ the beam filling factor assumed to be 1, and J(T) = hν/k (exp(hν/kT)−1)⁻¹ with Tex = 20 K, and T_bg = 2.73 K. The calculated values of T are presented in Table B.1.

Almost all values of τ are << 1, showing that the lines are indeed optically thin. Some values slightly larger than the unity (see cores 29 and 30 for the H₂CS and SO) could affect to the calculated column density in a factor about 2 (using the optical depth correction of τ/(1 −exp(−τ)), which does not change any result.

Appendix C Calculation of T_rot from the CH₃OH

Table C.1 presents the integrated emission from the CH₃OH 7(1,7)–6(1, 6)++ (at 335.582 GHz, E_u = 78.97 K) and 12(1,11)-12(0,12)−+ (at 336.865 GHz, E_u = 197.07 K) lines used to estimate temperatures. Such integrated intensities, obtained from Gaussian fittings (for an example see Fig. G4 available on Zenodo³), were used to construct rotational diagrams following Goldsmith & Langer (1999). By assuming LTE conditions, optically thin lines, and a beam filling factor equal to the unity, we can estimate the rotational temperatures (T_rot). This analysis is based on a derivation of the Boltzmann equation, $$\ln \left(\frac{{\text{N}}_{\text{u}}}{{\text{g}}_{\text{u}}}\right)=\ln \left(\frac{{\text{N}}_{\text{tot}}}{{\text{Q}}_{\text{rot}}}\right)-\frac{{\text{E}}_{\text{u}}}{{\text{kT}}_{\text{rot}}},$$(C.1)

where N_u represents the molecular column density of the upper level of the transition, g_u the total degeneracy of the upper level, N_tot the total column density of the molecule, Q_rot the rotational partition function, and k the Boltzmann constant.

Following Miao et al. (1995), for interferometric observations, the left-hand side of Eq. C.1 can be estimated from $$\ln \left(\frac{{\text{N}}_{\text{u}}}{{\text{g}}_{\text{u}}}\right)=\ln \left(\frac{2.04\times {10}^{20}}{{\theta }_{\text{a}}{\theta }_{\text{b}}}\frac{W}{{\text{g}}_{\text{k}}{\text{g}}_1{\nu }_0^3{\text{S}}_{\text{u}1}{\mu }_0^2}\right),$$(C.2)

where θ_a and θ_b are the major and minor axes of the clean beam, respectively, W is the integrated intensity, g_k and g_l are the K-level and reduced nuclear spin degeneracies, ν0 is the rest frequency of the line, Sul is the line strength of the transition, and µ₀ is the permanent dipole moment of the molecule. The free parameters, (N_tot/Q_rot) and T_rot were determined by a linear fitting to Eq. C.1. The obtained T_rot values are those included in Table D.2 as ${\text{T}}_{{\text{CH}}_3\text{OH}}$.

The T_rot was estimated by assuming optically thin lines (e.g. van der Tak et al. 2000). However, we consider the possibility of high optical depths in the used CH₃OH lines. Thus, following Goldsmith & Langer (1999) and Purcell et al. (2006), we compute the opacity correction factors for the N_u of both lines. We observe that in the cases that it makes sense to apply such a correction, it modifies the N_u value in a similar way for both transitions, and therefore it does not affect the T_rot estimation.

Appendix D Obtained parameters

Table D.1 displays the obtained column densities of the sulfur-bearing molecular species, and Table D.2 presents some measured and derived core parameters.

Appendix E Molecular abundances and measured ∆v FHWM relations

Figures E.1 and E.2 present the relations among the molecular abundances and the measured ∆v FHWM of the sulfur-bearing species, respectively.

Fig. E.1 Comparison of the molecular abundances of all analyzed molecules, presented in logarithmic scale. The lines are the result of linear fittings whose results are included at the top left corner of each panel. The colors correspond to the kinetic temperature measured in the cores, and the scale is presented in the bar at the right of each panel.

Fig. E.2 Comparison of the ∆v FWHM of the molecular lines. The points are color-coded based on the kinetic temperature measured in the cores. The dashed gray line represents the y = x locus. The Pearson correlation coefficient (ρ) is displayed in each panel.

References

All Tables

All Figures

	Fig. 1 Molecular abundances (X(x) = N(x)/N(H2)) (in logarithmic scale) vs. kinetic temperature. The lines are the result of linear fits whose results are included in the top left corner of each panel. The colors correspond to the volume density of the cores (in cm−3), and the scale is presented in the bar at the right of each panel.
In the text

	Fig. 2 Molecular abundances (X(x) = N(x)/N(H2)) vs. column density (N(x)) presented on a logarithmic scale. The lines are the result of linear fits whose results are included in the top left corner of each panel. The colors correspond to the kinetic temperature measured in the cores, and the scale is presented in the bar at the right of each panel.
In the text

	Fig. 3 Distribution of the measured line widths ∆v (FWHM) of the sample. The vertical dashed lines are the average values of ∆v of each molecule. In the case of SO and SO+, these averages are almost the same, and the dashed lines appear to overlap. The horizontal axis is truncated at 20 km s−1 for a better visualization.
In the text

	Fig. E.1 Comparison of the molecular abundances of all analyzed molecules, presented in logarithmic scale. The lines are the result of linear fittings whose results are included at the top left corner of each panel. The colors correspond to the kinetic temperature measured in the cores, and the scale is presented in the bar at the right of each panel.
In the text

	Fig. E.2 Comparison of the ∆v FWHM of the molecular lines. The points are color-coded based on the kinetic temperature measured in the cores. The dashed gray line represents the y = x locus. The Pearson correlation coefficient (ρ) is displayed in each panel.
In the text