Sulfurbearing molecules are commonly detected in dense cores within starforming regions but the total sulfur budget is significantly lower when compared to the interstellar medium value. The properties of sulfurbearing molecules are not well understood due to the absence of large sample studies with uniform observational configurations. To deepen our understanding of this subject we conducted a study using Atacama Large Millimetersubmillimeter Array 870 m observations of 11 massive protoclusters. By checking the spectra of 248 dense cores in 11 massive protoclusters a total of 10 sulfurbearing species CS SO HSUB2SUBCS NS SOSUB2SUB SUP33SUPSO SUP34SUPSOSUB2SUB SUP33SUPSOSUB2SUB SOSUP18SUPO and OCSUP34SUPS were identified. The parameters including systemic velocities line widths gas temperatures column densities and abundances were derived. Our results indicate that SO appears to be more easily detected in a wider range of physical environments than HSUB2SUBCS despite these two species showing similarities in gas distributions and abundances. Molecules SUP34SUPSOSUB2SUB and HSUB2SUBCS are good tracers of the temperature of sulfurbearing species in which HSUB2SUBCS traces the outer warm envelope and SUP34SUPSOSUB2SUB is associated with hightemperature central regions. Highmass starforming feedback outflow and other nonthermal motions significantly elevates the sulfurbearing molecular abundances and detection rates specifically for SOSUB2SUB and SO. A positive correlation between the SOSUB2SUB abundance increasing factor F and temperatures suggests that SOSUB2SUB could serve as a sulfur reservoir on the grain mantles of dense cores and then can be desorbed from dust to gas phase as the temperature rises. This work shows the importance of a large and unbiased survey to understand the sulfur depletion in dense cores.
['Molecular clouds', 'Star formation', 'Interstellar line emission', 'Submillimeter astronomy', '1072', '1569', '844', '1647', 'Astrophysics - Astrophysics of Galaxies']
A Survey of Sulfurbearing Molecular Lines toward the Dense Cores in 11 Massive Protoclusters
{'A survey of sulfur-bearing molecular lines toward the dense cores in eleven massive protoclusters': "Mengyao Tang , 1 Sheng-Li Qin , 2 Tie Liu , 3 Luis A. Zapata , 4 Xunchuan Liu , 3 Yaping Peng , 5 Fengwei Xu , 6, 7 Chao Zhang , 8 and Ken'ichi Tatematsu 9, 10  \n1 Institute of Astrophysics, School of Physics and Electronic Science, Chuxiong Normal University, Chuxiong 675000, PR China 2 School of Physics and Astronomy, Yunnan University, Kunming, 650091, PR China  \n3 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, PR China  \n4 Instituto de Radioastronom'ıa y Astrof'ısica, Universidad Nacional Aut'onoma de M'exico, P.O. Box 3-72, 58090, Morelia, Michoac'an, M'exico  \n- 5 Department of Physics, Faculty of Science, Kunming University of Science and Technology, Kunming 650500, PR China  \n6 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, PR China  \n7 Department of Astronomy, School of Physics, Peking University, Beijing, 100871, PR China  \n- 8 Institute of Astronomy and Astrophysics, School of Mathematics and Physics, Anqing Normal University, Anqing, PR China\n- 9 Nobeyama Radio Observatory, National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 462-2 Nobeyama, Minamimaki, Minamisaku, Nagano 384-1305, Japan  \n10 Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan", 'ABSTRACT': 'Sulfur-bearing molecules are commonly detected in dense cores within star-forming regions, but the total sulfur budget is significantly low, when compared to the interstellar medium (ISM) value. The properties of sulfur-bearing molecules are not well understood due to the absence of large sample studies with uniform observational configurations. To deepen our understanding of this subject, we conducted a study using ALMA 870 µ m observations of 11 massive protoclusters. By checking the spectra of 248 dense cores in 11 massive protoclusters, a total of 10 sulfur-bearing species (CS, SO, H 2 CS, NS, SO 2 , 33 SO, 34 SO 2 , 33 SO 2 , SO 18 O, OC 34 S) were identified. The parameters including systemic velocities, line widths, gas temperatures, column densities, and abundances were derived. Our results indicate that SO appears to be more easily detected in a wider range of physical environments than H 2 CS, despite these two species show similarities in gas distributions and abundances. 34 SO 2 and H 2 CS are good tracers of the temperature of sulfur-bearing species, in which H 2 CS traces the outer warm envelope and 34 SO 2 is associated with high-temperature central-regions. High-mass star-forming feedback (outflow and other non-thermal motions) significantly elevates the sulfur-bearing molecular abundances and detection rates specifically for SO 2 and SO. A positive correlation between the SO 2 abundance increasing factor ( F ) and temperatures suggests that SO 2 could serve as a sulfur reservoir on the grain mantles of dense cores and then can be desorbed from dust to gas phase as the temperature rises. This work shows the importance of a large and unbiased survey to understand the sulfur depletion in dense cores.  \nKeywords: Molecular clouds(1072) - Star formation(1569) - Interstellar line emission(844) - Submillimeter astronomy(1647)', '1. INTRODUCTION': "Sulfur is the tenth most abundant element in our galaxy, and sulfur-bearing (S-bearing, hereafter) molecules are widely detected in different environments  \nof star formation, e.g., circumstellar envelopes (van der Tak et al. 2003; Zapata et al. 2019), prestellar core (Vastel et al. 2018), hot core (Charnley 1997; Hatchell et al. 1998; Wilson & Rood 1994; Vidal & Wakelam 2018), and shock regions (Bachiller & P'erez Guti'errez 1997; Bachiller et al. 2001; Wakelam et al. 2005). Consequently, S-bearing molecules have attracted broad attention and have generally been adopted as molecular  \nprobes to investigate the properties of star-forming regions (Artur de la Villarmois et al. 2019; Oya et al. 2019; Tychoniec et al. 2021; Artur de la Villarmois et al. 2023).  \nHowever, the properties of S-bearing species are still unclear, a persistent unsolved problem in the research on sulfur is so called sulfur-depletion problem. Since the 1970s, the sulfur depletion problem became evident because it was found that the observed abundances of S-bearing molecules in dense regions generally lower than diffused regions (S/H ∼ 10 -5 ) (Penzias et al. 1971; Oppenheimer & Dalgarno 1974; Tieftrunk et al. 1994; Palumbo et al. 1995). Since then, sulfur depletion has always been a mystery in the research field of astronomy. Ruffle et al. (1999) suggested that sulfur predominantly existed as S + in translucent gas, frozen quickly onto negatively charged grains during cloud collapse. Chemical models suggest that hydrogenation is most effective on grains, and H 2 S is a main sulfur reservoir in grain mantles (Charnley 1997). However, H 2 S has not been identified in the solid phase, with its gas phase abundance estimated to be ∼ 10 -7 (e.g. Minh et al. 1990; van Dishoeck &Blake 1998), notably lower than sulfur levels in diffuse regions. Over the last two decades, various potential candidates for solid sulfur reservoirs have suggested, including hydrated sulphuric acid (H 2 SO 4 · n H 2 O), H 2 S n , sulfide minerals (FeS and MgS), sulfur chains (S n ), organo-sulfur species (Keller et al. 2002; Hony et al. 2002; Scappini et al. 2003; Wakelam et al. 2004; Jim'enezEscobar & Mu˜noz Caro 2011; Druard & Wakelam 2012; Vidal et al. 2017; Laas & Caselli 2019; Kama et al. 2019). Regrettably, detection in interstellar ices has only been successful for OCS (Geballe et al. 1985; Palumbo et al. 1995; McClure et al. 2023) and SO 2 (Boogert et al. 1997; Yang et al. 2022; McClure et al. 2023; Rocha et al. 2024),  \nleaving sulfur depletion as an ongoing enigma.  \nThe bulk of the sulfur reservoir is proposed to be locked onto grain mantels during the collapsing process (Ruffle et al. 1999; Woods et al. 2015; Vidal et al. 2017; Navarro-Almaida et al. 2020), and is subsequently released into gas phase when the temperature exceeds 100 K. Therefore, the collapsing regions could be ideal for investigating sulfur depletion. The first aim of this study is to investigate the sulfur depletion problem in massive star-forming regions. Therefore, we initiated a study using the data of an ALMA survey of 11 massive protoclusters in band 7. These 11 protoclusters, situated at distances ranging from 2.6 to 7.6 kpc, were recognized as infalling clumps with'blue-profile' signatures (Liu et al. 2016; Yue et al. 2021). Moreover, observations of other star-forming activities such as outflows (Baug et al. 2021) and H II regions (Xu et al. 2024) suggest that these clumps are active star-forming regions.  \nThe second main purpose of this study is to explore the similarities and differences among properties of Sbearing species in various dense cores. Despite numerous works conducted in past decades, most are performed by single-dish telescopes or from different observation projects with different observational frequency ranges, angular resolutions, and sensitivities. Recent observations of C 34 S and SO 2 in 10 sources (Artur de la Villarmois et al. 2019) and CS, SO, SO 2 , and 34 SO in 50 sources (Artur de la Villarmois et al. 2023) have revealed interesting statistical findings, suggesting that conducting statistical analyses on larger sample sizes could yield more valuable insights. However, their studies focused on low-mass star-forming regions like Ophiucus and Perseus with a limited number of detected species and samples. The lack of statistical results from extensive samples of dense cores in star-forming regions, particularly high-mass star-forming regions, poses a significant obstacle to advancing research on S-bearing species. Therefore, a statistical investigation involving several hundred dense cores spanning diverse physical environments and evolutionary stages could greatly enhance comprehension of S-bearing species.  \nIn our seuveyed 11 protoclusters, a total of 248 dense cores were identified before (Xu et al. 2024). These 248 dense cores cover diverse physical environments, including high-mass, low-mass, prestellar, protostellar, with/without outflows, and hot core/corinos (Qin et al. 2022; Xu et al. 2024; Liu et al. 2023). They serve as ideal subjects for statistically exploring the characteristics of S-bearing species. A comparative analysis of Sbearing species' parameters across various environments can unveil both similarities and differences among these species, advancing our knowledge of sulfur chemistry.  \nThe paper is structured as follows: Section 2 introduces the observational information and the samples of dense cores. Section 3 presents detected species, spectra, and derived parameters. In Section 4, statistical analyses and sulfur depletion problem are discussed. Section 5 summarizes the key highlights of the paper. The appendix contains additional information related to our results.", '2. OBSERVATIONS': "2.1. 870 µ m data  \nThe data utilized in this study are from ALMA Cycle 5 observations (PI: Tie Liu, Project ID: 2017.1.00545S), conducted between 18th and 20th May 2018 (UTC). The antenna configuration of C43-1 was adopted, employing a total of 48 12-meter antennas for the observations.  \nTable 1. Information of observations  \nNote -a : Dust temperatures listed in the table come from (Urquhart et al. 2018). b : Distance of each clump is estimated by Xu et al. (2024).  \nIn the project, a total of 11 massive protoclusters were observed. The basic information of these 11 protoclusters is listed in Table 1. To cover the entire regions of these sources, a mosaic observation mode was used to extend the field of view. Phase, amplitude, and bandpass calibration were conducted, the calibrators used are listed in Table 1. The flux calibration error of our observations is ∼ 10%, which is consistent with the flux calibration accuracy in Band 7 as stated in ALMA documentation (Braatz 2020). For the continuum data of the 11 sources, three rounds of phase self-calibration and one round of amplitude self-calibration were performed using CASA (CASA Team et al. 2022). In the 3 rounds of phase-only self-calibration, the solution intervals start with a single scan length and terminate with an integration times of 6.05s. For the amplitude self-calibration, the solution interval ranges were set as infinite. The detailed information of self-calibration is presented in Appendix A.  \nDuring the imaging process, the deconvolution was set as ' hogbom ' with a weighting parameter of ' briggs ', and the robust parameter of 0.5 was used to balance the sensitivity and resolution of the final images. A primary beam correction was applied. The resulting field of view is approximately 46 arcseconds with an angular resolution of approximately 0.8 arcseconds, corresponding to spatial scales ranging from 0.01 to 0.04 pc. The largest angular scales (LAS) of the observations are from 7.2 '' to 9.2 '' .  \nFour spectral windows, namely spw25, spw27, spw29, and spw31, were set to observe spectral line data. The center frequencies of spw25 and spw27 are 354.4 GHz and 356.7 GHz, respectively, with a spectral resolution of 0.244 MHz. The center frequencies of spw29 and spw31 are 345.1 GHz and 343.2 GHz, respectively, with a spectral resolution of 0.976 MHz. Each spectral window has a total of 1920 channels, resulting in a bandwidth  \nof 468.75 MHz for spw25 and spw27, and 1875 MHz for spw29 and spw31. The observational data used in our survey covers a total frequency bandwidth of 4.68 GHz. The spectra have been smoothed to a uniform spectral resolution of 0.976 MHz across all spectral windows. To improve the signal-to-noise ratio of the spectral images, the solutions obtained from continuum self-calibration were applied to the spectral line data.", '2.2. The sample of dense cores': 'By using the continuum data of the 11 protoclusters, Xu et al. (2024) identified 248 continuum dense cores. The physical parameters (such as core mass, peak flux, H 2 column density, number density, surface density, and core radius) were estimated and provided in their paper. The idenfitited cores are classified into 6 types in Xu et al. (2024). In this study, we used the dense core sample reported by them, but we classified all the samples into two categories: prestellar core candidates and protostellar cores based on specific criteria. Cores without any signs of star formation were classified as prestellar core candidates. Cores already associated with starforming activities or showing signs of warming up were all classified as protostellar cores. The spatial distributions of the dense cores are shown in the left-most panels of Figure 1 as blue ellipses, with labeled core names.', '3.1. Detections of S-bearing species': "We extracted spectra from the position of the brightest pixel of continuum emissions within the core's region. Line surveys were conducted to identify S-bearing species in dense cores based on the following criteria. Firstly, the line emission should have at least three continuous spectral channels with intensity exceeding three times the root-mean-square (rms) noise of spectrum. Secondly, the velocity of the observed transitions must  \nFigure 1. Illustration of dense core locations, continuum emissions, and S-bearing molecular gas distributions. The leftmost sequence of panels, numbered (1) to (11), showcases the continuum emissions, depicted with a red background and outlined in black contours. Contour levels are set at 6 σ , 12 σ , 24 σ , 48 σ , and continue until reaching the peak intensities. The rms noises σ are listed in Table 1, alongside the locations of dense cores, indicated by numbered blue ellipses. The small panels, such as (1a) to (1h), display the distributions of S-bearing species, represented by green contours. Contour levels range from 5% to 95% of peak values in 20% increments. The dust continuum emissions are superimposed as red backgrounds.. Gas distribution is generated by emissions with intensities 3 times higher than the rms noise. To ensure a clear representation of the gas distributions, the small panels are presented as magnified views. The red rectangles in the leftmost panels correspond to the areas that have been enlarged for detailed examination.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 1. -Continued  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 1. -Continued  \n<!-- image -->  \nalign with the systemic velocity of protocluster. By checking the spectra of 248 dense cores, a total of 10 S-bearing species (CS, H 2 CS, SO, 33 SO, SO 2 , 33 SO 2 , 34 SO 2 , SO 18 O, OC 34 S, and NS) were detected in our observations. The basic informations of detectable transitions is listed in Table B1, which are taken from CDMS (Muller et al. 2005). It should be noted that not all the transitions listed in the table can be detected in each core due to line blending effect and the different physical environments of the core.  \nFigure 2 presents the spectra of the I14498-7 core as an illustrative example. For 33 SO 2 , SO 18 O, and SO 2 , multiple detectable transitions fall within the observed frequency range, only one transition of each species that is less affected by line blending effects is shown in Figure 2. In the case of 34 SO 2 , all detectable transitions are covered by spw29, and most of them are less affected by the line blending effect. Therefore, we present all transitions of 34 SO 2 in Figure 2. A total of 8 detectable transitions of H 2 CS exist, with only 4 unblended transitions shown in the figure.  \nThe spectra of CS generally present obvious optically thick absorption dips and outflow wings, and the spectra of NS are significantly blended with CO outflow wings. We detected the CS and NS emissions, but the useful parameters are hardly derived from their spectra. Therefore, we presented the observed spectra of CS and NS without fitting results, and excluded them from further analyses. The S-bearing molecular spectra of all dense cores can be found in Appendix C as figure sets.  \nFigure 3 illustrates the detection rates of different molecules. The number of detected species for each core is listed in the last column of Table 2. The majority of dense cores (236 cores in total) detected at least one Sbearing species, these cores are labeled with number '1' to'10' in the last column of Table 2. A total of 12 dense cores do not have any detection of S-bearing species, as abeled with '0' in the last column of Table 2. There are 10 dense cores detected 10 S-bearing species, as labeled with '10' in the last column of Table 2. Artur de la Villarmois et al. (2023) reported the detection rates of 90% of CS, 78% of SO, and 44% of SO 2 in a sample of 50 Class 0/I protostellar sources in the Perseus molecular cloud. Our observations show a comparable detection rate for CS and SO in protostellar cores with theirs. The SO 2 detection rate is slightly lower because the SO 2 emissions in some cores are blended with HCO + (4-3), preventing us from claiming the detection of SO 2 in these cores.  \nThe integrated intensity maps (moment 0) for the detected species were generated and presented as green contours within the small panels of Figure 1. In the process of generating these maps, we have excluded emis-  \nsion that falls below three times the root mean square (rms) noise level. From Figure 1, gas distributions of H 2 CS and SO exhibit a notable similarity and are more spatially extended compared to those of other species. The distributions of the species, including SO 2 , OC 34 S, 33 SO, 34 SO 2 , 33 SO 2 , and SO 18 O, are found to be analogous to one another with compact morphologies. The distribution of SO 2 are similar to the isotopologues but with a more pronounced extension, which attributed to the partly blending with the HCO + line. The extended SO 2 features in the gas distribution are identified as contamination from the outflow wings of HCO + .", '3.2.1. Parameters derived from emissions': "We employed the SLIM tool from the MADCUBA package 1 (Mart'ın et al. 2019) to interactively estimate the physical parameters and their associated uncertainties under local thermodynamic equilibrium (LTE) assumption. The MADCUBA has garnered widespread application in the field of star formation (e.g. Mininni et al. 2023; Zeng et al. 2021; Jim'enez-Serra et al. 2021; Coletta et al. 2020; Rivilla et al. 2020). The parameters were derived by fitting an LTE model on the observed spectral emissions. The detailed fitting processes are presented in Appendix B. The LTE assumption is pivotal for parameter estimation. To ascertain the validity of the LTE assumption for our observations, we calculated the critical densities for SO, H 2 CS, and SO 2 using collisional data from the LAMDA database (van der Tak et al. 2020) at a uniform kinetic temperature of 30 K 2 . The critical densities were determined to be 3.7 × 10 3 , 8.3 × 10 4 , and 7.5 × 10 5 cm -3 , respectively. From Table 2, only one dense core, I15520-core32, exhibits a molecular hydrogen number density ( n H 2 ) below 10 5 cm -3 . Consequently, for SO, H 2 CS, and SO 2 , the LTE assumption is deemed reliable. We were unable to estimate the critical densities for 33 SO, 33 SO 2 , 34 SO 2 , OC 34 S, and SO 18 O due to the lack of collisional data. Goldsmith (2001) suggested that dust and gas are well coupled at molecular hydrogen densities exceeding 10 6 cm -3 . In our sample, almost all the dense cores 3 with S-bearing species detections have n H 2 higher than 10 6 cm -3 . Given these conditions, the LTE assumption in our study is considered reliable for all species in the majority of dense cores.  \nFigure 2. Spectra of the I14498-7 core are depicted. The observation is represented by the black stepped line in each panel. The spectrum exhibiting a non-Gaussian profile has been fitted with multiple velocity components, and the different components are shown as colored curves in each panel. The red stepped line represents the aggregated modeling of all velocity components. Molecular names and component numbers are indicated in the upper left corner of each panel. Multiple transitions of H 2 CS and 34 SO 2 are also plotted. The quantum numbers (QNs) of transitions for H 2 CS and 34 SO 2 are presented in the panels as well. The spectra of CS and NS were also presented with only the black stepped line, it is clear that the CS spectrum is double-peaked and with obvious outflow wings, and the NS spectrum is blended by CO outflow wings. Therefore, the spectra of CS and NS are excluded from the analyses.  \n<!-- image -->  \nH 2 CS and 34 SO 2 exhibit multiple unblended transitions that span a wide range of upper energy levels (see Table B1). We priorly conducted LTE fitting on the spectra of H 2 CS and 34 SO 2 to ascertain the gas temperatures. These specific transitions have also been utilized for temperature estimation in other studies (e.g. Liu et al. 2011; Chen et al. 2024). Out of the 115 cores that exhibited detectable H 2 CS transitions, 99 of them possessed a minimum of four transitions suitable for the derivation of gas temperatures. For the 34 SO 2 transitions, detections were made in 27 cores, and all were found to be adequate for calculating excitation temperatures.  \nIt should be noted that 4 transitions of H 2 CS less affected by line blending effect were used to estimate the rotational temperature. We marked these transitions with ' ∗ ' in Table B1. 34 SO 2 exhibits 12 detectable transitions in the observed data, and these transitions are not significantly affected by the line blending effect.  \nThis characteristic renders 34 SO 2 a reliable tracer of gas temperature for S-bearing species. However, the low detection rate of 13% suggests that it can only be utilized to trace the temperatures of a limited number of dense cores.  \nDue to the absence of multiple unblended transitions for the species SO, 33 SO 2 , OC 34 S, SO 2 , SO 18 O, and 33 SO 2 , their respective temperatures could not be determined through fitting processes. As elaborated in Section 3.2, the gas distributions of H 2 CS and SO are analogous, and a similar pattern is observed in their abundances, as detailed in Section 4.1. Therefore, we assume that H 2 CS and SO are likely to emanate from analogous environments, and may have similar gas temperatures. With this assumption, we utilized the gas temperature derived from H 2 CS as a fixed parameter in the LTE fitting for the SO spectra. Consequently, only the column density, systemic velocity, and line width  \nFigure 3. The detection rate of S-bearing species. The upper, middle, and bottom panels present detection rates of all cores, prestellar core candidates, and protostellar cores, respectively.  \n<!-- image -->  \nwere treated as free parameters in the fitting process for SO.  \nSimilarly, the gas distributions of SO 2 , 33 SO, OC 34 S, 33 SO 2 , and SO 18 O are analogous to those of 34 SO 2 . Based on this similarity, we hypothesize that these species also originate from comparable environments and exhibit similar temperatures. Accordingly, the temperature obtained from 34 SO 2 was adopted as a fixed parameter in the LTE fitting for their respective spectra, with only the column density, systemic velocity, and line width serving as free parameters for the fitting.  \nIt is important to note that the assumption of similar origins for the S-bearing species introduces additional uncertainties into the derived physical parameters. However, in the absence of superior alternatives, this approach is currently the most viable for estimating the parameters of SO, 33 SO, OC 34 S, SO 2 , SO 18 O, and 33 SO 2 . The dust temperatures of 11 massive protoclusters were given by Fa'undez et al. (2004) and Urquhart et al. (2018), but their angular resolutions are about 20 to 30 times larger than ours. They cannot be directly  \napplied to estimate the parameters of the species under investigation without causing significant uncertainty.  \nDuring the fitting processes, a Gaussian profile function is generally assumed. It is evident that the spectra of some detected species, SO in Figure 2 for instance, exhibit non-Gaussian profiles. This phenomenon could be attributed to three main reasons. Firstly, the presence of inhomogeneous/bulk motion of the gas envelope can lead to asymmetric non-Gaussian profiles. Secondly, fragmentation may occur on scales below our minimum spatial resolution, resulting in substructures with different systemic velocities. Thirdly, the interferometer is filtering-out extended structures larger than the LAS of the observations. Therefore, the observed spectra represent a stacking of Gaussian profiles from several velocity components. To properly obtain parameters from spectra, we fitted spectra with multiple velocity components in case that they exhibited obvious non-Gaussian profiles. Initially, a fitting with a single Gaussian component is carried out on a spectrum. If the fitting residual shows at least three consecutive spectral channels with intensities exceeding three times of rms noise, an additional Gaussian component is incorporated into the model. This iterative process continues until the fitting residual is subdued to below three times the rms noise threshold. Results of each individual component's fit are graphically represented by distinct colored curves in Figure 2. The aggregated fitting results of all components are shown with a red stepped line, giving a comprehensive overview of the spectral fitting.  \nThe parameters extracted from the SLIM fittings are compiled in the tables in Appendix B. It is noteworthy that the SLIM fitting accounts for optical depth in its analysis. During the fitting process, only the H 2 CS 10(0,10)-9(0,9) transition and the SO emissions in some massive cores exhibited optical depths exceeding the value of 1. Most line emissions are optically thin, which implies that the influence of optical depth on our results and subsequent analyses is negligible.", '3.2.2. Column density estimations for SO absorptions': 'The SO spectra of certain cores display discernible absorption features, as illustrated in Figure Set C1. Under the LTE assumption, the conventional approach to calculate the column density from SO absorption lines involves an integration over velocity with optical depth. This process is followed by the application of the subsequent equation (Mangum & Shirley 2015):  \nN = 3 h 8 π 3 Sµ 2 Q rot g u exp ( E u kT ex ) exp ( hν kT ex ) -1 × ∫ τ ν dv, (1)  \nwhere k and h are Boltzmann and Planck constants, respectively. The permanent dipole moment µ , line strength S are taken from CDMS (Muller et al. 2005), g u is level degeneracy. The partition function Q rot of SO can be estimated by an approximation for linear polyatomic molecules (McDowell 1988) as:  \nQ rot ≃ kT hB 0 exp ( hB 0 3 kT ) , (2)  \nwhere B 0 is rotational constant.  \nOptical depth for SO absorption can be estimated by the following approximation:  \nτ ν = -ln [ 1 -T R f ( J ( T cont ) -J ( T ex )) ] , (3)  \nwhere T R is radiation temperature of SO absorption, J ( T cont ) is Rayleigh-Jeans equivalent temperature of background continuum emissions, and J ( T ex ) is Rayleigh-Jeans equivalent temperature of T ex . The gas temperatures derived from H 2 CS were utilized to calculate the optical depth. In instances where the excitation temperature is not available, the dust temperature listed in Table 1 will be employed as a surrogate for these calculations. The column densities of SO absorptions are listed in Table B4.', '3.3. Abundances': 'Column densities of S-bearing species are estimated from the spectra extracted from pixels coinciding with  \nthe peak positions of continuum emissions within these dense cores. Consequently, the peak H 2 column density ( N H 2 ) from Xu et al. (2024) is adopted for calculating the relative abundances of S-bearing species to H 2 . The peak N H 2 and the column densities of S-bearing species were calculated from the same emitting areas. The resulting abundances of the S-bearing species for each core are detailed in Table 2. It is important to note that the abundances listed in Table 2 are computed from the integrated column density across all velocity components. As such, each species within a dense core is associated with a singular abundance value, irrespective of the number of velocity components detected.  \nThe estimation of N H 2 may be influenced by the optical depth of the continuum emission, given that the initial estimate was based on the assumption of an optically thin continuum emission (Xu et al. 2024). To ascertain whether the continuum emissions are optically thin or thick, we have calculated the optical depth of the continuum emission employing the following equation (Frau et al. 2010; Gieser et al. 2021):  \nτ cont = ln [1 -S ν Ω B ( T d ) ] (4)  \nIn the equation, S ν represents the flux density at the continuum peak position of the core, Ω denotes the solid angle of the beam, and B ( T d ) is the Planck function evaluated at the dust temperature T d . Utilizing the dust temperatures listed in Table 1, we estimated the optical depth from 10 -5 to 10 -2 . These values suggest that the estimations of N H 2 are in line with the optically thin assumption for the continuum emission.  \nTable 2 . Physical properties and molecular abundances of the dense cores  \nTable 2 (continued)Table 2 continued  \nTang et al.', 'Table 2 (continued)': 'Table 2 continued  \nTable 2 (continued)Table 2 continued  \nTang et al. Table 2 (continued)', '4. DISCUSSION': "The sample of 248 dense cores covers a wide range of densities ( N H 2 :2.7 × 10 22 ∼ 2.5 × 10 24 cm -2 , n H 2 :1.0 × 10 4 ∼ 3.7 × 10 8 cm -3 ), mass (0.1 ∼ 52.6 M ⊙ ), temperature ( T H 2 CS : 39.1 ∼ 134.7 K; T 34 SO 2 : 72.3 ∼ 158 K), molecular abundance (total abundance: 4.4 × 10 -12 ∼ 6.9 × 10 -7 ), and evolutionary stages (prestellar core candidates and protostellar cores). Statistical analyses of  \nthese data provide an opportunity to investigate the properties of S-bearing species in dense cores, which cannot be obtained through individual case studies alone. In this section, we will primarily conduct a statistical analysis of the parameters we have obtained to glean more information regarding the characteristics of Sbearing species, as well as discuss the similarities and differences in gas distributions, abundance, and gas temperatures among various S-bearing molecules, culminating in a discussion on the issue of sulfur depletion.  \nFigure 4. Two-dimensional KDE plots of molecular abundances, number of detected species, and temperatures versus low- and high-mass core categories. Prestellar core candidates and protostellar cores are analyzed separately and shown in blue and red colors, respectively. The vertical axis in each panel shows the mass categories of low- and high-mass cores, while the horizontal axis in each panel displays a specific parameter used for KDE analysis. Some species such as 33 SO, OC 34 S, SO 2 , 34 SO 2 , 33 SO 2 , and SO 18 O were not detected in prestellar core candidates, labeled as 'No detection' in the panels.  \n<!-- image -->", '4.1. Similarities and differences between H 2 CS and SO': "Kernel density estimation (KDE) is a non-parametric technique that employs kernel smoothing to infer the probability density function of a dataset. This method stands out as a valuable tool for visualizing the distribution of data. Unlike histograms, which may be subject to the arbitrariness of bin width selection, KDE provides a more fluid representation of the underlying data patterns. We categorized core samples into distinct groups: low-mass, high-mass prestellar core candidates and protostellar cores. For each category, we performed separate KDE analyses to elucidate the distribution characteristics. Figure 4 presents the two-dimensional KDE plots, which illustrates parameters of molecular abundances, total abundance, number of detected species, and excitation temperatures for H 2 CS and 34 SO 2 , respectively.  \nAmong the S-bearing species, a notable distinction is found: SO and H 2 CS are present in both prestellar core candidates and protostellar cores, whereas other species are exclusively detected in protostellar environments. As depicted in panels (1a) to (2b) of Figure 4, the abundances of SO and H 2 CS exhibit a common trend, with prestellar core candidates displaying lower abundances on the order of approximately 10 -9 , and protostellar cores exhibiting higher abundances around 10 -8 . This suggests a potential enhancement in both SO and H 2 CS abundances as dense cores evolve. The high detection rates of SO and H 2 CS indicate that these two species span various evolutionary stages from prestellar to protostellar cores (Chen et al. 2024). The gas distributions of SO and H 2 CS in Figure 1 demonstrate similar morphologies, underscoring a significant resemblance between SO and H 2 CS.  \nWhile SO and H 2 CS exhibit several similarities, they also display notable differences. One key distinction lies in their detection rates, as depicted in Figure 3. The higher detection rate of SO may suggest that it forms more easily in prestellar cores compared to H 2 CS. This enhanced formation could be due to the efficient production of SO through barrierless neutral-neutral reactions in the gas phase (Fuente et al. 2016; Semenov et al. 2018; Laas & Caselli 2019; Taquet et al. 2020; Booth et al. 2023). Furthermore, SO has been recognized as a tracer of diverse physical conditions, including protoplanetary disks, planet formation regions, shock zones, outflows, hot cores (e.g. Sakai et al. 2014, 2017; Tabone et al. 2017; Lee et al. 2018; Taquet et al. 2020; van Gelder et al. 2021; Tychoniec et al. 2021).  \nThe KDE plot of SO abundances reveals an extension to lower values, reaching approximately 10 -11 , while H 2 CS abundances are at higher levels, above 10 -10 . Our  \nresults support the notion that SO may be formed in more diffused regions compared to H 2 CS. The abundance range of H 2 CS reported in other similar studies, from 1.2 × 10 -10 to 5 × 10 -8 (e.g. Minh et al. 2011; Shimajiri et al. 2015; Minh et al. 2018; Luo et al. 2019; Rivi'ereMarichalar et al. 2019; Fuente et al. 2021; Moller et al. 2021; Bouscasse et al. 2022; Chen et al. 2024; Fontani et al. 2023), is consistent with the our results. There are also regions detected H 2 CS with abundance as low as ∼ 10 -11 in PDRs (Rivi'ere-Marichalar et al. 2019), where is quiet different with our sources.  \nIn summary, SO demonstrates a propensity for detection across a wider range of physical and chemical environments although numerous similarities in the properties of SO and H 2 CS are found.", '4.2. H 2 CS and 34 SO 2 trace environments with different temperatures': 'Panel (11a) to (12b) of Figure 4 present the twodimensional KDE plots of T (H 2 CS) and T ( 34 SO 2 ) in prestellar core candidates and protostellar cores. The analysis reveals that the kernel density peaks for T (H 2 CS) fall below 70 K across both low- and high-mass cores. In contrast, the peaks for T ( 34 SO 2 ) are distinctly situated at temperatures exceeding 100 K for cores of both mass categories.  \nThe KDE plot of H 2 CS shows an extension towards temperatures as high as 130 K in high-mass cores, which can be attributed to the presence of some hot cores in our samples (Liu et al. 2023; Chen et al. 2024; Xu et al. 2024), and they generally have gas temperatures exceeding 100 K. Comparing to the KDE plot of T (H 2 CS), the kernel density for T ( 34 SO 2 ) extends to even more elevated temperatures ( > 150 K) for both low- and highmass cores. This indicates that 34 SO 2 and H 2 CS trace the regions with different temperatures in both low- and high-mass cores. The Spearman correlation coefficient between T (H 2 CS) and T ( 34 SO 2 ) was calculated as 0.57, indicating a merely moderate positive correlation between them, also suggesting that H 2 CS and 34 SO 2 seem to trace gases with different temperatures.  \nFuente et al. (2021) investigated the hot core in Mon R2, noted a spatial offset between the distribution of H 2 CS gas and the continuum peak, suggesting that H 2 CS may not be an accurate tracer of the hot core. The gas distributions of H 2 CS and 34 SO 2 illustrated in Figure 1 demonstrate that the H 2 CS is more extended than 34 SO 2 . Additionally, the presence of H 2 CS is not invariably associated with hot cores/corinos 4 , but the  \nmajority of 34 SO 2 detections ( ∼ 89%) were indeed found in association with hot cores/corinos. Consequently, our results indicate that H 2 CS appears to trace the outer warm envelope of the core (Fuente et al. 2021; Tychoniec et al. 2021; Bouscasse et al. 2022), while 34 SO 2 is more representative of the high-temperature central regions.', '4.3. The differences between high-mass and low-mass cores': 'Panels (3a) to (8b) in Figure 4 present that 33 SO, OC 34 S, SO 2 , 33 SO, OC 34 S, 34 SO 2 , 33 SO 2 , and SO 18 O were exclusively found in protostellar cores, with their abundances notably differing between high-mass and low-mass cores. Specifically, 33 SO, 34 SO 2 , 33 SO 2 , and SO 18 O display a higher prevalence and abundances in high-mass cores, implying an enhancement of these species in high-mass star-forming regions. Additionally, within protostellar cores, both the diversity of detected species and the total molecular abundance are more pronounced in high-mass cores. Most strikingly, the dense cores that exhibit detections of all S-bearing species are predominantly high-mass, accounting for 80%. Our findings suggest that the process of high-mass star formation is associated with an enhancement in both the diversity of S-bearing species and the abundance.  \nIn the case of SO 2 , the kernel density peaks for both low-mass and high-mass cores are approximately 10 -7 . SO 2 has been observed in low-mass cores with abundances as low as 10 -10 , while in high-mass cores, the majority exhibit SO 2 abundances exceeding 10 -8 , with a kernel density peak centered around 10 -7 .  \nIt is widely recognized that SO and SO 2 are potential tracers of shocked gas, with their abundances potentially enhanced by the effects of shocks (Sakai et al. 2014, 2017; Artur de la Villarmois et al. 2019; Oya et al. 2019; van Gelder et al. 2021; Tychoniec et al. 2021). To ascertain whether the abundances of the observed S-bearing species were influenced by shocks or other non-thermal processes, we calculated the non-thermal velocity dispersion ( σ NT ) for identified species as following equation:  \nσ NT = [ σ 2 obs -kT ex m ] 1 2 , (5) √  \nwhere σ obs = v FWHM / 8 ln 2 is the observed total velocity dispersion, and v FWHM is line width from tables in Appendix B. k is Boltzmann constant, and m is the molecular mass. The relationship between σ NT and their abundances are investigated and shown in Figure 5. The figure presents each panel with the Spearman correlation coefficients, where high-mass cores and lowmass cores are denoted by orange and blue, respectively.  \nFor SO, the Spearman correlation ( r ) values of 0.37 and 0.51 in high-mass and low-mass cores suggest that  \nFigure 5. The correlation between the non-thermal velocity dispersions ( σ NT ) of SO (upper panel) and SO 2 (bottom panel) versus their abundances. In each panel, high-mass cores and low-mass cores are shown as orange and blue colors, respectively. The Spearman correlation coefficient ( r ) is displayed in the upper-right corner of each panel.  \n<!-- image -->  \nSO abundance is not significantly enhanced by nonthermal motions. Concerning SO 2 , with r values of 0.06 and 0.79 in high-mass and low-mass cores, indicating non-thermal motions appear to well correlated with SO 2 abundance in low-mass cores. The robust correlation of low-mass cores suggests that the abundances of SO 2 can be enhanced by non-thermal motions, and a significant subset of low-mass cores remains in relatively quiescent environments. In high-mass cores, no discernible correlation exists between SO 2 abundances and σ NT . However, these cores are clustered in the upperright corner of the panel, indicating a general trend of higher σ NT levels coincident with increased SO 2 abundance across all high-mass cores.  \nOur results suggest that high-mass cores likely exhibit more vigorous star-forming activities and potentially higher-velocity shocks, resulting in enhanced SO 2 abundances across all high-mass cores.  \nFigure 6. Violin plots for total abundance (upper panel) and the number of detected species (bottom panel). The regions with orange and blue colors represent the presence and absence of outflows, respectively. In each panel, the left plots are statistics of low-mass cores, and the right plots are statistics of high-mass cores. The dotted lines in the plots represent the quartiles.  \n<!-- image -->', '4.4. S-bearing molecular abundances enhanced by outflows': "It is well-established that outflows can enahnce Sbearing species in the gas phase (e.g. Pineau des Forets et al. 1993; Charnley 1997; Hatchell et al. 1998; Viti et al. 2001; Taquet et al. 2020; Artur de la Villarmois et al. 2023; Zhang et al. 2023). Based on the information provided by Xu et al. (2024), we categorized the protostellar cores into two groups: the presence and absence of outflows. Comparisons of abundances of each species were conducted. We only found that the total abundances and the number of detected species in high-mass cores with outflows exceeded those cores without outflows. But this trend was less pronounced in low-mass cores. A comparison of the total abundance and the number of detected S-bearing species between these two categories is presented in Figure 6. Our result suggests that outflows have indeed elevated the abundances and detection rates of S-bearing species in high-mass cores, underscoring the significant impact of high-mass starforming feedback on the formation of S-bearing species. This resuls align with our earlier discussion in Section 4.3.  \nFigure 7. Relationship between the abundance-increasing factor F of SO 2 and temperatures. The linear fitting is demonstrated as the blue dashed line. The fitting results are shown in legend. The definition of restoration factor can be found in Section 4.5. The F of other studies are estimated by the abundances and gas temperatures collected from literatures (Wang et al. 2013; Holdship et al. 2019; Rivi'ereMarichalar et al. 2019; Fuente et al. 2021; Moller et al. 2021; Bouscasse et al. 2022; Fontani et al. 2023).  \n<!-- image -->", '4.5. Sulfur depletion': "As stated in Section 1, the issue of sulfur depletion has been observed over the past few decades, but the reservoir of the sulfur element is still unknown. Our observed total S-bearing molecular abundance ranges from 4.4 × 10 -12 to 6.9 × 10 -7 . All of these values are lower than the sulfur abundance in the diffuse and highly ionized regions of ∼ 10 -5 (e.g. Savage & Sembach 1996; Mart'ın-Hern'andez et al. 2002; Howk et al. 2006; Garc'ıaRojas et al. 2006; Yamamoto 2017), indicating sulfur depletion is occurring in all dense cores. Possible solid sulfur reservoir candidates, for instance, H 2 S, OCS, SO 2 , organo-sulfur species, sulfur chains (S n , n =2 · · · 8), and sulfide minerals (FeS and MgS) were proposed (Keller et al. 2002; Hony et al. 2002; Scappini et al. 2003; Wakelam et al. 2004; Jim'enez-Escobar & Mu˜noz Caro 2011; Druard & Wakelam 2012; Vidal et al. 2017; Laas & Caselli 2019; Kama et al. 2019). Unfortunately, only OCS (Geballe et al. 1985; Palumbo et al. 1995; McClure et al. 2023) and SO 2 (Boogert et al. 1997; Yang et al. 2022; McClure et al. 2023; Rocha et al. 2024) have been detected in interstellar ices.  \nThe solid sulfur reservoir may sublimate from grains to gas phase as the temperature increases, making it a temperature-sensitive species. The temperature tracers H 2 CS and 34 SO 2 enable us to determine the gas  \ntemperatures of S-bearing species in dense cores. A positive correlation between molecular abundance and gas temperature maybe a indicator of a sulfur reservoir. However, as temperatures elevate, chemical reactions also intensify, potentially leading to increased gas-phase abundances of numerous species. Therefore, assessing the correlation between abundance and gas temperature alone may be inadequate to investigate sulfur depletion. This is because many species, not just sulfur-bearing ones, could exhibit positive correlations with temperature, complicating the interpretation of the depletion process.  \nSolid sulfur reservoirs, such as OCS and SO 2 , are anticipated to exhibit heightened temperature sensitivity compared to other S-bearing species due to they are locked onto ice grains. With the ascent of temperature, the abundances of these species are predicted to escalate at a more accelerated pace relative to other Sbearing species. This is attributed to their dual production pathways, which encompass not only gas-phase reactions but also the desorption mechanism from dust grains into the gas phase. In light of this, we have introduced the abundance increasing factor ( F ) to quantitatively evaluate the rate at which the abundance of a specific sulfur-bearing species increases in comparison to other S-bearing species.  \nF = f x ( f total /n ) , (6)  \nwhere f x represents the abundance of a certain S-bearing species, f total denotes the total S-bearing molecular abundance, and n is the number of detected S-bearing species in a dense core.  \nConsequently, we have calculated the abundance increasing factor ( F ) for SO 2 exclusively within dense cores. A value of F exceeding unity for SO 2 suggests that it is more abundant than other S-bearing species in a given dense core. An elevated F value denotes a relatively higher abundance of SO 2 in comparison to other S-bearing species. Thereby, the correlation between F and the T (H 2 CS) 5 serves as a metric for gauging the rate at which SO 2 abundance increases. The interrelation between F and gas temperature for SO 2 is illustrated in Figure 7. A linear correlation has been found  \nwith a slope of 0.04 and an intercept of -0.36. For comparison, we have collated abundances and gas temperatures from analogous studies to estimate F , and these literature-derived results are concurrently presented in the figure.  \nAlinear regression analysis was conducted on the samples within our study, which disclosed a positive correlation between the abundance increasing factor ( F ) and temperature. This result is corroborated by the results of previous studies, indicating a consistency in the observed correlation. The presence of this correlation suggests that SO 2 is likely to be significantly depleted onto grain mantels and may undergo desorption as the temperature elevates. This desorption process could accelerate the increase in SO 2 abundance in the gas phase relative to other sulfur-bearing species. Furthermore, this correlation potentially explains the relatively high abundance of SO 2 ( ∼ 10 -7 ) observed in star-forming regions. SO 2 can be produced through both gas-phase reactions and direct desorption from grain mantles. Additionally, as shown in Figure 1, the SO 2 emission is more compact and probably coming from the desorption of ice mantles located at higher temperatures, closer to the protostar.  \nHowever, our observations indicate that the abundance of SO 2 constitutes only approximately 1% of the total sulfur inventory. When accounting for other volatile species, such as H 2 S and OCS, their combined contribution remains minimal, at less than 4% (Geballe et al. 1985; Palumbo et al. 1995; Boogert et al. 1997; Kama et al. 2019). In areas influenced by shocks, the sulfur abundance may be elevated to approximately 10% of the total sulfur inventory (Anderson et al. 2013), yet approximately 90% of the sulfur remains unaccounted for. Our results, combined with analogous studies, imply that while SO 2 could serve as a sulfur reservoir on the grain mantles of dense cores, but the predominant solid sulfur reservoir is more likely to be composed of refractory elements. For instance, the sulfide mineral FeS, as proposed by Hony et al. (2002); Kama et al. (2019), is a strong candidate due to its high sublimation temperature of around 655 K (Larimer 1967; Lodders 2003; Kama et al. 2019).", '5. SUMMARY': "Based on the ALMA Band 7 (870 µ m) data, we have investigated the properties of S-bearing molecules in 248 dense cores from 11 massive star-forming clumps. The main results of this study are summarized as follows:  \n- 1. Ten S-bearing species of CS, SO, H 2 CS, NS, SO 2 , 33 SO, OC 34 S, 34 SO 2 , 33 SO 2 , and SO 18 O are detected. The parameters of them (except CS and NS) were estimated and presented.  \n- 2. Our analysis shows that H 2 CS and SO exhibit similar gas distribution. The 2D KDE analysis suggests they have comparable abundance properties, with levels of ∼ 10 -9 in prestellar core candidates and ∼ 10 -8 in protostellar cores. SO appears to be more easily detected in a broader range of physical and chemical environments than H 2 CS despite many similarities found between SO and H 2 CS.\n- 3. Temperatures were determined from multiple transitions of H 2 CS and 34 SO 2 . Gas distributions and 2D KDE analyses indicate that H 2 CS and 34 SO 2 trace regions with varying temperatures. H 2 CS appears to trace the outer warm envelope, while 34 SO 2 is associated with high-temperature central-regions.\n- 4. Abundances of S-bearing species and detection rates are higher in high-mass cores. SO and SO 2 show correlations between abundances and σ NT in low-mass cores, suggesting non-thermal motions may enhance these species. Particularly for SO 2 , its abundance and non-thermal velocity dispersion ( σ NT ) appear notably increased in high-mass cores, suggesting that feedback from high-mass star formation significantly boosts the abundance of SO 2 .\n- 5. The outflows elevate the total abundance and detection rates of S-bearing species in high-mass cores, underscoring the significant impact of highmass star-forming feedback on the formation of S-bearing species.\n- 6. Sulfur depletion was observed to be prevalent in dense cores. A positive correlation between abundance increasing factor ( F ) and temperature for SO 2 was discovered. The results from similar previous works also conform to this correlation, suggesting that a large amount of SO 2 may exist in the solid phase and desorb rapidly into a gaseous phase as the temperature rises. This result can also explain a high SO 2 abundance of ∼ 10 -7 and compact gas distributions in massive star-forming regions. However, SO 2 abundance only accounts for ∼ 1% of the total sulfur inventory, implying that the predominant solid sulfur reservoirs would be more likely to be composed of refractory species with higher sublimation temperatures.  \nlow-mass, high-mass, prestellar, and protostellar cores. The analysis, based on a large sample size, offers valuable insights of sulfur-depletion and properties of Sbearing species. However, the observed number of Sbearing species remains limited. Future investigations on S-bearing species could be enhanced by expanding observations to include as many as unblending transitions in larger sample sizes, leading to more accurate parameters estimations.  \nThe authors deeply appreciate the crucial comments from the anonymous reviewers, which have undoubtedly improved the quality of this paper. This work has been supported by the National Science Foundation of China (12203011, 12033005). MYT acknowledge the support by Yunnan province talent support program, Yunnan provincial Department of Science and Technology through grant No.202101BA070001261, Yunnan University Laboratory Open Project, and PhD research startup foundation of Chuxiong Normal University. S.-L.Q. thanks the Xinjiang Uygur Autonomous Region of China for their support through the Tianchi Program. T.L by the international partnership program of Chinese Academy of Sciences through grant No.114231KYSB20200009, Shanghai Pujiang Program 20PJ1415500, the science research grants from the China Manned Space Project with no. CMSCSST-2021-B06. L.A.Z. acknowledges financial support from CONACyT-280775, UNAM-PAPIIT IN110618, and IN112323 grants, M'exico. Y.P. Peng acknowledges support from NSFC through grant No. 12303028. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 This paper uses the following ALMA data: ADS/JAO.ALMA#2017.1.00545.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. 22 23 24 25 26 27 28 29  \nFacilities: ALMA(Atacama Large Milimeter/Submilimeter Array)  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018, 2022), MADCUBA(Mart'ın et al. 2019), CASA(CASA Team et al. 2022)  \nIn conclusion, our study presents statistical findings on S-bearing species across 249 dense cores, including  \nTable A1. Improvements of Signal-to-noise ratios (SNRs) of selfcalibrations  \nNote  \n-  \n- a : The SNR of without performing self-calibration.\n- b : The SNR after performing phase-only self-calibration.\n- c : The SNR after performing phase-only and amplitude self-calibration.\n- d : The improvement of SNR by comparing images before and after self-calibration.", 'A. SELF-CALIBRATION': 'Self-calibration could be a crucial technique for enhancing image quality when the signal-to-noise ratio (SNR) exceeds 20 (Brogan et al. 2018). As shown in Table A1, our data quality justifies the use of selfcalibration. We conducted three rounds of phase-only self-calibration followed by one round of amplitude selfcalibration.  \nThe phase-only self-calibration began with a solution interval of a scan length ( solint=inf ), which is ideal for the initial phase-only self-calibration (Richards et al. 2022). The second round was set with a solution interval of 6.05s ( solint=int ). If this setting results in over 10% failed solutions, intervals of 12.1s or 18.15s will be used. The third round of phase-only self-calibration utilized 6.05s to properly correct the phase over time. Table A1 presents the self-calibration results for each source, showing that the phase-only self-calibrations improved the SNRs of the images.  \nPhase-only self-calibrations can alter image fluxes, prompting us to conduct a subsequent round of amplitude self-calibration. The amplitude self-calibration has a greater potential to affect images than phase-only self-calibration. To avoid systemic errors, we set the solution interval to solint=inf since amplitude does not vary rapidly over time. We also enabled solnorm=True to maintain the ensemble amplitude correction around 1, which helps prevent an overall reduction in integrated flux (Brogan et al. 2018). Finally, we applied the amplitude solutions alongside the best phase solutions. The results in Table A1 indicate that amplitude self-calibrations did not induce systemic errors in the  \nimages. The SNRs improved by 5% to 47% compared to images without self-calibration.', 'B. PARAMETERS DERIVED FROM SPECTRA': "In this section, we present the spectroscopic parameters of detected species in Table B1. The parameters of systemic velocities, line widths, temperatures, and column densities of detected S-bearing species in each core are presented in Table B2 to Table B9. The dense cores without detection of S-bearing species are excluded from the tables.  \nThe derivation of these parameters was executed by SLIM from the MADCUBA package (Mart'ın et al. 2019). SLIM is capable of generating synthetic spectra under the LTE assumption, utilizing five input parameters: excitation temperature, column density, line width, velocity, and source size. The integrated 'auto-fit' feature within SLIM facilitates the best non-linear least-squares fit, employing the Levenberg-Marquardt algorithm to fit the LTE model on the observed spectra. The output of the fitting process encompasses parameters such as column density, temperature, velocity, line width, and source size.  \nThe spectral data for each core were extracted from pixels that coincide with the peak positions of the continuum flux. Continuum intensity has been subtracted from these spectra. The resultant reduced spectra were then imported into SLIM as observational data for the fitting procedure. During this process, the source size for each core was set as a fixed parameter, with the core sizes referenced from Xu et al. (2024). For molecules H 2 CS and 34 SO 2 , the column density, temperature, velocity, and line width were treated as free parameters during the fitting. For 33 SO, OC 34 S, SO 2 , 33 SO 2 , and SO 18 O, the fitting was performed with column density, velocity, and line width as the free parameters.  \nTable B1. Detected transitions and corresponding spectroscopic parameters  \n(  \na  \n)  \n:  \n- × : Transitions of CS and NS molecules are detected. The CS line is typically optically thick, displaying distinct absorption profiles and outflow wings. NS lines are heavily impacted by line blending;  \n∗  \n:  \nThe transitions adopted to estimate gas temperatures;  \nBlended: The transition is blended with other molecular transitons, and excluded from the presentation in Figure 2.  \nNote -The basic information of transitions are from CDMS(Muller et al. 2005).  \nTable B2 . H 2 CS parameters.  \nTable B2 continued  \nNote  \n-  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.", 'Tang et al.': 'Table B2 (continued)', 'Sulfur-bearing molecules in dense cores': 'Table B3 . 34 SO 2 parameters.  \nNote -  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.  \nTable B4 . SO parameters.Table B4 continued  \nTang et al.', 'Table B4 (continued)': 'Table B4 continued  \nTable B4 (continued)  \nTable B4 continued  \nNote  \n-  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.\n- absorption: The component is observed as absorption. Their parameters of V lsr and FHWM are not estimated. The column densities are calculated by equation (1).', 'Table B5 . 33 SO parameters.': 'Table B4 (continued)  \nNote  \n-  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.', 'Table B6 . OC 34 S parameters.': 'Table B7 continued  \nNote -  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.  \nTable B7 . SO 2 parameters.  \nTang et al.  \nTable B7 (continued)  \nNote -  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.  \nTable B8 . SO 18 O parameters.  \nNote -  \n- a : The number of velocity components defined in the spectra.  \nNote  \n-  \n- a : The number of velocity components defined in the spectra.\n- b : The column density are calculated by intergrating all different components.  \nFigure C1. SO spectra and fitting results.The black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (159 images) can be found in the online journal.  \n<!-- image -->', 'C. SPECTRA OF DETECTED SPECIES': 'In this section, the spectra of all detected S-bearing species are presented. In each plot, the observed spectra were illustrated as a black solid line, and the LTE fitting of velocity components was denoted as colored curves. The red lines represent the summed fitting results of all velocity components. It is worth noting that the column densities of SO absorptions are estimated by equation (1). Therefore, the LTE fitting results presented in this section are only from the emissions of detected species. The absorptions have been excluded from the fittings.  \nFigure C2. 33 SO spectra and fitting results. In the figure, the black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (36 images) can be found in the online journal.  \n<!-- image -->  \nTable B9 . 33 SO 2 parameters.  \nFigure C3. OC 34 S spectra and fitting results. In the figure, the black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (39 images) can be found in the online journal.  \n<!-- image -->  \nFigure C4. SO 18 Ospectra and fitting results. In the figure, the black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (13 images) can be found in the online journal.  \n<!-- image -->  \nFigure C5. 33 SO 2 spectra and fitting results. The black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (15 images) can be found in the online journal.  \n<!-- image -->  \nFigure C6. H 2 CS spectra and fitting results.The black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (115 images) can be found in the online journal.  \n<!-- image -->  \nFigure C7. 34 SO 2 spectra and fitting results. The black stepped line represents observed data, the colored curves depict fitting results of decomposed velocity components, and the red stepped line shows the summed spectral of all fitting components. The complete figure set (27 images) can be found in the online journal.  \n<!-- image -->', 'REFERENCES': "Xu, F., Wang, K., Liu, T., et al. 2024, ApJS, 270, 9, doi: 10.3847/1538-4365/acfee5  \nYamamoto, S. 2017, Introduction to Astrochemistry: Chemical Evolution from Interstellar Clouds to Star and Planet Formation, doi: 10.1007/978-4-431-54171-4  \nYang, Y.-L., Green, J. D., Pontoppidan, K. M., et al. 2022, ApJL, 941, L13, doi: 10.3847/2041-8213/aca289  \n- Yue, Y.-H., Qin, S.-L., Liu, T., et al. 2021, Research in Astronomy and Astrophysics, 21, 014,  \ndoi: 10.1088/1674-4527/21/1/14  \n- Zapata, L. A., Garay, G., Palau, A., et al. 2019, ApJ, 872, 176, doi: 10.3847/1538-4357/aafedf  \nZeng, S., Jim'enez-Serra, I., Rivilla, V. M., et al. 2021,  \nApJL, 920, L27, doi: 10.3847/2041-8213/ac2c7e  \n- Zhang, Z. E., Yang, Y.-l., Zhang, Y., et al. 2023, ApJ, 946, 113, doi: 10.3847/1538-4357/acbdf7"}