¹¹institutetext: Departamento de Astronomia, Instituto de Astronomia, Geofísica e Ciências Atmosféricas da USP, Cidade Universitária, 05508-090 São Paulo, SP, Brazil, [email protected] ²²institutetext: Astronomical Observatory of the Jagiellonian University, Orla 171, 30-244 Kraków, Poland ³³institutetext: National Center for Nuclear Research, ul. Pasteura 7, 02-093 Warsaw, Poland ⁴⁴institutetext: Scuola Internazionale Superiore di Studi Avanzati (SISSA), Via Bonomea 265, IT-34136, Trieste, Italy ⁵⁵institutetext: Institute for Fundamental Physics of the Universe (IFPU), Via Beirut 2, IT-34151, Trieste, Italy ⁶⁶institutetext: INAF - Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti 93/3, I-40127 Bologna, Italy, ⁷⁷institutetext: ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy ⁸⁸institutetext: Universidade Federal de Santa Catarina, Campus Universitário Reitor João David Ferreira Lima, 88040-900, Florianópolis, SC, Brazil ⁹⁹institutetext: European Southern Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile ¹⁰¹⁰institutetext: Joint ALMA Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile ¹¹¹¹institutetext: Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, E–44001 Teruel, Spain ¹²¹²institutetext: Centro de Astrobiología (CAB), CSIC-INTA, Carretera de Ajalvir km 4, Torrejón de Ardoz, 28850, Madrid, Spain ¹³¹³institutetext: National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA ¹⁴¹⁴institutetext: Max-Planck-Institut für Radioastronomie, Auf-dem-Hügel 69, 53121 Bonn, Germany ¹⁵¹⁵institutetext: Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 830011 Urumqi, China ¹⁶¹⁶institutetext: National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan ¹⁷¹⁷institutetext: Institute of Astronomy and Astrophysics, Academia Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4, Roosevelt Rd, Taipei 106319, Taiwan ¹⁸¹⁸institutetext: Department of Astronomy, School of Science, The Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo, 181-1855 Japan ¹⁹¹⁹institutetext: Institute of Astrophysics, Facultad de Ciencias Exactas, Universidad Andrés Bello, Sede Concepción, Talcahuano, Chile ²⁰²⁰institutetext: New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM, 87801, USA ²¹²¹institutetext: International Gemini Observatory/NSF NOIRLab, Casilla 603, La Serena, Chile ²²²²institutetext: Instituto de Astrofísica de La Plata, CONICET-UNLP, Paseo del Bosque s/n, B1900FWA, Argentina ²³²³institutetext: Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden ²⁴²⁴institutetext: Department of Physics, Faculty of Science and Technology, Keio University, 3-14-1 Hiyoshi, Yokohama, Kanagawa 223–8522 Japan ²⁵²⁵institutetext: Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan ²⁶²⁶institutetext: NOAO, 950 North Cherry Ave. Tucson, AZ 85719, United States ²⁷²⁷institutetext: GMTO Corporation, N. Halstead Street 465, Suite 250, Pasadena, CA 91107, United States

Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone

I. Studying the star formation in extragalactic giant molecular clouds

Pedro K. Humire 0000-0003-3537-4849 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Subhrata Dey 0000-0002-4679-0525 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Tommaso Ronconi 0000-0002-3515-6801 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Victor H. Sasse 0009-0008-8763-7050 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Roberto Cid Fernandes 0000-0001-9672-0296 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Sergio Martín 0000-0001-9281-2919 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Darko Donevski 0000-0001-5341-2162 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Katarzyna Małek 0000-0003-3080-9778 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Juan A. Fernández-Ontiveros 0000-0001-9490-899X Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Yiqing Song 0000-0002-3139-3041 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Mahmoud Hamed 0000-0001-9626-9642 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Jeffrey G. Mangum 0000-0003-1183-9293 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Christian Henkel 0000-0002-7495-4005 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Víctor M. Rivilla 0000-0002-2887-5859 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Laura Colzi 0000-0001-8064-6394 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone N. Harada 0000-0002-6824-6627 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Ricardo Demarco 0000-0003-3921-2177 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Arti Goyal 0000-0002-2224-6664 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone David S. Meier 0000-0001-9436-9471 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Swayamtrupta Panda 0000-0002-5854-7426 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneGemini Science FellowGemini Science Fellow Ângela C. Krabbe 0000-0003-4630-1311 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Yaoting Yan 0000-0001-5574-0549 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Amanda R. Lopes 0000-0002-6164-5051 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone K. Sakamoto 0000-0001-5187-2288 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone S. Muller 0000-0002-9931-1313 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone K. Tanaka 0000-0001-8153-1986 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Y. Yoshimura 0000-0003-0563-067X Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone K. Nakanishi 0000-0002-6939-0372 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Antonio Kanaan 0000-0002-2484-7551 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Tiago Ribeiro 0000-0002-0138-1365 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone William Schoenell 0000-0002-4064-7234 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone Claudia Mendes de Oliveira 0000-0002-5267-9065 Spatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular ZoneSpatially-resolved spectro-photometric SED Modeling of NGC 253’s Central Molecular Zone

(Received December 20, 2024; accepted december 25, 2024)

Abstract

Context. Studying the interstellar medium in nearby starbursts is essential for gaining insights into the physical mechanisms driving these extreme objects, thought to be analogs of young, primeval, star-forming galaxies. This task is now feasible due to deep spectro-photometric data enabled by rapid advancements in ground- and space-based facilities. To fully leverage this wealth of information, extracting insights from the spectral line properties, and the spectral energy distribution (SED) is imperative.

Aims. This study aims to produce and analyze the physical properties of the first spatially-resolved multi-wavelength SED of an extragalactic source that covers six decades in frequency (from near-ultraviolet to centimeter wavelengths) at an angular resolution of 3^′′ which corresponds to linear scale of $\sim$ 51 pc at the distance of NGC 253. We focus on the central molecular zone (CMZ) of this starburst galaxy, which contains giant molecular clouds (GMCs) responsible for half of the galaxy’s star formation.

Methods. We retrieve archival data from optical to centimeter wavelengths, covering six decades of spectral range. We compute SEDs to fit the observations using the GalaPy code and confronting results with the CIGALE code for validation. We also employ starlight library to analyze the stellar optical spectra of the GMCs.

Results. Our results reveal significant differences between central and external GMCs in terms of stellar and dust masses, star formation rates (SFRs), and bolometric luminosities, to provide a few. We have obtained tight relations between monochromatic stellar tracers and star-forming conditions obtained from panchromatic emission. We find that the best SFR tracers are radio continuum bands at 33 GHz, radio recombination lines (RRLs), and the total infrared (IR) luminosity range (L_IR; 8–1000 $\mu$ m) as well as the IR emission at 60 $\mu$ m. The emission line diagnostics based on the BPT and WHAN diagrams suggest that the nuclear region of NGC 253 exhibits shock signatures, placing it in the composite zone typically associated with AGN/star-forming region hybrids, while the AGN fraction from panchromatic emission is negligible ( $\leq$ 7.5%).

Conclusions. Our findings demonstrate, using an extensive dataset and through robust methods, the significant heterogeneity within the CMZ of NGC 253, with central GMCs exhibiting high densities, elevated SFRs, and greater dust masses compared to their external counterparts. We confirm the effectiveness of certain centimeter photometric bands as a reliable method to estimate the global SFR, in accordance with previous studies but this time at GMC scales.

Key Words.:

galaxies: starburst – galaxies: individual: NGC 253 – galaxies: star formation – galaxies: ISM

1 Introduction

The panchromatic spectral energy distribution (SED) of an astrophysical system encapsulates the intricate interplay among its various components, including stars in different evolutionary phases and their remnants; molecular, atomic, and ionized gas; dust, and compact objects such as supermassive black holes. It contains the imprints of baryonic processes that drive its formation and evolution across the cosmic times. Comparison of SEDs at different emission wavebands provides key insights into the origin and nature of emission and the factors that set their energy balance. Therefore, modeling a SED is one of the most effective tools, if not the most effective, for estimating the specific star formation rate (sSFR), from which the star formation rate (SFR) and the stellar mass ( $M_{\rm{\star}}$ ) can be derived. This method relies on reconstructing the stellar emission through star formation history (SFH) models, which can be either purely theoretical or non-parametric (the latter refers to models that do not assume a specific functional form and rely on templates to infer the SFH; see, e.g., Cid Fernandes et al. 2005; da Cunha et al. 2008; Conroy 2013a). SED fitting allows for a detailed analysis of the light emitted across different wavelengths, providing insights into the age, metallicity, and dust content of the stellar population (e.g., Bruzual & Charlot, 2003).

While SEDs have been extensively applied to high-redshift galaxies, there are inherent resolution limitations, particularly at far-infrared (FIR) wavelengths. For example, the Herschel SPIRE beam has a resolution of 36 $\arcsec$ at 500 $\mu$ m, corresponding to roughly 1 kpc at $z=0.05$ , meaning that the observed SEDs often conflate contributions from recent star formation, AGN, and/or stellar outflows, and sub-thermally excited gas, to mention a few. Considering also that a photometric aperture of 2.5 times the beam, assuming a Gaussian profile for the emission source, ensures the capture of extended emission, the problem becomes more critical. Given these constraints, to gain a deeper understanding of the relevant processes affecting galaxies, such as epochs with galaxy-galaxy interactions or times when their star formation began to decline (e.g., Hopkins & Beacom, 2006; Hopkins et al., 2010), it is imperative to focus on local analogs. In this regard, the extensively studied nearby starburst galaxy NGC 253 serves as an ideal target due to its proximity and intense nuclear star formation, which mirrors conditions seen in distant galaxies (Leroy et al., 2015). Focusing on spatially-resolved objects, such as giant molecular clouds (GMCs), allows us to take a significant step forward in producing SED fittings with unprecedented detail, overcoming the aforementioned limitations.

NGC 253 is a nearby starburst galaxy about 3.5 Mpc away with a steep inclination ranging from approximately 70^∘ to 79^∘ (Rekola et al., 2005; Pence, 1980; Iodice et al., 2014). Its inner kiloparsec (kpc) region is a hub of vigorous star formation, generating a rate of roughly 1.7 solar masses per year, which constitutes about half of the galaxy’s total star formation activity (Bendo et al., 2015; Leroy et al., 2015). This active star-forming core encompasses a central molecular zone (CMZ), spanning around 300 parsecs in length and 100 parsec in width projection (Sakamoto et al., 2011) and hosting ten giant molecular clouds (GMCs) discerned through molecular and continuum emissions at high resolution (Leroy et al., 2015). Notably, studies indicate that an active galactic nucleus (AGN) has minimal impact on the molecular gas within NGC 253’s CMZ (Fernández-Ontiveros et al., 2009; Müller-Sánchez et al., 2010; Pérez-Beaupuits et al., 2018), making it an excellent test-bed to study pure star-forming conditions.

The CMZ of NGC 253, along with its GMCs, has been the focus of the ALMA Comprehensive High-resolution Extragalactic Molecular Inventory (ALCHEMI) program (Martín et al., 2021). This project delves into chemical changes caused by shocks or cosmic ray flux density variations, some of them also aiming for a direct comparison with galaxy conditions, the discovery of new molecules, and the characterization of the GMCs that highlights differences between the very central part of the CMZ and its outskirts (Holdship et al., 2021; Harada et al., 2021; Martín et al., 2021; Behrens et al., 2022; Holdship et al., 2022; Humire et al., 2022; Harada et al., 2022; Butterworth et al., 2024; Huang et al., 2023; Behrens et al., 2024; Bouvier et al., 2024; Tanaka et al., 2024; Harada et al., 2024).

Previous studies of the NGC 253 spectral energy distribution (SED) have focused on specific wavelength regimes like optical/IR (Fernández-Ontiveros et al., 2009), mid-IR to (sub-)millimeter (Pérez-Beaupuits et al., 2018), far-UV to far-IR (Beck et al., 2022), the Rayleigh-Jeans tail (Martín et al., 2021), radio (1 GHz) to sub-mm (350 GHz) (Belfiori et al., submitted) and low-frequency radio (MHz to 11 GHz) regimes (Kapińska et al., 2017). However, there remains significant potential in extending these investigations across a broader wavelength regime. This extension can provide valuable insights into various physical processes, including opacity levels, the (s)SFR, and the mass-to-dust rates, to provide a few.

In the present article, we aim to create the first extragalactic panchromatic SED (near-ultraviolet to centimeter wavelengths) performed at a linear resolution of 51 pc or 3″, at the order of magnitude of a giant molecular cloud. For the specific case of this article, given the massive amount of information extracted, we will focus on the stellar processes derived per GMC, such as the SFRs and SFHs.

Along with the acquisition of our new results, we also highlight the importance of testing tracers or proxies that focus on a single wavelength range against the values provided by the full SED to understand the scope of the former in terms of luminosity, ages, object classifications, and calculated SFRs. This way, one can gauge their effectiveness on a larger spectro-photometric wavelength scale.

In the following sections, we will present the collected observations (Sect. 2), our SED and spectroscopic modeling (Sect. 3), our results with small discussions inside (Sect. 4), and a general discussion on common limitations of our methods and the placement of our findings in a high-z perspective (Sect. 5), before providing our conclusions in Sect. 6.

2 Observations

To ensure consistency in the physical scales analyzed, we extract fluxes using a uniform aperture size with a 3 $\arcsec$ diameter, despite the varying angular resolutions across the sample. The coordinates of all GMCs are listed in Table 1. A summary of the dataset used in this study is provided in Table 5. An over-view of the full information used in this work is presented in Fig. 1. Below, we provide a detailed description of the observations conducted with each facility.

2.1 S-PLUS

The Southern Photometric Local Universe Survey (S-PLUS; Mendes de Oliveira et al. 2019) will cover $\sim$ 9,300 deg² of the southern hemisphere (currently 80% complete), employing 7 narrow-band and 5 broad-band filters, encompassing a wavelength range from 3,533 Å ( $u$ filter) to 8,941 Å ( $z$ filter). The angular resolution is subject to weather conditions, with a typical seeing range of 0 $\aas@@fstack{\prime\prime}$ 8 to 2 $\aas@@fstack{\prime\prime}$ 0 and an average of 1 $\aas@@fstack{\prime\prime}$ 5.

For this work, we have used images in all 12 available filters from the S-PLUS fifth data release (available for the internal collaboration), maintaining a consistent angular resolution of 1 $\aas@@fstack{\prime\prime}$ 5 ( $\sim$ 25.4 pc). For further details on the S-PLUS survey, we refer the reader to Mendes de Oliveira et al. (2019) and Herpich et al. (2024).

Table 1: Giant Molecular Cloud (GMC) positions.

Position	$\alpha_{\rm ICRS}$	$\delta_{\rm ICRS}$	$v_{\rm LSR}$
	[00^h:47^m:–^s]	[-25^∘:17’:–”]	[km s^-1]
GMC 1	31.93	29.00	304
GMC 2	32.36	18.80	330
GMC 3	32.81	21.55	286
GMC 4	32.97	19.97	252
GMC 5	33.16	17.30	231
GMC 6	33.33	15.70	180
GMC 7	33.65	13.10	174
GMC 8	33.94	10.90	205
GMC 9	34.14	12.00	201
GMC 10	34.24	7.84	144

¹¹1GMC positions taken from Harada et al. (2024) except for GMCs 3, 4, and 10, that were obtained using CPROPS (Rosolowsky & Leroy, 2006) on ALCHEMI data at 1

\aas@@fstack{\prime\prime}

6 angular resolution. The Local Standard of Rest (LSR) velocities were obtained by fitting a local thermodynamic equilibrium model (LTE) to the ALCHEMI spectra and averaging between the A– and E–CH₃OH symmetry species (Humire et al., 2022). The ALCHEMI spectral resolution of

\sim

8–9 km s^-1 (see Subsect. 2.4) dominates the velocity uncertainty.

Refer to caption — Figure 1: Work overview. From left to right and top to bottom, the figures show: The S-PLUS footprint scan over the southern hemisphere, where each light blue square represents a single S-PLUS observation (FoV of $1.4\times 1.4\,\text{deg}^{2}$ ). A zoom-in into the S-PLUS field where NGC 253 is located is called field $-$ s20s09, where an inset provides a closer view into the nuclear regions of this galaxy. Inside the latter inset, there is another inset zooming in towards NGC 253’s central molecular zone that is shown in the second to last upper panel, where we indicate the H $\alpha$ (MUSE) emission as red contours and shaded areas, showing the strong stellar outflow observed at optical wavelengths. Complementing the latter, in the same previous to last upper panel, we show mainly in blue and white the integrated intensity of the CN 1–0 line at $\sim$ 113.4 GHz taken from ALCHEMI observations at 1 $\aas@@fstack{\prime\prime}$ 6 angular resolution. The leftmost bottom symbol of the beam in the same figure indicates this. The CN emission also presents features of molecular outflows, the most notorious one presumably coming from GMC 3, which is labeled –along with the other nine studied in this work– in the white framed numbers; GMC’s 3 $\arcsec$ apertures are indicated by circles and connected to the numbers by black arrows. The upper rightmost panel shows the SED of GMC 5 from GalaPy (see our Fig. 3 for more details), the one with the highest signal-to-noise ratio of the sample. Additionally, two zoomed-in insets are provided in the middle panels: the optical spectrum (obtained with MUSE) is shown in the middle-left panel, and the sub-millimeter spectrum (from ALCHEMI) is presented in the middle-right panel, while the magenta dots on top of both spectra correspond to the photometric observations listed in Table 5. Key emission lines are labeled in both spectra for reference. Finally, the bottom panel, and its legend in the rightmost position of the figure, summarizes the transmission curves of all the observations used in this work to perform the SED fits at 3 $\arcsec$ apertures and also to test their validity when accounting for bigger apertures, and FIR (Herschel) observations can be included (see Sect. 5.1).

2.2 VLT & HST

The study incorporated Adaptive Optics images obtained from both the Very Large Telescope/NaCo (VLT/NaCo) and VISIR (VLT/VISIR), in addition to Hubble Space Telescope (HST) images sourced from Fernández-Ontiveros et al. (2009), and we refer the reader to that work for further details.

The VLT/NaCo observations were conducted on December 2 and 4, 2005, utilizing the IR wavefront sensor (dichroic N90C10) for atmospheric correction. The observation included the following filters and integration times: J (600 s; 3% Strehl ratio²²2which is a measure of the quality of adaptive optics correction, indicating the ratio of the peak intensity of the observed point spread function to the ideal diffraction-limited one (SR)), Ks (500 s; 20% SR), L (4.375 s; 40% SR), and NB $\_$ 4.05 (8.75 s, centered on Br $\alpha$ ; 40% SR). The resulting images achieved different Full Width at Half Maximum (FWHM) resolutions: 0 $\aas@@fstack{\prime\prime}$ 29 (J), 0 $\aas@@fstack{\prime\prime}$ 24 (Ks), 0 $\aas@@fstack{\prime\prime}$ 13 (L), and 0 $\aas@@fstack{\prime\prime}$ 14 (NB $\_$ 4.05). Images from VLT/VISIR, captured on December 1, 2004, and October 9, 2005, encompassed the N (PAH2 $\_$ 2, $\lambda$ 11.88 $\mu$ m, $\delta\lambda$ 0.37 $\mu$ m; 2826 s) and Q (Q2, $\lambda$ 18.72 $\mu$ m, $\delta\lambda$ 0.88 $\mu$ m; 6237 s) bands. The achieved FWHM resolutions were 0 $\aas@@fstack{\prime\prime}$ 4 (N) and 0 $\aas@@fstack{\prime\prime}$ 74 (Q). Furthermore, HST images were utilized within the central $\sim$ 300 pc region of NGC 253, employing the F555W (V band), F656N (H $\alpha$ ), F675W (R band), F814W (I band), and F850LP filters.

