An ALCHEMI inspection of sulphur-bearing species towards the central molecular zone of NGC 253

M. Bouvier; S. Viti; E. Behrens; J. Butterworth; K.-Y. Huang; J. G. Mangum; N. Harada; S. Martín; V. M. Rivilla; S. Muller; K. Sakamoto; Y. Yoshimura; K. Tanaka; K. Nakanishi; R. Herrero-Illana; L. Colzi; M. D. Gorski; C. Henkel; P. K. Humire; D. S. Meier; P. P. van der Werf; Y. T. Yan

doi:10.1051/0004-6361/202449186

Home

All issues

Volume 689 (September 2024)

A&A, 689 (2024) A64

Full HTML

Open Access

Issue		A&A Volume 689, September 2024


Article Number		A64
Number of page(s)		30
Section		Extragalactic astronomy
DOI		https://doi.org/10.1051/0004-6361/202449186
Published online		02 September 2024

A&A, 689, A64 (2024)

An ALCHEMI inspection of sulphur-bearing species towards the central molecular zone of NGC 253

M. Bouvier¹^,⋆, S. Viti¹, E. Behrens², J. Butterworth¹, K.-Y. Huang¹, J. G. Mangum³, N. Harada⁴^,5, S. Martín⁶^,7, V. M. Rivilla⁸, S. Muller⁹, K. Sakamoto⁵, Y. Yoshimura¹⁰, K. Tanaka¹¹, K. Nakanishi⁴^,5, R. Herrero-Illana⁶, L. Colzi⁸, M. D. Gorski¹², C. Henkel¹³^,14^,15, P. K. Humire¹⁶, D. S. Meier¹⁷^,18, P. P. van der Werf¹ and Y. T. Yan¹³

¹ Leiden Observatory, Leiden University, PO Box 9513 23000 RA Leiden, The Netherlands
² Department of Astronomy, University of Virginia, PO Box 400325 530 McCormick Road, Charlottesville, VA 22904-4325, USA
³ National Radio Astronomy Observatory, 520 Edgemond Road, Charlottesville, VA 22903-2475, USA
⁴ National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
⁵ Department of Astronomy, School of Science, The Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-1855, Japan
⁶ European Southern Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile
⁷ Joint ALMA Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile
⁸ Centro de Astrobiología (CAB), INTA-CSIC, Carretera de Ajalvir km 4, Torrejón de Ardoz, 28850 Madrid, Spain
⁹ Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 43992 Onsala, Sweden
¹⁰ Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan
¹¹ Department of Physics, Faculty of Science and Technology, Keio University, 3-14-1 Hiyoshi, Yokohama, Kanagawa 223–8522, Japan
¹² Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
¹³ Max-Planck-Institut für Radioastronomie, Auf-dem-Hügel 69, 53121 Bonn, Germany
¹⁴ Astron. Dept., Faculty of Science, King Abdulaziz University, PO Box 80203 Jeddah 21589, Saudi Arabia
¹⁵ Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 830011 Urumqi, China
¹⁶ Departamento de Astronomia, Instituto de Astronomia, Geofísica e Ciências Atmosféricas da USP, Cidade Universitária, 05508-900 São Paulo, SP, Brazil
¹⁷ New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA
¹⁸ National Radio Astronomy Observatory, PO Box O 1003 Lopezville Road, Socorro, NM 87801, USA

Received: 8 January 2024
Accepted: 11 May 2024

Abstract

Context. Sulphur-bearing species are detected in various environments within Galactic star-forming regions and are particularly abundant in the gas phase of outflow and shocked regions in addition to photo-dissociation regions. Thanks to the powerful capabilities of millimetre interferometers, studying sulphur-bearing species and their region of emission in various extreme extra-galactic environments (e.g. starburst and active galactic nuclei) and at a high-angular resolution and sensitivity is now possible.

Aims. In this work, we aim to investigate the nature of the emission from the most common sulphur-bearing species observable at millimetre wavelengths towards the nuclear starburst of the nearby galaxy NGC 253. We intend to understand which type of regions are probed by sulphur-bearing species and which process(es) dominate(s) the release of sulphur into the gas phase.

Methods. We used the high-angular resolution (1.6″ or ∼27 pc) observations from the ALCHEMI ALMA Large Program to image several sulphur-bearing species towards the central molecular zone (CMZ) of NGC 253. We performed local thermodynamic equilibrium (LTE) and non-LTE large velocity gradient (LVG) analyses to derive the physical conditions of the gas where the sulphur-bearing species are emitted, and their abundance ratios across the CMZ. Finally, we compared our results with previous ALCHEMI studies and a few selected Galactic environments.

Results. To reproduce the observations, we modelled two gas components for most of the sulphur-bearing species investigated in this work. We found that not all sulphur-bearing species trace the same type of gas: strong evidence indicates that H₂S and part of the emission of OCS, H₂CS, and SO are tracing shocks, whilst part of SO and CS emission rather traces the dense molecular gas. For some species, such as CCS and SO₂, we could not firmly conclude on their origin of emission.

Conclusions. The present analysis indicates that the emission from most sulphur-bearing species throughout the CMZ is likely dominated by shocks associated with ongoing star formation. In the inner part of the CMZ where the presence of super star clusters was previously indicated, we could not distinguish between shocks or thermal evaporation as the main process releasing the S-bearing species.

Key words: astrochemistry / ISM: molecules / galaxies: abundances / galaxies: active / galaxies: ISM / galaxies: starburst

^⋆

Corresponding author; bouvier@strw.leidenuniv.nl.

© The Authors 2024

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Sulphur-bearing species are widely detected in various environments within Galactic star-forming regions such as hot cores (HCs) and hot corinos (e.g. Blake et al. 1996; Charnley 1997; van der Tak et al. 2003; Wakelam et al. 2004a; Esplugues et al. 2013; Crockett et al. 2014a; Li et al. 2015; Drozdovskaya et al. 2018; Luo et al. 2019; Codella et al. 2021; Bouscasse et al. 2022; Artur de la Villarmois et al. 2023), outflows and shocks (e.g. Bachiller et al. 2001; Codella et al. 2005; Lefloch et al. 2005; Podio et al. 2014; Holdship et al. 2016, 2019; Kwon et al. 2015; Ospina-Zamudio et al. 2019; Taquet et al. 2020; Feng et al. 2020; Tychoniec et al. 2021; Artur de la Villarmois et al. 2023), photo-dissociation regions (PDRs; e.g. Jansen et al. 1995; Goicoechea et al. 2006; Ginard et al. 2012; Goicoechea et al. 2021; Rivière-Marichalar et al. 2019), and molecular clouds and cores (e.g. Wakelam et al. 2004b; Agúndez & Wakelam 2013; Vastel et al. 2018; Navarro-Almaida et al. 2020; Hily-Blant et al. 2022; Spezzano et al. 2022; Fuente et al. 2023; Fontani et al. 2023).

In the dense gas phase of star-forming regions, most of the sulphur (S) is locked onto the icy grain mantles and its most efficient release into the gas phase can occur primarily via thermal evaporation in warm regions or non-thermal evaporation through shocks (e.g. Minh et al. 1990; Charnley 1997; Bachiller et al. 2001; Hatchell & Viti 2002; Esplugues et al. 2014; Woods et al. 2015). Once in the gas phase, S-bearing species are known to be extremely reactive, meaning their abundance greatly depends on the thermal and kinetic properties of the gas (Viti et al. 2004). S-bearing species are particularly useful for reconstructing the chemical history and dynamics of a variety of objects as they have been used as chemical clocks for the shocked gas and HCs (e.g. Pineau des Forêts et al. 1993; Charnley 1997; Codella & Bachiller 1999; Viti et al. 2001; Wakelam et al. 2004b; Li et al. 2015; Feng et al. 2020; Taquet et al. 2020; Codella et al. 2021; Bouscasse et al. 2022; Esplugues et al. 2023).

Thanks to the powerful capabilities of millimetre interferometers such as the Atacama Large Millimeter/(sub-)millimeter Array (ALMA) and the NOrthern Extended Millimetre Array (NOEMA), it is now possible to study many sulphur-bearing species in the extreme environments of external galaxies – for example starburst or active galactic nuclei (AGN) – compared to our own Galaxy. A variety of sulphur-bearing species have been previously detected in various types of galaxies. For example, SO and SO₂ likely trace HCs in the Large Magellanic Cloud (Sewiło et al. 2022), whilst H₂S, SO, SOs₂, and CS seem to probe shocks in AGNs (e.g. Scourfield et al. 2020; Kameno et al. 2023). In starburst galaxies, sulphur-bearing species were found to be associated with PDRs (e.g. Meier & Turner 2005; Bayet et al. 2009), dense molecular gas (e.g. Martín et al. 2005; Aladro et al. 2011; Ginard et al. 2015), shocks (e.g. Martín et al. 2003, 2005; Viti et al. 2014; Meier et al. 2015; Sato et al. 2022), or irradiated gas by forming stars (i.e. hot-core-like regions; Minh et al. 2007; Sato et al. 2022). However, these studies were limited by either a low angular resolution (Martín et al. 2005) or by a low number of available transitions (Minh et al. 2007) to allow the physical conditions of the gas where S-bearing species are emitted to be constrained, and hence reveal which environment they trace in starburst galaxies.

The ALMA Comprehensive High-resolution Extragalactic Molecular Inventory (ALCHEMI; Martín et al. 2021) large programme is the first and broader unbiased imaging molecular survey performed towards the starburst central molecular zone (CMZ) of the nearby galaxy NGC 253 beyond the previous single-dish unbiased surveys by Martín et al. (2006) and Aladro et al. (2015), giving us access to both high-angular resolution (1.6″ or ∼27 pc) and a large number of transitions for many molecular species, including sulphur-bearing species. Several ALCHEMI studies have already been published, which derived a high cosmic-ray ionisation rate (CRIR) towards the NGC 253 CMZ using various molecular species and abundance ratios (C₂H: Holdship et al. 2021; H₃O⁺/SO: Holdship et al. 2022; HCO⁺/HOC⁺: Harada et al. 2021, and HCN/HNC: Behrens et al. 2022), reported the first extragalactic detection of a phosphorus-bearing species, PN (Haasler et al. 2022), characterised various types of shocks throughout the CMZ (using HNCO and SiO: Huang et al. 2023; methanol masers: Humire et al. 2022; or HOCO⁺: Harada et al. 2022), explored the general physical properties of the CMZ of NGC 253 (Tanaka et al. 2024), and performed a full principal component analysis (PCA; Harada et al. 2024).

In this paper we aim to use the molecular richness and high-angular resolution of ALCHEMI to perform a complete investigation of the most abundant sulphur-bearing species in Galactic star-forming regions (CS, H₂S, OCS, SO, H₂CS, SO₂, and CCS; e.g. Hatchell et al. 1998; Li et al. 2015; Widicus Weaver et al. 2017; Humire et al. 2020) towards the starburst CMZ of NGC 253, where the detection of these species has been previously reported by Martín et al. (2021). We aim to determine the main process(es) dominating the release of sulphur-bearing species into the gas phase towards the CMZ of NGC 253, a prototypical starburst environment.

The paper is structured as follows: An overview of the NGC 253 CMZ structure is presented in Sect. 2 and the observations are in Sect. 3. Sect. 4 shows the emission distribution of the targeted S-bearing species and the regions chosen to perform the analyses. The results of the Gaussian fit and the LTE and non-LTE analyses are described in Sect. 5 and then are further discussed in Sect. 6. Finally, the conclusions are summarised in Sect. 7.

2. The structure of the central molecular zone of NGC 253

NGC 253 is a nearby starburst galaxy located at a distance of 3.5 Mpc (Rekola et al. 2005) with an inclination of 76° (McCormick et al. 2013). The inner kpc region exhibits intense star formation with a rate of ∼2 M_⊙ yr⁻¹, which represents about half of the global star formation activity of the galaxy (Bendo et al. 2015; Leroy et al. 2015). Within this starburst nucleus, the CMZ is a ∼300 pc × 100 pc region containing ten giant molecular clouds (GMCs) of a size of ∼30 pc identified in both molecular and continuum emission (Sakamoto et al. 2011; Leroy et al. 2015). The CMZ is particularly chemically rich (e.g. Martín et al. 2006, 2021; Aladro et al. 2015; Pérez-Beaupuits et al. 2018; Mangum et al. 2019; Krieger et al. 2020a) and relatively turbulent. Indeed, shocks have been identified in several GMCs lying close to the orbital intersections of the central bar hosted by NGC 253 (Harada et al. 2022; Huang et al. 2023; Humire et al. 2022), and have been claimed to explain its global chemical abundances (Martín et al. 2009a). Additionally, the presence of the starburst-driven large-scale outflow (Bolatto et al. 2013) and superbubbles (due to expanding HII regions of supernovae; Sakamoto et al. 2006; Gorski et al. 2017; Krieger et al. 2020b; Harada et al. 2021) could interact with the CMZ inducing streamers (Walter et al. 2017; Tanaka et al. 2024; Bao et al. 2024) and cloud-cloud collisions (Huang et al. 2023, Tanaka et al. 2024).

At a higher angular resolution, the GMCs located in the inner 160 pc of the CMZ can be resolved into several massive star-forming clumps called super proto- super star clusters (proto-SSCs; Ando et al. 2017; Leroy et al. 2018; Rico-Villas et al. 2020; Mills et al. 2021; Levy et al. 2021). These proto-SSCs are compact (≤1 pc), massive (M ≥ 10⁵ M_⊙) and young (≤3 Myr). Some of the proto-SSCs host super hot cores (SHCs; Rico-Villas et al. 2020), which are scaled up versions of Galactic HCs with high dust temperatures (∼200 − 300 K) and high densities (∼10⁶ cm⁻³), and show clear outflow signatures (Gorski et al. 2019; Levy et al. 2021). High kinetic temperatures (≥300 K) were also measured from H₂CO observations at scales ≤1″ (≤17 pc), very likely associated with dense gas heating sources (Mangum et al. 2019). Some studies suggested that the proto-SSCs formed following an inside-out formation (Rico-Villas et al. 2020; Mills et al. 2021) but this scenario seems to be contentious (Krieger et al. 2020a; Levy et al. 2022).

3. Observations, data processing, and line selection

The observations used in this work are part of the ALCHEMI Large Program (co-PIs.: S. Martín, N. Harada, and J. Mangum). We give here only a summary of the observations. The full details of the data acquisition, calibration and imaging are described by Martín et al. (2021). The observations were performed during Cycles 5 and 6, under the project codes 2017.1.00161.L and 2018.1.00162.S. The observations covered Bands 3 through 7, with a frequency ranging between 84.2 and 373.2 GHz and were centred at the position α(ICRS) = 00^h47^m33.26^s and δ(ICRS) = − 25° 17′17.7″; the common region imaged corresponds to a rectangular area of size 50″ × 20″ (850 × 340 pc) with a position angle of 65°, covering the CMZ of NGC 253. All the data cubes have a common beam size of 1.6″ × 1.6″ (∼27 pc), and a maximum recovered scale of ≥15″ (≥255 pc). The spectral resolution is ∼10 km s⁻¹. Finally, following Martín et al. (2021), we adopted a flux calibration error of 15%.

We used the continuum-subtracted FITS image cubes provided by ALCHEMI which contain the sulphur-bearing species and transitions of interest to this work. We processed the FITS image cubes using the packages MAPPING and CLASS from GILDAS¹ to produce the velocity-integrated maps, extract the spectra and perform a Gaussian fitting of the lines. We extracted the spectra from beam-sized regions (see Sect. 4) and converted them in brightness temperature using the equation:

$\begin{matrix} T (K) = 13.6 \times {(\frac{300 GHz}{ν})}^{2} (\frac{1^{″}}{θ_{\min}}) (\frac{1^{″}}{θ_{\max}}) S_{ν} (Jy) \end{matrix}$ $\begin{aligned} T(\mathrm{K})=13.6\times \left(\frac{300 \, \mathrm{GHz}}{\nu } \right)^2 \left(\frac{1^{\prime \prime }}{\theta _{\mathrm{min}}} \right)\left(\frac{1^{\prime \prime }}{\theta _{\mathrm{max}}} \right)S_\nu (\mathrm{Jy}) \end{aligned}$ (1)

where ν, θ_min, θ_max, and S_ν are the frequency, the minor and major axis of the beam, and the flux density, respectively. For each spectrum, we fitted a zeroth order polynomial to the line-free channels from the continuum-subtracted data cubes to retrieve the information about the spectral root-mean-square (rms).

We targeted the most abundant S-bearing species detected in Galactic star-forming regions, namely CS, H₂S, H₂CS, OCS, SO, SO₂ and CCS (e.g. Hatchell et al. 1998; Li et al. 2015; Widicus Weaver et al. 2017). All the species used in this work were previously detected in NGC 253 by Martín et al. (2021). We also extracted measurements of the detected isotopologues ³⁴SO and $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ . Isotopologues of CS are also detected but they will be the focus of a separate study so we do not consider them in the present paper. For the analysis, the threshold for line detection was set to 3σ for the integrated intensity. We did not consider lines that are contaminated by more than 15% by another transition, or lines that are too blended to allow a multiple Gaussian fit to be produced. To assess whether a line is contaminated, we used the existing line identification (Martín et al. 2021 and priv. comm.), and performed an additional check using the CLASS extension WEEDS (Maret et al. 2011). WEEDS uses the Jet Propulsion Laboratory (JPL; Pickett et al. 1998) and the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2001, 2005) databases. The list of all transitions used in this work, as well as their spectroscopic parameters, is presented in Table B.1.

4. Emission distribution and targeted regions

The velocity-integrated line intensity maps of the different sulphur-bearing species are shown in Figures 1, 2 and 3. The maps for the isotopologues $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ and ³⁴SO, are presented in Fig. A.1. For each species, we selected a sample of three transitions which are unblended and isolated or next to a line whose flux does not exceed the level of the flux calibration error of 15% (see Sect. 3) and integrated the emission over the full line profile ranging within [0;400] km s⁻¹. When possible, we selected lines covering different upper-level energies to show the variation of the emission distribution with different excitation levels reflecting the various excitation conditions. Before describing the emission distribution of each species, we give an explanation of how the regions where we extracted the spectra were chosen and how we named them.

Fig. 1.

Velocity-integrated maps of CS, H₂S, and SO, which show extended emission throughout the CMZ with their peak of intensity towards the inner CMZ. The positions of TH2 (Turner & Ho 1985) and the kinematic centre (Müller-Sánchez et al. 2010) are labelled by a filled magenta star and triangle, respectively. The regions where the spectra were extracted are shown by the magenta circles and are labelled in the left-most panel. The beam is depicted in the lower left corner of each plot. The scale bar of 100 pc corresponds to ∼6″. Top: velocity-integrated maps of CS (2–1), (4–3) and (6–5). Levels start at 3σ (1σ = 46, 150, and 148 mJy.beam⁻¹, respectively; white contours) with steps of 30, 30, and 50σ (black contours), respectively. Middle: velocity-integrated maps of H₂S (1_1, 0 − 1_0, 1), (2_2, 0 − 2_1, 1) and (3_3, 0 − 3_2, 1). Levels start at 3σ (1σ = 79, 50, and 81 mJy beam⁻¹, respectively; white contours) with steps of 15, 40, and 30σ (black contours), respectively. Bottom: velocity-integrated maps of SO (3₄ − 2₃) and (3₅ − 3₄). Levels start at 3σ (1σ = 50 and 28 mJy.beam⁻¹, respectively; white contours) with steps of 7 and 3σ (black contours), respectively.

Fig. 2.

Same as Fig. 1 but for OCS and H₂CS, which show their peak of emission towards each of the GMCs rather than the inner CMZ. Top: velocity-integrated maps of H₂CS (3_0, 3 − 2_0, 2), (7_1, 7 − 6_1, 6), and (11_1, 11 − 10_1, 10). Levels start at 3σ (1σ = 21, 50, and 27 mJy beam⁻¹, respectively; white contours) with steps of 4, 10, and 4σ (black contours), respectively. Bottom: velocity-integrated maps of OCS (8–7), (13–12) and (17–16). Levels start at 3σ (1σ = 14, 63.7 and 90 mJy beam⁻¹, respectively; white contours) with steps of 7, 4, and 7σ (black contours), respectively.

Fig. 3.

Same as Fig. 1 but for CCS and SO₂, which are mostly detected towards the inner part of the CMZ. Top: velocity-integrated maps of CCS (8₉ − 7₈), (11₁₂ − 10₁₁), and (18₁₉ − 17₁₈). Levels start at 3σ (1σ = 15.7, 30, and 33 mJy beam⁻¹, respectively; white contours) with steps of 7σ (black contours). Bottom: velocity-integrated maps of SO₂ (7_2, 6 − 6_1, 5), (12_1, 11 − 12_0, 12) and (19_4, 16 − 19_3, 17), from left to right. Levels start at 3σ (1σ = 18.4, 91.5 and 125 mJy beam⁻¹, respectively; white contours) with steps of 5, 3, and 3σ (black contours), respectively.

We selected the regions based on where the species show the most intense emission. Table 1 lists the selected regions and their associated coordinates. Each of the selected regions have been named after existing identifications. Regions GMC10, GMC7, and GMC6 correspond to the Giant Molecular Clouds 10, 7, and 6, respectively, as previously defined by Leroy et al. (2015). The S-bearing species have a peak intensity that is shifted from the GMC9, GMC8, GMC1 and GMC2 regions defined in Leroy et al. (2015), but which correspond to the newly defined regions GMC8a, GMC9a, GMC1a, GMC1b and GMC2b in Huang et al. (2023). We did not notice a particular peak of emission of the S-bearing species towards the position GMC2a and do not consider this position in this work, as we aim to focus on where the emission of the S-bearing species is the strongest. The positions SSC5 and SSC2 have their centre corresponding to the proto-SSCs 5 and 2, respectively, which were defined in Leroy et al. (2018). The emission of the S-bearing species peaks towards these proto-SSCs rather than GMC3 and GMC4. Hence, we used the coordinates of SSC5 and SSC2 to extract the spectra towards these two positions. It is important to note that as the ALCHEMI beam encompasses several proto-SSCs, the regions named SSC5 and SSC2 include also the SSC4 to SSC7 and the SSC1 to SSC3, respectively. Finally, we did not take into account the GMC5 region, which is close to the position of the strongest radio continuum called TH2 (Turner & Ho 1985) and the kinematic centre of NGC 253 (Müller-Sánchez et al. 2010). In previous studies, it was found that GMC5 shows particular features, such as absorption, due to the strong radio continuum emission at this location (Meier et al. 2015; Humire et al. 2022). The study of GMC5 is thus impossible to perform with the analysis proposed in this paper. Finally, the S-bearing species studied in this paper are not clearly detected in the streamers associated with the large-scale outflow, except for CS (Walter et al. 2017). We will therefore not investigate vertical structures or gas properties along these features.

Table 1.

Selected regions where the S-bearing species’ emission peaks, their coordinates, and their position within the CMZ.

Overall, we defined ten regions, shown by magenta circles in Figs. 1, 2 and 3. For the rest of the paper, we will refer to the part of the CMZ containing the regions GMC10, 9a, 8a, 1a, 1b, and 2b as the ‘outer CMZ’ (between 10″ and 16″ from the kinematic centre), and we will refer to the part of the CMZ encompassing the regions GMC7 to SSC2 (within ∼6″ around the kinematic centre) as the ‘inner CMZ’.

We find that the emission distribution varies depending on the species and the transitions. CS shows the most extended emission among the S-bearing species and is detected throughout the CMZ, with its emission peaking towards the inner CMZ. The emission distribution of H₂S and SO are relatively similar to that of CS (Fig. 1). For the isotopologues $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ and ³⁴SO, the former is detected only towards the inner CMZ whilst the latter is detected throughout the CMZ (Fig. A.1). However, as H₂S is less intense in the outer CMZ, the lack of detection of its isotopologue in these regions could be due to a lack of sensitivity. On the other hand, the emissions of OCS and H₂CS peak towards each of the GMCs/proto-SSCs throughout the CMZ rather than towards the GMCs/proto-SSCs of the inner CMZ (Fig. 2). Both species also peak towards SSC2 but not towards SSC5, as found in the case of the other species. Finally, CCS and SO₂ are the species which show most of their emission towards the inner CMZ, with SO₂, in particular, detected only towards the three innermost regions GMC6, SSC5, and SSC2 (Fig. 3). For all the species, lower-energy transitions (E_u ≤ 80 − 100 K for H₂S and OCS, E_u ≤ 50 − 60 K for H₂CS and SO, and E_u ≤ 25 K for CCS ) are extended throughout the CMZ whilst the higher-energy transitions are more concentrated towards the inner part of the CMZ.

5. Results

5.1. Gaussian fit and spectra

We extracted each spectrum from a circular region of a beam size (1.6″) centred at the coordinates given in Table 1. The spectra, the line-fitting results of the Gaussian fits and the rms for each species and region are publicly available on Zenodo². We centred each spectrum at the kinematic local standard of rest (LSRK) frequency of the associated transition before extracting the line parameters. To extract the integrated intensities (∫T_B dV), we used two methods depending on whether the lines are blended or not. In the case of non-blended lines, we measured the direct integration of the channel intensities. This is because most of the lines do not present a clear single-Gaussian profile. We nonetheless performed Gaussian fitting to extract the line widths (FWHM) and the peak velocity (V_lsr) to have a quantitative assessment towards these parameters. For the lines that are blended with another line, we performed a multi-Gaussian fit (with two Gaussian components) and we took the integrated intensity directly from the result of the fit, as well as the FWHM and V_lsr. If the blending with another line does not allow such a fit to be performed, the transition has been dismissed from the analysis. Finally, for the line profiles of CS towards GMC7 and GMC2b, two velocity components could be used for the fitting. However, for the rest of the analysis, we selected only the component with the highest peak intensity as this component is also detected in the other species (the other component is too weak or undetected, i.e. below the 3σ threshold). We thus discarded this component as we would not have a complete data set for a comparison with the other component and other regions. The results of the systemic velocities and line widths are discussed in Sects. 5.1.1 and 5.1.2 below.

5.1.1. Systemic velocities

Figure 4a shows the averaged systemic velocities derived for each species per region. The uncertainties derived from the Gaussian fit being usually smaller than the spectral resolution of 10 km s⁻¹, we thus adopted this value for the error bars. First, it appears that not all transitions of the same species peak at the same systemic velocity: the low upper-level energy transitions of CS (E_u < 20 K), H₂CS (E_u < 30 K), OCS (E_u < 100 K), ³⁴SO (E_u < 30 K) and CCS (E_u < 35 K) show a different peak compared to the higher upper-level transitions. However, depending on the species, this difference is not necessarily significant. Indeed, for CS the difference varies in the range ∼14 − 16 km s⁻¹ in most regions. The line widths vary between 65 and 81 km s⁻¹. The differences in peak velocity are thus less than 1/3 of the line widths and cannot be considered as significant. Similarly, the difference in peak velocities for OCS is at most 1/3 of the line widths. Additionally, and this is particularly visible with the CS lines, non-Gaussian line profiles could also lead to a shift in the derived systemic velocity. Therefore, for CS and OCS the difference in peak velocities between the low- and high-upper-energy level transitions is not significant. On the other hand, for H₂CS, ³⁴SO, and CCS, the difference in peak velocity between the low- and high upper-level energy transitions varies between 1/3 (CCS) and more than half (H₂CS and ³⁴SO) the line widths of these species and, hence, can be considered as significant.

Fig. 4.

Results of Gaussian fits. The listed sequence of clouds follows the layout of the GMCs of the CMZ as shown in Figs. 1–3. (a) Systemic velocities for each species as a function of the regions. If a strong difference is seen in the low- and high upper-level energy transitions of the species, the two velocity components are shown. (b) FWHM for each species as a function of the regions. If, for one species, a difference in line width is seen between the low- and high upper-level energy transitions, the two components are shown.

A second result is that, within each region, not all species peak at the same systemic velocity. Throughout the CMZ, a clear distinction is seen between the velocity peak of the low upper-energy levels of H₂CS and CCS, with that of the other species and transitions. On average, the difference in peak velocity is 10–20 km s⁻¹ compared to the low upper-level energy transitions of CS and OCS, and the high upper-level energy transitions of CCS. Here again, non-Gaussian line profiles could lead to a difference in the derived systemic velocities. Such small shifts in velocity might thus not be significant. A larger difference is seen compared to H₂S, SO, the low upper-level energy transitions of ³⁴SO, and the high upper-level energy transitions of CS and H₂CS with a shift of 30–40 km s⁻¹. Finally, the most important shift in peak velocity (40–50 km s⁻¹) is seen for SO₂ and the high upper-level energy transitions of ³⁴SO and OCS, when comparing with the systemic velocities of the other species and transitions. If we consider the line widths derived, all these peak differences are likely significant towards GMC10 where the peak difference is similar to or greater than the line width of the various transitions. Towards the regions SSC2 and GMC2b, where the line widths range between 20 and 60 km s⁻¹ (see Fig. 4b), peak velocity shifts higher than 30 km s⁻¹ can be considered as significant. Finally, in the other regions where most line widths are within the 40–80 km s⁻¹ range, shifts of at least 40 km s⁻¹ can be considered as significant. SSC5 is the only region where all species and transitions peak at almost the same velocity of ∼220 − 230 km s⁻¹, without any clear difference between transitions or species.

5.1.2. FWHMs

Figure 4b shows the FWHM analysis results. Towards GMC10, GMC2b, and SSC2 line widths are narrower compared to the other regions. The narrowest FWHMs are found towards GMC10, with all the species having FWHM < 40 km s⁻¹. For SSC2, the smaller line width (∼30 km s⁻¹) compared to most of the other regions is only seen for some species and transitions: ³⁴SO, SO₂ and OCS (E_u > 100 K). The rest of the species and transitions have line widths varying between 40 and 90 km s⁻¹. Then, within each region, there is some scatter in the FWHMs: in the outer CMZ, H₂S and CS have the largest FHWM, except towards GMC10 where only H₂S shows the larger FWHM. In the inner CMZ, CS has the largest line widths compared to the other species whilst SO₂ and ³⁴SO (E_u > 30 K) show the narrowest line widths. No general trend is found for the other species and transitions. Finally, for species where the low- (E_u < 20 − 30 K for CS, H₂CS, and ³⁴SO; E_u < 50 K for CCS; E_u < 100 K for OCS) and high upper-level energy transitions are peaking at different systemic velocities (see Sect. 5.1.1), there is no clear difference in the line widths, except for GMC1b where the line width of CS (E_u < 20 K) is larger than that of the other transitions by ∼22 km s⁻¹. Although some trends can be seen, it is important to keep in mind that optically thick lines and presence of multiple components (which cannot be disentangled at our current spectral resolution) can induce larger line widths or non-Gaussian line profiles. The optical depths are derived as part of the non-LTE LVG analysis (see Sect. 5.2.3) and the results of both line optical depths and whether the derived trends in line widths can be trusted will be discussed in Sect. 5.2.3.