Table 2: Extracted fluxes for sampled GMCs

Source	$\lambda$ [Å]	FWHM [ $\arcsec$ ]	$\theta_{\rm{LAS}}$ [ $\arcsec$ ]	GMC 1	GMC 2	GMC 3	GMC 4	GMC 5	GMC 6	GMC 7	GMC 8	GMC 9	GMC 10
SPLUS/uJAVA	3563	$1.5$	–	0.03	0.03	0.03	0.04	0.03	0.04	0.05	0.06	0.05	0.04
SPLUS/J0378	3770	$1.5$	–	0.04	0.04	0.04	0.04	0.05	0.04	0.08	0.09	0.08	0.06
SPLUS/J0395	3940	$1.5$	–	0.05	0.05	0.05	0.07	0.06	0.07	0.11	0.10	0.11	0.08
SPLUS/J0410	4094	$1.5$	–	0.07	0.07	0.08	0.10	0.10	0.12	0.16	0.19	0.18	0.13
SPLUS/J0430	4292	$1.5$	–	0.09	0.09	0.10	0.13	0.12	0.14	0.20	0.24	0.20	0.18
SPLUS/gSDSS	4751	$1.5$	–	0.13	0.15	0.24	0.31	0.26	0.30	0.40	0.44	0.38	0.31
SPLUS/J0515	5133	$1.5$	–	0.15	0.17	0.29	0.37	0.31	0.35	0.49	0.54	0.45	0.37
HST/F555W	5443	$0.1$	–	–	–	–	0.74	0.63	0.70	0.89	0.79	–	–
SPLUS/rSDSS	6258	$1.5$	–	0.35	0.42	1.07	1.56	1.22	1.24	1.43	1.37	1.10	0.87
SPLUS/J0660	6614	$1.5$	–	0.43	0.56	1.89	3.27	2.39	2.16	2.35	1.94	1.49	1.13
HST/F675W	6718	$0.1$	–	–	–	–	2.46	1.96	1.86	2.05	1.69	–	–
SPLUS/iSDSS	7690	$1.5$	–	0.65	0.84	2.56	3.94	3.08	2.90	2.80	2.53	1.91	1.53
HST/F814W	7996	$0.1$	–	–	–	–	6.28	4.81	4.15	3.78	2.98	–	–
SPLUS/J0861	8611	$1.5$	–	0.88	1.24	3.94	6.46	5.20	4.71	4.08	3.40	2.55	2.01
SPLUS/zSDSS	8831	$1.5$	–	1.07	1.52	5.11	9.07	7.08	6.07	4.97	4.06	2.93	2.35
HST/F850LP	9114	$0.13$	–	–	–	–	12.53	8.93	6.95	5.48	3.87	–	–
VLT/J	12650	$0.29$	–	–	–	22.01	45.75	31.77	20.08	9.88	–	–	–
VLT/K	21800	$0.24$	–	–	–	31.16	93.29	71.15	39.01	–	–	–	–
Spitzer/IRAC 3.6 $\mu$ m	36000	$1.66$	–	2.71	5.72	28.40	106.38	109.02	62.01	23.41	12.60	6.27	6.30
Spitzer/IRAC 4.5 $\mu$ m	45000	$1.72$	–	2.43	4.78	28.08	150.28	150.45	69.47	21.64	10.77	5.20	5.79
Spitzer/IRAC 5.8 $\mu$ m	58000	$1.88$	–	7.23	17.75	121.03	357.56	464.85	250.73	89.61	45.28	22.27	22.98
Spitzer/IRAC 8.0 $\mu$ m	80000	$1.98$	–	18.98	52.71	335.32	1016.10	1095.98	644.08	243.56	120.51	60.60	61.62
VLT/N	118800	$0.4$	–	–	–	287.88	3207.94	1273.06	477.52	121.49	–	–	–
VLT/Q	187200	$0.74$	–	–	–	1178.25	18901.04	9984.58	3340.52	654.91	355.37	–	–
ALMA/B7	9252854	$1.6$	15	11.48	9.63	188.92	289.28	343.56	255.51	51.86	18.37	18.35	17.51
ALMA/B6	12337138	$1.6$	15	2.58	3.28	68.74	122.77	165.38	103.49	17.11	5.82	5.17	6.45
ALMA/B5	16031682	$1.6$	15	2.85	1.91	35.92	72.02	110.39	58.05	8.34	3.07	2.86	2.64
ALMA/B4	20818921	$1.6$	15	1.29	1.57	24.13	56.96	101.44	40.57	5.76	2.37	1.69	2.50
ALMA/B3	29979246	$1.6$	15	–	–	18.41	50.46	102.15	33.61	3.60	–	–	1.27
EVLA Band Ka (DnC)	91330642	1.76 $\times$ 1.39	44	–	–	21.80	59.07	141.84	40.53	4.91	–	–	1.54
VLA Band K (B)	127245834	0.47 $\times$ 0.26	7.9	–	–	–	–	135.60	–	–	–	–	–
VLA Band Ku (B)	200665639	0.67 $\times$ 0.40	12	–	–	30.67	71.92	181.58	50.21	–	–	–	–
EVLA Band X (A)	299773911	0.46 $\times$ 0.19	5.3	–	–	–	98.03	264.56	81.33	–	–	–	–
VLA Band C (A)	616844217	0.69 $\times$ 0.39	8.9	–	–	41.34	104.45	318.02	70.16	9.08	–	–	–
VLA Band L (A)	2012164964	2.13 $\times$ 1.25	36	6.57	12.38	89.09	175.32	274.92	160.83	36.35	16.68	15.20	5.69

⁵⁵5 The fifth to last columns display the extracted fluxes [mJy] for each sampled Giant Molecular Cloud (GMC) within a 3

\arcsec

aperture diameter. A dash – indicates the lack of (enough) data filling the aperture. Uncertainties were assumed to be of the order of 10% for all the sample except for ALMA data, where we assume 15%. Additionally, observational parameters, namely, wavelength, FWHM of the beam/psf, and maximum recoverable scales in the case of interferometric observations are indicated in columns 2, 3, and 4, respectively, while the first column specifies the source (facility or survey) and filter, plus the array configuration mode provided parenthetically for the (E)VLA observations. The maximum recoverable scale, also known as the largest angular scale,

\theta_{\rm{LAS}}

, is pertinent for interferometric observations only, namely: ALMA (Martín et al., 2021), EVLA³³3https://science.nrao.edu/facilities/vla/docs/manuals/oss2012B/performance/referencemanual-all-pages, and VLA⁴⁴4https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/resolution.

2.3 Spitzer/IRAC

We retrieved Spitzer/IRAC mosaic images from the Spitzer Heritage Archive⁶⁶6https://irsa.ipac.caltech.edu/data/SPITZER/Enhanced/SEIP/. The mean FWHM of the point spread function (PSF) of the four IRAC bands at 3.6, 4.5, 5.8, and 8.0 $\mathrm{\mu m}$ are 1 $\aas@@fstack{\prime\prime}$ 66, 1 $\aas@@fstack{\prime\prime}$ 72, 1 $\aas@@fstack{\prime\prime}$ 88, and 1 $\aas@@fstack{\prime\prime}$ 98, respectively (Fazio et al., 2004, their Table 3). Thanks to the large field of view (FoV) of Spitzer, which covers the entire emission from NGC 253, we were able to use data from these four IRAC filter images in all our GMCs.

2.4 ALMA

ALMA observations of NGC 253 were conducted during Cycles 5 and 6 as part of the ALCHEMI large program. The observations consisted on an unbiased complete spectral scan across the ALMA bands 3 to 7 (84–373 GHz, $\lambda$ =3.6–0.8 mm) across the whole CMZ region (inner 500 pc; Sakamoto et al. 2006) of the galaxy. The observations were carried out using both 12m and 7m antenna arrays aiming to retrieve extended emission, and the data were imaged to common angular and spectral resolutions of 1 $\aas@@fstack{\prime\prime}$ 6 and $\sim$ 8–9 km s^-1, respectively. The flux density RMS noise ranges from 0.18 to 5.0 mJy beam^-1 in the 8–9 km s^-1 channels. In Fig. 2 we provide integrated intensity maps of the continuum emission (in mJy) of the CMZ of NGC 253 at ALMA Bands 3, 4, 6, and 7 (100, 144, 243, and 324 GHz, respectively).

The absolute flux density uncertainty was reported to be of 10 to 15%. We conservatively assume an uncertainty of 15% for our extracted ALCHEMI fluxes. However, the relative flux uncertainty within the individual ALCHEMI observations is $2-3\%$ after the alignment performed on the data enabled by the continuous frequency coverage. A full description of the dataset, calibration, and imaging is provided by Martín et al. (2021).

2.5 EVLA

We obtained Expanded Karl G. Jansky Very Large Array (EVLA) archival data in two different ways. One was by directly searching into the VLA archive. In this way, we obtained the reduced NGC 253’s CMZ image in the X Band at 3 cm (PI: N. Harada), which has been observed in A configuration mode, leading to a maximum recoverable scale (MRS) of 5 $\aas@@fstack{\prime\prime}$ 3 according to the EVLA documentation⁷⁷7https://science.nrao.edu/facilities/vla/docs/manuals/oss2012B/performance. We applied a cleaning to this data using tclean with a standard gridder and a hogbom deconvolver inside the Common Astronomy Software Applications (CASA) package⁸⁸8https://casa.nrao.edu/index.shtml. Additionally, we used 1 cm (33 GHz; Ka-Band) continuum image presented in Kepley et al. (2011). These observations were done at DnC hybrid configuration, which leads to an MRS of 44 $\arcsec$ ⁹⁹9according to the smaller array extension of the hybrid configuration, configuration C for this case. We used CASA imsmooth to convolve the X Band EVLA continuum image to a common circular beam of 1 $\aas@@fstack{\prime\prime}$ 6, matching that of ALCHEMI. In the same way, for the Ka-Band (26.5–40.0 GHz) image, we choose a slightly larger (1 $\aas@@fstack{\prime\prime}$ 8) beam due to an original beam of 1 $\aas@@fstack{\prime\prime}$ 76 $\times$ 1 $\aas@@fstack{\prime\prime}$ 38. Since the images are in Jy/beam, we converted them into Janskys by dividing their flux by the beam area in pixels [pix/beam] available in the CASA Viewer (Statistics/BeamArea) panel. We note that the largest angular scale, $\theta_{\rm{LAS}}$ (see Table 5), is smaller than the aperture size for the Ka-Band observations. However, we did not see significant deviations of these observations with respect to the rest of the dataset. In fact, the model adjusts well to the observations when they are available (GMCs 3–7), and any expected departure or uncertainty in flux is compensated for by the added uncertainties in our SED models, which is up to 30% for the case of GalaPy (Ronconi et al., 2024).

2.6 VLA

We retrieve Karl G. Jansky Very Large Array (VLA) continuum maps processed with the AIPS VLA pipeline by Lorant Sjouwerman (NRAO) from the webpage: https://www.vla.nrao.edu/astro/nvas/. Of those, VLA Band L (VLA A conf.) is the only one were all the 10 GMCs studied in this work are detectable. Details about wavelengths, fluxes per GMC, array configurations, and FWHMs, are summarized in Table 2. We find very good consistency between VLA and EVLA points, delivering similar fluxes within 1 or 2 $\sigma$ of the SED fitting. There is a noticeable missing flux in L-band for GMC 5 (see Fig. 3), which can be associated to the free-free absorption effect (e.g., Kellermann & Pauliny-Toth, 1969) observed in compact IR-bright starbursts (Condon et al., 1991; Clemens et al., 2010; Dey et al., 2022).

Since we are considering flux densities, we decided not to deconvolve our VLA continuum maps, as we noticed that it can produce artifacts that may show up in GMC positions without clear emission, leading to false positives. Upon examining the summed fluxes within 3 $\arcsec$ diameter apertures before and after deconvolution with CASA imsmooth, we observed differences of less than 10% (i.e., less than instrumental uncertainties).

2.7 GMC Positions

The selection of GMC positions to be analyzed in this study is driven by the ALMA data from the ALCHEMI project. Thus, we use the GMC positions determined by Harada et al. (2024), which incorporate both the continuum and molecular emission peaks, except for GMC 3, GMC 4, and GMC 10. For these three clouds, and given the similarity between the dataset used in Leroy et al. (2015) and the ALCHEMI data, A. K. Leroy (private communication) provided modified GMC coordinates tailored to the ALCHEMI data for the collaboration, which were retained for GMCs 3, 4, and 10. These three GMC positions were derived using CPROPS (Rosolowsky & Leroy, 2006) considering intensity peaks of the CS $J=2{-}1$ and H¹³CN $J=1{-}0$ transitions. Our region coordinates are very close to those used by Leroy et al. (2015).

The general reason for this selection is to avoid, as much as possible, overlapping between our apertures of 3 $\arcsec$ , while at the same time extracting most of the intensity in each aperture. We note that the variation in the GMC positions from different source extraction criteria is relatively marginal ( $\sim$ 0 $\aas@@fstack{\prime\prime}$ 2) compared to the apertures used in this paper.

Along the course of this work, we will compare our GMC characteristics with those of the GMCs in previous ALCHEMI works with the same numbering, unless it is explicitly mentioned.

3 Analysis: SED and Spectroscopic modeling

The main results from this study come from two methods used to explore in detail the spectroscopic and photometric information of the CMZ of NGC 253. We have employed GalaPy for the photometry information, from near-UV to centimeter wavelengths. To cross-validate the results obtained from GalAPy, we compared them with those from CIGALE. CIGALE is a robust SED modeling, continuously improved over the years by both its original developers and the broader astronomical community.

Additionally, to complement this study, we applied starlight for the spectral MUSE data plus photometric (S-PLUS) information in the near-UV to near-IR range. Below, we elaborate on our procedures. In Fig. 3 we provide our best-fit models to the extracted flux densities summarized in Table 5. This Figure primarily shows GalaPy results but also presents the CIGALE best-fit models (solid magenta lines). CIGALE best-fitting results are further provided in detail (i.e., with all the considered components) in Fig. 14.

3.1 GalaPy photometric SED fitting

We perform the SED fitting of our selected GMCs with GalaPy (Ronconi et al., 2024), an open-source application programming interface developed in Python/C $++$ tailored for this task across a range from X-ray to radio frequencies. GalaPy assumes a Chabrier initial mass function (IMF; Chabrier, 2003) and a flat $\Lambda$ CDM cosmology with approximate parameters: matter density around 0.3, baryonic density around 0.05, and a Hubble constant $H_{0}=$ 100 $h$ km s^-1 Mpc^-1, where $h$ is roughly 0.7 (Planck Collaboration et al., 2020). GalaPy is expected to include tools for optical spectroscopy and AGN component on top of its SED fitting algorithm. In the subsequent subsections, we detail the relevant configurations pertinent to the objectives of this study.

3.1.1 The In-Situ model

GalaPy allows the user to choose among a range of SFH models, both empirical and non-parametric. For our models, we use the default In-Situ SFH model, which has proven successful in predicting the emission from both late and early type galaxies in the local Universe as well as for young highly star-forming systems up to high redshift (Ronconi et al., 2024). The physically motivated In-Situ model, developed by Lapi et al. (2018), adopts a star formation rate ( $\psi(t)$ ) of:

\psi(t)\propto e^{-x}-e^{-s\gamma x},

(1)

with $x\equiv\tau/s\tau_{\star}$ , where $\tau$ is the galactic age, $\tau_{\star}$ is the characteristic star-formation timescale, $s\approx 3$ is related to the gas condensation, and $\gamma$ is related to the gas dilution, recycling, and stellar feedback (see Lapi et al., 2020, for more details).

In the in situ model, the evolution of the gas and dust masses and of the gas and stellar metallicities can be followed analytically as a function of the galactic age, and self-consistently with respect to the evolution of the SFR. This means that an isolated system is assumed for each molecular cloud (MC), implying an interdependence among parameters. For example, a correspondence between metallicities and ages. This approach guarantees that the derived physical properties that directly result from the evolution of the stellar population (e.g. the components’ masses and metallicities), are consistently propagated to the other models contributing to the overall emission of the object.

3.1.2 Simple Stellar Populations

Among the available Simple Stellar Population (SSP) libraries in GalaPy, we have selected the “refined” version of $parsec22.NTL$ (Bressan et al., 2012; Yan et al., 2019; Ronconi et al., 2024) as our preferred choice. This library ensures a dense wavelength grid with a minimum of 128 values per order of magnitude within its 2189-point grid. It incorporates emission from dusty Asymptotic Giant Branch (AGB) stars and accounts for nebular emission and free-free continuum in addition to the stellar continuum and non-thermal synchrotron radiation from core-collapse supernovae.

For the purposes of our study, the $parsec22.NTL$ SSP library was modified to correctly process the photometric points obtained from the narrow bands of the VLT/VISIR instrument (see Sect. 2.2). This adjustment will be available in the next release of GalaPy.

3.1.3 Dust model

GalaPy implements an age-dependent two component dust model which is constructed in order not to assume an attenuation curve but to derive it instead from structural parameters (e.g. density and extension) which can be obtained by tuning the model to best represent an observed dataset.

The two different components comprise the contribution of (1) a hot molecular cloud phase (normally referred in the literature as a “birth-molecular cloud”; see, e.g., Charlot & Fall 2000) of new stellar populations and (2) a diffuse and extended dusty medium which further attenuates the stellar emission. Both components also contribute to the IR emission which is then the combination of two separate modified grey-bodies.

Among the several possible free parameters on which this dust model depends, we chose to fix some known values and ranges, in order to both ensure physically plausible results and speed up the computing process. This consists in setting a range of reliable physical terms, such as the total number of molecular clouds to $N_{\text{MC}}=1$ , as in this work we are assuming isolated GMCs, and their dusty and molecular radii, which were set to be within a 10–100 pc range, although they will likely be closer to the lower limit considering the 30 pc diameter inferred from their sub-mm emission (Leroy et al., 2015). A list of all common parameter ranges and individual values used to model the SEDs in the ten GMCs is provided in Table 3. There, the redshift was taken from the SIMBAD Astronomical Database (Wenger et al., 2000). These ranges are fixed for all GMCs. On the other hand, the parameters that were fine-tuned for specific GMCs along with their corresponding adjusted ranges are presented in Table 4. The parameters listed in this second table vary across different GMCs, with their respective ranges for each GMC. For example, the range of the molecular cloud radius, R_CM, was fine-tuned and is set between 0.0 and 1.7 for GMCs 2 and 8, while for GMCs 1, 3, 4, 5, 6, 7, 9, and 10, it remains between 0.0 and 2.0 (in log10 scale).

We use different initial conditions for unknown parameter ranges such as the GMC ages (age), the time of maximum star-formation (sfh.tau_star), the time that stars need to escape from its parental molecular cloud (ism.tau_esc), the maximum SFR (sfh.psi_max), and the average radius of the molecular cloud (ism.R_MC). We mainly modified the R_MC, ism.R_MC, and sfh.psi_max in cases where the results were bimodal, in order to achieve a single result in the parameter space distribution. For example, the R_MC was fixed to not exceed 50 pc in GMCs 2 and 8, motivated by Table 3 in Leroy et al. (2015), where their molecular radii are expected to be around 52 and 44 pc, respectively. Interestingly, these two GMCs are the largest in the sample, according to the mentioned study.

We recommend the reader to see Table B.3 in Ronconi et al. (2024), also described in its $ReadtheDocs$ web-page¹⁰¹⁰10https://galapy.readthedocs.io/en/latest/general/free_parameters.html, to find a complete description of all GalaPy parameters, beyond the ones described in the present work, which assumes and is restricted to the In-Situ model.