5.2. Derivation of physical parameters

To derive the physical parameters of each species and for each region, we used two different methods. As a first step, we performed rotation diagrams (RDs) without beam correction. As we cover a large range of E_u for most of the species, it can allow us to identify whether different components with different excitation conditions (temperature and density) are present within each region studied. In this case, we expect to clearly see different slopes (hence different excitation temperatures) in the RDs. The main assumptions of the RD method are local thermodynamic equilibrium (LTE) conditions and optically thin emission of the lines. Although we do not expect these two assumptions to hold for all transitions and species, the results we obtain can be used as a first indication of the rotation temperature and beam-averaged column densities of each species. As examples of the results we derived from the LTE analysis, we will present below only those for the species for which we could not perform a non-LTE analysis as a second step.

We could not perform the non-LTE analysis for CCS, SO₂, and OCS. In the case of CCS, no collisional excitation rates are available in the online databases. For SO₂, the results of the RD predicted temperatures in the range 55 − 79 K for the first component. The LAMDA database (Schöier et al. 2005) provides two sets of collisional excitation rates; one for temperatures lower than 50 K and the other one for temperatures larger than 100 K. None of the two sets were thus usable to model the physical conditions of the first component of SO₂. For the second component, although we could use the set for which the rates were calculated for temperatures above 100 K, the non-LTE code we use cannot converge towards a solution if the first 200 levels of a molecule are populated, which is the case here. Finally, in the case of OCS, the collisional coefficient rates were initially calculated by Green & Chapman (1978) for the first thirteen levels and for a temperature up to 100 K. An extrapolation of these rates is available on the LAMDA database for transition levels up to J = 99 and for a temperature up to 500 K (Schöier et al. 2005). However, no matter which range of physical conditions (T, n) we run, the code does not converge and is unable to find a solution³. We could find a solution by using the rates available in the BASECOL database (Dubernet et al. 2013) but the available rates only concern the 13 first levels which corresponds to the first five transitions that are detected in this work. Consequently, we could not model the second component on the regions GMC6, SSC5, and SSC2. In this section, we will thus present the LTE results for OCS, SO₂, and CCS. The LTE results for the other species are publicly available on Zenodo⁴.

5.2.1. LTE analysis: OCS, CCS, and SO₂

Thanks to a large enough (≥10 for most species and regions) number of transitions detected, as well as the wide range of upper-level energy transitions available per species, we could perform rotation diagram (RD) analyses. The beam filling factor was set to 1 as (1) we do not know the size of emission of the species and (2) the emission is more extended than the beam size in most cases. We included the calibration error of 15% and the spectral RMS in the integrated flux (see Sect. 3). The best-fit was calculated by minimising the reduced chi-square, $χ_{r}^{2}$ $\chi^2_{\rm r}$ . Depending on the regions and species, we could fit more than one rotation temperature. This behaviour can be due to the presence of multiple components with different physical conditions, or to non-LTE and opacity effects (Goldsmith & Langer 1999). For CCS and OCS, the break in the temperature found at ∼30 K and ∼100 K for CCS and OCS, respectively, is consistent with the difference in systemic velocity seen in Sect. 5.1.1, which could favour the multiple component scenario. For SO₂, no clear shift in systemic velocity was found. The various possibilities which could cause the presence of more than one rotation temperature cannot be distinguished with rotation diagrams only but this will be investigated during the non-LTE analysis, except for CCS, SO₂.

Figures 5, 6, and 7 show the rotation diagrams OCS, CCS, and SO₂, respectively. The rotation temperatures (T_rot) and beam-averaged column densities (N_beam) derived are reported on Zenodo (see footnote).

Fig. 5.

Rotation diagrams of OCS. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.

Fig. 6.

Rotation diagrams of CCS. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.

Fig. 7.

Rotation diagrams of SO₂. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.

For OCS, two rotation temperature components were needed in the three innermost regions, GMC6, SSC5, and SSC2. The first temperature component includes transitions up to E_u = 111 K and ranges from T_rot = 10 K in GMC10 to T_rot = 65 K in SSC5. The range of the second temperature component is T_rot = 124 − 312 K. The beam-averaged column densities range between (0.2 and 11)×10¹⁴ cm⁻².

For CCS, several transitions with upper-level energy ranging between 20 and 230 K were detected. For the regions GMC6, SSC5, and SSC2, we found the fit was better if we considered three rotation temperature components ( $χ_{r}^{2}$ $\chi^2_{\rm r}$ = 1.5 vs. $χ_{r}^{2}$ $\chi^2_{\rm r}$ = 0.8). As mentioned above, the first break in rotation temperature corresponds to the shift seen in the systemic velocity for transitions at E_u ∼ 30 K (see Sect. 5.1.1). For GMC6 and SSC5, the second break in temperature corresponds to E_u ∼ 50 K. From the Gaussian fit, another shift in systemic velocity at this upper-level energy seem to be occur but the difference being of the order of the spectral resolution, we could not conclude whether this shift is significant, and we thus did not report and discuss it in Sect. 5.1.1. In the outer CMZ, only two transitions were detected. This resulted in relatively large errors, in particular for the beam-averaged column densities. Overall, the rotation temperatures range from ≤10 K in the outer CMZ to 72 K in the inner CMZ (towards GMC6) and the beam-averaged column densities are ≤3 × 10¹⁴ cm⁻².

Finally, SO₂ was detected towards only the three innermost regions of the CMZ. Two rotation temperatures components were needed to fit the observations. The first and second components have T_rot = (55 − 79) K and (115 − 298) K, respectively. The beam-averaged column densities are (1.4 − 6)×10¹⁴ cm⁻².

5.2.2. Methodology of the LVG analysis

As a second step, we performed a non-LTE analysis using the large velocity gradient (LVG) code grelvg⁵. The collisional excitation rates used for the non-LTE analysis, as well as the corresponding references and temperature ranges for which the coefficients are calculated, can be found in Table 2. The collision excitation rates were taken from the BASECOL database, except for H₂CS and SO₂ for which the rates were taken from the LAMDA database. For each region and species, we ran a large grid of models (≥10 000) varying the column density of the species, the gas density and the gas temperature. The range of values chosen are shown in Table 2 and are large enough to cover the physical conditions in molecular clouds, outflows/shocks and HCs, as we expect the emission from S-bearing species to be associated with either of these types of environments. For the gas temperature, we used the range of temperature for which the collisional rates were calculated for each species.

Table 2.

References for the collision rates and range of physical parameters.

As a result of the LTE analysis (see Sect. 5.2.1 and the online material on Zenodo), we identified between one and three temperature (T_rot) components, depending on the species and the region. For each species, we tried to fit all the transitions together to understand whether all the transitions could probe one same gas component (and thus the multiple rotation temperature we identified in Sect. 5.2.1 are due to the presence of non-LTE and opacity effects) but we could not find models reproducing simultaneously the low- and high-upper-level energy transitions. We, therefore, performed the non-LTE analysis for each of these components individually. We left as free parameters the column density, N_tot, the gas temperature, T_gas, and the gas density, n_gas. For components with at least 4 lines available, we left the source size, θ_s, as a free parameter. In grelvg the sizes can vary between 0.1″ and 150″. Given that the beam size of the observations (θ_b) is 1.6″, leaving the source size as a free parameter is equivalent of varying the beam filling factor defined as $θ_{s}^{2} / (θ_{s}^{2} + θ_{b}^{2})$ $\theta_{\rm s}^2/(\theta_{\rm s}^2+\theta_{\rm b}^2)$ from 0.004 to 1. We could leave the size as a free parameter in most cases, except when three or fewer transitions were available. In these cases, we fixed the source size as a first step before varying it around the fixed value to see the range of sizes for which the solution stays the same (same temperature and density range, and same χ²). The choice of the initial fixed size is indicated in Appendix C for the corresponding species and components. In some cases, the size could be left as a free parameter, but it did not result in constrained column densities and the best-fitted column density provided is in fact the product θ_s × N_tot. This can happen if multiple gas components are still mixed together and, thus, our two component assumption was not suitable. Another possibility is if all the lines are optically thin (τ < 1). In these cases, we ran again the best-fitting procedure, this time by fixing the source size to the best-fitted value obtained initially (i.e. the size associated with the best χ² value obtained when leaving the size as a free parameter). We then varied the source size around its best-fit value until the degeneracy disappears, hence when the θ_s × N_tot product does not give the same chi square. We had to use this method for H₂CS, and SO.

Finally, we assumed the geometry to be a semi-infinite slab (Scoville & Solomon 1974; de Jong et al. 1975). Choosing this geometry seemed relevant in case S-bearing species are tracing mostly shocks, as in our Galaxy. We nonetheless performed some modelling using the geometry of a uniform sphere which led to similar results (not shown). The model takes into account the cosmic background radiation (CBR) field in the level population excitation, with a temperature of 2.7 K. The ortho-to-para ratio (OPR) for H₂ is assumed to be set to the statistical value of 3. The assumed line widths are the average ones measured from the spectral lines for each species (and component in the case of multiple component identified) in each region. We included a calibration error of 15% in the observed intensities (see Sect. 3). The lines optical depths τ_L, are also an output of the LVG analysis. Further information concerning the LVG modelling of each species (e.g. ³²S/³⁴S ratio, OPR of H₂S and H₂CS) are presented in Appendix C.

5.2.3. Results of the LVG analysis

We present the results for individual species, followed by a comparison between the species across regions. Figure 8 shows an overview of the derived gas temperature, gas density, and column densities for each species analysed in this section, and for each region. Some regions were found to be relatively similar in terms of physical conditions derived for each species. In the outer CMZ, GMC2b shows different results from the other outer regions. In the inner CMZ, SSC2 shows fairly different results from the other inner regions. Therefore, Figure 9 and Table 3 summarise the results for GMC9a and GMC2b as representative regions for the outer CMZ, and GMC6 and SSC2, as representative regions for the inner CMZ. The results for the other regions are shown in Appendix C.

Fig. 8.

Results of the LVG analysis for CS (square markers), H₂CS (starred markers), H₂S (circle marker), OCS (crossed marker), and SO (diamond marker) as a function of the regions. The first and second components for CS, H₂CS, OCS, and SO are labelled “_1” and “_2”, respectively. Lower limits are indicated by upward arrows. The listed sequence of clouds follows the layout of the GMCs of the CMZ as shown in Figs. 1–3. From top to bottom, the results for the gas temperature, gas density, and column densities, are shown.

Table 3.

Best-fit results and 1σ confidence level range from the non-LTE LVG analysis for GMC9a, GMC6, SSC2, and GMC2b, respectively (from top to bottom).

CS. Based on the RD analysis, we modelled two different components for CS with the LVG code, the first one comprising the N = 2 − 1 to N = 4 − 3 levels, and the second one comprising the levels N = 5 − 4 to N = 7 − 6. The column densities range between 2 × 10¹³ and 1.7 × 10¹⁶ cm⁻². The lowest values are found towards GMC10, whilst the highest are found towards GMC6 and SSC5. For these regions, only a lower limit for the column density of the first component could be derived. This is because all the observed lines are optically thick, and the emission becomes that of a black body. In this case, the range of optical depths derived (τ ∼ 4 − 7) are likely lower limits to the true optical depths. The first component is very likely optically thick with high τ_L values derived for all the regions whilst the second component is likely not (τ_L ≤ 0.8 in all the regions). Therefore, the possible trend seen in FWHM in Sect. 5.1.2 with the low upper-level energy transitions of CS having larger line widths than the other species is likely due to the line’s optical depth and could thus not be related to the kinematics of the gas. The high upper-level energy transitions seem to have larger line widths compared to other species in some of the regions. The output of the LVG analysis indicates that these transitions are optically thin. However, if the emission size were to be smaller, the optical depths could be underestimated. Therefore, we cannot conclude firmly whether the higher line width for these transitions is due to the kinematics of the gas or to the line’s optical depth. The gas temperature derived for both components is relatively constant throughout the CMZ, with values between T_gas = 5 − 14 K, and between T_gas = 15 − 110 K for the first and second components, respectively. The highest values are derived towards GMC6 and SSC2 for the first component. For the second component, the lower values are derived towards GMC6 and SSC5, where the temperature does not exceed 30 K. The gas density traced by CS varies depending on the component and the region. Indeed, for the first component, gas densities traced are n_gas > 10⁵ − 10⁶ cm⁻³ in the outer CMZ, and n_gas < 10⁶ cm⁻³ in the inner CMZ. For the second component, the gas densities traced are n_gas > 10⁵ cm⁻³ for all the regions. We could not leave the size as a free parameter (see Appendix C), which means that other solutions might be possible. Here, we consider the very likely extension of CS to arise on scales of GMCs (∼30 pc). For the first component, the lower limit found for the size is ∼1 − 1.5″ for the entire CMZ whilst for the second component, a lower limit can be as low as ∼0.3 − 0.6″ (∼5 − 10 pc) indicating that this component could arise from a more compact region. Overall, due to the low number of transitions available for each of the component, the results for CS are among the less constrained and should be taken with caution. In particular, it is not really possible to distinguish between low temperature/high density and high temperature/low density solutions for the high upper-level energy transitions which is a classical problem occurring for linear molecules.

H₂S. We could derive the physical conditions of H₂S only in the inner CMZ and GMC10. The derived total column densities range between 8 × 10¹⁵ and 3 × 10¹⁷ cm⁻² in the inner CMZ, and between (0.6 and 3)×10¹⁵ cm⁻² for GMC10. The gas temperatures and densities derived are constant in all the regions with T_gas = 30 − 159 K and n_gas ≥ 10⁷ cm⁻³ (≥2 × 10⁶ for GMC7), respectively. Concerning the emission size, the LVG model predicts relatively compact emissions with values ranging between 0.17 and 0.5″ (∼3 − 8 pc). Finally, the lines are optically thin in GMC10 and GMC7 (τ_L = 0.01 − 0.4), whilst they are optically thick in GMC6, SSC5, and SSC2. In these regions, the $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ (1_1, 0 − 1_0, 1) transition is, however, optically thin. As we could not run the LVG code for the regions GMC8a, GMC1a, and GMC1b, we do not have the information on the line optical depths. Hence, we cannot conclude whether in these regions, the fact that H₂S seems to have a larger line width compared to the other species is a significant feature or not.

H₂CS. Based on the rotation diagrams, we performed the LVG analysis for two components of H₂CS in the inner CMZ. The two components are comprised of transitions with E_u < 50 and > 50 K, respectively. The derived column densities range between 5 × 10¹³ and 7 × 10¹⁷ cm⁻². The highest values are found towards GMC6 and SSC2 for both components. For the gas temperature, the derived temperature lies in the range 10–70 K for the first component and ≥35 K for the second component. The gas densities are ≥3 × 10³ cm⁻³ and between 7 × 10⁴ and 10⁷ cm⁻³, for the first and second components, respectively. We could leave the size as a free parameter most of the time and derived sizes of emission of ≤0.70″ (∼12 pc) for most regions. In the inner CMZ, where a second component of H₂CS has been identified, the sizes derived are more compact than for the first component, with an average best-fit value of 0.20″ (∼3 pc). Finally, for the first component, the lines are usually optically (moderately) thick for all the regions but GMC2b and GMC9a, for which the lines were found to be optically thin within the range of size of emission derived. For the second component the lines are optically thin towards SSC5 and GMC7 and likely optically thick towards GMC6 and SSC2.

OCS. As explained in Sect. 5.2.2, we could only use the first five transitions of OCS to derive the physical conditions of the gas emitting the first component of OCS. We assume that the results associated with these first five transitions are representative enough for the entire first component of OCS, which includes transitions up to the level N = 19 − 18 (E_u = 111 K). As we derive gas temperature around T_gas ∼ 14 K in the outer CMZ, and around 30 − 50 K in the inner CMZ, which correspond to the excitation temperatures found in the rotational diagrams, our assumption holds and we are likely not underestimating the temperature for this component of OCS. The total column densities derived range between 2 × 10¹⁵ and 5 × 10¹⁷ cm⁻². The highest values are found towards the regions of the outer CMZ. For the gas density, only a lower limit of ∼2 × 10⁴ cm⁻³ for all regions except for GMC9a (lower limit of 10⁵ cm⁻³) could be derived. The size was left as a free parameter, and the best-fit value of ∼0.17″ (∼3 pc) on average is the same for the entire CMZ. Across the range of sizes derived, the lines could be optically thick as the range of optical depths lies between τ_L = 0.1 − 3.0, except for SSC5 for which the lines are always optically thin (τ_L = 0.1 − 0.2) over the range of size derived.

SO. Both the RDs of SO and ³⁴SO show two components in the inner CMZ. We thus ran separately the LVG code for the two components, this first one including transitions with E_u < 40 K. For the first component, the derived total column densities in the range (0.1 − 5)×10¹⁶ cm⁻² are highest towards the outer CMZ (except GMC10) compared to the inner CMZ, where derived values lie in the range (0.3 − 5)×10¹⁵ cm⁻². The gas temperatures are also poorly constrained with lower limits of ∼30 K for most of the regions. SSC2 and GMC7 have the lower range of derived gas temperature with T_gas = 15 − 40 K and T_gas = 11 − 15 K, respectively. For the second component, the gas temperature of ≤60 K is found towards GMC6 and SSC2, whilst for GMC7 and SSC5, T_gas ≥ 50 K are found. In the outer CMZ, the gas densities all lie around n_gas = 10⁴ − 10⁵ cm⁻³ whilst for the inner CMZ, a lower limit for the two components of n_gas = 10⁵ cm⁻³ is found. For the emission size, two cases were identified between the outer and inner CMZ for the first component. In the outer CMZ, relatively compact sizes of emission were found, with θ_s ≤ 0.50″ (∼8.5 pc), whilst in the inner CMZ, the emission is more extended with θ_s ≥ 0.8″ (∼13.5 pc). The second component of the inner CMZ is more compact with best-fit values in the range θ_s = 0.20″ (∼3 pc) on average. Finally, for the first component, the emission is optically thin towards the inner CMZ and GMC 10 (τ_L ≤ 0.1) whilst potentially optically thick in the outer CMZ where τ_L ranges between 0.1 and 1.6. For the second component from the inner CMZ, the emission could be (moderately) optically thick with τ_L = 0.1 − 1.4 except for GMC7, for which the emission is optically thin (τ_L ≤ 0.2). The ³⁴SO transitions are always optically thin.

Different physical conditions were found for the gas emitting the various S-bearing species. In the outer CMZ, the first component of CS shows the coldest temperature with T_gas < 10 K, followed by OCS with T_gas < 20 K, and H₂CS with T_gas < 50 K, except in GMC2b where the gas temperature goes up to 70 K. Such cold temperatures for CS were already derived in a previous study performing an LVG analysis but with lower angular resolutions (Bayet et al. 2009). The second component of CS and the SO show the warmest gas temperature with T_gas ≥ 30 K. In terms of gas densities, within each region, the lower ones are usually derived for H₂CS with n_gas > 10³ − 10⁴ cm⁻³ except for GMC2b, for which the lower values are found for SO (n_gas = [1.5 − 2]×10⁴ cm⁻³). Then, OCS has larger density values derived with lower limits around 2 × 10⁴ cm⁻³, followed by CS with n_gas ≥ 10⁵ − 10⁶ cm⁻³).

In the inner CMZ, the situation is different. First, comparing the gas densities, the lowest values are found for OCS and the first component of CS and of H₂CS with n_gas ≥ 10³ − 10⁴ cm⁻³ on average. H₂S shows the largest gas density derived with n_gas ≥ 10⁷ cm⁻³ (n_gas ≥ 10⁶ cm⁻³ in GMC7). Within the inner CMZ, smaller densities are found towards GMC7 for H₂S and the first component of CS and H₂CS. A striking result concerns the gas temperatures derived with values below ∼50 K (≤100 K for the second component of CS) towards the region SSC2. In the remaining inner regions, the first components of CS and H₂CS have the lowest gas temperatures (T_gas ≤ 40 K) whilst H₂S and the second components of SO and H₂CS have the highest ones (T_gas > 45 K).

Comparing with the outer CMZ, OCS, H₂S, H₂CS, and the second component of CS show relatively constant gas density throughout the CMZ, whilst the first component of CS is emitted from a less dense gas in the inner CMZ compared to the outer CMZ but this result could be a result of fixing the emission size and, hence, possibly underestimating the density of CS in the regions. Concerning the gas temperature, the first components of CS and H₂CS, and SO in the case of GMC7, show the coldest gas temperatures within each region. Finally, when comparing the temperatures derived in the inner CMZ to those from the outer CMZ (see Tables 3 and C.1), the derived gas temperatures are higher in the inner CMZ for CS (first component) and OCS, the same for CS (second component) and H₂CS, and lower for SO in the case of SSC2 and GMC7.

6. Discussion

For the purpose of the following discussion and to ease the reading, we will refer to the ‘low-E_u’ component when we discuss the component corresponding to the low upper-level energy transitions of the species, and the ‘high-E_u’ component when we discuss the component associated with the high upper-level energy transitions of the species. For CCS where we identified 3 components, we will refer as the ‘low-’, ‘mid-’, and ’high-E_u’ component, respectively. Table 4 summarises the range of upper-level energies associated with each component of each species, as well as the derived physical conditions (T_gas, n_gas) across the CMZ.

Table 4.

Range of physical conditions derived towards the CMZ of NGC 253 for each species and component.

6.1. Comparison with previous ALCHEMI studies

A first step towards understanding what the S-bearing species trace is to compare the physical conditions of the gas from which they emit with the conditions derived from other species. Previous ALCHEMI studies derived physical conditions of the gas emitting C₂H (Holdship et al. 2021), H₃O⁺ and SO (Holdship et al. 2022), HCN and HNC (Behrens et al. 2022), and SiO and HNCO (Huang et al. 2023). From these studies, it was found that C₂H, SO, H₃O⁺ and the HCN/HNC abundance ratio were all sensitive to cosmic rays (CRs), C₂H, SO and H₃O⁺ being enhanced with the CRIR, and the HCN/HNC abundance ratio being anti-correlated with the CRIR. On the other hand, Huang et al. (2023) studied the low- and high-velocity shock tracers, HNCO and SiO, respectively, and found that shocks were present in most of the GMCs. To compare this with the S-bearing species, we reported the derived gas density and temperature of C₂H, H₃O⁺ and SO, HCN and HNC, and HNCO in Figures 9 and C.1. The physical conditions of the gas from which SiO is emitted could not be compared with our results, as the gas density is usually lower (∼10² − 10⁴ cm⁻³) and the gas temperature much higher (∼500 − 600 K) than what we derived here for the S-bearing species.

Fig. 9.

Density and temperature contour plots derived from the LVG analysis for representative regions in the outer CMZ (GMC9a and GMC2b) and for the inner CMZ ( GMC6 and SSC2). The contours show the 1σ solutions obtained for the minimum value of the χ² in the column density parameter derived for each species (and components) and region. Derived gas density and parameters for other species (HCN, HNC; HNCO; SO, H₃O⁺; CCH) studied in previous ALCHEMI studies were plotted for comparison: [1] Behrens et al. (2022), [2] Huang et al. (2023), [3] Holdship et al. (2022), [4] Holdship et al. (2021).

In the outer CMZ, the derived gas density and temperature for SO usually correlates well with that derived for HCN and HNC. For GMC 10, the gas density traced by SO is reaching that of HCN and HNC at the very high range of the gas temperature (T_gas > 200 K). No clear trend is seen between the other S-bearing species and HNCO. For C₂H and H₃O⁺/SO,, studies were performed only for GMC3 to GMC7.

In the inner CMZ, the gas properties traced by the S-bearing species are quite different between regions. First, OCS and H₂S trace the same gas conditions than HNCO, except for SSC5 in the case of OCS. However, the result towards SSC5 could have been impacted by the fact that we had to fix the emission size due to the low number of OCS lines available (because of blending) in this region. The low-E_u component of CS and H₂CS also show the same gas physical conditions across the inner CMZ and do not seem to consistently trace the same type of gas of any of the tracers from the previous ALCHEMI studies. Then, towards GMC6 and SSC5, the high-E_u component of CS does not trace the same gas as any other tracer whilst same physical conditions with multiple tracers are found towards SSC2 and GMC7 (e.g. C₂H, HNCO). The high-E_u component of H₂CS and SO also shows different behaviour as they probe the same gas conditions as HNCO in SSC2, SSC5, and GMC6 but not towards GMC7. Finally all the S-bearing species and components in SSC2 show similar physical conditions to those of HNCO, which could indicate that this region is largely governed by slow shocks. The region seems to be indeed located at the position of the south-west (SW) streamer associated with the large-scale outflow of the galaxy (e.g. Walter et al. 2017; Zschaechner et al. 2018; Krieger et al. 2019). Enhancement of shocks at the position of SSC2 are likely occurring coherently with the location of the base of the SW streamer (e.g. Bao et al. 2024).

We note the discrepancy between the physical parameters we derive for SO towards GMC7 and SSC2 and what was previously found in Holdship et al. (2022). This could result from the fact that we performed different LVG analyses (SO versus H₃O⁺/SO) and that we did not take into account the SO transitions with A_ij ∼ 10⁻⁶ s⁻¹ as they could not be well reproduced by the LVG models (see Sect. C.3).

Overall, from comparing the physical conditions of the gas traced by the various S-bearing species with that of previous ALCHEMI studies, we can get a first hint of the type of gas traced by some of the S-bearing species. Indeed, in the outer CMZ, SO traces the same type of gas (same gas density and temperature) as HCN and HNC which are well probing the GMCs of NGC 253 (Leroy et al. 2015). In the inner CMZ, SSC2 could be a region where shocks are relatively prominent compared to the other regions, leading to the release/formation of the sulphur-bearing species mostly via non-thermal desorption due to the shocks. Throughout the inner CMZ, H₂S and OCS show similar gas density and temperature as the low-velocity shock tracer HNCO. This is also the case for the high-E_u component of H₂CS and SO (also the low-E_u component of SO but the solutions are unconstrained in GMC6 and SSC5), except in GMC7. As none of the S-bearing species trace the same gas as SiO, it is likely that strong and high-velocity shocks are not at the origin of the release of S-bearing species but rather low-velocity and milder shocks (traced by HNCO; Huang et al. 2023). None of the S-bearing species seems to trace consistently the same gas as C₂H.

6.2. Comparison with Galactic environments

A second step towards understanding the origin of emission of the S-bearing species in the CMZ of NGC 253 is to compare our results with what is found in different Galactic environments. In this work, we focused on the most abundant S-bearing species typically found in Galactic star-forming regions (e.g. Hatchell et al. 1998; Li et al. 2015; Widicus Weaver et al. 2017; Humire et al. 2020). As other S-bearing species than those investigated in this work are detected towards NGC 253 (e.g. NS, SO⁺; Martín et al. 2021), we cannot determine the total S-budget of NGC 253. Therefore, we cannot accurately determine the level of S depletion in NGC 253 and assess how it compares to what is found in various Galactic environments. In order to remove this uncertainty, we compare various S-bearing abundance ratios.

We selected the abundance ratios depending on the region of emission of the S-bearing species constrained by the LVG analysis in Sect. 5.2.3. For H₂S, OCS (low-E_u component), H₂CS (both low- and high-E_u components), and SO (both low- and high-E_u components, in the outer and inner CMZ, respectively), we derived a range for the emission size relatively similar, in the range (∼0.1″ − 0.5″). To compare with the Galactic environments, we investigated all the possible abundance ratios involving these six species or their low-/high-E_u components and selected those which showed the most interesting trends helping us understand their region of emission. The selected ratios are: [H₂S]/[OCS], [H₂S]/[SO], [H₂S]/[H₂CS], [OCS]/[H₂CS], [OCS]/[SO], and [H₂CS]/[SO]. For the other species, CS (both low- and high-E_u components), SO (low-E_u component in the inner CMZ), CCS, SO₂, and OCS (high-E_u component), the emission size were either unconstrained, or assumed to be extended (for the RDs). For these species, after investigating the various abundance ratios available, we selected the following ones: [CS]/[CCS], [SO]/[SO₂], [CS]/[SO], [OCS]/[SO₂], and [OCS/CS]. The range of values for each abundance ratios and for each region of NGC 253 are shown in Table D.1.

For the choice of Galactic regions, we selected those most relevant to compare with the starburst CMZ of NGC 253. This includes regions dominated by shocks and outflows, HCs, PDRs, and molecular clouds. We would like to raise a word of caution as we are comparing GMC-size scales in NGC 253 with parsec and sub-parsec scales in the various Galactic environments. Nonetheless, by performing this comparison, we aim to perform a first order comparison of the abundance ratios found towards the CMZ of NGC 253 and those derived in different Galactic environments. We selected the regions for which the column densities or abundances of several S-bearing species relevant to this work were available.

For the shocks’ and outflows’ dominated regions, we selected the prototypical L1157-B1 bow shock, a well-known shocked region within the protostellar outflow of L1157-mm, and the Orion-KL (Kleinmann-Low) Plateau regions located in the well-studied Orion high-mass star-forming region, which is a mixture of shocks, outflows, and interactions with the ambient cloud (Blake et al. 1987; Wright et al. 1996; Lerate et al. 2008; Esplugues et al. 2013). These two different shocked/outflow regions are representative of two different shock scenarios: the L1157-B1 bow shock could be representative of sporadic shocks due to star formation, whilst the Orion-KL plateau could be representative of the interaction between outflows and an ambient cloud, hence faster (up to 100 km s⁻¹ for Orion Plateau; e.g. O’Dell 2001; Goddi et al. 2009 vs. 20–40 km s⁻¹ for L1157-B1; e.g., Lefloch et al. 2010; Viti et al. 2011; Benedettini et al. 2013, 2021) and more turbulent shocks. These two scenarios are part of the three shock scenarios Huang et al. (2023) proposed to explain the SiO and HNCO emission in the CMZ of NGC 253. If the abundance of S-bearing species varies between these two types of shocked regions, we could put constraints on which of the shocked scenarios proposed by Huang et al. (2023) is the more likely to happen in the CMZ of NGC 253.