Parameter	Value/Range	Brief description
redshift	0.000864
sfh.tau_quench	1e+20	Star formation quenching (yrs.)
ism.f_MC	([0.0, 1.0], lin)	MCs’ fraction into the ISM
ism.norm_MC	100.0	MCs’ normalization factor
ism.N_MC	1.0	number of MCs
ism.dMClow	1.3	Extinction index $<$ 100 $\mu$
ism.dMCupp	1.6	Extinction index $\gtrsim$ 100 $\mu$
ism.norm_DD	1.0	Diffuse dust norm. factor
ism.Rdust	([0.0, 2.0], log)	Radius of the diffuse dust (DD, pc)
ism.f_PAH	([0.0, 1.0], lin)	DD fraction radiated by PAH
ism.dDDlow	0.7	DD extinction index $<$ 100 $\mu$
ism.dDDupp	2.0	DD extinction index $\gtrsim$ 100 $\mu$
syn.alpha_syn	0.75	Spectral index
syn.nu_self_syn	0.2	Self-absorption frequency (GHz)
f_cal	([-5.0, 0.0], log)	Calibration uncertainty

Table 3: GalaPy input parameter values or ranges common to the ten GMCs studied in this work. The terms “log” and “lin” next to the ranges indicate logarithmic (log10) and linear values, respectively. For a detailed explanation, see Sect. 3.

Parameter GMC(s) with varied parameters Common parameter ranges Brief description age GMC2: ([6.0, 12.0]) ([6.0, 11.0]) Age of the MC sfh.tau_star GMC2: ([6.0, 12.0]) ([6.0, 11.0]) Characteristic timescale ism.R_MC GMC2: ([0.0, 1.7]), GMC8: ([0.0, 1.7]) ([0.0, 2.0]) Average radius of a MC ism.tau_esc GMC2: ([6.0, 12.0]) ([6, 11.0]) Stars’ escape time from MC sfh.psi_max GMC3: ([-2.0, 1.5]), GMC5: ([-2.0, 1.5]) ([-3.0, 1.5]) Maximum SFR, at ism.tau_star

Table 4: Parameter ranges set to GalaPy that very depending on the modeled GMC. All values are given in log10 scale. For a detailed explanation, see Sect. 3.

3.2 CIGALE photometric SED fitting

In addition to the main results obtained from GalaPy, for comparison we also provide simultaneous UV-to-radio SED modeling with the recently updated CIGALE 2022.1 version (Yang et al., 2022). This software operates on the principle of energy balance, where energy attenuation in the UV-to-optical range is re-emitted in the infrared in a self-consistent manner. Its modular and parallelized design, along with its user-friendly interface, has made it a popular choice in the literature for estimating the astrophysical characteristics of various targets. Below, we provide a brief overview of the models employed to simulate the galaxy emission.

To construct the SEDs with CIGALE, we utilized the full photometric coverage available for our GMCs. The modeling was based on the widely-used Bruzual & Charlot (2003) SSPs, adopting the IMF of Chabrier (2003) and assuming a solar-like metallicity. The stellar populations evolved using a delayed SFH with an additional burst component (Ciesla et al., 2015). This SFH approach is more realistic compared to the simple delayed model and has been demonstrated to minimize biases in the estimation of SFR and stellar mass (Schreiber et al., 2018). Additionally, the inclusion of the burst is consistent with the ALMA detections (Hamed et al., 2021, 2023). For the SFH parameters, we maintained flexibility in the age and the e-folding time of the main stellar population to avoid introducing artificial biases in the efficiency of the initial burst.
Dust attenuation was modeled following the Calzetti et al. (2000) attenuation law, allowing the color excess of the young stellar population to vary as a free parameter. Dust emission was a critical component of our SED modeling, given the extensive IR coverage, particularly from ALMA. We used the models provided by Draine et al. (2014). To ensure a good fit, we explored the parameter space for PAH fractions and the radiation field, U_min, aiming to minimize $\chi^{2}$ while keeping the models physically realistic. Considering the rich radio data available for the GMCs, we extended our UV-IR modeling to radio regime taking into account the PL nature of the synchrotron spectrum and the ratio of the FIR/radio (q_IR) correlation in the CIGALE fitting.
CIGALE minimizes the $\chi^{2}$ by measuring the difference between observed data and the model predictions, weighted by the uncertainties in the observations. The goal of minimization is to find the best-fitting model parameters that reduce this difference (Boquien et al., 2019). The CIGALE modeling was performed in the rest frame, which is directly applicable to NGC 253, as no corrections were necessary. Additional details, including the addition of AGN-related parameters, are provided in Appendix A.

3.3 Stellar population fitting with starlight

GalaPy and CIGALE are used in the full SED fits (all wavelengths from UV to Radio), while starlight will be used for the optical data only to fit stellar population models. starlight has been successfully applied to fit the stellar component of a diverse range of galaxy types, such as low-luminosity AGN, Seyfert nuclei, and normal galaxies.

While traditionally used for optical spectra, we applied the latest version of the code, described by Werle et al. (2019), to fit simultaneously the combined MUSE and S-PLUS data (see Fig. 4), thus covering the 3,500–10,000 Å range. The combination of these datasets provides a comprehensive view of the galaxy’s stellar populations, spanning a wide range of wavelengths. In particular, the S-PLUS bands densely cover the highly informative region around the 4,000 Å break which falls outside the range of MUSE spectra that starts at 4,600Å.

starlight uses a chi-squared minimization technique to find the best-fit non-parametric combination of stellar populations of different ages and metallicities from a pre-defined library. This analysis allows us to estimate properties, such as mean age, metallicity, and extinction, as well as the velocity dispersion of the stellar component. The base is built from SSP models from the Charlot & Bruzual 2019 library models (Martínez-Paredes et al., 2023) using PARSEC evolutionary tracks (Chen et al., 2015) and a Chabrier IMF. A base of composite stellar population is built from the SSP models assuming constant SFRs in 16 logarithmically spaced age bins with 5 different metallicities. Extinction is modeled with two components: one applied to all populations ( $A_{\rm{V}}$ ) and an extra one applied exclusively to $\leq 10$ Myr stars ( $\delta A_{V}$ ).

The list of parameters to fit includes the flux (or mass) fractions of the different populations, global and selective extinction, the velocity shift and velocity dispersion.

We have adjusted the weight scheme for the code used to extract residuals for H $\alpha$ , [N II], [O III], and H $\beta$ by augmenting them around those lines. This modification allows a detailed examination on the contribution of stellar absorption features to these emission lines, compared to a general fit to the MUSE+SPLUS spectra/photometry. These lines will subsequently serve as input for diagnostic diagrams and to estimate the extinction affecting the emission line regions.

4 Results

In the following, we supplement our SED analysis with existing literature. Since the GMC selection was based on ALMA observations, our primary information source is the ALCHEMI results (Martín et al., 2021). After detailing the characteristics of the ten GMCs by dividing them in two main groups, we will compare our panchromatic results – derived from SED fitting – on star formation rates (SFRs) and histories (SFHs) with those obtained from pure-optical and pure-sub-mm/cm observations, that is, from monochromatic tracers such as the H $\alpha$ and H40 $\alpha$ emission lines and different infrared and radio continuum regimes or bands. This approach will help us understand the key limitations to consider when working with mono-wavelength datasets. For a more detailed description of our results in light of previous literature information and for each of the ten GMCs, we refer the reader to our Appendix B.

4.1 SED fitting

The wealth of archival data described in Sect. 2 does not always cover the ten GMCs. Only S-PLUS, Spitzer/IRAC, ALCHEMI, VLA L (1.5 GHz), and Ku (15 GHz) Bands fully cover the NGC253’s CMZ (see Table 2). EVLA Band X does not detect emission outside GMCs 3 to 6, similar to the case in VLA K (23.6 GHz) band. Also, the VLA C (4.9 GHz) Band does not detect emission in GMCs 1 and 2, and our aperture of 3″ is not fully covered by emission in GMCs 8–10, hence we discarded the latter contribution. From the above, we can say that the best-studied GMCs are GMCs 3 to 6, corresponding to the core of NGC253 and representing most of the star formation produced there (Bendo et al., 2015). It is worth saying that, in addition to the possible lack of continuum emission, the not imaged GMCs may fall below the detection limit of our ALMA Band 3 observations and the different centimeter bands mentioned above, the detection limit at 3 $\sigma$ of EVLA Band X, and VLA Bands K and C are 2.4 $\times$ 10^-2, 5.4 $\times$ 10^-3, and 6.6 $\times$ 10^-4 Jy, respectively.

In general, the GMC’s SED shapes, which can be seen in Fig. 3, are similar to what is found in local and high-redshift starbursts (e.g., Swinbank et al., 2010). These SEDs exhibit significant absorption in the near-UV and optical wavelength regimes ( $\sim 3\times 10^{3}$ – $10^{4}$ Å), caused by an exceptionally large amount of dust. This leads to a substantial difference between the predicted and observed emissions (i.e., stellar emission with and without extinction), illustrated by the red and green lines in Fig. 3. At the mentioned wavelengths, the flux difference can reach up to four orders of magnitude (e.g., from $\sim 10^{3}$ mJy down to $\sim 10^{-1}$ mJy in GMCs 4–6 at 10⁴ $\AA$ ). The total extinction ( $A_{v}$ ) derived from the stellar component exceeds 5.0 magnitudes in GMCs 3–6, based on analyses from GalaPy, CIGALE, and the Balmer lines (H $\alpha$ and H $\beta$ ; see Sect. 4.2.2).

Tables 5 and 6 respectively contain GalaPy and CIGALE results related to star-formation processes, such as the SFR, the stellar mass ( $M_{\bigstar}$ ), the characteristic timescale ( $\tau_{\bigstar}$ ), or the stellar age (Age). In addition, results on the gas and dust characteristics, the attenuation, the radiation field, the time for stars to escape from their parental MC (ism.tau_esc), and uncertainties of the model, to mention a few, are given in Tables 9 and 10 for GalaPy and CIGALE, respectively.

We summarize below the main characteristics of the molecular clouds studied here by dividing them basically in two different media. The central GMCs, numbered from 3 to 6 and the external ones: GMCs 1,2 and 7–10. We find strong differences among them in terms of stellar and dust masses, and star formation rates. Appendix B provides detailed information of individual GMCs.

GMC SFR $M_{\bigstar}$ $\tau_{\bigstar}$ $Z_{\bigstar}$ $L_{\rm{bol}}$ Age [ $\times 10^{-4}M_{\odot}$ yr^-1] [log₁₀ ( $M_{\odot}$ )] [log₁₀ (yrs)] [ $\times$ 10^-3] [log₁₀ ( $L_{\odot}$ )] [log₁₀ (yrs)] 1 50 ${}_{-10}^{+10}$ 6.535 ${}_{-0.097}^{+0.159}$ 10.88 ${}_{-0.21}^{+0.08}$ 0.35 ${}_{-0.08}^{+0.18}$ 8.065 ${}_{-0.073}^{+0.052}$ 9.30 ${}_{-0.13}^{+0.15}$ 2 80 ${}_{-10}^{+10}$ 7.534 ${}_{-0.086}^{+0.106}$ 9.21 ${}_{-0.14}^{+0.14}$ 3.88 ${}_{-0.15}^{+0.19}$ 8.391 ${}_{-0.062}^{+-0.002}$ 9.66 ${}_{-0.12}^{+0.13}$ 3 870 ${}_{-90}^{+110}$ 8.573 ${}_{-0.100}^{+0.108}$ 9.07 ${}_{-0.46}^{+0.28}$ 7.88 ${}_{-0.94}^{+2.27}$ 9.354 ${}_{-0.042}^{+0.071}$ 9.62 ${}_{-0.23}^{+0.15}$ 4 2740 ${}_{-280}^{+390}$ 8.850 ${}_{-0.136}^{+0.251}$ 7.48 ${}_{-0.10}^{+0.13}$ 20.84 ${}_{-1.47}^{+2.31}$ 10.132 ${}_{-0.114}^{+0.050}$ 8.54 ${}_{-0.08}^{+0.16}$ 5 6450 ${}_{-590}^{+600}$ 8.620 ${}_{-0.105}^{+0.124}$ 7.57 ${}_{-0.11}^{+0.12}$ 17.59 ${}_{-1.14}^{+1.17}$ 10.184 ${}_{-0.065}^{+0.016}$ 8.40 ${}_{-0.03}^{+0.13}$ 6 1870 ${}_{-170}^{+180}$ 8.605 ${}_{-0.095}^{+0.072}$ 8.17 ${}_{-0.12}^{+0.10}$ 12.87 ${}_{-0.63}^{+0.82}$ 9.602 ${}_{-0.045}^{+0.061}$ 8.56 ${}_{-0.04}^{+0.06}$ 7 220 ${}_{-20}^{+30}$ 8.156 ${}_{-0.087}^{+0.096}$ 8.50 ${}_{-0.18}^{+0.31}$ 8.12 ${}_{-1.18}^{+0.87}$ 8.903 ${}_{-0.080}^{+0.041}$ 9.39 ${}_{-0.14}^{+0.23}$ 8 110 ${}_{-10}^{+10}$ 7.837 ${}_{-0.127}^{+0.098}$ 9.24 ${}_{-0.24}^{+0.24}$ 4.45 ${}_{-0.40}^{+0.71}$ 8.533 ${}_{-0.048}^{+0.045}$ 9.76 ${}_{-0.16}^{+0.15}$ 9 90 ${}_{-10}^{+10}$ 7.752 ${}_{-0.099}^{+0.106}$ 9.12 ${}_{-0.18}^{+0.24}$ 4.53 ${}_{-0.53}^{+0.62}$ 8.605 ${}_{-0.212}^{+0.102}$ 9.73 ${}_{-0.10}^{+0.12}$ 10 50 ${}_{-10}^{+10}$ 7.581 ${}_{-0.150}^{+0.095}$ 9.55 ${}_{-0.25}^{+0.21}$ 3.31 ${}_{-0.32}^{+0.37}$ 8.507 ${}_{-0.090}^{+0.075}$ 9.64 ${}_{-0.12}^{+0.12}$

Table 5: GalaPy results for star-formation-related properties. A complete list of results is provided in Table 9.

GMC SFR_INST SFR_100Myr SFR_10Myr $M_{\bigstar}$ $\tau_{\rm{burst}}$ $M_{\bigstar,\text{old}}$ $M_{\bigstar,\text{young}}$ Age [ $\times 10^{-4}\,M_{\odot}\,\text{yr}^{-1}$ ] [ $\times 10^{-4}\,M_{\odot}\,\text{yr}^{-1}$ ] [ $\times 10^{-4}\,M_{\odot}\,\text{yr}^{-1}$ ] [log₁₀( $M_{\odot}$ )] [log₁₀(yrs)] [log₁₀( $M_{\odot}$ )] [log₁₀( $M_{\odot}$ )] [log₁₀(yrs)] 1 20.76 $\pm$ 3.19 23.27 $\pm$ 3.77 20.95 $\pm$ 3.20 6.70 $\pm$ 0.07 8.15 $\pm$ 0.30 6.70 $\pm$ 0.03 7.87 $\pm$ 0.11 9.41 $\pm$ 8.88 2 31.02 $\pm$ 3.67 26.37 $\pm$ 7.05 31.80 $\pm$ 3.62 7.17 $\pm$ 0.06 8.32 $\pm$ 0.28 7.17 $\pm$ 0.00 8.10 $\pm$ 0.12 9.52 $\pm$ 8.63 3 262.01 $\pm$ 68.63 292.05 $\pm$ 53.19 264.88 $\pm$ 70.81 7.82 $\pm$ 0.07 8.10 $\pm$ 0.30 7.82 $\pm$ 0.05 7.74 $\pm$ 0.04 9.36 $\pm$ 8.85 4 642.50 $\pm$ 32.13 446.83 $\pm$ 49.84 655.24 $\pm$ 32.76 8.50 $\pm$ 0.02 8.39 $\pm$ 0.30 8.50 $\pm$ 0.00 8.16 $\pm$ 0.00 9.54 $\pm$ 8.55 5 1956.59 $\pm$ 129.59 204.02 $\pm$ 13.25 2024.70 $\pm$ 128.51 8.17 $\pm$ 0.02 8.26 $\pm$ 0.30 8.17 $\pm$ 0.01 7.92 $\pm$ 0.11 9.44 $\pm$ 8.65 6 476.20 $\pm$ 49.31 509.41 $\pm$ 61.56 481.12 $\pm$ 51.15 7.94 $\pm$ 0.06 8.23 $\pm$ 0.30 7.94 $\pm$ 0.00 7.93 $\pm$ 0.00 9.44 $\pm$ 8.85 7 191.88 $\pm$ 32.85 214.98 $\pm$ 34.37 193.54 $\pm$ 32.60 7.68 $\pm$ 0.06 8.03 $\pm$ 0.30 7.68 $\pm$ 0.01 7.68 $\pm$ 0.09 9.33 $\pm$ 8.85 8 47.72 $\pm$ 11.49 46.49 $\pm$ 9.52 49.32 $\pm$ 11.65 7.51 $\pm$ 0.07 8.02 $\pm$ 0.30 7.51 $\pm$ 0.01 7.75 $\pm$ 0.01 9.36 $\pm$ 8.76 9 35.29 $\pm$ 5.96 41.35 $\pm$ 7.97 35.72 $\pm$ 5.93 7.18 $\pm$ 0.06 8.15 $\pm$ 0.28 7.18 $\pm$ 0.00 7.71 $\pm$ 0.00 9.35 $\pm$ 8.75 10 32.41 $\pm$ 5.03 38.11 $\pm$ 8.36 32.98 $\pm$ 5.10 7.17 $\pm$ 0.07 8.22 $\pm$ 0.30 7.17 $\pm$ 0.00 7.82 $\pm$ 0.00 9.40 $\pm$ 8.86

Table 6: CIGALE results for star-formation-related properties, including 1

\sigma

uncertainties. Additional results are listed in Tables 11, 12, and 13.

4.1.1 GMC characteristics

The Central Molecular Zone (CMZ) of NGC 253 can mainly be divided into two primary groups based on the physical and chemical properties found in their giant molecular clouds (GMCs) listed in Table 1: internal and external GMCs. The internal GMCs, corresponding to GMCs 3, 4, 5, and 6, are located in the very nuclear region of the CMZ, about 120 pc extension (see, e.g., the scale bar in Fig. 13), and are characterized by high densities, elevated temperatures, fast shocks likely triggered by a tremendous star formation activity (Bouvier et al., 2024; Bao et al., 2024). Conversely, the external GMCs encompass GMCs 1, 2, 7, 8, 9, and 10, located near (GMC 7) or fully immersed in the periphery of the internal bar. There, orbital intersections such as inner Lindblad resonances (Iodice et al., 2014), bar/spiral arm interactions –also known as $x_{1}/x_{2}$ interactions (Kim et al., 2012)–, and cloud-cloud collisions (Ellingsen et al., 2017, see their Sect. 4.3) take place (see, e.g., ellipses in Fig. 2). These GMCs are dominated by a lower-density molecular gas, as compared to the internal GMCs, and exhibit signatures of slow shocks (up to $\sim$ 30 km s^-1) rather than intense stellar feedback. This division reflects significant dynamical and chemical variations across the CMZ, previously noted by ALCHEMI studies (e.g., Tanaka et al., 2024; Behrens et al., 2024; Harada et al., 2024).

The study by Harada et al. (2024) utilized ALMA observations to conduct a detailed survey of molecular lines, applying principal component analysis (PCA) to distinguish patterns of chemical excitation and physical conditions across the CMZ. Their results indicate that the internal GMCs (GMCs 3–6) exhibit strong emissions from dense gas tracers such as HCN and HC₃N, along with radio recombination lines (RRLs), which are associated with active star formation and high-excitation regions (e.g., Bendo et al., 2015). In contrast, the external GMCs (GMCs 1, 2, and 7 to 9 in the work of Harada et al. (2024), which did not account for GMC 10) are characterized by tracers of slow shocks, such as CH₃OH and an increment of HNCO over SiO (Meier et al., 2015; Humire et al., 2022), reflecting less energetic dynamical processes. The PCA revealed that molecular emission correlates strongly with dynamical features, with high-excitation species dominating the central regions, while shock tracers are more prevalent in the outskirts.