Concerning the HCs, we selected the Orion KL (Feng et al. 2015; Luo et al. 2019), Sagittarius B2 N core (Neill et al. 2014), and Cyg X-N12 (El Akel et al. 2022). For Cyg X-N12, we used the column densities derived for the warm chemistry (Table 5 in El Akel et al. 2022). Cyg X-N12 is close to another HC, Cyg X-N30, for which abundances of sulphured species were also derived. However, the authors classified Cyg-X N30 as a “non-traditional” HC as the molecular richness towards this source are likely the result of the combination of processes rather than only the HC chemistry (van der Walt et al. 2021). We thus did not include this source in our comparison.

For PDRs, we included only highly UV-irradiated PDRs. These types of PDRs are particularly relevant for our study as we expect high-UV irradiation to occur towards the CMZ of NGC 253 due to the intense star formation. The presence of strong PDRs was deduced by Harada et al. (2021) using the HCO⁺/HOC⁺ abundance ratios. We thus selected the Orion Bar PDR (Jansen et al. 1995; Leurini et al. 2006) and Monoceros R2 (Mon R2; Ginard et al. 2012) PDRs.

Finally, for the molecular clouds, we selected two different types of clouds. On the one hand, we selected prototypical dark clouds such as L1544 (Vastel et al. 2018), L483 (Lattanzi et al. 2020), and B1-b (Fuente et al. 2016). On the other hand, we compare with Galactic Centre clouds (towards Sgr B2 and A molecular complexes; Nummelin et al. 2000; Martín et al. 2005 and references therein, Armijos-Abendaño et al. 2015). Choosing these two types of clouds allows us to compare our results with cold (T_kin ∼ 10 K) and dense (n_gas ≥ 10⁵ cm⁻³) quiescent clouds (Dark clouds) and with Galactic Centre clouds which are much warmer (T = 50 − 200 K) with an average density of n_gas ∼ 10⁴ cm⁻³ (Güsten & Philipp 2004).

Figure 10 shows the various abundance ratios derived for the S-bearing species in NGC 253 as a function of the selected regions and for various Galactic environments described above. For CCS and SO₂, as we could not confirm whether the multiple temperature component is correct are not due to non-LTE or optically thick effects, we used the lowest and highest value for the beam-averaged column density derived across the components to calculate the abundance ratios. We listed below the main results from this comparison:

The [H₂S]/[OCS_1] abundance ratios of the inner regions are similar to those found in HCs, PDRs, and the two outflows/shocks regions. The ratio seems to increase towards SSC5.
The [H₂S]/[SO_2] and [H₂S]/[H₂CS_2] abundance ratios are similar to what is found in the L1157-B1 shock, HCs, and clouds (only the case for [H₂S]/[H₂CS_2]). Whilst [H₂S]/[SO_2] remains constant across the inner regions, [H₂S]/[H₂CS_2] is the highest towards SSC5 and the lowest towards GMC7.
The [OCS_1]/[H₂CS_1] ratio seems to increase towards the inner CMZ. In the outer CMZ, the ratios derived agree with what is found in most Galactic regions. In the inner CMZ, the ratio is different to what is found in GC clouds in GMC7, whilst it is different to what is found in the L1157-B1 shock towards SSC5.
The [OCS_1]/[SO_1] ratio mostly agrees with what is found in HCs and GC clouds.
The [H₂CS_1]/[SO_2] ratio decreases towards SSC5, where the ratio is similar to what is found in the Orion KL Plateau region. In GMC6 and SSC2, the ratio is similar to those found in the Galactic regions, except the Orion KL Plateau region. In GMC7, the ratio is higher than what is found in the Galactic environments.
The [CS_2]/[CCS] ratio is unconstrained in the outer regions due to the lower limit in beam-average column density derived for CCS (see Sect. 5.2.1). However, in the inner regions, this ratio is constant and agrees well with what is derived in the L1157-B1 shock and the molecular cloud regions.
The [SO_1]/[SO₂] ratio in NGC 253 is similar to most of the Galactic regions, except for the dark clouds and the Orion KL Plateau region.
The [CS_1]/[SO_1] abundance ratio agrees well with what is found in PDRs and GC clouds.
The [OCS_2]/[SO₂] and [OCS_2]/[CS_2] ratios are both constant in the three innermost regions of NGC 253. The [OCS_2]/[SO₂] ratio agrees with all Galactic regions, whilst the [OCS_2]/[CS_2] ratio mostly agrees with dark clouds and the L1157-B1 shock (only towards GMC6).

Fig. 10.

Abundance ratio of the S-bearing species in each region of NGC 253 (red shaded area), and in various Galactic environments. For NGC 253, the listed sequence of clouds follows the layout of the GMCs of the CMZ as shown in Figs. 1–3. The ratios are between species for which the emission size is similar, and is constrained (top panel) or unconstrained/extended (≥1.6″; bottom panel). The underscore numbers 1 and 2 refer to the low- and high-E_u component, respectively. For the outflows and shocks (blue shaded area), we compare with the prototypical shock L1157-B1 (Holdship et al. 2019, and references therein) and the Orion KL Plateau (Persson et al. 2007; Tercero et al. 2010; Esplugues et al. 2013; Crockett et al. 2014a,b; Gong et al. 2015). The hot cores (HCs; yellow shaded area) are Orion KL (Feng et al. 2015; Luo et al. 2019), Sgr B2 (N) (Neill et al. 2014), and Cyg X-N12 (El Akel et al. 2022). The PDRs (green shaded region) are the Orion Bar (Jansen et al. 1995; Leurini et al. 2006) and Mon R2 (Ginard et al. 2012), two high-UV irradiated PDRs. Finally, the molecular clouds category (purple shaded area) comprises dark clouds on the one hand (L1544, L483, and B1-b; Fuente et al. 2016; Vastel et al. 2018; Lattanzi et al. 2020), and Galactic Centre (GC) Clouds on the other hand (Nummelin et al. 2000; Martín et al. 2005 and references therein, Armijos-Abendaño et al. 2015). Lower and upper limits are marked with filled triangles. The legend on the bottom part of the plot is associated with NGC 253. For Galactic environments, the same colour code applies.

6.3. Assembling the puzzle: Regions traced by S-bearing species

After comparing the physical conditions of the gas emitting the S-bearing species in Sect. 6.1 and comparing their abundance ratios with prototypical Galactic regions in Sect. 6.2, we can attempt to understand what S-bearing species trace in the CMZ of NGC 253.

6.3.1. Primary shock tracers

First, from Sect. 6.2, we found that the [H₂S]/[OCS_1] and [H₂S]/[SO_2] abundance ratios all agree with the abundance ratios found in the L1157-B1 shock and HCs. This suggests that H₂S, the low-E_u component of OCS, and the high-E_u component of SO emit from either a shocked gas due to star formation or a warm gas due to the presence of proto-SSCs in the inner CMZ.

In HCs, we expect temperatures of a few hundred Kelvins, much higher than the temperatures derived in this study for H₂S and the low-E_u component of OCS. Additionally, the gas temperature derived for these species is lower than their desorption temperature (T_des) (T_des > 100 K for H₂S and T_des ∼ 150 K for OCS; e.g. Collings et al. 2004; Jiménez-Escobar & Caro 2011; Jiménez-Escobar et al. 2014; Cazaux et al. 2022) in most regions. If we assume that H₂S and the low-E_u component of OCS trace the same gas component in all the regions of the CMZ, then we can conclude that thermal desorption is not the main process explaining the H₂S emission in the CMZ of NGC 253. On the other hand, the physical conditions of H₂S and the low-E_u component of OCS are similar to those of HNCO (only in the inner CMZ for OCS) and the derived emission size of these two species, which is a few parsec (3-8 pc) scale, are consistent with a shock origin. This scenario can be further strengthened for OCS: the derived gas temperatures are lower in the outer CMZ (T_gas < 20 K) compared to the inner CMZ (T_gas = 30 − 50 K), indicating that the shocks traced by OCS could be older in the outer CMZ as the gas needs time to cool down after the shock passage, consistent with the shock timescale derived by Huang et al. (2023), where older shocks are present towards the outer CMZ and younger shocks towards the inner CMZ. Finally, although the physical conditions of the gas emitting OCS and HNCO are different in the outer CMZ, this discrepancy could be explained by a difference in the cooling timescale of the two species. Indeed, since the gas densities derived for HNCO are, on average, lower than those derived for OCS (n_gas ≤ 2 × 10⁴ cm⁻³ vs. ≥2 × 10⁴ cm⁻³ for OCS), the gas traced by HNCO would cool slower than the gas traced by OCS. Hence, HNCO would thus trace a warmer gas than OCS, as in the case of the outer CMZ. Another scenario to explain the difference in gas physical conditions of OCS and HNCO in the outer CMZ could be that the two species probe different types of shocks (e.g. sporadic shocks, outflow-induced shocks, cloud-cloud collision), as proposed by Huang et al. (2023).

Then, we can also conclude that there is a shock origin for the high-E_u component of SO. The LVG analysis showed that this SO component has a high density (10⁵ − 10⁶ cm⁻³). The gas temperature is either lower than ∼60 K (GMC 6 and SSC2) or higher than 60 K (GMC7 and SSC5). Finally, the species arise from a compact region (∼3 pc on average) and the physical conditions are usually similar to those derived for HNCO in most regions. SO is known to be a post-shocked gas tracer as it is a product of the subsequent reaction of H₂S with OH in the gas phase (e.g. Pineau des Forêts et al. 1993; Charnley 1997; Viti et al. 2001). More recently, Taquet et al. (2020) suggested that SO could be abundant both in the front shock and behind the shock. The absence of this component in the outer CMZ would indicate that SO is present in the front shock and has not yet been fully converted to SO₂. Additionally, SO is thought to co-desorb with H₂O from the ice grain mantles (e.g. Viti et al. 2004) so thermal desorption of SO is unlikely to occur in GMC6 and SSC2 where the gas temperature is lower than 100 K. If we assume SO is tracing the same feature throughout the inner CMZ, then this SO component is likely tracing shocks.

Finally, we also propose that the low-E_u component of H₂CS traces shocks. As OCS, H₂S, and the high-E_u component of SO, the low-E_u component of H₂CS emits on ∼3 pc scales. In the outer CMZ, the physical conditions traced by H₂CS are close to that of OCS, although H₂CS shows a lower density most of the time (≥10⁴ cm⁻³ for OCS and ∼10³ − 10⁵ cm⁻³ for H₂CS) and a larger range of gas temperature (T_gas = 10 − 70 K for H₂CS versus T_gas = 12 − 16 K for OCS). However, whilst OCS and H₂S are known to trace mainly front shock events as they would directly come off the ice mantles of grains (e.g. Hatchell et al. 1998; Podio et al. 2014; Taquet et al. 2020), H₂CS rather traces post-shocked gas as it would be formed in the gas phase after the release of H₂S and OCS (e.g. Millar & Herbst 1990; Charnley 1997; Hatchell et al. 1998; Codella et al. 2005). The abundance ratios [OCS_1]/[H₂CS_1] and [H₂CS_1]/[SO_2] do not give much constrain as the ratios are similar to those found in most Galactic regions, but they indicate that the amount of low-E_u component of H₂CS decreases towards the inner CMZ, where shocks are younger (Huang et al. 2023), which is consistent with H₂CS tracing the post-shocked gas rather than the front shocked gas. The difference seen in the systemic velocities of the low-E_u component of H₂CS compared to the rest of the species (see Sect. 5.1.1), also strengthens this scenario.

In summary, H₂S, the low-E_u component of OCS and H₂CS and the high-E_u component of SO likely trace shocks, although they do not seem to trace the same shock component or event based on the information we could gather. However, none of these S-bearing species are probing fast and, thus, strong shocks since they do not trace the same gas as SiO, a well-known fast shock tracer (e.g. Schilke & Walmsley 1997; Gusdorf et al. 2008; Kelly et al. 2017), and as most of the abundance ratios investigated here agree on average with what is found in the L1157-B1 shock rather than in the Orion KL Plateau. In the region SSC5, however, some of the ratios agree better with what is found towards the Orion KL Plateau which could indicate the presence of stronger shocks (likely due to a more intense star formation activity) compared to the rest of the regions. Overall, our results seem to agree with a scenario where H₂S, the low-E_u components of OCS, H₂CS, and the high-E_u component of SO probe shocks which are likely associated with star formation in both the outer and inner CMZ, although star-formation is known to be particularly intense in the inner CMZ where most of radio continuum sources are located (Ulvestad & Antonucci 1997). Our results are consistent with previous low-angular resolution studies of S-bearing species which concluded that S-bearing species are associated with low-velocity shocks similar to those found in the Sgr B2 molecular cloud complex, where star-formation is ongoing (Martín et al. 2003, 2005). Nonetheless, performing chemical modelling will be crucial to confirm these results and assess more accurately the nature of the shocks traced by these species.

6.3.2. Primary molecular cloud tracers

The [OCS_1]/[SO_1] abundance ratio derived in the outer CMZ is constant and agrees consistently with what is found in HCs and GC clouds. However, we just determined that the low-E_u component of OCS arises likely from shocks, indicating that the low-E_u components of OCS and SO do not trace the same gas. The physical conditions of the gas emitting the low-E_u component of SO correlate well with that traced by HCN and HNC (Behrens et al. 2022), which are well-known dense gas tracers of the GMCs of NGC 253 (Leroy et al. 2015). From the LVG analysis, this component of SO emits from a warm (≥40 K) and moderately dense (n_gas ∼ 10⁴ − 10⁵ cm⁻³) gas in most regions. The physical conditions derived are similar to that of the high-density component of the NGC 253 GMCs (T_gas ∼ 85 K and n_gas ∼ 10⁴ cm⁻³) derived by Tanaka et al. (2024) suggesting that SO traces the dense clouds, in particular the compact regions (≤0.3 − 0.4″ or ≤5 − 7 pc) as derived with the LVG analysis in the outer CMZ. This conclusion agrees with a previous study of SO towards NGC 253 by Holdship et al. (2022), who concluded that the SO abundance was unlikely driven by PDRs or shocks. Finally, the physical conditions of the low-E_u component of SO in the inner CMZ being different than those in the outer CMZ (cooler gas temperature and larger emission region) could either indicate a different origin of emission of SO in the inner CMZ or could result from the influence of the higher concentration of CRs in the inner towards the inner CMZ (Mangum et al. 2019; Behrens et al. 2022), which are known to destroy SO (Bayet et al. 2011; Holdship et al. 2022).

The low-E_u component of CS also very likely traces the dense gas from molecular clouds. The [CS_1]/[SO_1] abundance ratio is similar to that found in the PDRs and GC clouds. The low-E_u component of CS emits from a cold (< 20 K) and dense (≥10⁴ cm⁻³) gas on GMC scales (θ_s ≥ 1″ or ∼17 pc) which would favour the molecular cloud origin rather than the PDR one. These results needs to be taken with caution as they originate from the fact that we initially fixed the emission size to 1.6″ (see Sect. 5.2.3). If the emission region were smaller, the physical conditions derived could differ. Nonetheless, our results suggest that the low-E_u component of CS traces the dense GMCs of NGC 253, thus agreeing with previous observations (e.g. Bayet et al. 2009; Martín et al. 2009b; Aladro et al. 2011; Leroy et al. 2015; Holdship et al. 2022).