In the current analysis, the differentiation between internal and external GMCs is further supported by variations in star formation rate (SFR), stellar mass, and dust mass. The internal GMCs exhibit SFRs above 0.087 $M_{\odot}$ yr^-1, as found by GalaPy, which provides an instantaneous SFR listed in the second column of Table 5, as well as greater stellar (third column; Table 5) and dust masses ( $M_{\rm{dust}}$ ; Table 9), consistent with the high-density, high-excitation environments described by Harada et al. (2024). For the sake of clarity, we provide a plot showing the dust mass, $M_{\rm{dust}}$ versus stellar mass $M_{\rm{star}}$ in Fig. 5, they were obtained from GalaPy. In the latter figure, we label the ranges for these quantities for nuclear and external regions or GMCs. The internal or nuclear GMCs (GMCs 3-6) display stellar masses in the 3.7–7.1 $\times$ 10 ${}^{8}\leavevmode\nobreak\ M_{\odot}$ range, and dust masses in the 1.6–5.8 $\times$ 10 ${}^{5}\,M_{\odot}$ range, in strong contrast to the external GMCs, where the stellar and dust masses range between 0.034–1.4 $\times$ 10 ${}^{8}\leavevmode\nobreak\ M_{\odot}$ and 2.6–5.1 $\times$ 10 ${}^{4}\,M_{\odot}$ , respectively.

The above indicates that the dust mass is the clearest distinction between the two primary groups, while the stellar mass, although enhanced in the internal GMCs, is only a factor of two larger in GMC 6 (internal) with respect to GMC 7 (external).

The external GMCs are also characterized by a lower SFR than the inner ones, with values in the 0.005–0.022 $M_{\odot}$ yr^-1 range. These findings align with the chemical and dynamical distinctions previously observed in multiple ALCHEMI studies (e.g., Tanaka et al., 2024; Behrens et al., 2024; Harada et al., 2024), reinforcing the conclusion that the physical properties of the CMZ are tightly coupled with its molecular and dynamical characteristics.

4.2 Optical spectral fitting

We fit archival MUSE data (ID:0102.B-0078(A), PI: Laura Zschaechner) in the 4600–9350Å range with SPLUS photometric data using starlight (Cid Fernandes et al., 2005) in its latest version which is capable to fit spectra and photometry simultaneously (López Fernández et al., 2016; Werle et al., 2019).

To correct for the mismatch between spectroscopic and photometric fluxes, we scale MUSE spectra to match the iSDSS band flux of S-PLUS. The S-PLUS bands used for the fits are uJAVA, J0395, J0410 and J0430, iSDSS. the filter F0378 was not used to avoid contamination by [O II] 3727 emission. Of these, the blue ones are the most important, because they extend the MUSE coverage to the age sensitive 4,000 Å region.

The starlight fits to the MUSE spectra and the fitted S-PLUS photometric points for the ten GMCs are shown by the blue lines and the cyan points, respectively, in Fig. 4. Some important stellar absorption features are labeled in the top panels and indicated by vertical dashed green lines for all the GMCs. The outputs of the fitting are summarized in Table 10.

As is commonly known (see, e.g., Conroy, 2013a), low-mass stars dominate both the total mass and the number of stars in a galaxy, but they contribute only a small fraction to the overall light output, or bolometric luminosity, since young stars outshine older stars. This can be seen in the star formation histories per GMC, which are plotted in Fig. 6. Independent of the molecular cloud, light fraction is always shifted towards young stars, while the mass fraction is concentrated in old stars. For example, in GMC 5, while nearly 100% was already formed 10 Gyr ago $\sim$ 50% of the light comes from stars $\leq$ 30 Myr.

Also in Fig. 6, for GMC 5, we have over-plotted (in magenta) the normalized skewed Gaussian derived from the age-metallicity relation (AMR; see Appendix C), this same distribution was adapted to be cumulative (cAMR) and plotted (in yellow) in the same panel for GMC 5 of Fig. 6. GMCs 3 to 6 also show the estimated ages from super star clusters derived by Butterworth et al. (2024) in vertical dashed blue lines.

GalaPy (black) provides older stellar age estimations than averaged starlight values in GMC 10 and this happens because, although GalaPy is affected by the newest stellar production in the IR (the IR bump), it is not limited by ages of stellar templates, which is the case for starlight. Nevertheless, for all the other GMCs (GMCs 1–9), we observe that the stellar ages derived from GalaPy fall in between the old and new stellar populations that best fit the stellar contribution in the optical using starlight. This is more clearly seen in Fig. 19, where we plot, for all ten GMCs studied in this work, the average old and young stellar population ages from starlight, and the global stellar ages from GalaPy. Both codes agree in that GMC 4 is the youngest in the sample, with younger stellar populations in GMCs 3 to 7 and older ones located in GMCs 1, 2, and 8 to 10. In addition, in the same plot we added results from CIGALE, whose output delivers youngest ages in general with the youngest stellar burst in GMC 5, very close to the youngest burst in GMC 4 derived by GalaPy.

4.2.1 Emission line diagnostic diagrams

The fluxes of emission lines from the MUSE spectra were measured in the pure nebular spectra, obtained after the subtraction of the stellar population continuum from the observed spectra. The contribution of the stellar population was determined using the stellar population synthesis code starlight (Cid Fernandes et al., 2005).

These estimates were used to analyze the dominant ionization mechanism of the emitting gas in the galaxy. To do this, we used two diagnostic diagrams:

(a)

The [O III] $\lambda 5007$ /H $\beta$ versus [N II] $\lambda 6584$ /H $\alpha$ diagnostic diagram proposed by Baldwin et al. (1981), commonly known as the BPT diagram. This diagram includes the theoretical separation between Hii-like and AGN-like objects proposed by Kewley et al. (2001) [K01], and the empirical star-forming limit proposed by Kauffmann et al. (2003) [K03]. The region between the theoretical and empirical limits is commonly referred to as the “composite region,” although there are claims against this designation (see below).
(b)

The WHAN diagram, namely, $\log(\mathrm{EW\leavevmode\nobreak\ H}\alpha)$ versus $\log([\mathrm{N\leavevmode\nobreak\ II}]\lambda 6584)$ (Cid Fernandes et al., 2010), which is used to differentiate the nature of the ionization sources, classifying objects as star-forming, strong AGNs, weak AGNs, or retired galaxies. The line at $\log(\mathrm{EW\leavevmode\nobreak\ H}\alpha)=6\AA\,$ represents the limit between weak and strong AGN emission. We note that in this diagram, he chosen log10([N II] $\lambda 6584$ /H $\alpha$ ; x-axis) limit to disentangle between SF and AGN gaseous ionization origin is normally taken from Stasińska et al. (2006) [S06] in the literature, although for NGC 253’s CMZ we can rely on K03 as there are no GMCs between S06 and K03 limits.

The results are shown in Fig. 7 where we extract the emission lines in our common 3 $\arcsec$ diameter aperture using archival MUSE observations (ID:0102.B-0078(A), PI: Laura Zschaechner). Interestingly, our GMCs do not occupy the regions in these diagrams that are typically associated with star-forming regions. Instead, they are located in what is often referred to as the “composite” zone, lying between star-forming and AGN-dominated areas. According to some studies, the latter behavior can happen in strong starburst, when stellar shocks are fast (e.g., Kewley et al., 2001).

However, given the lack of evidence for an actual AGN in NGC 253, and the fact that these line-ratios are observed well away from the nucleus, our GMCs are most likely not starburst $+$ AGN composites. Instead, we prefer to describe them as “starburst $+$ shocks” composites, where the observed line ratios are a mixture of star-formation and shocks from stellar winds and supernovae. Shocks are indeed known to be present in NGC 253 (Meier et al., 2015; Holdship et al., 2022; Humire et al., 2022; Behrens et al., 2022; Harada et al., 2022; Huang et al., 2023; Behrens et al., 2024; Bouvier et al., 2024; Bao et al., 2024), whereas there is extensive literature arguing against the presence of an AGN (e.g., Fernández-Ontiveros et al., 2009; Brunthaler et al., 2009; Müller-Sánchez et al., 2010).

From the SED perspective, a clear argument against the existence of an AGN is provided in Appendix A, where the AGN fraction, $f_{\rm{AGN}}$ in Table 13, is negligible ( $<$ 7.5%) and also does not peak at GMC 5 as one might expect in the presence of an AGN. Although this behavior could potentially suggest the presence of a low-luminosity AGN (LLAGN), it is more plausibly interpreted as an artifact of CIGALE’s attempt to accurately fit the AGN-free SED. Further arguments specifically against the presence of an LLAGN can be found, for instance, in Mangum et al. (2019, their Sect. 7.1).

4.2.2 Attenuation estimations

The attenuation from dust in Hii regions can be estimated using the Balmer lines H $\alpha$ and H $\beta$ by contrasting their observed ratios with their expected ratios without attenuation. The intrinsic H $\alpha$ /H $\beta$ flux ratio is 2.86 under typical conditions of electron density ( $n_{e}\sim 100$ cm^-3) and temperature ( $T_{e}\sim 10,000$ K), although variations in these conditions can change the ratio by up to 4% (Osterbrock & Ferland, 2006). The so called Balmer extinction is given by

A_{\rm{V}}^{\rm Balmer}=\frac{2.5}{\frac{A_{H\beta}}{A_{\rm{V}}}-\frac{A_{H\alpha}}{A_{\rm{V}}}}\log_{10}\left(\frac{(F_{H\alpha}/F_{H\beta})_{\text{obs}}}{2.86}\right),

(2)

where $A_{\mathrm{H}\beta}/A_{\rm{V}}=1.14$ and $A_{\mathrm{H}\alpha}/A_{\rm{V}}=0.82$ for a Calzetti et al. (2000) reddening law. In the following, we will apply this method to the GMCs studied in this work by measuring $(F_{\mathrm{H}\alpha}/F_{\mathrm{H}\beta})_{\text{obs}}$ for each region and estimating their respective attenuations.

For comparison, we also considered the attenuation model proposed by Charlot & Fall (2000), a widely used approach for estimating dust extinction in stellar populations. This method distinguishes between the contributions of young and old stellar populations, assigning different power-law attenuation slopes to each. The attenuation for young populations, attributed to birth clouds, is represented by $\alpha_{\rm{BCs}}$ , while that for older populations, associated with the ISM, is denoted as $\alpha_{\rm{ISM}}$ , as described in Eq. 4. According to this framework, stars are initially formed within interstellar birth clouds (BCs) and migrate to the ISM after approximately $10^{7}$ years. The transition between these two attenuation regimes occurs at a wavelength of 5,500 Å. In CIGALE, this model can be adopted using the dustatt_2powerlaws module. The functional form of the law is commonly expressed in the literature as:

A(\lambda)=\begin{cases}A_{\lambda}(\text{BC})+A_{\lambda}(\text{ISM}),&\text{for young stars, age}<10^{7}\,\text{years},\\ A_{\lambda}(\text{ISM}),&\text{for old stars, age}>10^{7}\,\text{years}.\end{cases}

(3)

where $A_{\lambda}$ is the wavelength-dependent attenuation expressed as $A_{\lambda}=A_{V}\left(\frac{\lambda}{\lambda_{V}}\right)^{\delta}$ , with $\lambda_{V}=5,500\,\text{Å}$ , and $A_{V}$ represents the attenuation in the V band. In practice, this translates to the following expression:

A(\lambda)=A_{V}\left[\left(\frac{\lambda}{\lambda_{0}}\right)^{\alpha_{\rm{BCs}}}(\lambda<5500\,\text{\AA})+\left(\frac{\lambda}{\lambda_{0}}\right)^{\alpha_{\rm{ISM}}}(\lambda\geq 5500\,\text{\AA})\right],

(4)

where $A_{V}$ is the overall attenuation at 5,500 Å that we assume to be unity for simplicity. The exponents $\alpha_{\rm{BCs}}$ and $\alpha_{\rm{ISM}}$ describe the slopes for the birth clouds and the ISM, respectively. In the implementation, the values for $\alpha_{\rm{BCs}}$ range from -0.69 to -0.9, while $\alpha_{\rm{ISM}}$ ranges from -0.75 to -0.3. We find these values starting from those commonly used in Charlot & Fall (2000) and also being influenced by other works that assume a -0.7 (an original assumption from C&F 2000) value for both stellar populations. Then, looking at our GalaPy-derived curves from the SED, we set values that better reproduce the observations.

We compute the attenuation curves over a wide wavelength range, from the near-UV to the near-infrared, using these parameter ranges. The curves are normalized at 5,500 Å to allow a consistent comparison. As shown in Fig. 8, the model’s attenuation curves are sensitive to the chosen values of $R_{\rm{V}}$ for the Calzetti et al. (2000) attenuation law, which assumes a single attenuation slope for all stellar populations, and the values of $\alpha_{\rm{BCs}}$ and $\alpha_{\rm{ISM}}$ in the case of the Charlot & Fall (2000) models.

Table 7 lists the measured Balmer decrement, the corresponding $A_{\rm{V}}^{\mathrm{Balmer}}$ , as well as the $A_{\rm{V}}$ values obtained with CIGALE and starlight. As expected from the extremely red shape of the optical spectra (see Fig. 4), we obtain large values of $A_{\rm{V}}^{\mathrm{Balmer}}$ : from 1 to 6 mag. The largest values are found for GMCs 3, 4, 5, and 6, which are closer to the galactic bar and in dense star-forming areas. We note, however, the very substantial uncertainties in $(F_{\mathrm{H}\alpha}/F_{\mathrm{H}\beta})_{\text{obs}}$ (ranging from 10 to nearly 90%). These happen because of the weak $\mathrm{H}\beta$ emission, which is barely visible in most cases in Fig. 4.

Table 7 further lists the starlight results for the extinction. In the case of starlight the table gives $\tilde{A}_{\rm{V}}$ , which we define as the average of the extinctions affecting all populations ( $A_{\rm{V}}$ ) and that affecting only those younger than 10 Myr ( $A_{\rm{V}}+\delta A_{\rm{V}}$ ), weighting them by the fluxes of the respective components as derived from the fits. Again, the extinction values are large (1.8–5.7 mag), as expected from the redness of the optical continuum. The derived values are close to, but in general somewhat smaller than those derived from the Balmer decrement.

Table 7: Updated attenuation estimates for GMCs derived from the Balmer decrement, CIGALE, and effective

A_{\rm{V}}

for starlight.

GMC	$(F_{H\alpha}/F_{H\beta})_{\text{obs}}$	$A_{\rm{V}}^{\text{Balmer}}$	$A_{\rm{V}}^{\text{CIGALE}}$	$\tilde{A_{\rm{V}}}^{\textsc{starlight}}$
		[mag]	[mag]	[mag]
1	$3.81\pm 0.32$	$1.08\pm 0.11$	$3.57\pm 0.02$	$2.12$
2	$8.32\pm 1.70$	$2.89\pm 0.30$	$3.58\pm 0.02$	$2.42$
3	$29.79\pm 16.90$	$5.06\pm 0.56$	$6.04\pm 0.11$	$4.44$
4	$45.69\pm 37.94$	$5.87\pm 0.67$	$6.06\pm 0.00$	$5.67$
5	$43.43\pm 37.37$	$5.81\pm 0.65$	$6.24\pm 0.02$	$4.66$
6	$29.14\pm 16.45$	$5.04\pm 0.54$	$6.08\pm 0.01$	$4.10$
7	$21.21\pm 10.11$	$4.38\pm 0.47$	$5.07\pm 0.02$	$2.97$
8	$19.18\pm 10.04$	$4.18\pm 0.45$	$3.53\pm 0.01$	$2.28$
9	$17.01\pm 7.84$	$3.96\pm 0.42$	$3.54\pm 0.01$	$1.84$
10	$11.13\pm 5.05$	$3.36\pm 0.37$	$3.54\pm 0.12$	$2.11$

In Table 7, we provide our results. For starlight, we list both the value of $A_{\rm{V}}$ affecting all populations $\tilde{A}_{V}$ and the flux-weighted extinction considering both $A_{\rm{V}}$ and the additional extinction ( $\delta A_{\rm{V}}$ ) applied only to populations younger than 10 Myr to account for the dust in their birth clouds. Table 7 also lists the CIGALE values of $A_{\rm{V}}$ , which differ from those derived from the Balmer lines, $A_{\rm{V}}^{\rm Balmer}$ , by a 16% in median. On the other hand, starlight-based values also differ from those of $A_{\rm{V}}^{\rm Balmer}$ by a median value of 27%.

The only GMC where $\tilde{A}_{V}$ exceeds $A_{\rm{V}}^{\rm Balmer}$ is GMC 1. This discrepancy may suggest additional sources of attenuation, potentially linked to warm or dense dust. One possible explanation is the influence of strong shocks, which could sputter dust grains and enhance local dust densities. Supporting evidence for the presence of shocks in GMC 1 comes from its elevated SiO $J=$ 5–4/HNCO (10_0,10–9_0,9) transition ratios observed at about 220 GHz (ALMA Band 6). This region shows the highest ratio values (about 2), as noted in the right panel of Fig. 11 in Humire et al. (2022). These elevated ratios are well-known tracers of strongly shocked gas, indicating a prevalence of strong (traced by SiO) over weak/mild shocks (traced by HNCO). Additionally, GMC 1 hosts methanol masers with the strongest deviations from LTE expectations among the 10 GMCs, as illustrated in the left panel of Fig. 13 in the same study. While these observations strongly point to the presence of shocks, they are not direct indicators of warm or dense dust. However, shocks may indirectly affect the dust environment by altering its spatial distribution or properties. Future studies are needed to assess whether gas and dust are tightly coupled in this region and to determine the extent to which shocks influence the observed attenuation.

We note that, even though the $R_{\rm{V}}$ parameter was set to 3.1 in CIGALE, which corresponds to the typical value of the dust extinction law for the Milky Way (Calzetti et al., 2000), the differences between the resulting $A_{\rm{V}}$ and any other attenuation-related parameter in Table 11 do not differ by that 3.1 factor from $A_{\rm{V}}$ . The $R_{\rm{V}}$ parameter controls the relationship between the total and the selective extinctions of starlight by dust. A higher value indicates that, as the amount of dust increases, it less affects the attenuation, something that may be clearer seen in the following expression: $A_{\rm{V}}/R_{\rm{V}}=E(B-V)$ .

GalaPy derived attenuation curves are shown in Fig. 8, where curves indicate the attenuation in the visible band and legends show its value at 5500Å, on a linear scale, to be compared with starlight outputs (Table 10). The results from this code are normalize at that wavelength value. $R_{V}$ values following Calzetti et al. (2000) are also included for comparison. These values describe the relationship between the total amount of dust extinction and the amount of reddening it causes in the context of the Milky Way’s dust extinction curve and corresponds to the ratio between the total extinction $A_{V}$ and the color excess $E(B-V)$ . Considering the ten GMCs studied, most $R_{V}$ values range between 6 and 9, far beyond the Galactic value of 3.1. Surprisingly, extreme values come from clouds at the boundaries of the CMZ, with the highest $R_{V}$ lying in GMC 1, 9, and 10, and the lowest value in GMC 2.

As a final remark, when comparing the different attenuation estimates, it is important to recognize that the gas attenuation (from the Balmer decrement) and stellar attenuation (from GalaPy, CIGALE, and STARLIGHT in our case) have distinct origins. Recent studies suggest that these two types of attenuation should become increasingly similar as the observed region ages (Li & Li, 2024). Although, we do not see a clear trend for attenuation differences between the ones derived from the Balmer decrement and those obtained from CIGALE, namely $\Delta A_{\rm{V}}=A_{\rm{V}}^{\rm{Balmer}}-A_{\rm{V}}^{\rm{CIGALE}}$ , taken from Table 7, versus derived stellar ages from either GalaPy or CIGALE (last columns of Tables 5 and 6, respectively), we do see a trend when comparing Attenuation differences between Balmer lines and the gas attenuation derived from the stellar continuum by starlight versus stellar ages derived by CIGALE. For an easy visualization of that finding, in Fig. 9 we present such a relation.

4.3 Star Formation Rates

4.3.1 Kennicutt relation

There are many ways to obtain the SFR from monochromatic emission. One of the most commonly used is the relation found by Kennicutt (1998) regarding H $\alpha$ :

\rm{SFR}\leavevmode\nobreak\ [\mathrm{M_{\odot}\leavevmode\nobreak\ year^{-1}}]=7.9\times 10^{-42}L(\mathrm{H\alpha})\leavevmode\nobreak\ [\mathrm{erg\leavevmode\nobreak\ s^{-1}}].

(5)

where $L(\mathrm{H\alpha})$ is corrected for extinction using $A_{H\alpha}=0.82\times A_{\rm{V}}^{\rm Balmer}$ .

Fig. 12 shows a comparison between the SFR obtained from the above considerations and the one obtained from the SED fitting using GalaPy and CIGALE, while the dashed-red line indicates a one-to-one proportion. To create this plot we considered the attenuation of the H $\alpha$ emission line derived SFR considering the $A_{\rm{V}}^{\rm{Balmer}}$ , since it is expected to be the most accurate method: starlight and SED derived attenuation rely on the continuum level.

We notice that SFRs derived by GalaPy (red symbols in Fig. 12) tend to be larger than expected by the Kennicutt relation, this is an expected result as the H $\alpha$ ’s derived SFR corresponds to the stellar formation averaged over the last 10 Myr (see, e.g., Stroe et al., 2015). Interestingly, however, considering the 10 Myr and 100 Myr CIGALE’s averaged SFRs (grey symbols in Fig. 12), the Kennicutt relation is followed at a higher degree when considering 100 and not 10 Myr values (see Table 6), and is less dispersed in relation to the instantaneous SFR derived by CIGALE (SFR_INST; blue symbols in Fig. 12).

4.3.2 Radio recombination lines

There are several ways to estimate the star formation rate of an astronomical source. One of the most accurate methods not affected by dust extinction involves radio recombination lines (RRLs). Following the approach described by Bendo et al. (2015), where the H $40\alpha$ line at 99.02 GHz was utilized, we determined the SFR for the central nuclear regions (specifically GMCs 3 to 6). As shown in Figure 13, the SFR density per spaxel aligns with values obtained by GalaPy for the same regions where RRLs are detected. Our calculation, focusing on GMCs 3 to 6, within a 3″ aperture, yielded a value of 1.852 M_⊙ yr^-1 (refer to Table 9), which is consistent with the SFR of 1.73 $\pm$ 0.12 reported by Bendo et al. (2015), considering uncertainties. We highlight that the SFR obtained by GalaPy corresponds to the instantaneous (i.e., current) SFR while RRLs are expected to also trace recent SF of the last $\sim$ 5–10 Myr (Scoville & Murchikova, 2013; Emig et al., 2020). It is worth noting that due to our 3″ aperture extraction, there is some overlap in the fluxes from GMCs 5 and 6. This overlap likely contributes to our estimations being at the upper limit of RRL estimations. However, our photometric extraction does not include fluxes between GMCs, as illustrated in Fig. 13. Considering a 9 $\arcsec$ -aperture-radius SED involving some of the mentioned inner GMCs (GMCs 4 to 6), as it can be seen in the bottom inset of the left panel of Fig. 20 and considering Herschel PACS observations, we obtain a SFR of 1.620 ${}_{-0.156}^{+0.201}$ M_⊙ yr^-1. After adding the SFR from GMC 3, which is not covered by the 9 $\arcsec$ aperture and is 0.087 ${}_{-0.009}^{+0.011}$ M_⊙ yr^-1, we obtain a SFR for GMCs 3 to 6 of 1.707 ${}_{-0.156}^{+0.201}$ , also consistent with the results of Bendo et al. (2015). However, the latter should be taken as a lower limit because, although devoid of any slight overlapping issue between GMCs 5 and 6, the 9 $\arcsec$ aperture does not fully encompass GMC 4, as shown in the lower inset of Fig. 20’s left panel.

4.3.3 Radio cm continuum and FIR correlation

As summarized in Murphy et al. (2011), there are several SFR tracers in the wavelength range considered in this work. According to these authors, the 33GHz continuum emission is a reliable indicator of current star formation, as it primarily arises from thermal free-free emission. Although this frequency can sometimes be affected by rotating dust, this is not the case for NGC253 at the studied scales, as shown in the SEDs. We indeed find a strong correlation between the emission in the EVLA Ka Band, as originally published in Kepley et al. (2011), and the total SFR computed by GalaPy. Furthermore, accounting for scaling offsets, we also observe generally good correlations—close to a 1:1 correspondence—between cm tracers and the GalaPy results, as seen in Fig.10. The latter was obtained using the same approach as for the 33GHz, namely, by performing a linear regression in logarithmic space with all retrieved (E)VLA continuum images that detect at least four GMCs (see Table 5), aiming to obtain the best-fit assembly valid for the GMCs in the starburst galaxy NGC 253. A caveat should be noted, as different galaxy types (e.g., MW-like, dwarf galaxies, AGN hosts) may present different relationships. In this regard, an immediate prospect for the current work is to evaluate the different relations between SFR tracers.

Accounting for IR points as SFR tracers, shown in Fig. 11, we generally find a lower dispersion and higher quality of the fit in terms of mean R in relation to the radio bands. However, that is an artifact since we are not accounting for the error bars in any of the plots, which are larger in the x-axis, as can be seen in the Fig. Given that our FIR points are taken from the SED fit itself, they are not independent to the derived SFR, although we note larger dispersion when altering the SED points in other bands (not shown). The general conclusion that we can still obtain from the IR tracers is that the best ones correspond to the integrated IR emission between 8 and 1000 $\mu$ m, normally referred as total IR emission and the IR emission at 60 $\mu$ m.

We note that, contrary that for the case of optical lines, we do not assume a relation factor a priory, we only convert our mJy measurements to erg s Hz for an easy comparison between this study and the literature. For the mJy erg s Hz conversion we have assumed a distance of 3.5 Mpc to NGC 253 (Rekola et al., 2005).

5 Discussion

5.1 Reliability of our SED model along the FIR window

One of the critical, if not the most critical, challenges that astronomers have to face when attempting to fit a panchromatic SED at high angular resolution is the impossibility to acquire far infrared (FIR) observations at resolutions below 6 $\arcsec$ . Current facilities in the regime, such as the James Webb (JWST) and Euclid telescopes, cover mainly the near-IR and mid-IR part of the electromagnetic spectrum, with the largest wavelength window ranging from 0.5 (Euclid) to 28 $\mu$ m (JWST). Indeed, only the JWST can keep angular resolutions below $\sim$ 0.8 $\arcsec$ at 28 $\mu$ m. However, the FIR bump primarily occurs at wavelengths larger than those, mainly from about 50 to 300 $\mu$ m.

While we made use of our well covered sub-mm regime, with ALCHEMI data, to fine-tune the Rayleigh-Jeans tail, at the emergence of the FIR bump on the radio side, and used VLT/NACO data at 11 and 18.7 $\mu$ m on the MIR side for GMCs 3–7, we know the caveat that this is likely not enough to properly ensure the fit on the IR bump itself, relying on our SED model extrapolations. This is why we further retrieved archival Herschel PACS observations selecting level 2.5 PACS images at 70, 100, and 160 $\mu$ m. While these FIR observations are not able to resolve the targeted GMCs, they are sufficient to fit emission from the inner part of the CMZ, covering at most GMCs 3 to 6 in the worst-case scenario (at 11.3 $\arcsec$ ), and partially covering GMCs 4 to 6 at 6 and 9 $\arcsec$ (see Appendix E for more information on this latter).

In Fig. 20 it can be noted that the far infrared peak agrees between the 3 $\arcsec$ aperture SED, devoid of FIR photometric points and the SED at 9 $\arcsec$ , whose diffuse dust temperature, $T_{\rm{dd}}$ , of 68.833 ${}_{-13.916}^{+7.808}$ K is in agreement to the ones of the covered GMCs (GMCs 4–6), as noted in Table 5. Specifically, we get $T_{\rm{dd}}$ estimations of 84.532 ${}_{-12.551}^{+6.791}$ , 67.929 ${}_{-3.896}^{+4.023}$ , and 63.318 ${}_{-2.800}^{+3.233}$ K for GMCs 4, 5, and 6, respectively. All in all, this arguments in favor of our 3 $\arcsec$ apertures as the FIR bump and the SED shape in general peaks at a very similar point. It is worth noting that the SED shape is responsible for specific star formation rate estimates and, thus, for the stellar mass and SFR estimates as well (Conroy, 2013a).

Given the lack of FIR photometric points in the FIR regime at 3 $\arcsec$ , our SED models experienced considerable freedom in this specific wavelength range. We found that GalaPy produced SED shapes that were closer, compared to CIGALE’s models, to those observed at 9 $\arcsec$ , where Herschel PACS observations are available. Nevertheless, it is worth emphasizing that CIGALE performs well in fitting the SED of GMC 5 at 9 $\arcsec$ (which also includes GMCs 4 and 6), both with and without an AGN component, as discussed in Appendix A.

Table 8: Comparison of specific star formation rates (sSFR) in log₁₀ scale derived from GalaPy and CIGALE for our ten selected GMCs.

GMC	sSFR (GalaPy)	sSFR (CIGALE)
	[log₁₀ yr^-1]	[log₁₀ yr^-1]
1	$-8.84\pm 0.09$	$-9.38\pm 0.07$
2	$-9.63\pm 0.05$	$-9.68\pm 0.05$
3	$-9.63\pm 0.05$	$-9.40\pm 0.11$
4	$-9.41\pm 0.04$	$-9.69\pm 0.02$
5	$-8.81\pm 0.04$	$-8.88\pm 0.03$
6	$-9.33\pm 0.04$	$-9.26\pm 0.05$
7	$-9.81\pm 0.04$	$-9.40\pm 0.07$
8	$-9.80\pm 0.04$	$-9.83\pm 0.11$
9	$-9.80\pm 0.05$	$-9.63\pm 0.07$
10	$-9.88\pm 0.09$	$-9.66\pm 0.07$

5.2 Comparison to high-z dusty starbursts

The question of whether local starbursts can serve as analogues of primeval dusty galaxies at high redshifts ( $z>4-8$ ) remains a topic of active debate. A useful parameter in this context is the dust-to-stellar mass ratio, which provides insight into the balance between dust production and destructive processes such as supernova-driven shocks and astration (Nanni et al. 2020, Donevski et al. 2020, Kokorev et al. 2021, Donevski et al. 2023, Witstok et al. 2023, Palla et al. 2024, Sawant et al. 2024).

To place the ISM conditions of NGC 253 into perspective against those of distant dusty galaxies, we compare the dust-to-stellar mass ratios ( $M_{\rm dust}/M_{\star}$ ) derived in this work to those inferred from ALMA observations of high- $z$ starbursts (Donevski et al. 2020, Salak et al. 2024, Sawant et al. 2024). Regardless of the SED fitting code used, all GMCs in NGC 253 exhibit relatively high dust-to-stellar mass ratios, with a median value of $\log(M_{\rm dust}/M_{\star})=-2.32\pm 0.29$ (values from CIGALE), or $\log(M_{\rm dust}/M_{\star})=-2.89\pm 0.45$ , if values from GalaPy are considered. These ratios are close to those observed in young dusty starburst galaxies in the distant Universe. Specifically, the average dust-to-stellar mass ratio is $\log(M_{\rm dust}/M_{\star})=-2.13\pm 0.3$ at $z\sim 2-5$ (Donevski et al. 2020) and $\log(M_{\rm dust}/M_{\star})=-2.58\pm 0.4$ at $z\sim 4-6$ (Sawant et al. 2024). A caveat in comparing the dust-to-stellar mass ratio found here with literature is in the method used to derive galaxy properties. Both in Donevski et al. (2020) and Sawant et al. (2024), the SED fitting is performed using CIGALE, which results in average lower values with respect to GalaPy, as also reported above.

Moreover, inferred attenuation values in the GMCs of NGC 253 are comparable to those observed in high- $z$ dusty starbursts but also significantly higher than those typically found in JWST-detected galaxies at $z>6-8$ (e.g., Nanni et al. 2020, Álvarez-Márquez et al. 2023, Markov et al. 2024). It is important to note that, despite their comparable ”dustiness,” GMCs in NGC 253 have $M_{\star}$ , SFR and sSFR several orders of magnitude lower than their high- $z$ counterparts. In such low stellar mass regimes, dust production is generally expected to be dominated solely by stellar sources, such as Type II supernovae (Galliano et al. 2018, Sommovigo et al. 2020).

However, recent models suggest that dust production in young starbursts may be significantly diminished due to ISM removal in the presence of strong stellar outflows, similar to those observed in NGC 253. To reconcile the high dust-to-stellar mass ratios, these models propose efficient growth of large dust grains through ISM accretion in rapidly cooling gas within stellar outflows (Hirashita & Chen 2023, Donevski et al. 2023, Romano et al. 2024, Palla et al. 2024). While a detailed exploration of this intriguing possibility lies beyond the scope of this work, we note that the high dust-to-stellar mass ratios and attenuations observed in the GMCs may signal that NGC 253 is undergoing a shift in the dominant dust production mechanisms, as proposed by recent simulations (see e.g., Narayanan et al. 2024, also Schneider & Maiolino 2024 for a review).

6 Conclusions and Prospects

This study represents a significant step forward in understanding the CMZ of NGC 253, utilizing the first panchromatic SED analysis at a 51-pc resolution towards extragalactic regions. The integration of diverse modeling approaches, from GalaPy to CIGALE and starlight on the stellar continuum features, enabled a thorough characterization of physical properties such as star formation rates, dust and stellar masses, and attenuation effects across ten GMCs using both parametric and non-parametric models.

The integration of multi-wavelength SED modeling tools, such as GalaPy and CIGALE plus the integration of mono-wavelength tracers in the optical using starlight and the sub-mm using CASSIS, provided detailed insights into star formation histories, attenuation, and stellar population properties at different timescales. High extinction values (e.g., A_V exceeding 5 magnitudes in internal GMCs) were robustly derived, aligning well with other methods and showcasing the importance of panchromatic approaches.

We summarize our main findings below.

1.
Dynamical and Chemical Diversity:
1. •
  
  The Central Molecular Zone (CMZ) of NGC 253 exhibits significant dynamical and chemical distinctions between internal (nuclear) and external GMCs.
2. •
  
  Internal GMCs (e.g., GMCs 3–6) are characterized by higher star formation rates (SFR), stellar masses, and dust masses, consistent with active star formation, high-excitation regions, and higher metallicities.
3. •
  
  External GMCs (e.g., GMCs 1, 2, 7–10) display lower density, slower shocks, and chemical signatures indicative of less energetic processes as compared to the internal GMCs.
2.
Comparison of Internal vs. External GMCs:
1. •
  
  Internal GMCs exhibit stellar masses in the range of $3.7-7.1\times 10^{8}M_{\odot}$ and dust masses in the range of $1.6-5.8\times 10^{5}M_{\odot}$ , both greatly surpassing (doubling) maximum values observed in external GMCs.
2. •
  
  Specific star formation rates (sSFR) further illustrate these differences, with significantly higher activity in the nuclear regions.
3.
Insights from Age-Metallicity Relations:
1. •
  
  Differences in stellar ages and metallicities highlight the varied evolutionary states of the studied GMCs.
4.
Emission Line Diagnostics:
1. •
  
  The BPT and WHAN diagrams indicate shock signatures in the nuclear region of NGC 253, placing it in the composite zone. However, these shocks likely result from stellar winds and supernovae in the starburst environment, not from AGN activity.
5.
Star Formation, Panchromatic vs. Monochromatic Results:
1. •
  
  We have found a strong 1:1 correlation between cm tracers and SED derived SFR values, being the 33 GHz (Ka) band the best among the available dataset. This is in line with previous works based on entire galaxy properties (e.g., Murphy et al., 2011), and is now confirmed at GMC scales. Additionally, we validate the effectiveness of RRLs (H40 $\alpha$ ) as a SF tracer, although limited to the brightest sources (GMCs 3–6).

The current study provides a robust framework for leveraging spatially-resolved multi-wavelength datasets to investigate the interplay between star formation, gas dynamics, and attenuation in starburst galaxies. In addition, this work not only advances our knowledge of starburst environments but also serves as a benchmark for studying high-redshift galaxies, where similar extreme conditions prevail.

In the future, we plan to focus on several key directions to broaden the applicability of these findings. These include utilizing JWST instead of Spitzer observations, achieving aperture diameters as fine as 1″ in the mid-IR, enabling more detailed analyses of star-forming regions and/or augmenting the redshift regime for our potential targets. We plan to apply the same approach to other well-known starburst like M 82. Alternatively, we want to improve the robustness of our fit by restricting ourselves to very nearby galaxies where FIR data from Herschel can still resolve GMC-scale structures. Additionally, this framework has the potential to be extended to diverse systems, including Milky Way like galaxies, LINERs, ULIRGs, and AGN, to explore differences in star formation processes across various galactic environments. Lastly, focused studies of GMCs in our Galaxy will offer a critical baseline for contrasting extragalactic GMC properties, contributing to a thorough understanding of star formation under varying conditions.

Table 9: Gas, star-formation (SF) related parameters and dust-related parameters for all SED results derived from GalaPy. For main SF related parameters, see Table 5.

GMC	$T_{\rm{MC}}$	$M_{\rm{gas}}$	ism.R_MC	$Z_{\rm{gas}}$	sfh.psi_max	ism.f_MC	ism.tau_esc
	[K]	[log₁₀ ( $M_{\odot}$ )]	[log₁₀ (pc)]	[ $\times 10^{-3}$ ]	[log₁₀ ( $M_{\odot}$ yr^-1)]		[log₁₀ (yrs)]