6.3.3. Shocks, proto-SSCs or molecular clouds: Unconstrained origin

For some of the species or components, we cannot constrain the region of emission. However, it is still possible to propose some hypotheses based on our results.

For the high-E_u component of H₂CS and OCS, and SO₂, we derived temperatures high enough to be consistent with thermal desorption (T_des ≥ 100 K⁶; e.g. Collings et al. 2004; Viti et al. 2004; Kaňuchová et al. 2017; Vidal & Wakelam 2018). These hot components are only present towards GMC6, SSC2, and SSC5 (also GMC7 in the case of H₂CS) and could thus arise from the proto-SSCs. This hypothesis could be strengthened by the fact that such high temperatures for H₂CS and SO₂ were found in higher angular resolution observations towards the proto-SSCs of NGC 253 by Krieger et al. (2020a). However, we cannot rule out a possible shock origin. Young shocks create dense and warm environments by compressing and heating the surrounding gas before cooling (e.g. Frank et al. 2014). In the example of the L1157 outflow, the younger shock position called B0 is much warmer (T > 120 K) than the older shock positions B1 and B2 (T < 80 K; e.g. Lefloch et al. 2012; Feng et al. 2022). If both the low-E_u and high-E_u components of OCS trace shocks in the inner CMZ, then the high-E_u component would trace hot young shocks in the inner CMZ whilst the low-E_u component traces older cooled shocks. In the case of H₂CS, this species is efficiently formed on the grains at high densities (e.g. Jiménez-Escobar & Caro 2011; Drozdovskaya et al. 2018; Laas & Caselli 2019) and strong enough shocks could directly sputter H₂CS into the gas phase. The higher concentration of the fast shock tracer SiO towards the inner CMZ (Huang et al. 2023) could indicate the presence of stronger shocks in these regions. In Sect. 5.1.1, we found a significant difference between the systemic velocity of the low- and high upper-level energy transitions of H₂CS, suggesting the two components being emitted from two different gas or from different shock events. Finally, SO₂ is a well-known post-shocked gas tracer due to the reaction between SO and OH in the gas phase, which takes place within 10⁴ yr (e.g. Pineau des Forêts et al. 1993). The sole detection of SO₂ in the inner CMZ could be explained by the fact that for timescales ≥3 × 10⁴ yr, SO₂ is converted to atomic sulphur (e.g. Pineau des Forêts et al. 1993; Hatchell et al. 1998; Wakelam et al. 2004b; Feng et al. 2015). Since Huang et al. (2023) derived that the shocks in the outer CMZ have ∼10⁵ yr time scales. SO₂ could have been already converted to atomic sulphur. The [H₂S]/[H₂CS_2] and [OCS_2]/[SO₂] abundance ratios do not give further constrain on the possible region of emission of these species as the values are similar to that of most Galactic environments.

Finally, for the high-E_u component of CS and CCS, we cannot distinguish between a molecular cloud origin or a shock/outflow origin. The [OCS_2]/[CS_2] abundance ratio is consistently similar to those found in dark clouds throughout the three innermost regions of the CMZ. We mentioned that the high-E_u component of OCS could either arise from shocked regions or proto-SSCs. The abundance ratio indicates that the high-E_u component of CS and OCS are not arising from the same gas or the same shock event, if they both trace shocks. If CS traces shocks, it would likely trace a post-shocked gas following the destruction of OCS and H₂S (e.g. Charnley 1997; Wakelam et al. 2004b, 2005; Codella et al. 2005; Viti et al. 2014; Gómez-Ruiz et al. 2015; Taquet et al. 2020). On the other hand, the physical conditions of the high-E_u CS component (n ≥ 10⁵ cm⁻³ and T_gas ∼ 20 − 100 K) could indicate that it traces more compact and dense regions of the GMCs than the low-E_u transitions. As mentioned already, the physical conditions of CS need to be taken with caution as we had to fix the source emission size for the LVG analysis.

The [CS_2]/[CCS] abundance ratio agrees well with the ratios found in the L1157-B1 shock and the molecular cloud regions. Since we mentioned the high-E_u component of CS traces either dense molecular gas or post-shock gas, CCS could also trace these types of regions: On the one hand, CCS is widely detected in the galactic young and cold dark clouds (e.g. Suzuki et al. 1992; Rathborne et al. 2008; Hirota et al. 2009; Roy et al. 2011; Marka et al. 2012; Tatematsu et al. 2014; Seo et al. 2019; Koley 2022) and the excitation temperatures derived for the various components (T_ex ∼ 10 − 72 K) could correspond to cold dark clouds but also to GMCs (for which the temperature can be several tens of Kelvins). On the other hand, CCS could originate from shocked gas. In the L1157-B1 shock, Holdship et al. (2019) also derived two different excitation temperatures (∼6 and 50 K) which are within the range of what we derived in this work. Additionally, the low-E_u component of CCS seem to be well correlated with the systemic velocity and line width of that of the low-E_u component of H₂CS (see Sect. 5.1) which could strengthen the shock origin of the emission of this species. Additionally, if CCS mainly traces shocks, its concentration towards the inner CMZ could be explained by the fact that its lifetime is thought to be around 10⁴ yr (de Gregorio-Monsalvo et al. 2006) which is lower than the age of the shock induced star formation derived by Huang et al. (2023) towards the outer CMZ (∼10⁵ yr). It is also possible that the various components of CCS are emitted from different gas. This would agree with the difference in systemic velocity seen between the low- and high- upper-level energy transitions in Sect. 5.1.1.

In summary, we found that most of the S-bearing species could be associated with shocks due to star formation activity within the CMZ of NGC 253. Figure 11 shows a schematic of the proposed scenario for the outer and inner CMZ, as well as a summary of the average physical conditions traced by each species. In the inner CMZ, the presence of stronger and younger shocks due to an intense star forming region would lead to the emission of most S-bearing species, although we cannot exclude that the presence of proto-SSCs contributes to thermally desorb some of the S-bearing species. In the outer CMZ, shocks could also account for the emission of most of S-bearing species, although they are likely to be less strong or older (due to less intense/older star formation) than those from the inner CMZ. Performing chemical modelling will be the next crucial step to disentangle the various possible scenarios for these species (Bouvier et al., in prep.).

Fig. 11.

Schematic (not to scale) showing the proposed origin of emission of the S-bearing species towards the CMZ of NGC 253 (Sect. 6.3) and summarising the main results for the physical properties (T_gas, n_gas) of the gas emitting the S-bearing species (Sect. 5.2). We illustrate two scenarios, for a giant molecular cloud in the outer (left) CMZ and inner (right) CMZ. The timescale for the shocks induced by star-formation are from Huang et al. (2023). The temperatures and densities reported are averaged over the regions and variations could be present in some of the regions.

6.4. Other considerations from the general morphology of the CMZ

We mainly discussed the emission and origin of sulphur-bearing species in the context of ongoing star formation, since these species are known to be enhanced in these regions in the Milky Way. However, the CMZ of NGC 253 is relatively turbulent and other physical effects could be involved in the enhancement of S-bearing species.

First, as mentioned in Sect. 2, shocks could be caused by large-scale cloud-cloud collisions, induced by orbital intersections between the central bar and the spiral arms (Harada et al. 2022; Humire et al. 2022; Huang et al. 2023), the streamers from the large-scale outflow (Walter et al. 2017; Tanaka et al. 2024; Bao et al. 2024), or the superbubles from supernovae remnants (e.g., Sakamoto et al. 2006; Gorski et al. 2017; Krieger et al. 2020b). Orbital intersections are thought to happen at the north-east and SW part of the CMZ, that is, close to the positions we labelled as GMC 1a, 1b, 8a, 9a, and 10 (see e.g., Sorai et al. 2000; Das et al. 2001; Levy et al. 2022 for more details of the various bar orbits). Previous ALCHEMI studies explained the enhancement of HOCO⁺ and HNCO (Harada et al. 2022; Huang et al. 2023), and the presence of CH₃OH masers (Humire et al. 2022) with potential cloud-cloud collisions happening in these regions. Among the S-bearing species analysed in this study, OCS is the species for which the emission distribution is the most resembling that of the above-mentioned species. However, our analysis showed that in the outer CMZ, hence where these cloud-cloud collisions occur, OCS and HNCO do not trace the same gas (as we also mentioned in Section 6.3.1), which could indicate that these two species are likely not tracing the same shock event or type. Further investigations about the type of shock using chemical modelling (Bouvier et al., in prep.) are needed to conclude whether cloud-cloud collisions are responsible for the OCS emission in the outer CMZ.

Then, studies performed on the large-scale outflow of NGC 253 and its link with the CMZ showed the presence of molecular streamers emerging from the centre of the CMZ (Bolatto et al. 2013; Walter et al. 2017; Krieger et al. 2019), one of them, in the SW being particularly prominent. These streamers are powered by the ongoing star formation occurring inside the inner GMCs and the presence of strong shocks is detected in these GMCs, as well as at the base of the streamers, near the GMCs 7, 6, and 3 (corresponding to our region SSC2) (e.g. Walter et al. 2017; Krieger et al. 2019; Huang et al. 2023; Bao et al. 2024). The present analysis of the S-bearing also indicates the possibility of stronger shocks happening in the inner CMZ (see Sects. 6.3.1 and 6.3.3, and Fig. 11.), due to more enhanced ongoing star formation in the inner CMZ.

Finally, in the inner CMZ, other types of sources associated with star formation are present, due to the inflow of gas along the bar which triggers star formation inside the GMCs. Several supernovae explosions, supernovae remnants (SNRs), and HII regions were identified in the inner CMZ (Ulvestad & Antonucci 1997), which could contribute to shocking and heating the molecular gas. Whilst GMC6 is bright in all the S-bearing species, differences can be seen between the positions SSC5 and SSC2. Indeed, the emission of CS, H₂S, and CCS seems to be slightly more intense towards SSC5 compared to SSC2 (see Figs. 1 and 3). Behrens et al. (2022) showed that a higher number of HII regions was located towards SSC5 (GMC 4 in their study) compared to the other inner GMCs. Enhancement of CS, H₂S, and CCS due to the presence of these HII regions is however unlikely, as we showed in Sect. 6.3 that the emission of S-bearing species is not associated with PDRs, powered by HII regions. On the other hand, Behrens et al. (2022) showed that SNRs are more numerous towards GMC 6 compared to the other GMCs studied in this paper, but there is no distinct patterns seen in the gas properties or column densities of the S-bearing species that could indicate a specific enhancement in GMC 6 due to the presence of SNRs.

7. Conclusions

We investigated the origin of the emission of the most abundant sulphur-bearing species (CS, H₂S, SO, SO₂, H₂CS, OCS, and CCS) detected in Galactic star-forming regions towards the CMZ of the nearby starburst galaxy NGC 253. Thanks to the high sensitivity and broad frequency range of ALCHEMI, we had access to a variety of transitions for each species, allowing us to perform rotational diagram and LVG analyses. We then first compared our results with other various species from previous ALCHEMI studies, and we then compared several molecular abundance ratios found in NGC 253 with those found in various Galactic regions. We summarise the main results below.

CS, H₂S, SO, SO₂, H₂CS, OCS, and CCS show different emission distributions and peaks. The low upper-level energy transitions of CS and H₂S show the most extended emission throughout the CMZ. On the other hand, SO₂ is concentrated within the 6″ region around the kinematic centre of the CMZ. OCS and H₂CS emissions peak towards each of the GMCs rather than towards the centre of the CMZ as the other species.
We performed Gaussian fits to derive the line widths and systemic velocities. We found significant differences in the systemic velocity peak of the low and high upper-level energy transitions of H₂CS, ³⁴SO, and CCS, indicating differences in the emission which are likely due to excitation conditions, and between the low upper-level energy transitions of H₂CS and CCS with the other species and their higher upper-level energy transitions. The narrowest line widths for each species were found towards GMC 10.
We performed both LTE and non-LTE LVG analyses to derive the physical conditions of the gas emitting the sulphur-bearing species. We found that towards the inner CMZ, we need at least two gas phase components (associated with low and high upper-level energy transitions) to reproduce the observed line intensities indicating the presence of different excitation conditions between the inner and the outer CMZ
Across the CMZ, the low upper-level energy transitions of CS and H₂CS trace the coldest gas (T_gas < 50 K). In the outer CMZ, the highest gas temperatures were derived for SO and the high J transitions of CS; whereas, in the inner CMZ H₂S and the high upper-level energy transitions of SO and H₂CS show the highest gas temperature. Towards SSC2, the overall gas temperature is colder (T_gas < 50 K) than in the other inner regions. The highest density is traced by H₂S (n_gas ≥ 10⁶ − 10⁷ cm⁻³) throughout the CMZ. SO and the low upper-energy level transitions of CS trace higher densities in the inner CMZ; whereas, for the other species, the density stays relatively constant.
We found that most of the sulphur-bearing species investigated in this work are likely tracing shocks throughout the CMZ. We also found that they are not likely probing strong shocks as SiO would do. Only the low upper-level energy transitions of SO and CS are tracing dense molecular gas. In the inner CMZ, where the proto-SSCs are located, we could not always disentangle whether the S-bearing species are released through shocks or thermal evaporation. Chemical modelling will be a crucial next step to distinguish between the two possibilities (Bouvier et al., in prep.).

We investigated the region of emission of various S-bearing species assuming that they are emitting from the same gas throughout the CMZ, or distinguishing between the inner and outer part of the CMZ. However, as seen with the LVG analysis (Sect. 5.2.3), the physical condition of a species could differ in some regions compared to the others, like in the case of SSC2 (see Sect. 6.1). Harada et al. (2024) also concluded that the various GMCs of the CMZ were distinct following their PCA analysis. Therefore, our conclusions could be not entirely valid and further investigation would be needed for these specific regions.

¹

http://www.iram.fr/IRAMFR/GILDAS

²

https://zenodo.org/records/11351042

³

The collisional rate for a same temperature and transition differs by one order of magnitude between the LAMDA and BASECOL files, even though both databases are using the Green & Chapman (1978) paper. We notified the LAMDA database team about this issue.

⁴

https://zenodo.org/records/11351042

⁵

Initially developed by Ceccarelli et al. (2003).

⁶

New calculations of the binding energies of SO₂ (Perrero et al. 2022) indicate that SO₂ can sublimate at a lower temperature (Bianchi et al. 2023).

Acknowledgments

We thank the anonymous referee for their comments that helped improve the paper. This work is founded by the European Research Council (ERC) Advance Grant MOPPEX 833460. L.C. and V.M.R. acknowledge support from the grants No. PID2019-105552RB-C41 and PID2022-136814NB-I00 by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033 and by ERDF, UE. VMR also acknowledges support from the grant number RYC2020-029387-I funded by MICIU/AEI/10.13039/501100011033 and by “ESF, Investing in your future”, and from the Consejo Superior de Investigaciones Científicas (CSIC) and the Centro de Astrobiología (CAB) through the project 20225AT015 (Proyectos intramurales especiales del CSIC). N.H. acknowledges support from JSPS KAKENHI grant No. JP21K03634. MG gratefully acknowledges Swedish Research Council grant 621-2011-4143. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00161.L and ADS/JAO.ALMA#2018.1.00162.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.