1	50.365 ${}_{-6.532}^{+11.259}$	8.573 ${}_{-0.173}^{+0.095}$	0.90 ${}_{-0.48}^{+0.60}$	0.53 ${}_{-0.12}^{+0.27}$	-1.69 ${}_{-0.24}^{+0.23}$	0.01 ${}_{-0.01}^{+0.01}$	7.75 ${}_{-0.43}^{+0.34}$
2	67.295 ${}_{-4.749}^{+6.536}$	7.126 ${}_{-0.119}^{+0.129}$	1.02 ${}_{-0.13}^{+0.13}$	4.11 ${}_{-0.16}^{+0.20}$	-1.67 ${}_{-0.06}^{+0.07}$	0.02 ${}_{-0.01}^{+0.01}$	9.19 ${}_{-0.26}^{+0.17}$
3	66.739 ${}_{-6.464}^{+7.394}$	8.010 ${}_{-0.439}^{+0.260}$	0.79 ${}_{-0.18}^{+0.20}$	8.77 ${}_{-0.94}^{+2.50}$	-0.54 ${}_{-0.17}^{+0.38}$	0.04 ${}_{-0.02}^{+0.09}$	8.83 ${}_{-1.96}^{+0.29}$
4	61.899 ${}_{-10.657}^{+22.586}$	6.925 ${}_{-0.124}^{+0.135}$	1.27 ${}_{-0.40}^{+0.41}$	25.38 ${}_{-2.09}^{+3.44}$	1.15 ${}_{-0.14}^{+0.22}$	0.43 ${}_{-0.23}^{+0.13}$	7.98 ${}_{-0.45}^{+0.25}$
5	54.833 ${}_{-5.765}^{+6.791}$	7.386 ${}_{-0.134}^{+0.133}$	1.60 ${}_{-0.35}^{+0.28}$	21.14 ${}_{-1.45}^{+1.57}$	0.84 ${}_{-0.12}^{+0.12}$	0.20 ${}_{-0.10}^{+0.17}$	6.98 ${}_{-0.68}^{+0.71}$
6	35.814 ${}_{-5.174}^{+13.118}$	7.441 ${}_{-0.121}^{+0.095}$	1.53 ${}_{-0.60}^{+0.34}$	14.70 ${}_{-1.08}^{+1.09}$	0.26 ${}_{-0.12}^{+0.11}$	0.35 ${}_{-0.12}^{+0.13}$	6.18 ${}_{-0.13}^{+0.40}$
7	40.327 ${}_{-7.569}^{+10.562}$	6.833 ${}_{-0.164}^{+0.288}$	1.41 ${}_{-0.51}^{+0.37}$	8.77 ${}_{-1.32}^{+1.00}$	-0.54 ${}_{-0.25}^{+0.17}$	0.38 ${}_{-0.19}^{+0.16}$	7.64 ${}_{-1.39}^{+1.10}$
8	59.476 ${}_{-2.716}^{+3.305}$	7.320 ${}_{-0.163}^{+0.147}$	1.05 ${}_{-0.11}^{+0.13}$	3.04 ${}_{-0.32}^{+0.41}$	-0.51 ${}_{-0.09}^{+0.10}$	0.02 ${}_{-0.01}^{+0.01}$	8.12 ${}_{-0.28}^{+0.25}$
9	64.423 ${}_{-8.143}^{+16.924}$	7.079 ${}_{-0.175}^{+0.234}$	0.54 ${}_{-0.33}^{+0.38}$	4.77 ${}_{-0.54}^{+0.65}$	-1.45 ${}_{-0.18}^{+0.19}$	0.04 ${}_{-0.02}^{+0.03}$	9.22 ${}_{-0.15}^{+0.12}$
10	41.421 ${}_{-2.957}^{+3.811}$	7.250 ${}_{-0.247}^{+0.213}$	1.53 ${}_{-0.65}^{+0.31}$	3.48 ${}_{-0.31}^{+0.37}$	-1.92 ${}_{-0.14}^{+0.15}$	0.11 ${}_{-0.05}^{+0.09}$	9.02 ${}_{-0.56}^{+0.40}$

GMC	$T_{\rm{dd}}$	$M_{\rm{dust}}$	ism.Rdust	ism.f_PAH	noise.f_cal	$\chi^{2}_{red}$
	[K]	[log₁₀ ( $M_{\odot}$ )]	[log₁₀ (pc)]
1	34.643 ${}_{-3.512}^{+2.957}$	4.410 ${}_{-0.195}^{+0.210}$	1.08 ${}_{-0.09}^{+0.10}$	0.17 ${}_{-0.06}^{+0.10}$	-0.64 ${}_{-0.11}^{+0.11}$	2.27
2	33.528 ${}_{-2.875}^{+4.437}$	4.435 ${}_{-0.105}^{+0.115}$	1.08 ${}_{-0.07}^{+0.07}$	0.32 ${}_{-0.07}^{+0.06}$	-2.25 ${}_{-0.14}^{+0.14}$	1.80
3	34.110 ${}_{-3.781}^{+14.707}$	5.766 ${}_{-0.315}^{+0.204}$	1.64 ${}_{-0.20}^{+0.11}$	0.12 ${}_{-0.07}^{+0.04}$	-0.61 ${}_{-0.06}^{+0.07}$	1.47
4	84.532 ${}_{-12.551}^{+6.791}$	5.195 ${}_{-0.108}^{+0.153}$	1.17 ${}_{-0.08}^{+0.15}$	0.04 ${}_{-0.02}^{+0.02}$	-0.55 ${}_{-0.07}^{+0.07}$	1.32
5	67.929 ${}_{-3.896}^{+4.023}$	5.566 ${}_{-0.107}^{+0.133}$	1.43 ${}_{-0.09}^{+0.08}$	0.03 ${}_{-0.01}^{+0.02}$	-0.62 ${}_{-0.08}^{+0.08}$	1.14
6	63.318 ${}_{-2.800}^{+3.233}$	5.453 ${}_{-0.102}^{+0.079}$	1.35 ${}_{-0.08}^{+0.06}$	0.04 ${}_{-0.01}^{+0.01}$	-0.73 ${}_{-0.08}^{+0.10}$	1.35
7	62.911 ${}_{-10.806}^{+4.631}$	4.474 ${}_{-0.135}^{+0.145}$	1.08 ${}_{-0.10}^{+0.12}$	0.06 ${}_{-0.03}^{+0.02}$	-0.94 ${}_{-0.09}^{+0.10}$	1.38
8	41.419 ${}_{-2.629}^{+4.510}$	4.492 ${}_{-0.127}^{+0.147}$	1.21 ${}_{-0.09}^{+0.08}$	0.17 ${}_{-0.05}^{+0.07}$	-1.38 ${}_{-0.10}^{+0.11}$	1.35
9	43.383 ${}_{-4.396}^{+5.060}$	4.679 ${}_{-0.124}^{+0.133}$	1.17 ${}_{-0.09}^{+0.08}$	0.15 ${}_{-0.04}^{+0.07}$	-1.11 ${}_{-0.08}^{+0.09}$	1.29
10	41.532 ${}_{-3.377}^{+5.290}$	4.574 ${}_{-0.117}^{+0.128}$	1.20 ${}_{-0.08}^{+0.08}$	0.19 ${}_{-0.05}^{+0.06}$	-1.25 ${}_{-0.09}^{+0.10}$	1.31

GMC adev Lum_tot $\log M_{*}$ v0 $\sigma_{*}$ $A_{\rm{V}}$ $\delta A_{V}$ x( $\delta A_{V}$ $>$ 0) $\tilde{A_{\rm{V}}}$ $\langle\log t_{*}\rangle_{L}$ $\langle\log t_{*}\rangle_{M}$ $\langle\log Z/Z_{\odot}\rangle_{L}$ $\langle\log Z/Z_{\odot}\rangle_{M}$ 1 1.85 9.82 4.94 -2.41 107.86 2.12 3.85 0.00 2.12 8.95 9.79 0.06 0.15 2 2 234.68 5.27 -11.02 113.8 2.22 7.24 2.76 2.42 9.15 8.77 -0.07 -0.07 3 4.3 160.71 6.28 11.63 154.93 4.44 0 23.78 4.44 8.72 9.86 0.02 0.13 4 3.56 6241.17 6.23 1.37 83.82 3.35 5.77 40.23 5.67 7.33 7.22 0.22 0.17 5 3.46 229.04 6.35 -18.47 126.71 4.66 0 4.75 4.66 8.51 9.89 -0.09 0.26 6 3.03 251.04 6.05 -39.4 125.21 3.59 2.04 24.95 4.1 8.66 9.88 -0.11 -0.25 7 2.73 310.66 5.8 -62.53 133.2 1.85 3.94 28.56 2.97 8.57 9.82 0.46 0.54 8 2.77 40.41 5.63 -59.2 106.07 2.28 4.46 0.00 2.28 9.69 9.95 -0.4 -0.51 9 2.65 22.74 5.35 -56.9 116.32 1.84 4.22 0.00 1.84 9.57 9.78 -0.12 -0.18 10 2.14 74.93 5.01 -69.57 124.54 1.43 3.61 18.89 2.12 8.6 9.08 0.55 0.55

Table 10: Best-fitting results derived for each studied GMC by starlight, described in Sect. 4.2.

Acknowledgements

P.K.H. gratefully acknowledges the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) for the support grant 2023/14272-4. P.K.H acknowledges L. Barcos, E. Caux, and G. Ortiz-León for their help regarding the initial aims of this project, the FIR cross-check methodology, and help with EVLA X-Band data reduction, respectively. We thank A. Kepley for providing us with the 33 GHz (Ka-band) EVLA image. S.D has been supported by the Polish National Science Center project UMO-2023/51/D/ST9/00147 TR is supported by the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations). D.D. acknowledges support from the National Science Center (NCN) grant SONATA (UMO-2020/39/D/ST9/00720). JAFO acknowledges financial support by the Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033), by “ERDF A way of making Europe” and by “European Union NextGenerationEU/PRTR” through the grants PID2021-124918NB-C44 and CNS2023-145339; MCIN and the European Union – NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA. M.H. acknowledges the support by the National Science Centre, Poland (UMO-2022/45/N/ST9/01336). V.M.R. and L.C. acknowledge support from the grant 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. V.M.R also acknowledges support from the grant 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); and from the grant CNS2023-144464 funded by MICIU/AEI/10.13039/501100011033 and by “European Union NextGenerationEU/PRTR”. A.R.L. acknowledges financial support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Agencia I+D+i (PICT 2019–03299) and Universidad Nacional de La Plata (Argentina). A.C.K thanks Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) for the support grant 2024/05467-9. SP is supported by the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America. R.D. gratefully acknowledges support by the ANID BASAL project FB210003

For the near-UV, optical, and near-IR regimes, this work utilizes data from the S-PLUS collaboration (Mendes de Oliveira et al., 2019), as well as VLT and HST observations reported by Fernández-Ontiveros et al. (2009).

Additionally, this study is partially based on archival data obtained with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA.

For radio observations, this work includes 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. This study also made use of (extended) Karl G. Jansky Very Large Array archival observations provided by the National Radio Astronomy Observatory (NRAO), a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

In constructing the SED fitting at 9 $\arcsec$ resolution, this work incorporates observations from the ESA Herschel Space Observatory (Pilbratt et al., 2010), specifically using the PACS instrument (Poglitsch et al., 2010).

References

Allard et al. (2006) Allard, E. L., Knapen, J. H., Peletier, R. F., & Sarzi, M. 2006, MNRAS, 371, 1087
Allard et al. (2005) Allard, E. L., Peletier, R. F., & Knapen, J. H. 2005, ApJ, 633, L25
Álvarez-Márquez et al. (2023) Álvarez-Márquez, J., Crespo Gómez, A., Colina, L., et al. 2023, A&A, 671, A105
Anantharamaiah & Goss (1996) Anantharamaiah, K. R. & Goss, W. M. 1996, ApJ, 466, L13
Ando et al. (2017) Ando, R., Nakanishi, K., Kohno, K., et al. 2017, ApJ, 849, 81
Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22
Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
Bao et al. (2024) Bao, M., Harada, N., Kohno, K., et al. 2024, A&A, 687, A43
Beck et al. (2022) Beck, A., Lebouteiller, V., Madden, S. C., et al. 2022, A&A, 665, A85
Behrens et al. (2022) Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119
Behrens et al. (2024) Behrens, E., Mangum, J. G., Viti, S., et al. 2024, ApJ, 977, 38
Belloche et al. (2013) Belloche, A., Müller, H. S. P., Menten, K. M., Schilke, P., & Comito, C. 2013, A&A, 559, A47
Bendo et al. (2015) Bendo, G. J., Beswick, R. J., D’Cruze, M. J., et al. 2015, MNRAS, 450, L80
Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
Bouvier et al. (2024) Bouvier, M., Viti, S., Behrens, E., et al. 2024, A&A, 689, A64
Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
Brunthaler et al. (2009) Brunthaler, A., Castangia, P., Tarchi, A., et al. 2009, A&A, 497, 103
Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
Bulatek et al. (2023) Bulatek, A., Ginsburg, A., Darling, J., Henkel, C., & Menten, K. M. 2023, ApJ, 956, 78
Butterworth et al. (2024) Butterworth, J., Viti, S., Van der Werf, P. P., et al. 2024, A&A, 686, A31
Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
Charlot & Fall (2000) Charlot, S. & Fall, S. M. 2000, ApJ, 539, 718
Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
Cid Fernandes et al. (2005) Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS, 358, 363
Cid Fernandes et al. (2010) Cid Fernandes, R., Stasińska, G., Schlickmann, M. S., et al. 2010, MNRAS, 403, 1036
Ciesla et al. (2015) Ciesla, L., Charmandaris, V., Georgakakis, A., et al. 2015, A&A, 576, A10
Clemens et al. (2010) Clemens, M. S., Scaife, A., Vega, O., & Bressan, A. 2010, MNRAS, 405, 887
Cohen et al. (2020) Cohen, D. P., Turner, J. L., & Consiglio, S. M. 2020, MNRAS, 493, 627
Condon et al. (1991) Condon, J. J., Huang, Z. P., Yin, Q. F., & Thuan, T. X. 1991, ApJ, 378, 65
Conroy (2013a) Conroy, C. 2013a, ARA&A, 51, 393
Conroy (2013b) Conroy, C. 2013b, ARA&A, 51, 393
da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
Dey et al. (2024) Dey, S., Goyal, A., Małek, K., & Díaz-Santos, T. 2024, ApJ, 966, 61
Dey et al. (2022) Dey, S., Goyal, A., Małek, K., et al. 2022, ApJ, 938, 152
Donevski et al. (2023) Donevski, D., Damjanov, I., Nanni, A., et al. 2023, A&A, 678, A35
Donevski et al. (2020) Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144
Draine et al. (2014) Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172
Ellingsen et al. (2017) Ellingsen, S. P., Chen, X., Breen, S. L., & Qiao, H. H. 2017, MNRAS, 472, 604
Emig et al. (2020) Emig, K. L., Bolatto, A. D., Leroy, A. K., et al. 2020, ApJ, 903, 50
Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
Fernández-Ontiveros et al. (2009) Fernández-Ontiveros, J. A., Prieto, M. A., & Acosta-Pulido, J. A. 2009, MNRAS, 392, L16
Galliano et al. (2018) Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673
Galvin et al. (2018) Galvin, T. J., Seymour, N., Marvil, J., et al. 2018, MNRAS, 474, 779
Gorski et al. (2019) Gorski, M. D., Ott, J., Rand, R., et al. 2019, MNRAS, 483, 5434
Hamed et al. (2021) Hamed, M., Ciesla, L., Béthermin, M., et al. 2021, A&A, 646, A127
Hamed et al. (2023) Hamed, M., Małek, K., Buat, V., et al. 2023, A&A, 674, A99
Harada et al. (2022) Harada, N., Martín, S., Mangum, J. G., et al. 2022, ApJ, 938, 80
Harada et al. (2021) Harada, N., Martín, S., Mangum, J. G., et al. 2021, ApJ, 923, 24
Harada et al. (2024) Harada, N., Meier, D. S., Martín, S., et al. 2024, ApJS, 271, 38
Hernández-Gómez et al. (2019) Hernández-Gómez, A., Sahnoun, E., Caux, E., et al. 2019, MNRAS, 483, 2014
Herpich et al. (2024) Herpich, F. R., Almeida-Fernandes, F., Oliveira Schwarz, G. B., et al. 2024, A&A, 689, A249
Hirashita & Chen (2023) Hirashita, H. & Chen, C.-C. 2023, MNRAS, 526, 4710
Holdship et al. (2022) Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89
Holdship et al. (2021) Holdship, J., Viti, S., Martín, S., et al. 2021, A&A, 654, A55
Hopkins & Beacom (2006) Hopkins, A. M. & Beacom, J. F. 2006, ApJ, 651, 142
Hopkins et al. (2010) Hopkins, P. F., Bundy, K., Croton, D., et al. 2010, ApJ, 715, 202
Huang et al. (2023) Huang, K. Y., Viti, S., Holdship, J., et al. 2023, A&A, 675, A151
Humire et al. (2022) Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33
Humire et al. (2024) Humire, P. K., Ortiz-León, G. N., Hernández-Gómez, A., et al. 2024, A&A, 688, L1
Humire et al. (2020) Humire, P. K., Thiel, V., Henkel, C., et al. 2020, A&A, 642, A222
Iodice et al. (2014) Iodice, E., Arnaboldi, M., Rejkuba, M., et al. 2014, A&A, 567, A86
Kapińska et al. (2017) Kapińska, A. D., Staveley-Smith, L., Crocker, R., et al. 2017, ApJ, 838, 68
Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
Kellermann & Pauliny-Toth (1969) Kellermann, K. I. & Pauliny-Toth, I. I. K. 1969, ApJ, 155, L71
Kennicutt (1998) Kennicutt, Robert C., J. 1998, ARA&A, 36, 189
Kepley et al. (2011) Kepley, A. A., Chomiuk, L., Johnson, K. E., et al. 2011, ApJ, 739, L24
Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
Kewley et al. (2001) Kewley, L. J., Heisler, C. A., Dopita, M. A., & Lumsden, S. 2001, The Astrophysical Journal Supplement Series, 132, 37
Kim et al. (2012) Kim, W.-T., Seo, W.-Y., & Kim, Y. 2012, ApJ, 758, 14
Kobayashi et al. (2020) Kobayashi, C., Karakas, A. I., & Lugaro, M. 2020, ApJ, 900, 179
Kobayashi et al. (2011) Kobayashi, C., Karakas, A. I., & Umeda, H. 2011, MNRAS, 414, 3231
Kokorev et al. (2021) Kokorev, V. I., Magdis, G. E., Davidzon, I., et al. 2021, ApJ, 921, 40
Lapi et al. (2020) Lapi, A., Pantoni, L., Boco, L., & Danese, L. 2020, ApJ, 897, 81
Lapi et al. (2018) Lapi, A., Pantoni, L., Zanisi, L., et al. 2018, ApJ, 857, 22
Leroy et al. (2015) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25
Li & Li (2024) Li, N. & Li, C. 2024, ApJ, 975, 234
López Fernández et al. (2016) López Fernández, R., Cid Fernandes, R., González Delgado, R. M., et al. 2016, Monthly Notices of the Royal Astronomical Society, 458, 184
Mangum et al. (2019) Mangum, J. G., Ginsburg, A. G., Henkel, C., et al. 2019, ApJ, 871, 170
Markov et al. (2024) Markov, V., Gallerani, S., Ferrara, A., et al. 2024, arXiv e-prints, arXiv:2402.05996
Martín et al. (2021) Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46
Martínez-Paredes et al. (2023) Martínez-Paredes, M., Bruzual, G., Morisset, C., et al. 2023, MNRAS, 525, 2916
Meier et al. (2015) Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ, 801, 63
Mendes de Oliveira et al. (2019) Mendes de Oliveira, C., Ribeiro, T., Schoenell, W., et al. 2019, MNRAS, 489, 241
Müller et al. (2005) Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, Journal of Molecular Structure, 742, 215
Müller-Sánchez et al. (2010) 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
Murphy et al. (2011) Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67
Nanni et al. (2020) Nanni, A., Burgarella, D., Theulé, P., Côté, B., & Hirashita, H. 2020, A&A, 641, A168
Narayanan et al. (2024) Narayanan, D., Stark, D. P., Finkelstein, S. L., et al. 2024, arXiv e-prints, arXiv:2408.13312
Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Sausalito, CA: University Science Books)
Palla et al. (2024) Palla, M., De Looze, I., Relaño, M., et al. 2024, MNRAS, 528, 2407
Pence (1980) Pence, W. D. 1980, ApJ, 239, 54
Pérez-Beaupuits et al. (2018) Pérez-Beaupuits, J. P., Güsten, R., Harris, A., et al. 2018, ApJ, 860, 23
Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2
Rekola et al. (2005) Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330
Rico-Villas et al. (2020) Rico-Villas, F., Martín-Pintado, J., González-Alfonso, E., Martín, S., & Rivilla, V. M. 2020, MNRAS, 491, 4573
Romano et al. (2024) Romano, M., Donevski, D., Junais, et al. 2024, A&A, 683, L9
Ronconi et al. (2024) Ronconi, T., Lapi, A., Torsello, M., et al. 2024, A&A, 685, A161
Rosolowsky & Leroy (2006) Rosolowsky, E. & Leroy, A. 2006, PASP, 118, 590
Sakamoto et al. (2006) Sakamoto, K., Ho, P. T. P., Iono, D., et al. 2006, ApJ, 636, 685
Sakamoto et al. (2011) Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, ApJ, 735, 19
Salak et al. (2024) Salak, D., Hashimoto, T., Inoue, A. K., et al. 2024, ApJ, 962, 1
Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
Sawant et al. (2024) Sawant, P., Nanni, A., Romano, M., et al. 2024, arXiv e-prints, arXiv:2412.02505
Schneider & Maiolino (2024) Schneider, R. & Maiolino, R. 2024, A&A Rev., 32, 2
Schreiber et al. (2018) Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85
Schwörer et al. (2019) Schwörer, A., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 628, A6
Scoville & Murchikova (2013) Scoville, N. & Murchikova, L. 2013, ApJ, 779, 75
Sommovigo et al. (2020) Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 497, 956
Stalevski et al. (2012) Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756
Stalevski et al. (2016) Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288
Stasińska et al. (2006) Stasińska, G., Cid Fernandes, R., Mateus, A., Sodré, L., & Asari, N. V. 2006, MNRAS, 371, 972
Stroe et al. (2015) Stroe, A., Oosterloo, T., Röttgering, H. J. A., et al. 2015, MNRAS, 452, 2731
Swinbank et al. (2010) Swinbank, A. M., Smail, I., Longmore, S., et al. 2010, Nature, 464, 733
Tanaka et al. (2024) Tanaka, K., Mangum, J. G., Viti, S., et al. 2024, ApJ, 961, 18
Turner & Ho (1985) Turner, J. L. & Ho, P. T. P. 1985, ApJ, 299, L77
van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
Vastel et al. (2015) Vastel, C., Bottinelli, S., Caux, E., Glorian, J. M., & Boiziot, M. 2015, in SF2A-2015: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 313–316
Walmsley (1989) Walmsley, C. 1989, in Interstellar Dust, ed. L. J. Allamandola & A. G. G. M. Tielens, Vol. 135, 263
Walter et al. (2017) Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265
Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
Werle et al. (2019) Werle, A., Cid Fernandes, R., Vale Asari, N., et al. 2019, MNRAS, 483, 2382
Wirström et al. (2010) Wirström, E. S., Bergman, P., Black, J. H., et al. 2010, A&A, 522, A19
Witstok et al. (2023) Witstok, J., Jones, G. C., Maiolino, R., Smit, R., & Schneider, R. 2023, MNRAS, 523, 3119
Xiang & Rix (2022) Xiang, M. & Rix, H.-W. 2022, Nature, 603, 599
Yan et al. (2019) Yan, Y. T., Zhang, J. S., Henkel, C., et al. 2019, ApJ, 877, 154
Yang et al. (2022) Yang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192
Zeng et al. (2020) Zeng, S., Zhang, Q., Jiménez-Serra, I., et al. 2020, MNRAS, 497, 4896
Zhao et al. (2023) Zhao, J. Y., Zhang, J. S., Wang, Y. X., et al. 2023, ApJS, 266, 29

Appendix A CIGALE SED Modeling with AGN Component

To assess potential AGN contributions, we incorporated the SKIRTOR module (Stalevski et al. 2012, 2016). SKIRTOR is based on the 3D radiative transfer code SKIRT (Baes et al. 2011). It includes obscuration by dusty torus and obscuration by dust settled along with the polar directions. The 2022.1 version of CIGALE offers significant advancements over its predecessors, particularly through its enhanced radio module, which includes both thermal (free-free emission from nebular gas) and nonthermal (synchrotron emission from star formation and AGN activity) contributions. The updated module also incorporates a new AGN component, enabling us to estimate the radio loudness, $R_{\rm{AGN}}$ , defined as the ratio of AGN luminosities measured at 5 GHz and 2500Å. Additionally, it provides the slope of the AGN power-law (PL) radiation, which is assumed to be isotropic (Yang et al. 2022).
CIGALE 2022.1 models the radio luminosities from star formation and AGN components using the FIR/radio correlation coefficient (q_IR) and radio loudness (R_AGN) parameters in scaling relations. At 1.4 GHz, the radio luminosity from star formation is predominantly nonthermal synchrotron emission, normalized using the q_IR parameter. The q_IR determines the normalization at this frequency, while a single power-law describes the synchrotron emission from star formation. Variations in the q_IR parameter contribute to uncertainties in modeling the radio component with CIGALE. The nonthermal synchrotron emission from the AGN is also described using a single PL form governed by the R_AGN parameter. However, the galaxy-averaged radio SEDs of starburst galaxies are rarely represented by single-PL models, showing low-frequency turnover due to free-free absorption (e.g., Galvin et al. 2018; Dey et al. 2022, 2024). It would be interesting to compare the physical conditions responsible for radio emission on the galaxy-wide scale to GMC scales, and that will be the aim of future studies using spatially-resolved radio SEDs of NGC 253. The list of input parameters to CIGALE analysis without and with AGN component is given in Table 11. For the CIGALE modeling including the AGN component, we only added the AGN module to the existing input parameters. These CIGALE inputs have yielded a mean reduced $\chi^{2}$ of 2.2 for the ten GMCs studied in this work. The resulting SEDs are shown in Figs. 14 (without an AGN component) and 15 (considering an AGN component). These Figs. also include, in the bottom middle panel each, the SED fitting centered at GMC 5 at 9 $\arcsec$ aperture diameter, which allows to include Herschel PACS observations (see Appendix E). The results on the main physical properties for each GMC are shown in Table 6, those regarding attenuation effects are also in Table 11 while additional results can be found in Tables 12, and 13.

Table 11: List of input parameters for CIGALE modeling

Parameters		Values
delayed star formation history + additional burst (Ciesla et al. 2015)
e-folding time of the main stellar population model [Myr]	$\tau_{\rm{main}}$	300-5000 by a bin of 300
e-folding time of the late starburst population model [Myr]	$\tau_{\rm{burst}}$	50-500 by a bin of 50
mass fraction of the late burst population	$f_{\rm{burst}}$	0.01, 0.03, 0.05, 0.1, 0.3, 0.6, 0.9
Age of the main stellar population in the galaxy [Myr]	age	1000, 2000, 3000, 4500, 5000, 6500, 10000, 12000
Age of the late burst [Myr]	age_burst	10-350 by a bin of 50
stellar synthesis population (Bruzual & Charlot 2003)
Initial mass function	IMF	(Salpeter 1955)
Metallicity	$Z$	0.02
Separation age		1 Myr
dust attenuation laws (Calzetti et al. 2000)
Colour excess of young stars	E(B-V)	0.1-2 by a bin of 0.2
Reduction factor⁽ⁱⁱⁱ⁾	$f_{\rm{att}}$	0.3, 0.44, 0.6, 0.7
dust emission (Draine et al. 2014)
Mass fraction of PAH	$q_{\rm{PAH}}$	0.47, 1.77, 2.50, 3.19, 4.58, 5.26, 6.63, 7.32
Minimum radiation field	$U_{\rm{min}}$	0.4, 0.8, 1.2, 1.6, 2.0, 2.5, 4, 8, 15
power law index of the radiation	$\alpha$	2
fraction illuminated from $U_{\mathrm{min}}$ to $U_{\mathrm{max}}$	$\gamma$	0.02, 0.06, 0.1, 0.15, 0.2
active nucleus model; Skirtor (Stalevski et al. 2012, 2016)
optical depth at 9.7 $\leavevmode\nobreak\ \mu$ m	$\tau_{9.7}$	3.0, 7.0
torus density radial parameter	pl	1.0
torus density angular parameter	q	1.0
angle between the equatorial plan and edge of the torus [deg]	oa	40
ratio of outer to inner radius	R	20
fraction of total dust mass inside clumps [%]	Mcl	97
inclination (viewing angle) [deg]	i	30 (type 1), 70 (type 2)
AGN fraction		0.0-0.4 by a bin of 0.05
extinction law of polar dust		SMC
E(B-V) of polar dust		0.01-0.7 by a bin of 0.5
Temperature of the polar dust	K	100
Emissivity index of the polar dust		1.6
Radio
SF FIR/radio parameter^∗	$q_{\rm{IR}}$	1.3 - 2.9 by a bin of 0.5
SF Power-law slope (Flux $\propto$ Frequency ${}^{\alpha_{synch}}$ )	$\alpha_{\rm{SF}}$	$-$ 2.0 to $-$ 0.2 by a bin of 0.5
Radio-Loudness parameter^∗∗	$R_{\rm{AGN}}$	0.01, 0.05, 0.1, 0.5… 1000, 5000
AGN Power-law slope	$\alpha_{\rm{AGN}}$	$-$ 2.0 to $-$ 0.2 by a bin of 0.5

^∗ computed as $\log_{10}L_{\rm{IR}(8-1000\mu\rm m)}$ - $\log_{10}L_{1.4\,\rm{GHz}}$ , where $L_{1.4\,\rm{GHz}}$ is the radio luminosity at 1.4 GHz.
^∗∗ defined as a ratio of AGN luminosities measured at 5 GHz and 2500Å.

Table 12: The extinction-related results obtained from SED modeling with CIGALE, including the AGN component.

GMC V_B90 B_B90 E_BV_factor E_BV_lines E_BVs [mag] [mag] [mag] [mag] [mag] 1 $2.87\pm 0.04$ $3.58\pm 0.05$ $0.10\pm 0.00$ $7.00\pm 0.09$ $0.70\pm 0.01$ 2 $3.00\pm 0.17$ $3.70\pm 0.21$ $0.33\pm 0.33$ $5.00\pm 2.80$ $0.73\pm 0.05$ 3 $6.14\pm 0.04$ $6.20\pm 0.02$ $0.40\pm 0.00$ $3.00\pm 0.00$ $1.20\pm 0.00$ 4 $6.10\pm 0.01$ $7.60\pm 0.01$ $0.45\pm 0.09$ $3.50\pm 0.90$ $1.50\pm 0.00$ 5 $6.10\pm 0.01$ $7.60\pm 0.01$ $0.45\pm 0.09$ $3.50\pm 0.90$ $1.50\pm 0.00$ 6 $4.90\pm 0.01$ $6.10\pm 0.01$ $0.40\pm 0.00$ $3.00\pm 0.00$ $1.20\pm 0.00$ 7 $3.20\pm 0.01$ $4.00\pm 0.01$ $0.80\pm 0.01$ $1.00\pm 0.05$ $0.80\pm 0.01$ 8 $3.00\pm 0.17$ $3.70\pm 0.21$ $0.29\pm 0.31$ $5.40\pm 2.70$ $0.73\pm 0.04$ 9 $2.80\pm 0.11$ $3.50\pm 0.14$ $0.10\pm 0.02$ $7.00\pm 0.34$ $0.70\pm 0.03$ 10 $2.90\pm 0.03$ $3.60\pm 0.04$ $0.10\pm 0.05$ $7.00\pm 0.41$ $0.70\pm 0.01$ ¹¹¹¹11The columns represent the GMC number, $A_{\rm{V}}^{\rm{Balmer}}$ as presented in Table 7, the slope of the attenuation curve, the scaling factor between $E(B-V)$ for lines and continuum, the color excess for nebular emission lines, and the color excess for the stellar continuum.

Table 13: Dust and SFH related results obtained from SED modeling with CIGALE, including AGN component.

GMC $M_{\rm dust}$ U_min $\gamma$ $q_{\rm PAH}$ age_burst $f_{\rm burst}$ $\tau_{\rm burst}$ $\tau_{\rm main}$ [log₁₀ ( $M_{\odot}$ )] [Myr] [Myr] [Myr] 1 $4.47\pm 0.09$ $4.39\pm 1.51$ $0.025\pm 0.033$ $2.236\pm 0.590$ $7.87\pm 99.45$ $0.03\pm 0.02$ $142.23\pm 85.13$ $1170.23\pm 645.23$ 2 $4.49\pm 0.09$ $3.06\pm 1.84$ $0.057\pm 0.067$ $4.821\pm 1.218$ $8.10\pm 32.48$ $0.02\pm 0.01$ $210.26\pm 75.10$ $290.93\pm 141.61$ 3 $5.89\pm 0.07$ $3.25\pm 0.72$ $0.011\pm 0.010$ $1.118\pm 0.036$ $7.74\pm 98.75$ $0.03\pm 0.02$ $126.62\pm 81.91$ $864.58\pm 527.05$ 4 $6.19\pm 0.04$ $0.68\pm 0.16$ $0.089\pm 0.038$ $0.470\pm 0.001$ $8.16\pm 8.01$ $0.01\pm 0.00$ $246.94\pm 60.90$ $203.13\pm 30.50$ 5 $6.23\pm 0.09$ $0.90\pm 0.39$ $0.078\pm 0.025$ $0.471\pm 0.026$ $7.92\pm 0.01$ $0.01\pm 0.02$ $180.44\pm 84.63$ $213.33\pm 61.83$ 6 $6.08\pm 0.06$ $0.58\pm 0.34$ $0.189\pm 0.033$ $1.036\pm 0.312$ $7.93\pm 101.79$ $0.03\pm 0.02$ $170.89\pm 86.57$ $1520.01\pm 683.09$ 7 $5.13\pm 0.11$ $9.03\pm 1.47$ $0.192\pm 0.029$ $3.116\pm 0.871$ $7.68\pm 93.06$ $0.03\pm 0.02$ $134.24\pm 83.65$ $844.18\pm 544.03$ 8 $4.79\pm 0.21$ $2.65\pm 1.02$ $0.160\pm 0.053$ $3.384\pm 1.132$ $7.75\pm 88.29$ $0.02\pm 0.01$ $176.04\pm 86.66$ $258.82\pm 125.10$ 9 $4.73\pm 0.23$ $5.37\pm 1.72$ $0.012\pm 0.014$ $3.105\pm 0.838$ $7.71\pm 93.91$ $0.03\pm 0.02$ $137.08\pm 83.70$ $592.04\pm 305.75$ 10 $4.83\pm 0.19$ $1.87\pm 0.71$ $0.051\pm 0.042$ $3.412\pm 0.918$ $7.82\pm 103.06$ $0.02\pm 0.02$ $170.82\pm 86.23$ $589.30\pm 352.26$

¹²¹²12The columns represent the GMC number, dust mass, minimum radiation field, mass fraction of PAH, age of late burst, mass fraction od late burst population, e-folding time of the late starburst population model, e-folding time of the main stellar population model

Table 14: Results of AGN and radio SED modeling using CIGALE.

GMC $f_{\rm{AGN}}$ accretion power $L_{\rm disk}$ polar dust luminosity $R_{\rm AGN}$ $\alpha_{\rm AGN}$ $\alpha_{\rm SF}$ $q_{\rm IR}$ [%] [log₁₀ erg s^-1] [log₁₀ erg s^-1] [log₁₀ erg s^-1] 1 $2.52\pm 1.37$ $32.70\pm 0.24$ $30.40\pm 0.24$ $31.97\pm 0.24$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.56\pm 0.11$ $2.20\pm 0.03$ 2 $3.86\pm 2.56$ $33.01\pm 0.30$ $30.70\pm 0.31$ $32.28\pm 0.30$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.59\pm 0.12$ $2.20\pm 0.00$ 3 $2.00\pm 0.01$ $34.01\pm 0.02$ $31.56\pm 0.10$ $33.28\pm 0.02$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.51\pm 0.04$ $2.21\pm 0.05$ 4 $2.05\pm 0.46$ $34.48\pm 0.10$ $32.01\pm 0.10$ $33.75\pm 0.10$ $0.01\pm 0.01$ $0.70\pm 0.00$ $0.50\pm 0.00$ $2.20\pm 0.00$ 5 $5.34\pm 1.48$ $34.87\pm 0.12$ $32.40\pm 0.12$ $34.13\pm 0.12$ $0.01\pm 0.01$ $0.70\pm 0.00$ $0.50\pm 0.00$ $2.20\pm 0.00$ 6 $4.66\pm 2.68$ $34.31\pm 0.26$ $31.92\pm 0.24$ $33.57\pm 0.26$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.50\pm 0.00$ $2.20\pm 0.01$ 7 $7.54\pm 2.85$ $33.79\pm 0.18$ $31.47\pm 0.25$ $33.06\pm 0.18$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.69\pm 0.11$ $2.64\pm 0.08$ 8 $2.56\pm 1.54$ $33.18\pm 0.24$ $30.85\pm 0.25$ $32.45\pm 0.24$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.73\pm 0.12$ $2.22\pm 0.09$ 9 $2.17\pm 0.93$ $32.95\pm 0.14$ $30.62\pm 0.22$ $32.22\pm 0.14$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.79\pm 0.11$ $2.20\pm 0.01$ 10 $2.51\pm 1.42$ $32.97\pm 0.26$ $30.64\pm 0.25$ $32.23\pm 0.26$ $0.04\pm 0.04$ $0.70\pm 0.00$ $0.50\pm 0.00$ $2.58\pm 0.01$ ¹³¹³13The columns represents the GMC number, AGN fraction, accretion power, disk luminosity and polar dust luminosity of AGN, radio loudness, AGN power-law slope, SF power-law slope, and radio-IR correlation par.

Appendix B Individual GMC properties

GMC 1: Located at the boundaries of the CMZ, where x₁/x₂ interactions are present, GMC 1 traces weak shocks at low-frequency transitions (Humire et al. 2022; Harada et al. 2022; Huang et al. 2023) and strong shocks at high- $J$ frequency transitions. Indeed, the strongest SiO(5–4)/HNCO(10 ${}_{0,10}-9_{0,9}$ ) ratios are observed at this position (Humire et al. 2022). The whole picture may also be that weak and strong shocks are actually mixed over the entire CMZ (Bao et al. 2024). The shocked environment results in some of the more prominent Class I methanol masers (MMCIs) in the $J_{-1}\rightarrow(J-$ 1) ${}_{0}-E$ and $J_{0}\rightarrow(J-$ 1) ${}_{1}-A$ line series (Ellingsen et al. 2017; Gorski et al. 2019; Humire et al. 2022). Devoid of a strong continuum, but showing the strongest maser over LTE emission (Humire et al. 2022, their Fig. 13), GMC 1 appears to be an extragalactic example of what is seen in the Milky Way object GMC G $+$ 0.693 (Zeng et al. 2020). In line with the latter, this cloud possesses the lowest amount of dust and stellar masses in the sample, according to GalaPy (see Table 9).

Contrary to the others GMCs in our sample, GMC 1 is the only one whose SED fitting could also deliver ages of the order of 10⁵ years. However, since we observe CO emission, a molecule that needs a few million years to be produced (e.g., Walmsley 1989), a permitted range of 10^6-11 or 10^6-12 years was imposed to all GMCs. This restriction yielded SED shapes that are closer to the observed ones compared to those with an age of $\sim$ 10⁵ years, and significantly improved the results at millimeter and centimeter wavelengths by better reproducing the observed flux densities.

GMC 2 This molecular cloud also overlaps with x₁/x₂ interactions and presents physical and chemical conditions that are similar to GMC 1, although GMC 2 is more massive in terms of stellar and dust mass (Table 9; see also Humire et al. 2022 and Bouvier et al. 2024), older (Table 9), and shows weaker MMCIs exclusively in the $J_{-1}\rightarrow(J-$ 1) ${}_{0}-E$ family (Humire et al. 2022). With slightly lower values than GMC 1, GMc 2 presents the lowest dusty radius, as calculated by GalaPy, and the lowest stellar mass, according to starlight, as indicated in Tables 9 and 10, respectively.