References

Agúndez, M., & Wakelam, V. 2013, Chem. Rev., 113, 8710 [Google Scholar]
Aladro, R., Martín-Pintado, J., Martín, S., Mauersberger, R., & Bayet, E. 2011, A&A, 525, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Aladro, R., Martín, S., Riquelme, D., et al. 2015, A&A, 579, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ando, R., Nakanishi, K., Kohno, K., et al. 2017, ApJ, 849, 81 [CrossRef] [Google Scholar]
Armijos-Abendaño, J., Martín-Pintado, J., Requena-Torres, M. A., Martín, S., & Rodríguez-Franco, A. 2015, MNRAS, 446, 3842 [CrossRef] [Google Scholar]
Artur de la Villarmois, E., Guzman, V. V., Yang, Y.-L., Zhang, Y., & Sakai, N. 2023, A&A, 678, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bachiller, R., Pérez Gutiérrez, M., Kumar, M. S. N., & Tafalla, M. 2001, A&A, 372, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bao, M., Harada, N., Kohno, K., et al. 2024, A&A, 687, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bayet, E., Aladro, R., Martin, S., Viti, S., & Martin-Pintado, J. 2009, ApJ, 707, 126 [NASA ADS] [CrossRef] [Google Scholar]
Bayet, E., Williams, D. A., Hartquist, T. W., & Viti, S. 2011, MNRAS, 414, 1583 [NASA ADS] [CrossRef] [Google Scholar]
Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119 [NASA ADS] [CrossRef] [Google Scholar]
Bendo, G. J., Beswick, R. J., D’Cruze, M. J., et al. 2015, MNRAS, 450, L80 [NASA ADS] [CrossRef] [Google Scholar]
Benedettini, M., Viti, S., Codella, C., et al. 2013, MNRAS, 436, 179 [Google Scholar]
Benedettini, M., Viti, S., Codella, C., et al. 2021, A&A, 645, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bianchi, E., López-Sepulcre, A., Ceccarelli, C., et al. 2023, Faraday Discuss., 245, 164 [NASA ADS] [CrossRef] [Google Scholar]
Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621 [Google Scholar]
Blake, G. A., Mundy, L. G., Carlstrom, J. E., et al. 1996, ApJ, 472, L49 [NASA ADS] [CrossRef] [Google Scholar]
Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450 [Google Scholar]
Bouscasse, L., Csengeri, T., Belloche, A., et al. 2022, A&A, 662, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cazaux, S., Carrascosa, H., Caro, G. M. M., et al. 2022, A&A, 657, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ceccarelli, C., Maret, S., Tielens, A. G. G. M., Castets, A., & Caux, E. 2003, A&A, 410, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Charnley, S. B. 1997, ApJ, 481, 396 [NASA ADS] [CrossRef] [Google Scholar]
Codella, C., & Bachiller, R. 1999, A&A, 350, 659 [NASA ADS] [Google Scholar]
Codella, C., Bachiller, R., Benedettini, M., et al. 2005, MNRAS, 361, 244 [Google Scholar]
Codella, C., Bianchi, E., Podio, L., et al. 2021, A&A, 654, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Collings, M. P., Anderson, M., Chen, R., et al. 2004, MNRAS, 354, 1133 [NASA ADS] [CrossRef] [Google Scholar]
Crockett, N. R., Bergin, E. A., Neill, J. L., et al. 2014a, ApJ, 781, 114 [NASA ADS] [CrossRef] [Google Scholar]
Crockett, N. R., Bergin, E. A., Neill, J. L., et al. 2014b, ApJ, 787, 112 [Google Scholar]
Dagdigian, P. J. 2020a, J. Chem. Phys., 152, 074307 [NASA ADS] [CrossRef] [Google Scholar]
Dagdigian, P. J. 2020b, MNRAS, 494, 5239 [NASA ADS] [CrossRef] [Google Scholar]
Das, M., Anantharamaiah, K. R., & Yun, M. S. 2001, ApJ, 549, 896 [NASA ADS] [CrossRef] [Google Scholar]
de Gregorio-Monsalvo, I., Gomez, J. F., Suarez, O., et al. 2006, ApJ, 642, 319 [CrossRef] [Google Scholar]
de Jong, T., Chu, S., & Dalgarno, A. 1975, ApJ, 199, 69 [NASA ADS] [CrossRef] [Google Scholar]
Denis-Alpizar, O., Stoecklin, T., Halvick, P., Dubernet, M.-L., & Marinakis, S. 2012, J. Comput. Phys., 137, 234301 [Google Scholar]
Denis-Alpizar, O., Stoecklin, T., Guilloteau, S., & Dutrey, A. 2018, MNRAS, 478, 1811 [Google Scholar]
Drozdovskaya, M. N., van Dishoeck, E. F., Jørgensen, J. K., et al. 2018, MNRAS, 476, 4949 [Google Scholar]
Dubernet, M.-L., Alexander, M. H., Ba, Y. A., et al. 2013, A&A, 553, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
El Akel, M., Kristensen, L. E., Le Gal, R., et al. 2022, A&A, 659, A100 [CrossRef] [EDP Sciences] [Google Scholar]
Esplugues, G. B., Tercero, B., Cernicharo, J., et al. 2013, A&A, 556, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Esplugues, G. B., Viti, S., Goicoechea, J. R., & Cernicharo, J. 2014, A&A, 567, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Esplugues, G., Rodríguez-Baras, M., Andrés, D. S., et al. 2023, A&A, 678, A199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Feng, S., Beuther, H., Henning, T., et al. 2015, A&A, 581, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Feng, S., Codella, C., Ceccarelli, C., et al. 2020, ApJ, 896, 37 [NASA ADS] [CrossRef] [Google Scholar]
Feng, S., Liu, H. B., Caselli, P., et al. 2022, ApJ, 933, L35 [CrossRef] [Google Scholar]
Fontani, F., Roueff, E., Colzi, L., & Caselli, P. 2023, A&A, 680, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 451 [Google Scholar]
Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fuente, A., Rivière-Marichalar, P., Beitia-Antero, L., et al. 2023, A&A, 670, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ginard, D., González-García, M., Fuente, A., et al. 2012, A&A, 543, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ginard, D., Fuente, A., García-Burillo, S., et al. 2015, A&A, 578, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Goddi, C., Greenhill, L. J., Humphreys, E. M. L., et al. 2009, ApJ, 691, 1254 [NASA ADS] [CrossRef] [Google Scholar]
Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Goicoechea, J. R., Aguado, A., Cuadrado, S., et al. 2021, A&A, 647, A10 [EDP Sciences] [Google Scholar]
Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209 [Google Scholar]
Gómez-Ruiz, A. I., Codella, C., Lefloch, B., et al. 2015, MNRAS, 446, 3346 [Google Scholar]
Gong, Y., Henkel, C., Thorwirth, S., et al. 2015, A&A, 581, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gorski, M., Ott, J., Rand, R., et al. 2017, ApJ, 842, 124 [CrossRef] [Google Scholar]
Gorski, M. D., Ott, J., Rand, R., et al. 2019, MNRAS, 483, 5434 [NASA ADS] [CrossRef] [Google Scholar]
Green, S., & Chapman, S. 1978, ApJS, 37, 169 [Google Scholar]
Gusdorf, A., Cabrit, S., Flower, D. R., & Pineau Des Forêts, G. 2008, A&A, 482, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Güsten, R., & Philipp, S. D. 2004, Springer Proceedings in Physics, 91, 253 [CrossRef] [Google Scholar]
Haasler, D., Rivilla, V. M., Martín, S., et al. 2022, A&A, 659, A158 [CrossRef] [EDP Sciences] [Google Scholar]
Harada, N., Martín, S., Mangum, J. G., et al. 2021, ApJ, 923, 24 [NASA ADS] [CrossRef] [Google Scholar]
Harada, N., Martin, S., Mangum, J., et al. 2022, ApJ, 938, 80 [NASA ADS] [CrossRef] [Google Scholar]
Harada, N., Meier, D. S., Martín, S., et al. 2024, ApJS, 271, 38 [NASA ADS] [CrossRef] [Google Scholar]
Hatchell, J., & Viti, S. 2002, A&A, 381, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hatchell, J., Thompson, M. A., Millar, T. J., & Macdonald, G. H. 1998, A&AS, 133, 29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hily-Blant, P., Forêts, G. P. D., Faure, A., & Lique, F. 2022, A&A, 658, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hirota, T., Ohishi, M., & Yamamoto, S. 2009, ApJ, 699, 585 [Google Scholar]
Holdship, J., Viti, S., Jimenez-Serra, I., et al. 2016, MNRAS, 463, 802 [NASA ADS] [CrossRef] [Google Scholar]
Holdship, J., Jimenez-Serra, I., Viti, S., et al. 2019, ApJ, 878, 64 [NASA ADS] [CrossRef] [Google Scholar]
Holdship, J., Viti, S., Martín, S., et al. 2021, A&A, 654, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89 [NASA ADS] [CrossRef] [Google Scholar]
Huang, K.-Y., Viti, S., Holdship, J., et al. 2023, A&A, 675, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Humire, P. K., Thiel, V., Henkel, C., et al. 2020, A&A, 642, A222 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Jansen, D. J., Spaans, M., Hogerheijde, M. R., & van Dishoeck, E. F. 1995, A&A, 303, 541 [NASA ADS] [Google Scholar]
Jiménez-Escobar, A., & Caro, G. M. M. 2011, A&A, 536, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Jiménez-Escobar, A., Muñoz Caro, G. M., & Chen, Y.-J. 2014, MNRAS, 443, 343 [CrossRef] [Google Scholar]
Kameno, S., Sawada-Satoh, S., Impellizzeri, C. M. V., et al. 2023, ApJ, 944, 156 [NASA ADS] [CrossRef] [Google Scholar]
Kaňuchová, Z., Boduch, P., Domaracka, A., et al. 2017, A&A, 604, A68 [CrossRef] [EDP Sciences] [Google Scholar]
Kelly, G., Viti, S., García-Burillo, S., et al. 2017, A&A, 597, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Koley, A. 2022, MNRAS, 516, 185 [NASA ADS] [CrossRef] [Google Scholar]
Krieger, N., Bolatto, A. D., Walter, F., et al. 2019, ApJ, 881, 43 [NASA ADS] [CrossRef] [Google Scholar]
Krieger, N., Bolatto, A. D., Leroy, A. K., et al. 2020a, ApJ, 897, 176 [Google Scholar]
Krieger, N., Bolatto, A. D., Koch, E. W., et al. 2020b, ApJ, 899, 158 [Google Scholar]
Kwon, W., Fernández-López, M., Stephens, I. W., & Looney, L. W. 2015, ApJ, 814, 43 [NASA ADS] [CrossRef] [Google Scholar]
Laas, J. C., & Caselli, P. 2019, A&A, 624, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lattanzi, V., Bizzocchi, L., Vasyunin, A. I., et al. 2020, A&A, 633, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lefloch, B., Cernicharo, J., Cabrit, S., & Cesarsky, D. 2005, A&A, 433, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lefloch, B., Cabrit, S., Codella, C., et al. 2010, A&A, 518, L113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lefloch, B., Cabrit, S., Busquet, G., et al. 2012, ApJ, 757, L25 [Google Scholar]
Lerate, M. R., Yates, J., Viti, S., et al. 2008, MNRAS, 387, 1660 [NASA ADS] [CrossRef] [Google Scholar]
Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25 [NASA ADS] [CrossRef] [Google Scholar]
Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2018, ApJ, 869, 126 [NASA ADS] [CrossRef] [Google Scholar]
Leurini, S., Rolffs, R., Thorwirth, S., et al. 2006, A&A, 454, L47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2021, ApJ, 912, 4 [NASA ADS] [CrossRef] [Google Scholar]
Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2022, ApJ, 935, 19 [NASA ADS] [CrossRef] [Google Scholar]
Li, J., Wang, J., Zhu, Q., Zhang, J., & Li, D. 2015, ApJ, 802, 40 [Google Scholar]
Luo, G., Feng, S., Li, D., et al. 2019, ApJ, 885, 82 [NASA ADS] [CrossRef] [Google Scholar]
Mangum, J. G., Ginsburg, A. G., Henkel, C., et al. 2019, ApJ, 871, 170 [NASA ADS] [CrossRef] [Google Scholar]
Maret, S., Hily-Blant, P., Pety, J., Bardeau, S., & Reynier, E. 2011, A&A, 526, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Marka, C., Schreyer, K., Launhardt, R., Semenov, D. A., & Henning, T. 2012, A&A, 537, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Martín, S., Mauersberger, R., Martin-Pintado, J., Garcia-Burillo, S., & Henkel, C. 2003, A&A, 411, L465 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Martín, S., Martín-Pintado, J., Mauersberger, R., Henkel, C., & García-Burillo, S. 2005, ApJ, 620, 210 [CrossRef] [Google Scholar]
Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, ApJS, 164, 450 [Google Scholar]
Martín, S., Martín-Pintado, J., & Mauersberger, R. 2009a, ApJ, 694, 610 [CrossRef] [Google Scholar]
Martín, S., Martín-Pintado, J., & Viti, S. 2009b, ApJ, 706, 1323 [Google Scholar]
Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
McCormick, A., Veilleux, S., & Rupke, D. S. N. 2013, ApJ, 774, 126 [NASA ADS] [CrossRef] [Google Scholar]
Meier, D. S., & Turner, J. L. 2005, ApJ, 618, 259 [NASA ADS] [CrossRef] [Google Scholar]
Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ, 801, 63 [NASA ADS] [CrossRef] [Google Scholar]
Millar, T. J., & Herbst, E. 1990, A&A, 231, 466 [NASA ADS] [Google Scholar]
Mills, E. A. C., Gorski, M., Emig, K. L., et al. 2021, ApJ, 919, 105 [NASA ADS] [CrossRef] [Google Scholar]
Minh, Y. C., Ziurys, L. M., Irvine, W. M., & McGonagle, D. 1990, ApJ, 360, 136 [NASA ADS] [CrossRef] [Google Scholar]
Minh, Y. C., Muller, S., Liu, S.-Y., & Yoon, T. S. 2007, ApJ, 661, L135 [NASA ADS] [CrossRef] [Google Scholar]
Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49 [Google Scholar]
Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215 [Google Scholar]
Müller-Sánchez, F., González-Martín, O., Fernández-Ontiveros, J. A., Acosta-Pulido, J. A., & Prieto, M. A. 2010, ApJ, 716, 1166 [CrossRef] [Google Scholar]
Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Neill, J. L., Bergin, E. A., Lis, D. C., et al. 2014, ApJ, 789, 8 [Google Scholar]
Nummelin, A., Bergman, P., Hjalmarson, A., et al. 2000, ApJS, 128, 213 [NASA ADS] [CrossRef] [Google Scholar]
O’Dell, C. R. 2001, ARA&A, 39, 99 [Google Scholar]
Ospina-Zamudio, J., Lefloch, B., Favre, C., et al. 2019, MNRAS, 490, 2679 [NASA ADS] [CrossRef] [Google Scholar]
Pérez-Beaupuits, J. P., Güsten, R., Harris, A., et al. 2018, ApJ, 860, 23 [Google Scholar]
Perrero, J., Enrique-Romero, J., Ferrero, S., et al. 2022, ApJ, 938, 158 [NASA ADS] [CrossRef] [Google Scholar]
Persson, C. M., Olofsson, A. O. H., Koning, N., et al. 2007, A&A, 476, 807 [EDP Sciences] [Google Scholar]
Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spectr. Rad. Transf., 60, 883 [Google Scholar]
Pineau des Forêts, G., Roueff, E., Schilke, P., & Flower, D. R. 1993, MNRAS, 262, 915 [CrossRef] [Google Scholar]
Podio, L., Lefloch, B., Ceccarelli, C., Codella, C., & Bachiller, R. 2014, A&A, 565, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Price, T. J., Forrey, R. C., Yang, B., & Stancil, P. C. 2021, J. Chem. Phys., 154, 034301P [NASA ADS] [CrossRef] [Google Scholar]
Rathborne, J. M., Lada, C. J., Muench, A. A., Alves, J. F., & Lombardi, M. 2008, ApJS, 174, 396 [NASA ADS] [CrossRef] [Google Scholar]
Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330 [NASA ADS] [CrossRef] [Google Scholar]
Rico-Villas, F., Martin-Pintado, J., Gonzalez-Alfonso, E., Martin, S., & Rivilla, V. M. 2020, MNRAS, 491, 4573 [NASA ADS] [CrossRef] [Google Scholar]
Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16 [Google Scholar]
Roy, N., Datta, A., Momjian, E., & Sarma, A. P. 2011, ApJ, 739, L4 [NASA ADS] [CrossRef] [Google Scholar]
Sakamoto, K., Ho, P. T. P., Iono, D., et al. 2006, ApJ, 636, 685 [NASA ADS] [CrossRef] [Google Scholar]
Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, ApJ, 735, 19 [CrossRef] [Google Scholar]
Sato, M. T., Aalto, S., Kohno, K., et al. 2022, A&A, 660, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schilke, P., Walmsley, C. M., Pineau des Forêts, G., & Flower, D. R., 1997, A&A, 321, 293 [Google Scholar]
Schöier, F., van der Tak, F., van Dishoeck, E. F., & Black, J. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Scourfield, M., Viti, S., Garcia-Burillo, S., et al. 2020, MNRAS, 496, 5308 [NASA ADS] [CrossRef] [Google Scholar]
Scoville, N. Z., & Solomon, P. M. 1974, ApJ, 187, L67 [NASA ADS] [CrossRef] [Google Scholar]
Seo, Y. M., Majumdar, L., Goldsmith, P. F., et al. 2019, ApJ, 871, 134 [Google Scholar]
Sewiło, M., Cordiner, M., Charnley, S. B., et al. 2022, ApJ, 931, 102 [CrossRef] [Google Scholar]
Sorai, K., Nakai, N., Kuno, N., Nishiyama, K., & Hasegawa, T. 2000, PASJ, 52, 785 [NASA ADS] [CrossRef] [Google Scholar]
Spezzano, S., Sipilä, O., Caselli, P., et al. 2022, A&A, 661, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Suzuki, H., Yamamoto, S., Ohishi, M., et al. 1992, ApJ, 392, 551 [Google Scholar]
Tanaka, K., Mangum, J. G., Viti, S., et al. 2024, ApJ, 961, 18 [NASA ADS] [CrossRef] [Google Scholar]
Taquet, V., Codella, C., de Simone, M., et al. 2020, A&A, 637, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tatematsu, K., Hirota, T., Ohashi, S., et al. 2014, ApJ, 789, 83 [NASA ADS] [CrossRef] [Google Scholar]
Tercero, B., Cernicharo, J., Pardo, J. R., & Goicoechea, J. R. 2010, A&A, 517, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Troscompt, N., Faure, A., Wiesenfeld, L., Ceccarelli, C., & Valiron, P. 2009, A&A, 493, 687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Turner, J. L., & Ho, P. T. P. 1985, ApJ, 299, L77 [NASA ADS] [CrossRef] [Google Scholar]
Tychoniec, L., van Dishoeck, E. F., van’t Hoff, M. L. R., et al. 2021, A&A, 655, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ulvestad, J. S., & Antonucci, R. R. J. 1997, ApJ, 488, 621 [NASA ADS] [CrossRef] [Google Scholar]
van der Tak, F. F. S., Boonman, A. M. S., Braakman, R., & van Dishoeck, E. F. 2003, A&A, 412, 133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
van der Walt, S. J., Kristensen, L. E., Jørgensen, J., et al. 2021, A&A, 655, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Vastel, C., Quénard, D., Gal, R. L., et al. 2018, MNRAS, 478, 5514 [NASA ADS] [CrossRef] [Google Scholar]
Vidal, T. H. G., & Wakelam, V. 2018, MNRAS, 474, 5575 [Google Scholar]
Viti, S., Caselli, P., Hartquist, T. W., & Williams, D. A. 2001, A&A, 370, 1017 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Viti, S., Collings, M. P., Dever, J. W., McCoustra, M. R. S., & Williams, D. A. 2004, MNRAS, 354, 1141 [NASA ADS] [CrossRef] [Google Scholar]
Viti, S., Jimenez-Serra, I., Yates, J. A., et al. 2011, ApJ, 740, L3 [NASA ADS] [CrossRef] [Google Scholar]
Viti, S., García-Burillo, S., Fuente, A., et al. 2014, A&A, 570, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wakelam, V., Caselli, P., Ceccarelli, C., Herbst, E., & Castets, A. 2004a, A&A, 422, 159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wakelam, V., Castets, A., Ceccarelli, C., et al. 2004b, A&A, 413, 609 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wakelam, V., Ceccarelli, C., Castets, A., et al. 2005, A&A, 437, 149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265 [NASA ADS] [CrossRef] [Google Scholar]
Widicus Weaver, S. L., Laas, J. C., Zou, L., et al. 2017, ApJS, 232, 3 [NASA ADS] [CrossRef] [Google Scholar]
Wiesenfeld, L., & Faure, A. 2013, MNRAS, 432, 2573 [NASA ADS] [CrossRef] [Google Scholar]
Woods, P. M., Occhiogrosso, A., Viti, S., et al. 2015, MNRAS, 450, 1256 [Google Scholar]
Wright, M. C. H., Plambeck, R. L., & Wilner, D. J. 1996, ApJ, 469, 216 [CrossRef] [Google Scholar]
Yang, B. H., Zhang, P., Qu, C., et al. 2020, Chem. Phys., 532, 110695 [NASA ADS] [CrossRef] [Google Scholar]
Zschaechner, L. K., Bolatto, A. D., Walter, F., et al. 2018, ApJ, 867, 111 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Velocity-integrated maps

In this section, the velocity integrated maps of $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ and ³⁴SO are shown.

Fig. A.1.

Same as Fig. 1 but for $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ and ³⁴SO. Top: Velocity-integrated maps of ³⁴SO (2₃-1₂). Levels start at 3σ (1σ = 50, 28 and 10.4 mJy.beam⁻¹, respectively; white contours) with steps of 7, 3, and 5σ (black contours), respectively. Bottom: Velocity-integrated maps of $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ (1_1, 0-1_0, 1). Level starts at 3 σ (1σ = 56 mJy.beam⁻¹; white contours) with steps of 7σ (black contours).

Appendix B: Targeted species and transitions

Table B.1.

Species and transitions detected and used in this work, as well as their spectroscopic parameters.

Appendix C: LVG modelling and results

We give in this section detailed information about the LVG modelling of each species. The list of collision coefficients used to perform the analysis are presented in Table 2 with their corresponding references. The results of the non-LTE analysis for the regions GMC10, GMC8a, GMC7, SSC5, GMC1a, and GMC1b are displayed in Figures C.1 and in Table C.1.

Fig. C.1.

Density and temperature contour plots for GMC10, GMC8a, GMC7, SSC5, GMC1a, and GMC1b. The contours show the 1σ solutions obtained for the minimum value of the χ² in the column density parameter derived for each species (and components) and region. Derived gas density and parameters for other species (HCN, HNC; HNCO; SO, H₃O⁺; CCH) studied in previous ALCHEMI studies were plotted for comparison: [1] Behrens et al. (2022), [2] Huang et al. (2023), [3] Holdship et al. (2022), [4] Holdship et al. (2021).

Table C.1.

Best-fit results and 1σ confidence level range from the non-LTE LVG analysis for GMC10, GMC7, SSC5, GMC1a, and GMC1b, respectively (from top to bottom).

C.1. CS

From the LTE analysis, two components were found, the first one including the transitions J = 2–1 to J = 4–3 and the second one including the transitions J = 5–4 to J = 7–6. For each component, only three lines were available. Hence, we fixed the size of the emission to be able to constraint the other parameters (N_tot, T_gas, n_gas). As we do not have information on where the CS could be emitted from we chose to fix the size of the emission to the beam size of 1.6″, which corresponds in having a beam filling factor of 0.5. The resulting reduced χ² were good enough ( $χ_{red}^{2} ~ 0.1 - 1.9$ $\chi^2_{{\rm red}}\sim 0.1-1.9$ ) to be considered as possible solutions. It is important to note, however, that other solutions with emission size more compact could be possible.

C.2. H₂S + $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$

In the case of H₂, we also included the line of $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$ at 167 GHz that we detect in several regions. We assumed a value of 10 for the ³²S/³⁴S ratio, as derived in Martin et al. 2021. We could perform the LVG analysis for the inner CMZ and GMC10 where we detect at least two transitions. Except for GMC10 where only two transitions were available, we left the size as a free parameter. For GMC10, we took the size found for GMC7 (0.34″), assuming the emission comes for the same component in the two regions. Finally, as we detected both ortho- and para- transitions of H₂S, we had to assume an ortho-to-para ratio (OPR). The value usually taken in the literature is the equilibrium value of three. However, as some studies found that the OPR could be lower in the case of H₂S, and in particular in Orion where a value of 1.7±0.8 has been measured in Orion (Crockett et al. 2014), we ran two types of models, one with an OPR of three and the other with an OPR of one. We found that the best OPR reproducing the observations in all the regions is one.

C.3. SO+³⁴SO

From the rotation diagram analysis, we found one component in the inner CMZ, and two components in the outer CMZ for SO and ³⁴SO. We thus modelled the two components separately. As for H₂S, we used the ³²S/³⁴S ratio of 10 derived in Martin et al. 2021 to include the ³⁴SO lines in the model. This ratio was reproducing well the observations for the first component in the inner CMZ but for the second component, we found that a lower ratio of 8 was necessary to reproduce the observations. In the outer CMZ, only the 99 GHz line of ³⁴SO is detected. We found that including this line was driving the fit to very low temperatures. We do not have the main isotopologue of this line and its E_u of 9K is lower than that of the other detected transitions. Hence, this transition could trace a different component and we thus chose to not include it in the LVG modelling.

For most of the regions we had enough lines to leave the size as a free parameter, except for the second component of GMC7 where we had only three lines. In this case, we fixed the size to 0.2″ assuming the emission comes from the same component as for the other regions in the inner CMZ (for which we derived sizes between 0.15 and 0.3 ″). On the other hand, for some components of regions, the best fit provided by the LVG analysis provided the product θ_s × N_tot because of a degeneracy between the size and column density. This happened because the all the lines of SO are optically thin (τ < < 1). We used the method described in Sec. 5.2.2 to constrained the column density. Finally, we could not include the SO lines at 100, 136 and 254 GHz with A_ij ∼ 10⁻⁶ s⁻¹, as for the LTE analysis, as these lines are not well reproduced by the LVG model when including them with the other lines. These lines are probably tracing a third component that we could not model due to the poor number of lines.

C.4. H₂CS

From the LTE analysis, up to two components were found in the inner CMZ. We thus modelled the two components separately in these regions. As for H₂S, H₂CS is detected in both ortho- and para- forms. We assumed an OPR ratio of 3, which reproduces well the observations. For most regions, enough lines were available to leave the size as a free parameter but the column density was not always constrained. In these cases, we used the method described in Sec.5.2.2 to better constrain the column density. For the region GMC10, and the second component in GMC7, there are only three lines available, preventing us from leaving the size as a free parameter. In the case of GMC10, we fixed the size to both extended emission and 0.2″ as these two values were found as best-fit values in the other regions. We adopted the solution with the best $χ_{red}^{2}$ $\chi^2_{{\rm red}}$ , which is with a size of 0.2″. For GMC7, we fixed it to 0.12″, which is the average size found in the best-fit of the other regions for this component.

C.5. OCS

From the rotation diagram analysis, we found up to two components in the inner CMZ. In all the regions but SSC5 and SSC2, we left the source as a free parameter. For SSC5 not enough lines were available whilst for SSC2, no χ² minimum was found. For these sources, we then fixed the size to 0.17″, which corresponds to the average best-fit found in the other sources.

C.6. Results

The physical parameters derived from the non-LTE LVG modelling for the regions GMC10, GMC8a, GMC7, SSC2, GMC1a, and GMC1b are displayed in Table C.1 and their temperature-density contour plots are shown in Figure C.1. The abundances derived for these regions are also indicated in Table C.1.

Appendix D: Molecular abundance ratios

The abundances ratios used in Fig. 10 for NGC 253 are shown in Table D.1.

Table D.1.

Abundance ratios derived in each region.

All Tables

Table 1.

Selected regions where the S-bearing species’ emission peaks, their coordinates, and their position within the CMZ.

	Fig. 5. Rotation diagrams of OCS. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.
In the text

	Fig. 6. Rotation diagrams of CCS. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.
In the text

	Fig. 7. Rotation diagrams of SO₂. The parameters N_u, g_u, and E_u are the column density, degeneracy and level energy (with respect to the ground state level) of the upper level, respectively. The error bars on ln(N_u/g_u) include the 15% calibration error. The blue line represent the best fits and the dashed grey lines are the extrapolations of the fit for the full range of E_u.
In the text

An ALCHEMI inspection of sulphur-bearing species towards the central molecular zone of NGC 253

1. Introduction

2. The structure of the central molecular zone of NGC 253

3. Observations, data processing, and line selection

4. Emission distribution and targeted regions

5. Results

5.1. Gaussian fit and spectra

5.1.1. Systemic velocities

5.1.2. FWHMs

5.2. Derivation of physical parameters

5.2.1. LTE analysis: OCS, CCS, and SO2

5.2.2. Methodology of the LVG analysis

5.2.3. Results of the LVG analysis

6. Discussion

6.1. Comparison with previous ALCHEMI studies

6.2. Comparison with Galactic environments

6.3. Assembling the puzzle: Regions traced by S-bearing species

6.3.1. Primary shock tracers

6.3.2. Primary molecular cloud tracers

6.3.3. Shocks, proto-SSCs or molecular clouds: Unconstrained origin

6.4. Other considerations from the general morphology of the CMZ

7. Conclusions

Acknowledgments

References

Appendix A: Velocity-integrated maps

Appendix B: Targeted species and transitions

Appendix C: LVG modelling and results

C.1. CS

C.2. H2S + H 2 34 S

C.3. SO+34SO

C.4. H2CS

C.5. OCS

C.6. Results

Appendix D: Molecular abundance ratios

All Tables

All Figures

5.2.1. LTE analysis: OCS, CCS, and SO₂

C.2. H₂S + $H_{2}^{34} S$ $\rm{H}_2^{34}\rm{S}$

C.3. SO+³⁴SO

C.4. H₂CS