GMC 3 Appears to be the source of the stellar outflow best viewed by high-density tracers such as CN, HCN (Walter et al. 2017), or an increment in the CO/¹³CO $J=$ 1–0 ratios owed to a more transparent environment (Bao et al. 2024). Although dense molecular gas possess a velocity dispersion of $\sim$ 40 km s^-1 in the streamer (Walter et al. 2017), the stellar velocity dispersion is the lowest among the sample (see Table 10) at the GMC position. This corresponds to the depiction of cold/low-dispersion gas actively producing stars at the convergence of the ring and the bar, as observed in M 100 (Allard et al. 2005). A more evolved and therefore, more easily observed stellar stream is in line with this position being the oldest among the nuclear regions, as deduced from our SED fitting in Table 9. In the same table, and also in Figure5, it can be seen than GMC 3 presents the largest dust reservoir and extension.

GMC 4 exhibits the highest SiO(2–1)/HNCO(4 ${}_{0,4}-3_{0,3}$ ) ratio in the sample (Meier et al. 2015; Humire et al. 2022), indicative of strong shocks. It is also the second youngest molecular cloud in the sample according to GalaPy and the youngest according to starlight fits, respectively (see Table 5 and Figs. 6 and 19). From CIGALE and GalaPy, it also possess the second largest SFR. This intense stellar formation, despite the absence of a distinct molecular stellar streamer, may be a consequence of its youthful nature or simply that radio emission is not strong enough to trace it (although the stellar streamers originating from the central GMCs is clearly seen in H $\alpha$ ; see second to last upper panel of Fig. 1). Additionally, and from GalaPy’s results, GMC 4 has the largest $M_{\rm{dust}}/M_{\rm{gas}}$ ratio in the sample ( $\log_{10}\left(\frac{M_{\rm{dust}}}{M_{\rm{gas}}}\right)=-1.730$ ; see Table 9), which could suggest that material is being rapidly funneled toward the nucleus via the bar. Similarly, GMCs 5 and 6 are also among the youngest in the sample, with the second and third largest dust-to-gas mass ratios ( $\log_{10}\left(\frac{M_{\rm{dust}}}{M_{\rm{gas}}}\right)=-1.820$ and -1.988, respectively). A possible explanation is that material is being concentrated in these GMCs due to inner-inner Lindblad resonances caused by the nuclear bar (Anantharamaiah & Goss 1996; Cohen et al. 2020). Indeed, looking at the left panel of Fig. 11 in Humire et al. (2022), and after considering our 3 $\arcsec$ photometric aperture and modified GMC positions, we can clearly see a peak in the SiO(2–1)/HNCO(4 ${}_{0,4}-3_{0,3}$ ) ratios around GMC 5, this latter showing a gap due to SiO(2–1) self-absorption.

Our mid-IR data from NACO and VISIR at the VLT show the highest fluxes in the sample at this position (see Fig. 3). GMCs 2 and 4 are the only giant molecular clouds where 10–36% of the stellar mass was produced in the last $\sim$ 2 million years. In contrast, the bulk of stellar mass in the remaining GMCs was produced approximately 1 to 10 billion years ago (see Fig. 6).

GMC 5 Among our sample, GMC 5 has the best signal across the different wavelengths, being the most clearly observed at centimeter wavelengths, with its emission covering almost the entire aperture. It is closest to the strongest radio continuum source, TH2 (Turner & Ho 1985), which also produces some (self-)absorption features in methanol (Humire et al. 2022), in contrast to the outermost GMCs, where methanol absorption is primarily attributed to anti-inversion against the cosmic microwave background (CMB; Bulatek et al. 2023). The strong continuum emission also leads to self-absorption in other molecules, such as SiO $J$ =2–1, as inferred from Humire et al. (2022, their Fig. 11, left panel), although this has not been investigated in detail. According to GalaPy and CIGALE results, GMC 5 exhibits by far the largest instantaneous and time-averaged (over 10 and 100 Myr) SFRs in our sample ( $\sim$ 0.7 and 0.2 $M_{\odot}$ yr^-1, according to GalaPy and CIGALE, respectively; see Tables 5 and 6), creating about three times more stars than GMC 4, which shows the closest characteristics among the ten GMCs studied.

GMC 6 Most species peak at this location, including CO isotopologues (Harada et al. 2024). Clear exceptions are the sub-mm continuum emissions at Bands 3-7, RRLs, CN, HCN, and CO ( $J=$ 2–1, and $J=$ 3–2), whose emission is strongest in GMC 5. It shows the highest column densities in both methanol symmetric types and presents the second highest SiO(2–1)/HNCO(4 ${}_{0,4}-3_{0,3}$ ) ratios after GMC 4, indicative of a strongly perturbed environment due to the presence of strong shocks (Meier et al. 2015; Humire et al. 2022; Huang et al. 2023).

GMC 7 GMC 7 hosts the strongest methanol emission at 84.53 and 132.89 GHz, which is enhanced by maser emission. It also presents the strongest difference between E- and A-CH₃OH from a single LTE component analysis, with a E/A ratio of 3.16 ${}^{+0.17}_{-0.18}$ , which can only be related to shocks since theoretical values are not expected to surpass 1.0, unless methanol is not formed from CO on cold grain surfaces in NGC 253 (Wirström et al. 2010). We note that E/A ratios of up to 6.0 have been seen in massive star-forming regions in our galaxy (see, e.g., Zhao et al. 2023, their Fig. 10(a)).

GMC 8 Using a two-component LTE analysis, this GMC shows one of the largest E/A-CH₃OH ratios in the cold layer (3.14 $\pm$ 0.58). This may be due to low temperatures, where E-CH₃OH is more populated as their upper level energies ( $E_{\rm{up}}$ ) are in general lower than its A-CH₃OH counterparts (Humire et al. 2022, 2024). The above is in line with cold dust emission increments, as can be inferred from enhanced 100 GHz continuum emission, in contrast to its spurious signal in GMCs 1 and 2, as accounted for in our Fig. 2. A similar conclusion can also be obtained from the distribution of sulfurated molecules seen in ALMA Band 3 (84–116 GHz), whose emission is also enhanced toward the northeast segment of NGC 253’s CMZ (Meier et al. 2015).

GMC 9 Presents a mild stellar velocity dispersion ( $\sigma_{\star}$ ; see Table 10) that may account for a relative low formation of stars (Allard et al. 2006) and be a result of its position in between the spiral arms and the inner bar, where (x₁/x₂) interactions and inner Lindblad resonances take place (see dashed gray ellipses in Fig. 2). A similar origin for the velocity dispersion can be inferred for GMCs 1,, 2, and 8–9, in contrast to the larger $\sigma_{\star}$ measured at GMCs 3, 5, 6, and 7, that are possibly due to a vigorous stellar formation instead (Bouvier et al. 2024; Bao et al. 2024).

GMC 10 we do not see much high-excitation methanol lines in this region, that shows the lowest excitation temperatures and number of methanol lines, in addition to the lowest column density and line-widths (FWHM $\sim$ 23 km s^-1) of both methanol symmetric types among the ten GMCs studied in this work. Using high density sulfured molecules, Bouvier et al. (2024) also determined the lowest line-widths (FWHM $\sim$ 20 km s^-1) and rest velocities (100–150 km s^-1) among the GMCs in the sample.

Appendix C Sulfur ratios and Age-metallicity relation in GMC 5

While many simple abundant molecules are prone to undergo strong opacity effects, impeding a confident derivation of line ratios, isotopologues of less abundant elements are normally optically thin. Moreover, line ratios can be used to derive metallicities such as [Fe/H], defined concerning the Solar metallicity by the following relation:

[\rm{Fe/H}]=log_{10}([N_{\rm{Fe}/N_{\rm{H}}}])_{\rm{star}}-log_{10}([N_{Fe}/N_{\rm{H}}])_{\rm{Sun}}.

(6)

One way to obtain [Fe/H] through optically thin molecules is using the relation in Table 3 of Kobayashi et al. (2011), well defined across the -0.5 $\leq$ [Fe/H] $\leq$ 0.0 range:

\frac{\it{N}\rm{(^{32}S})}{\it{N}\rm{(^{34}S)}}=-19.8\times[\rm{Fe/H}]+23.2.

(7)

This equation relates the column densities, $N$ , of sulfur isotopologues with the metallicities. With the purpose of obtaining the required column density ratio, we employed the ¹³CS and ¹³C³⁴S isotopologues, of which ¹³C³⁴S is only detectable in GMC 5 from our sub-mm data. Within the ALCHEMI frequency coverage (84.4–373.3 GHz), ¹³C³⁴S is always close but not fully blended with HC₃N (see Fig. 16, bottom panels). Its transitions at $J_{u}=$ 2 to 8, when present, fall $\sim$ 160 km s^-1 away (redshifed) from HC₃N at $J_{u}=$ 2 $x$ to 8 $x$ , with $x$ being equal to 5 and $J_{u}$ being always the unity above the lower $J$ ; for example $J_{u}=4$ being equal to the $J=$ 4–3 transition.

Assuming LTE conditions, we fitted the observed spectra with a single component to account for the very dense regions that CS isotopologues trace. To this end, we employed the Centre d’Analyse Scientifique de Spectres Instrumentaux et Synthétiques (CASSIS), a software developed by CESR/IRAP¹⁴¹⁴14http://cassis.irap.omp.eu/?page=cassis (see, e.g., Vastel et al. 2015) and designed for the analysis and synthesis of molecular spectra in astrophysical research. It enables spectral modeling, line identification, and data interpretation for astronomical observations, with the capability of performing both local thermodynamic equilibrium (LTE) and non-LTE modeling. When performing non-LTE modeling, CASSIS functions as a wrapper for the RADEX code (van der Tak et al. 2007), facilitating its utilization for users. CASSIS is not restricted to the radio regime (e.g., Hernández-Gómez et al. 2019) as it has also been used to distinguish between LTE and non-LTE emission, enabling the discovery of new maser transitions (Humire et al. 2024). While originally designed to process single spectra, the code can be adapted to handle data cubes (Harada et al. 2022)

To produce the synthetic spectrum that better reproduces the CS-related lines, we have utilized the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005) for spectroscopic data and restricted to LTE conditions, since these molecules are expected to trace dense gas. Our best-results are shown in Fig. 16.

This way, we obtained different solutions ranging from ³²S/³⁴S rates much lower than expected in the Milky Way ( $\sim$ 5) to more reasonable values ( $\sim$ 17). We also find that the line intensities can be reproduced under a large range of excitation temperatures. Due to this degeneracy, we decided to limit the ³²S/³⁴S range to values fulfilling the criterion: [Fe/H] $\leq$ 0.5, motivated by the metallicity values obtained from a sample of $\sim$ 247,000 sub-giant stars in the Milky Way (see Xiang & Rix (2022) and Fig. 17). This way, we expect ³²S/³⁴S isotopologue ratios no lower than 13.3, namely, within the range of values measured in most of the galaxy (see Humire et al. 2020, their Fig. 5). Considering that Sgr B2(N) shows properties comparable to a mini-starburst system (Belloche et al. 2013; Schwörer et al. 2019), and that the best ³²S/³⁴S column density ratio measured there is 16.3 ${}^{+2.1}_{-1.7}$ (Humire et al. 2020), our assumption of ³²S/³⁴S $\geq$ 13.3 is well justified.

Our best-fit LTE modeling parameters from the spectra of ¹³CS and ¹³C³⁴S are presented in Table 15. The resulting ³²S/³⁴S column density ratio is 14.26 ${}^{7.53}_{-4.72}$ , and we set it to be 14.26 ${}^{7.53}_{-0.96}$ considering our lowest threshold. Using Eq. 7, the resulting [Fe/H] metallicity values are 0.45 ${}^{+0.05}_{-0.38}$ , namely, between 0.07 and 0.5. Comparing this latter result with the age-metallicity relation presented in Fig. 2a of Xiang & Rix (2022), reproduced in the left panel of Fig. 17, we can estimate an age range. To this aim, we utilized a skewed Gaussian distribution model due to the uni-modal and asymmetric distribution of ages within our estimated metallicity range, as can be seen in the right panel of Fig. 17. The fitting process involved optimizing the parameters of the skewed Gaussian function, including amplitude $A$ , mean $\mu$ , standard deviation $\sigma$ , and skewness $\alpha$ , using a least-squares approach. The skewed Gaussian distribution function is defined as:

f(x)=\frac{A}{\sigma\sqrt{\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}\left(1+\text{erf}\left(\alpha\frac{x-\mu}{\sigma\sqrt{2}}\right)\right),

(8)

where the error function, $\text{erf}(z)$ , is defined as:

\text{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}\,dt.

(9)

Our analysis revealed a negative skewness parameter $\alpha$ , indicating a non-symmetric age distribution with a higher density toward young ( $\sim$ 2.5–7.5 Gyr) sub-giant stars, withing the age-metallicity relation distribution used here (Fig. 17). The resulting estimated age, namely, the distribution peak and its FWHM, is 7.85 ${}^{+0.18}_{-5.39}$ Gyrs, reflecting a stellar population 31 times older than the one traced by our SED fitting in GMC 5, summarized in Table 9, whose value is of 10^8.40 or $\sim$ 251 Myr.

Table 15: Best fit parameters from our LTE plus non-LTE modeling for CARMA-7.

Species	$N$ (Sp)	$T_{\rm ex}$	FWHM	$V_{\rm LSR}$	Size
	[ $\times$ 10¹⁴cm^-2]	[K]	[km s^-1]	[km s^-1]	[ $\arcsec$ ]
¹³CS	29.67 ${}^{+4.12}_{-9.63}$	36.17 ${}^{+46.95}_{-10.14}$	83.87 ${}^{+86.64}_{-22.81}$	198.89 ${}^{+39.20}_{-12.99}$	1.08 ${}^{+1.91}_{-0.25}$
¹³C³⁴S	2.08 ${}^{+0.02}_{-0.53}$	60.97 ${}^{+101.06}_{-23.88}$	39.83 ${}^{+0.16}_{-14.66}$	195.00 ${}^{+1.00}_{-1.00}$	2.93 ${}^{+0.07}_{-1.96}$

¹⁵¹⁵15Best-fitting parameters obtained through our LTE modeling of sulfur isotopologues. From left to right columns correspond to: Name of the species,

N

(Sp) stands for column density of species,

T_{\rm{ex}}

is the excitation temperature, the line-width measured from its FWHM, the local standard of rest velocity,

V_{\rm{LSR}}

, and the best-fit source size in arcseconds. Uncertainties correspond to a 3

\sigma

level.

Although magnetohydrodynamic models for sulfur may not yet be well constrained for the bulge of galaxies, the most likely explanation for sulfur ratios in the very center of our Milky Way is that its large-scale bar is funneling material from its edges toward the galaxy’s nucleus or bulge (Humire et al. 2020). Following this reasoning, the ³²S/³⁴S abundances can still be modeled according to Kobayashi et al. (2011) (C. Kobayashi, priv. comm.), by assuming that the origin of ³²S/³⁴S abundances occurs at a certain level in the disk of galaxies.

In a more recent paper, Kobayashi et al. (2020) presented an updated version of their age-metallicity relation (their Fig. 2b), which, as the SFH models used in GalaPy and CIGALE, also considers a single burst in its SFH model. This updated model integrates a bulge outflow component that is more suited for the NGC 253 central molecular zone (CMZ), where stars drive material outside the bar through significant stellar outflows, likely originating from supernovae. Applying the 0.0–0.5 iron metallicity range to the bulge outflow model yields ages starting from one Gyr, while the bulge component alone provides a more restricted range of one to five Gyr, peaking at 3 Gyr (or 10^9.48 years). Both results are in good agreement with the AMR-derived estimations.

The discrepancy arises because metallicities measured by molecular emitters provide a broader perspective, reflecting the average contributions of all stellar populations that have expelled material to the interstellar medium (ISM). In contrast, values derived from SED fitting are predominantly influenced by recent star formation, which emits most of its radiation in the mid-infrared (MIR; 1.5–50 $\mu$ m). Consequently, we anticipate deriving an average age estimate from the sub-millimeter regime and a younger age estimate from SED fitting. Additionally, employing molecular tracers at higher vibrational modes ( $v=1$ or 2) and angular resolutions allows for tracing the youngest stellar population nicely tracing the last burst in the light-weighted SFHs of Fig. 6, representing up to 43% of the total stellar light emitted in the optical regime (see dashed blue vertical lines in Fig. 6) (Rico-Villas et al. 2020; Butterworth et al. 2024).

However, it might be that the MIR emission we account for in the SED fitting mainly originates from low-mass/old stellar populations. If that is the case, we do not have a clear answer for this discrepancy.

Appendix D Comparison among different methodologies

Given the different SED modeling approaches, that differ in terms of star formation histories and several other assumptions, here we provide a few comparison plots. Depending directly on the SED shape, we will obtain a sSFR (Conroy 2013b). Then, each model assumes a certain mass-to-light ratio to derive the stellar mass and a specific SFH to derive the stellar ages, combining both results, mass and time, to obtain the SFR. Given the variation between GalaPy and CIGALE modeled SED shapes, as seen in Fig. 3, we expect dispersion in their resulting sSFRs, but a relative tendency to be a one-to-one relation. Such a tendency is seen in the bottom right panel of our Fig. 18, meaning that, at this level of comparison, both codes share similarities. From that same Figure it can be noted that GalaPy tends to produce larger SFRs and stellar masses, as also previously noted by Ronconi et al. (2024). The dust mass is mostly consistent between the GalaPy and CIGALE models within 1 $\sigma$ uncertainties for the less massive GMCs but systematically shifts toward the lower values of CIGALE compared to the higher values of GalaPy, as evidenced by the 1:1 red line in the bottom left panel of Fig. 18. This trend becomes more pronounced in the more nuclear and massive regions, specifically GMCs 4–6 and also GMC 7.

In Fig. 19, we compare the stellar population ages, burst ages, SF characteristic timescales, and ages weighted by light and mass, as provided by the different models utilized in this work. Specifically, CIGALE provides the main stellar population age, the late starburst age, the burst age, and the average age of the young and old stellar populations. From starlight, we include the averaged stellar ages weighted by light and by mass, which primarily correspond to the young and old stellar populations, respectively. Additionally, we present the stellar ages derived by GalaPy, along with the characteristic star-formation timescale, $\tau_{\star}$ , obtained from the same code.

Appendix E Testing the FIR bump accuracy

We have performed an SED fitting with both GalaPy and CIGALE considering a 9″ aperture diameter to include Herschel PACS data in the far infrared. This approach ensures that the infrared bump more accentuated in GalaPy than in CIGALE models (see Fig. 3 for a clear comparison) at around 100 $\mu$ m is properly fit without extrapolations. Additionally, it allows us to test whether our extrapolation was consistent with real observations, which cannot be verified at 3 $\arcsec$ due to instrumental limitations. We found a similar dust temperature, which is a result driven by the FIR SED peak, in both cases, at 3 and 9 $\arcsec$ , as already pointed out in Sect. 5.1. Overall, the SED model used in GalaPy for GMC 5 correlates well with the one derived for the larger 9″ aperture centered at the same GMC.

In contrast, while CIGALE models underpredict the FIR bump in several GMCs compared to GalaPy, incorporating Herschel photometric points at 70, 100, and 160 $\mu$ m enables the code to fit all SED points accurately, regardless of whether an AGN component is assumed. The 9 $\arcsec$ SED fitting using CIGALE is shown in the bottom middle panel of Figs. 14 and 15 for the cases without and with an AGN component, respectively. Interestingly, the SFRs differ by less than 5% between these models, reaffirming that the presence of an AGN is irrelevant from the SED perspective, as previously discussed in Sect. 4.2.1.

Finally, the 3 $\arcsec$ and 9 $\arcsec$ SED models from GalaPy are shown in the right panel of Fig. 20. The 3 $\arcsec$ aperture SED corresponds to what was already presented in Fig. 3 for GMC 5.