© The Authors 2026

1 Introduction

Astrochemistry is inherently a multidisciplinary field that involves physics, chemistry, biology, and geology to investigate the formation and evolution of molecules in the interstellar medium (ISM). Unveiling the chemistry in a particular environment provides crucial insight into the emergence of molecular complexity, the dominant chemical processes, and how these are inherited at different evolutionary stages. Furthermore, astro-chemical studies provide constraints on the physical conditions and lifetimes of these regions (Jørgensen et al. 2020). At its core, astrochemistry studies the potential precursors of life, and together with astrobiology, addresses one of the most fundamental open questions in science: the origin of life.

As of November 2025, approximately 340 molecules have been detected in the ISM or in circumstellar shells (Müller 2025). They include the so-called complex organic molecules (COMs), which are defined as carbon-bearing species with at least six atoms (Herbst & van Dishoeck 2009). COMs play a central role as they represent key steps towards increasing chemical complexity, potentially leading to the formation of prebiotic species. The latter are regarded as building blocks of life, as they share structural elements with biomolecules found in living organisms and are thought to play a role in prebiotic chemical processes (Caselli & Ceccarelli 2012; Ceccarelli et al. 2022).

Although the formation of COMs remains under debate, two main and potentially complementary pathways have been proposed: (i) surface reactions on interstellar dust grains (Simons et al. 2020; Garrod et al. 2022), and (ii) gas-phase chemical reactions (Balucani et al. 2015; Skouteris et al. 2018). During the early stages of star formation, cold molecular clouds are characterised by temperatures of ∼10–20 K and densities of ∼10³–10⁴ cm⁻³, which increase as the cloud evolves and collapses. With these conditions, atoms and simple species such as oxygen, nitrogen, and hydrogen accrete onto dust grains, forming icy mantles. In this scenario, grain-surface chemistry is dominated by hydrogenation reactions, leading to a chemical inventory dominated by simple species such as water (H₂O), methane (CH₄), or ammonia (NH₃; Tielens 2013; Boogert et al. 2015; Linnartz et al. 2015). As collapse proceeds and densities reach ∼10⁴–10⁶ cm⁻³, significant CO freeze-out occurs, enabling CO hydrogenation and the formation of species such as the formyl radical (HCO), formaldehyde (H₂CO), and eventually, methanol (CH₃OH), the simplest and most abundant COM (Watanabe & Kouchi 2002; Fuchs et al. 2009; Santos et al. 2022). Moreover, ultraviolet photons and cosmic rays enable the gradual build-up of molecular complexity (Öberg 2016; Padovani et al. 2020). As material accretes onto the protostar, the resulting accretion shocks significantly increase the central luminosity and raise the envelope temperature from 10 K to 100–300 K. This thermal processing triggers the sublimation of icy mantles from dust grains, releasing species into the gas phase (Caselli & Ceccarelli 2012). When they are in the gas phase, COMs can be detected with radio telescopes via their characteristic rotational transitions.

The COMs are commonly observed in high-mass star-forming regions at different evolutionary stages due to high gas densities, intense heating occurring during the hot molecular core (HMC) phase, and efficient desorption processes that enhance their detectability and abundance (Baek et al. 2022; Coletta et al. 2020). HMC regions have the richest chemistry in the Galaxy (Bisschop et al. 2007; Belloche et al. 2013; Rivilla et al. 2017). In addition to tracing chemical complexity, COMs probe a wide range of physical conditions (temperature and density) and provide valuable insights into the physical properties, kinematics, and formation mechanisms of high-mass star-forming regions. In this context, heavy COMs (those with ≥ 8 atoms) and prebiotic species are typically less abundant and more optically thin than standard HMC tracers such as CH₃CN, allowing them to better probe the innermost densest regions of high-mass cores (Rivilla et al. 2017).

However, it remains challenging to detect COMs because their emission is generally weaker than that of simpler and more abundant species. In chemically rich regions, spectra can become so densely populated with molecular transitions that individual lines overlap and cannot be uniquely identified, reaching the line-confusion limit, where increased sensitivity no longer enables secure detections of weaker species. This effect is particularly severe in turbulent environments. Consequently, a careful and detailed line identification is required to reliably assign molecular transitions and derive the physical properties of less abundant species.

Motivated by these considerations, we performed an exhaustive line identification towards core 9 of the AG318.9477−0.1960 clump (hereafter AG318−c9; Fig. 1), which is one of the most chemically rich cores in the ALMA evolutionary study of high-mass protocluster formation in the Galaxy (ALMA-GAL) large project (2019.1.00195.L; Molinari et al. 2025; Sánchez-Monge et al. 2025). ALMAGAL is a survey of high-mass star-forming regions in the Galactic plane that are observed with the Atacama Large Millimetre Array (hereafter ALMA; Wootten & Thompson 2009) in Band 6 (Kerr et al. 2004). The survey was selected from the Herschel Infrared Galactic Plane Survey (Hi-GAL, Molinari et al. 2016; Elia et al. 2018, 2021) and Red MSX Source survey (RMS; Hoare et al. 2005; Urquhart et al. 2007; Lumsden et al. 2013) catalogues. By observing 1013 dense clumps in various evolutionary stages and Galactic environments, ALMAGAL provides a statistically significant and uniformly observed sample to investigate the physical and chemical evolution of high-mass star-forming regions. Detailed descriptions of the ALMAGAL survey characterisation (Molinari et al. 2025), dataset and data reduction pipeline (Sánchez-Monge et al. 2025), fragmentation properties of the compact sources (Coletta et al. 2025), and early scientific results (Wells et al. 2024; Mininni et al. 2025; Elia et al. 2026; Schisano et al. 2026; Wallace et al. 2026) are available.

In a forthcoming study (Allande et al., in prep.), we will analyse the emission of three potentially chemically related COMs across the full ALMAGAL sample: GA, EG, and MF. Laboratory experiments and astrochemical models indicate that GA and EG are chemically linked and may share common formation pathways in the ISM (Sorrell 2001; Fedoseev et al. 2015; Coutens et al. 2018; Garrod et al. 2022). MF, while not directly connected through the same synthetic routes, is the most abundant isomer of GA and provides complementary constraints on the formation of GA and EG, as all three species are expected to arise under similar physical conditions (Simons et al. 2020).

As a necessary first step towards this survey-wide analysis, we present a detailed line identification of the ALMAGAL spectral band coverage (216.98–218.86 GHz and 219.06–220.94 GHz) towards AG318−c9. We identify all detectable molecular species contributing to the emission in this frequency range and perform an initial analysis of the physical properties of GA, EG, and MF, with particular emphasis on assessing line blending and selecting the least affected transitions for further study. In addition, we explore the feasibility of deriving spatially resolved physical parameter maps (column density and excitation temperature) for these species, evaluating their potential as tracers of the physical conditions of dense gas (Rivilla et al. 2017; Nazari 2026).

This study therefore serves multiple purposes within the framework of the ALMAGAL survey. First, it provides a detailed characterisation of the chemical inventory of one of the most chemically rich HMCs in the sample. Second, it provides a robust and reproducible template for molecular line identification and analysis that can be consistently applied to the full ALMAGAL dataset. Third, by carefully assessing line blending effects affecting GA, EG, and MF, this work ensures a reliable derivation of their physical parameters and supports forthcoming survey-wide analyses. Finally, the spatially resolved analysis of these species enables us to evaluate their ability to trace the innermost warm and dense regions of HMCs, where COM chemistry is expected to be most active.

Fig. 1 ALMA map of the continuum emission at 1.38 mm of the AG318.9477−0.1960 clump (left) and a zoom-in on core 9 (right). The white contours correspond to 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ, with σ ∼ 0.28 mJy beam−1. The cyan ellipses surrounded by black show the position and size of each core. In particular, the size of core 9 is 0.″66 × 0.″57, with a position angle of PA=42°. The synthesised beam (0.″46 × 0.″37, PA=7°) is shown in the bottom left corner.

2 The source

2.1 AG318.9477−00.1960 clump

AG318.9477−0.1960 is a high-mass clump located at a heliocentric distance of D_HEL ∼ 10.4 ± 0.4 kpc and a Galactocentric distance of D_GAL ∼ 6.84 kpc. It is centred at α_clump(ICRS) = 15^h 00^m 55.s293 and δ_clump(ICRS) = −58° 58′ 52.″511. The clump has a mass of M_clump ∼ 5870 M_⊙, a luminosity of L_clump ∼ 2.7 × 10⁵ L_⊙, an L/M ∼46.5 L_⊙/M_⊙, and a systemic velocity of V_lsr = −34.4 km s⁻¹ (Molinari et al. 2025).

The source is associated with compact mid-IR and far-IR emission, indicating a dense and embedded environment. The radio continuum emission observed with the Australia Telescope Compact Array (ATCA) at 5.5–22.8 GHz detected unresolved free-free emission with a rising spectral index of α = 0.69±0.22, which was interpreted as arising from a partially optically thick conical non-collimated radio jet (Purser et al. 2016). The possibility that the free-free emission is associated with an HII region was dismissed based on the low radio luminosity measured at 22.8 GHz, which is only 1.4% of that expected for an optically thin HII region. Additional evidence of jet activity includes aligned knots of H₂ emission (Lee et al. 2001; De Buizer 2003), a bipolar SiO outflow (De Buizer et al. 2009), and CH₃OH masers at 36 GHz and 44 GHz (Voronkov et al. 2014).

2.2 Core 9

The compact cores embedded in the ALMAGAL clumps were identified with the structure identification algorithm CUrvature Thresholding EXtractor (CuTEx; Molinari et al. 2011) based on their integrated dust continuum flux. This algorithm has allowed us to identify 6348 cores in 844 of the 1013 ALMAGAL clumps (Coletta et al. 2025). By fitting Gaussians to each extracted core, Coletta et al. (2025) have estimated the dust continuum emission from the line-free channels and the size of the continuum emission. Within clump AG318.9477−0.1960, 19 compact cores were identified (Fig. 1), among which core 9 stands out as a particularly chemically rich source. We discuss the criteria used to choose this source among all the others in the ALMAGAL sample in 4.1 in detail. The properties of AG318−c9 are described in Table 1.

3 Observations

The ALMAGAL large project observed a total of 1013 targets using ALMA. The observations were carried out with the main 12 m array in two different configurations: an extended configuration (C-5 and C-6), and a compact configuration (C-2 and C-3), hereafter referred to as TM1 and TM2, respectively; (Escoffier et al. 2007). The observations were also carried out with the 7 m Atacama Compact Array (Iguchi et al. (2009); hereafter 7M). We used data that combine all three configurations (7M + TM1 + TM2), yielding a synthesised beam of 0.″46 × 0.″37 with a PA=7 °. For AG318−c9, this corresponds to physical scales of ∼4800 au × 3900 au along the major and minor axes, respectively.

The observations were conducted in ALMA Band 6 at ∼ 1.38 mm (∼217 GHz). The spectral setup comprised four spectral windows (spw). We used spw0 and spw1, each covering a bandwidth of 1.87 GHz with 3840 channels of 488 kHz, corresponding to a spectral resolution of ≈1.4 km s⁻¹. These spectral windows cover the frequency ranges 216.98–218.86 GHz and 219.06–220.94 GHz, respectively. The remaining two spectral windows (spw2 and spw3), which cover narrower bandwidths of 0.468 GHz, were not considered here as they correspond to higher spectral resolution subwindows within the frequency ranges already covered by spw0 and spw1 and do not provide additional constraints for line identification, especially for weak species. The RMS noise of the continuum map is σ ≈ 0.28 mJy beam⁻¹ (Coletta et al. 2025). A detailed description of the ALMAGAL data-processing workflow, including calibration, data staging, continuum determination, self-calibration, imaging, and data products, is provided in Sánchez-Monge et al. (2025).

4 Method

To study the richness of the spectral line emission of each core, we extracted the spectra in the two broad spectral windows (spw0 and spw1) by integrating the emission over the ellipse that describes the size of each continuum core. The core identification procedure is described in detail in the ALMAGAL core catalogue (Coletta et al. 2025).

4.1 Identification of the richest chemical core

We sorted the spectra of all the ALMAGAL cores using the intensity of the K = 2 component of methyl cyanide (CH₃CN), where K denotes the projection of the total rotational angular momentum onto the principal axis of the molecule, at 220.730 GHz. This transition is one of the brightest in the spectral setup, is minimally affected by line blending, and is a well-established tracer of HMCs (e.g. Cesaroni 2005; Cesaroni et al. 2011). After sorting the spectra of the ALMAGAL cores, we found that sources with weak emission in the low-K components of CH₃CN exhibit little or no complex chemistry, whereas sources with strong low-K CH₃CN emission display chemically rich spectra, consistent with previous studies (e.g. G31.41+0.31; Beltrán et al. 2005; Mininni et al. 2020). We therefore selected sources with a very good signal-to-noise ratio (S/N) emission in the CH₃CN K = 2 transition because weaker and less abundant species would otherwise not be reliably detected.

Since the forthcoming study will focus on GA, EG, and MF, we further considered the intensity of one of the brightest MF transitions in the observed band, the 17_3,14 → 16_3,13 at 218.298 GHz. The two selection criteria consistently place AG318−c9 among the cores with the highest S/N in the ALMA-GAL sample, supporting its use as a template source for detailed chemical analysis.

Moreover, we selected AG318−c9 as a representative chemically rich core within the ALMAGAL sample based on additional criteria. First, it exhibits a high line density of approximately ∼80 lines GHz⁻¹ (Sánchez-Monge et al., in prep.). This line density is comparable to that observed in the most chemically rich cores of Sgr B2, one of the chemically richest regions in the Galaxy (Belloche et al. 2025). Second, AG318−c9 does not show strong multiple velocity components, which greatly facilitates reliable line identification and detailed spectral analysis. Third, AG318−c9 is the third most massive core in the sample. The two more massive cores were excluded because the most massive core shows prominent self-absorption features, while the second core exhibits highly noisy spectra.

4.2 Line identification

For the line identification in AG318−c9, we searched for molecular species with an S/N ≥ 5 (corresponding to ∼1.35 K) to ensure reliable detections. We first fitted species expected to be bright and abundant in HMCs, such as H₂CO, CH₃CN, or CH₃OH, which simultaneously account for multiple of the brightest transitions in the observed frequency range. For the remaining unidentified lines, we used Splatalogue¹ to search for additional molecular species with transitions consistent with the observed frequencies, drawing spectroscopic data from the CDMS² (Müller et al. 2001, 2005; Endres et al. 2016) and the JPL³ (Pickett et al. 1998) catalogues.

To confirm or discard the identification of the different candidate species, we used the spectral analysis software Madrid data cube analysis (MADCUBA⁴; Martín et al. 2019). This software allowed us to identify and fit different species using the tool called spectral line identification and modeling (SLIM), which generates a synthetic spectrum under the assumption of local thermodynamic equilibrium (LTE) while accounting for line opacity. Given the high gas density of AG318−c9, $n_{{\text{H}}_2}=1.72\times {10}^8{\text{cm}}^{-3}$ (Coletta et al. 2025), the assumption of LTE is well justified. The H₂ column density, $N_{{\text{H}}_2}=7.35\times {10}^{24}{\text{cm}}^{-2}$ (Table 1), was derived using $N_{{\text{H}}_2}=\frac{4M_{\text{core}}}{\pi \mu m_{\text{H}}{D}_{{\text{core}}^2}}$, where M_core is the mass of the core, and D_core = 6396 au is the geometric mean size of AG318−c9, as reported by Coletta et al. (2025). Here, m_H = 1.67 × 10⁻²⁴ g is the hydrogen atom mass, and µ = 2.8 is the mean molecular weight per particle. This $N_{{\text{H}}_2}$ was then used to compute molecular abundances relative to H₂ ($X=N_{\text{species}}/N_{{\text{H}}_2}$).

The synthetic spectra were generated and controlled by input physical parameters: column density (N), excitation temperature (T_ex), radial velocity of the source (V), full width at half maximum (FWHM), and source size (θ_s). Best-fit solutions were obtained using the AUTOFIT routine, which performs a non-linear least-squares minimisation based on the Levenberg-Marquardt algorithm (Levenberg 1944; Marquardt 1963).

All detectable molecular species towards AG318−c9 were initially identified and fitted manually using MADCUBA. As we aimed to establish a framework for the chemical analysis of the full ALMAGAL sample, AG318−c9 was used as a testbed to develop a script that automates the line identification and fitting procedures (Appendix A). This automation facilitates the evaluation of individual molecular species in the survey, improving the robustness of the derived physical parameters and their associated uncertainties, particularly for fainter species, such as EG and GA. Although we attempted to derive all the physical parameters for each species, in some cases, AUTOFIT did not converge, and certain parameters had to be fixed. In particular, T_ex was fixed to the value derived for MF ${\text{C}}_2{\text{H}}_5^{13}\text{CN}$ when only one transition was available or when the detected transitions did not span a sufficiently wide range of energies to converge in a solution. The full procedure is described in Appendix A.

The species for which we left T_ex free generally span a wide range of energies, with minimum values typically of a few tens of kelvin and maximum values reaching several hundred kelvin. The widest E_up range was found for NH₂CN (579 K), and the smallest range was found for CH₃COOH (57 K).

Since EG, GA, and MF have clearly been detected in AG318−c9, we also generated N and T_ex maps (Sect. 6.4). Pixel-by-pixel fitting of selected transitions was performed (Appendix B), allowing us to better constrain the physical conditions and chemical structure of the core.

5 Spectral line analysis

We present the complete chemical inventory analysed in AG318−c9 derived from the ALMAGAL frequency setup. To place these results, we compare them with those of the well-studied HMC G31.41+0.31.

5.1 Molecular inventory

The analysis of the spectra reveals molecular emission from a total of 65 different species (Table 2). This underlines the chemically rich environment in AG318−c9. Of these species, 31 correspond to rotational transitions of main species, and the remaining 34 species arise from isotopologues and vibrationally excited species. Figures 2 and 3 show the observed spw0 and spw1 spectra together with the global LTE fit, which accounts for the contribution of all identified species. Zoomed-in versions of the spectra, presented in Appendix C, illustrate the quality of the fits for weaker transitions. Despite the large number of detected species, not all observed spectral features could be reliable identified. A total of 91 lines remained unassigned and are labelled unidentified (U) in the spectra. Although the most recent CDMS and JPL catalogues were used, some of these unidentified features may correspond to transitions that are not yet included in current laboratory spectroscopic databases.

Overall, the identified molecular inventory reveals a chemically rich environment dominated by oxygen-bearing species (44), followed by nitrogen-bearing (28), sulphur-bearing (8), and silicon-bearing (2), with COMs accounting for more than half of the detections (33). In particular, GA, EG, and MF were identified through multiple transitions with good S/N ratios (≥ 5), of which 3, 8, and 14 transitions, respectively, were unblended. Additional heavy COMs, such as acetone (CH₃COCH₃) and ethyl formate (C₂H₅OCHO), are also present.

Remarkably, this rich chemical inventory is revealed by covering a narrow spectral setup of 4 GHz and with short integration times. The number of different species detected in AG318−c9 is comparable with those detected towards typical chemically rich HMCs such as G31.41+0.31, Sgr B2, and IRAS 16293−2422. However, the inventory of species in these well-known HMCs has been built through large unbiased spectral surveys: e.g., the GUAPOS survey (G31.41+0.31 Unbiased ALMA sPectral Observational Survey; Mininni et al. 2020) and the EMoCA survey (Exploring molecular complexity in the Galactic Centre with ALMA; Belloche 2017), which both span nearly the entire ALMA Band 3 (∼32 GHz), and the PILS survey (The ALMA Protostellar Interferometric Line Survey; Jørgensen et al. 2016) from 329 to 363 GHz (∼34 GHz), allowing a more complete census of molecular species. All this suggests that AG318−c9 is a chemically reference HMC within the ALMAGAL survey that might be as rich as other well-known and studied sources.

The detected species span a wide range of molecular complexity, from simple diatomic molecules such as the isotope of carbon monoxide ¹³CO to molecules with 11 atoms such as ethyl formate (C₂H₅OCHO). The complete list of identified species and their correspondent transitions is available in a repository table⁵. Together, this inventory supports the classification of AG318−c9 as a chemically rich HMC.

For all identified species, we derived the column density, excitation temperature, velocity, line width, and abundance (Table 3), following the procedure described in Appendix A. The excitation temperatures of the detected species are in the range 51-250 K, while the column densities are 10¹³–10¹⁸ cm⁻², and their relative abundances range from 10⁻¹¹ to 10⁻⁷. For species for which isotopologues were also detected (H₂CO, SO₂, SO, HC₃N, HNCO, HC(O)NH₂, CH₃OH, CH₃CN, t − HCOOH, CH₃OCH₃, CH₃CHO, and C₂H₅CN), the column densities of the main isotopologues were obtained by scaling those of the isotopologues by the measured isotopic ratios. These ratios were calculated taking the Galactocentric distance (D_GAL) of AG318−c9 into account, which is 6.84 kpc (calculated from the D_HEL of 10.4 kpc; Molinari et al. 2025). Specifically, we adopted ¹²C/¹³C = 46.6 ± 10 (Yan et al. 2019), ³²S/³³S = 88.9 ± 7.7, and ³²S/³⁴S = 20.7 ± 1.2 (Yan et al. 2023).

Overall, the global fit reproduces the observed spectra well in most frequency ranges, with the exception of some transitions of CH₃CN, H₂CO, and ¹³CO, which are not fully reproduced by the synthetic spectrum and present overestimated emission. This behaviour arises because the CH₃CN transitions with K = 0, 1, 2 are optically thick and show absorption features, as is also the case for the H₂CO 3_0,3 → 2_0,2 transition at 218.222 GHz. The CH₃CN,K = 3 ladder belongs to the A symmetry species associated with the internal rotation of the CH₃ group, whereas most other K ladders are of E symmetry. Since A and E species do not interconvert, differences in their populations or traced physical conditions prevent a single LTE model from reproducing all K ladders simultaneously, leading to an overestimation of the K = 3. In addition, the ¹³CO emission is spatially extended and is therefore filtered out by the interferometer, resulting in missing flux.

The detection of multiple COMs, and in particular, of EG and GA, constitutes a key result of this study, placing this source among the list of cores exhibiting signatures of COMs with potential prebiotic relevance. GA has been reported in only a limited number of sources to date (e.g. Hollis et al. 2000; Beltrán et al. 2009; Jørgensen et al. 2012; Chen et al. 2023), and its secure detection towards AG318−c9 reinforces the view that this core hosts physical and chemical conditions that are favourable for prebiotic chemistry, including grain-surface processes and warm-up timescales enabling ice sublimation.

Fig. 2 Spectra of spw0 towards core 9 of the AG318.9477−00.1960 clump. The black histogram and its grey shadow are the observational spectrum, and the red curve indicates the LTE fit after taking into account the contribution of all the species (see Table 3 and its complete version in Appendix D.1). The horizontal dashed green line indicates the 5σ ≈ 1.35 K level, with σ being the noise (σ ∼ 0.28 K). Only transitions contributing at ≥ 5σ are labelled. The vertical dashed black lines indicate the position where ν was converged. U labels indicate unidentified transitions.

Fig. 3 Spectra of spw1 towards core 9 of the AG318.9477−00.1960 clump. The black histogram and its grey shadow are the observational spectrum, and the red curve indicates the LTE fit after taking into account the contribution of all the species (see Table 3 and its complete version in Appendix D.1). The horizontal dashed green line indicates the 5σ ≈ 1.35 K level, with σ being the noise (σ ∼ 0.28 K). Only transitions contributing at ≥5σ are labelled. The vertical dashed black lines indicate the position where ν was converged. U labels indicate unidentified transitions.

5.2 Comparison with the hot core G31.41+0.31

To place our results in context, we compared the chemical reservoir of AG318−c9 with that of the well-studied HMC G31.41+0.31 (G31 hereafter). G31 is the first core in which GA was detected outside the Galactic centre (Beltrán et al. 2009), and it is one of the first cores in which EG was found (Rivilla et al. 2017). This HMC was also the target of the GUA-POS project, a spectral survey of the entire ALMA band 3 (84.05–115.91 GHz; Mininni et al. 2020). This HMC is located at a D_HEL of 3.75 kpc (Immer et al. 2019) with a L = 4.4 × 10⁴ L_⊙ (Osorio et al. 2009) and a gas mass M ∼ 70 M_⊙ (Cesaroni 2019). By contrast, AG318−c9 is significantly more distant, at a heliocentric distance of ∼10.4 kpc, and its estimated gas mass is M ∼ 120 M_⊙ (Coletta et al. 2025). However, this mass is likely an overestimate because it was obtained by assuming a uniform core temperature of 49 K, whereas our T_ex maps reveal a temperature gradient in the core with temperatures significantly >49 K (Sect. 6.4), implying a lower total gas mass. Since a core-specific luminosity is unavailable, we estimated it based on the host clump luminosity (2.7 × 10⁵ L_⊙), which is fragmented into 19 cores (Coletta et al. 2025). We apportioned the clump luminosity to AG318−c9 by scaling it by the ratio of its peak flux $(F_{{\text{peak}}_{\text{core}}})$ to the sum of all the core peak fluxes $(F_{{\text{peak}}_{\text{all\hspace{0.17em}cores}}})$, lead ing to L_core ≈ 1.2 × 10⁵ L_⊙, which is almost three times higher than the luminosity of G31. Although we would like to interpret the evolutionary stages of the two HMCs based on their L/M ratios, this approach is subject to significant uncertainties. We note that the interpretation of the evolutionary stage based on global indicators such as the luminosity-to-mass ratio is not straightforward in AG318−c9. The core mass is likely overestimated, while the luminosity is indirectly inferred from the host clump, introducing additional uncertainties. On the other hand, the mass of G31 should be considered as a lower limit (Cesaroni 2019). Additionally, in high-mass star-forming regions, L/M may not provide an unambiguous evolutionary indicator due to the effects of fragmentation, temperature gradients, and internal source structure (Motte et al. 2018).

The GUAPOS survey (Mininni et al. 2020; Colzi et al. 2021; Mininni et al. 2023; Fontani et al. 2024; López-Gallifa et al. 2024, 2025) explored the chemical richness of G31, detecting numerous O- and N-bearing COMs as formamide (HC(O)NH₂), methyl isocyanate (CH₃NCO), and ethyl iso-cyanide (C₂H₅NC). A one-to-one comparison of the column density for species detected in the two sources is presented in Fig. 4, where the spectra are extracted over comparable spatial scales (∼6400 au for AG318−c9 and ∼4500 au for G31), indicating that the sampled regions are broadly similar. A clear trend emerges: although G31 has a lower luminosity and mass, the column densities of the species detected towards this HMC are systematically higher than those estimated towards AG318−c9 for most of the species. This discrepancy is particularly pronounced for O-bearing COMs, such as CH₃COOH, CH₃COCH₃, and H₂C¹⁷O, which show significantly higher N values in G31. Only two species present enhanced N in AG318−c9, which are t − HCOOH, and C¹⁸O. N-bearing species also show higher values in G31, but not as high as O-bearing species. S- and Si-bearing species show column densities that are more comparable, clustering closer to the 1:1 correlation line (i.e. within ∼ 1 dex). We found 11 exceptions to the general trend, in which species present enhanced N in AG318−c9: CH₃CHO, t-HCOOH, C¹⁸O, H₂CNH, H¹⁵NCO, NH₂CN, SO₂, SO, ³⁴SO₂, SiO, and CH₃OCH₃.

The observed chemical segregation, characterised by a pronounced excess of O-bearing COMs in G31, strongly suggests a more efficient thermal sublimation of the ice mantles in G31 (Mininni et al. 2023). This implies that G31 is either intrinsically warmer or at a more advanced evolutionary stage than AG318−c9, having had more time to release and accumulate a larger fraction of these species into the gas phase.

N-bearing species have more diverse formation pathways, including active gas-phase chemistry, and require a longer timescale to be synthesised in the cold cloud (Lee et al. 2020). Their abundances are also sensitive to environmental factors such as the UV radiation field, shocks, cosmic-ray ionisation rate, and local N/O elemental ratios. This may explain why their column densities are more comparable between the two sources. Among the exceptions, the higher N of the prebiotic species cyanamide (NH₂CN) or the methanimine (H₂CNH) in AG318−c9 is particularly noteworthy, suggesting that the physical conditions or precursor availability in AG318−c9 are uniquely favourable for their production. Since D-bearing species are thought to be formed more efficiently, at the low temperatures (< 20K) of the preceding cold core phase, their detection in the hot gas of AG318−c9 implies that they were formed on the ice grains and subsequently sublimated (Caselli & Ceccarelli 2012). However, G31 also presents some of the D-bearing species found in AG318−c9 (V. Rivilla, priv. comm.). This suggests that the two cores share a similar chemical heritage from their cold pre-stellar phase.

The emission of Si- and S-bearing species is usually associated with shocks, as these species would be released from the dust mantles by sputtering. The fact that the column density of these species in both HMCs is comparable suggests that both cores are associated with shocks, likely produced by molecular outflows. For G31, Beltrán et al. (2022) has detected at least six outflows in the core, while for the AG318.9477−00.1960 clump, a preliminary analysis suggests that this core is associated with at least two outflows (M. Benedettini, priv. comm.).

In Fig. 5, we show a comparative histogram of the molecular abundances $(X=N_{\text{species}}/N_{{\text{H}}_2})$ for the two sources. We also included the species detected exclusively in G31 and those that are detected in AG318−c9, but were not reported for G31 to date. To ensure a fair comparison, we restricted this selection to species covered by the spectral setups of both surveys. For this reason, we excluded the following molecular species: NS, OCS, NH₃, NH₂D, N₂H+, HCN, CO, CS, CH₃CCH, CH₃NC, s-C₂H₅CHO, and ¹³CO. In addition, a direct comparison with the GUAPOS survey is not yet possible because their full deuterated inventory has not yet been published, and for this reason, DCN, HDCS, CH₃OD, CH₂DOH, CHD₂OH, and DC(O)NH₂ are not included in Figs. 4 and 5.

In addition to the enhanced column densities reported above, the comparison based on molecular abundances (Fig. 5) indicates that C₂H₅CN is comparatively more abundant in AG318−c9. Moreover, a distinct chemical divergence can be noted: G31 features complex N-bearing COMs such as acetamide (CH₃CONH₂) and methylamine (CH₃NH₂). The higher N of acetamide is particularly significant because it contains the peptide-like linkage (–NH–C=O), a key component for prebiotic chemistry (Colzi et al. 2021). The molecular inventory of AG318−c9 is distinguished by the presence of O-bearing species that have not been reported in G31 so far, most notably, ethyl formate (C₂H₅OOCH), along with others such as HOCH₂CN. The detection of ethyl formate in AG318−c9, together with its apparent absence in the G31 survey suggests a chemistry dominated by the thermal sublimation of O-rich ice mantles, releasing these O-bearing COMs into the gas phase (Belloche et al. 2009).

This comparison suggests that the two sources might correspond to different evolutionary stages. AG318−c9 might represent a younger HMC, where the chemistry is still dominated by the recent sublimation of ice mantles rich in ester-bearing species. In contrast, G31 appears to be a more evolved environment in which this initial chemical reservoir has likely been processed by warm gas-phase reactions.

Although the overall lower abundances in AG318−c9 are consistent with a less chemically processed source, abundance differences alone do not uniquely constrain the relative evolutionary stage of the two cores. Variations in temperature structure, radiation field, source geometry, or chemical pathways might also contribute to the observed differences (Ossenkopf et al. 2013; Heays et al. 2017). Therefore, the scenario in which AG318−c9 is younger than G31 should be regarded as plausible rather than definitive. These differences indicate distinct chemical evolutionary pathways, although detailed chemical modelling is required to draw firm conclusions.

Fig. 4 Column density comparison for all common species between AG318−c9 (x-axis) and G31 (y-axis). The solid black line indicates the 1:1 relation. The dark and light orange shaded regions represent deviations within 0.5 dex and 1.0 dex from equality, respectively.

Fig. 5 Top: comparison of molecular abundances (X = N/H2) towards G31 (blue bars) and AG318−c9 (orange bars), sorted by molecular groups and by the abundance of AG318−c9 inside each group. The species are grouped on the x-axis according to their detection status. Left: detected in both sources. Middle: detected in G31 alone. The orange triangles represent the corresponding upper limits for AG318−c9. Right: reported in AG318−c9. Bottom: logarithm of the ratio of the abundances of the two sources. HC3N* represents the group composed of HC3N, HC3N (v6 = 1), HC3N (v7 = 1), and HC3N (v7 = 2).

6 Glycolaldehyde, ethylene glycol, and methyl formate

We focused on the identification of GA, EG, and MF emission towards AG318−c9, paying particular attention to line-blending effects that might affect their detection. GA stands out as a key prebiotic species, being the simplest sugar-like monosaccharide molecule that plays a fundamental role in sugar-formation pathways, leading to the formation of essential biomolecules such as glyceraldehyde and ribose, the backbone of nucleic acids. GA has been detected in different star-forming regions in the ISM, including the first detection in Sgr B2(N) in the Galactic centre (Hollis et al. 2000), the HMC G31.41+0.31 (Beltrán et al. 2009), and the solar-type protostar IRAS 16293-2422 (Jørgensen et al. 2012).

The heavy COM EG, with ten atoms, and the fully reduced form of GA where the aldehyde group has been reduced to an alcohol, although not directly prebiotic, is chemically linked to GA, and they may both share the same formation pathways in the ISM (Sorrell 2001; Fedoseev et al. 2015; Coutens et al. 2018; Garrod et al. 2022). EG was first found in the ISM towards the Galactic centre (Hollis et al. 2002) in HMCs such as G31.41 + 0.31 (Rivilla et al. 2017) or G34.3 + 0.1 (Brouillet et al. 2015) and in hot corinos such as NGC 1333-IRAS 2A (Maury et al. 2014).

Finally, MF, although not chemically linked through the same formation pathways, can offer complementary constraints on the chemistry of GA and EG. All three species are expected to trace similar physical conditions (Simons et al. 2020).

6.1 Glycolaldehyde

In Fig. 6 we show all 19 transitions of GA, clearly above S/N ≥ 5. The transition 20_2,18 → 19_3,17 is the brightest and has the highest opacity (τ = 0.061), confirming that all the GA transitions are optically thin. Unlike for EG and MF (Sects. 6.2 and 6.3), in the available ALMAGAL frequency coverage, there are few unblended transitions of GA. There are three unblended transitions: 9_5,5 → 8_4,4 at 217626.13 MHz, 31_4,27 → 31_3,28 at 217271.61 MHz, and 20_3,17 → 19_4,16 at 218260.59 MHz. Most of the blending of GA is with CH₃COCH₃, $(X=N_{\text{species}}/N_{{\text{H}}_2})$, CH₃OCH₃, and CH₃NCO transitions, whose contribution is included in the global red fit in Fig. 6, and it explains the remaining features of the observed spectrum. The fitted transitions span a broad range of energies (E_up = 26–316 cm⁻¹ or 37–450 K).

Although SLIM converged to a solution for GA without the need to fix any parameter, we fixed the FWHM to 4.0 km s⁻¹, as the line width obtained from the unconstrained fit might be visually refined to match the observed spectra better. The resulting fit of the detected transitions of GA indicates T_ex = 168 ± 13 K, N = (1.45 ± 0.12) × 10¹⁶ cm⁻² and V = −36.89 ± 0.1 km s⁻¹.

6.2 Ethylene glycol

Ethylene glycol is a heavy COM and the fully reduced alcohol form of GA. The presence of both species, which are thought to be chemically linked through common formation pathways (Coutens et al. 2018; Simons et al. 2020; Garrod et al. 2022), confirms the complexity of the chemistry in AG318−c9. EG can adopt different conformations because torsions around its C–C bond and its two C–O bonds lead to four forms with an anti-arrangement (a-EG) of the hydroxyl groups (–OH), and six with a gauche-arrangement (g-EG) (Kazerouni et al. 1997; Rivilla et al. 2017). We focused on the a-EG conformer because it shows a significantly larger number of transitions within the ALMA-GAL spectral setup compared to the g-EG conformer. The latter shows many blended lines, and the few unblended transitions have similar energies, E_up ≈ 130 K, so that it is difficult to properly constrain the physical parameters.

Figure 7 shows selected lines corresponding to a-EG, eight of which are unblended and strong. Other transitions are more blended, mainly with MF, CH₃NCO, HNCO, HN¹³CO, C₂H₅OH, and CHD₂OH. The 24_1,24 → 23_1,23 transition at 217450 MHz is the brightest and the most optically thick with τ = 0.046. To fit the EG emission, we took 37 different transitions into account, which span a range of energies of E_up = 79–389 cm⁻¹ or 114–580 K. Fig. 7 shows that the LTE fit, with the joint contribution of all blending species, is able to properly explain the observed emission.

For a-EG, SLIM converged to a solution for a-EG without the need to fix any parameters. The resulting fit of the detected transitions of a-EG indicates T_ex = 240 ± 30 K, N = (5.1 ± 0.6) × 10¹⁶ cm⁻², FWHM = 6.19 ± 0.19 km s⁻¹, and V = −35.9 ± 0.08 km s⁻¹.

Fig. 6 Selection of GA transitions observed and fit (see Sect. 6.1). The black histogram and its grey shadow are the observational spectrum, and the red curve is the cumulative LTE fit considering all detected species. The blue curve is the fit of the individual transitions. The horizontal dashed green line indicates the 5σ level, corresponding to ∼1.35 K. Only transitions contributing at ≥5σ are labelled. The plots are sorted by decreasing intensity of the GA transitions (blue line). The upper left corner of each panel indicates the EUP of the corresponding GA transition. When multiple GA transitions are shown in the same panel, the EUP values are listed from top to bottom following the left-to-right order of the transitions.

6.3 Methyl formate

Methyl formate shows the highest number of clean transitions of the three COMs, with a total of 14 unblended transitions. In Fig. 8, we show selected lines corresponding to MF. To fit the MF emission, we took a total of 41 transitions into account (see Appendix A), with 17_3,14 → 16_3,13 at 218297.89 MHz being the brightest. The same transition also has the largest τ = 0.185, being the transition with highest opacity among all the fitted transitions of EG, GA, and MF. In any case, we can still consider all MF transitions as optically thin. The fitted MF transitions cover a wide energy range E_up = 69–318 cm⁻¹ or 100–457 K.

The estimated excitation temperature, T_ex = 167 ± 3 K, is consistent with that derived for EG and GA. The column density is one order of magnitude higher N = (2.2 ± 0.4) × 10¹⁷ cm⁻² than those of EG and GA, which indicates that MF is more abundant, as already observed in other works (Hollis et al. 2001; Mininni et al. 2020). In the case of MF, we did not fix any parameter. We obtained a value for the FWHM = 5.59 ± 0.08 km s⁻¹ and V = −35.22 ± 0.04 km s⁻¹, which are similar to those of EG and GA.

Fig. 7 Same as Fig. 6, but for a-EG, as described in Sect. 6.2.

6.4 Integrated intensity, T_ex, and N maps of EG, GA, and MF

A key aspect of our analysis is to understand how EG, GA, and MF trace the physical environment of AG318−c9. To this end, we used MADCUBA to compute their integrated intensity (moment 0), N and T_ex maps (Appendix B), together with radial profiles of N and T_ex obtained by averaging the values within annuli of 0.″1. The resulting maps and profiles are presented in Figs. 9–11, and the corresponding relative uncertainty maps are shown in Appendix B.1.

To generate the maps, we initially attempted to use the same transitions as we adopted to derive the physical parameters presented in Sect. 5.1. However, for EG and GA, a reduced set of transitions was selected, prioritising the least blended transitions while preserving a broad energy coverage. This choice was necessary to ensure convergence of the fits, particularly in the outer core regions with a lower S/N, where the inclusion of blended transitions prevented reliable solutions. The transitions used to produce the maps are listed in Appendix B.2. For EG, the selected transitions cover the same energy range as used in the global fit (E_up = 114–580 K), while for GA, they span a similar range (E_up = 37–410 K).

The moment 0 maps shown in Fig. 9 were constructed using representative transitions for each molecule: the 35_10,26 → 35_9,27 transition at 220.055 GHz for GA (E_up = 410 K), the 24_1,24 → 23_1,23 transition at 217.45 GHz for EG (E_up = 117 K), and the 17_3,14 → 16_3,13 transition at 218.298 GHz for MF (E_up = 100 K). While the EG and MF moment 0 maps were generated with comparable upper-level energies, the GA map relied on a significantly higher E_up transition. This difference requires caution when interpreting relative spatial extents based solely on moment 0 emission. Although the moment 0 maps were required to automate the generation of the N and T_ex maps, they were not used to derive these quantities directly. Instead, they were used solely to estimate the noise level, which was then applied to assess the reliability of the pixel-by-pixel spectral fits and to exclude solutions below a 3σ threshold (Appendix B). The N and T_ex maps are therefore based on multi-transition fits and are not biased by the specific energy of the moment 0 transition.

To generate N and T_ex maps, only the target species was considered in each case, without considering the contribution of other molecules (Appendix B). As a result, the core-averaged physical parameters reported in Tables 3 and D.1, which account for all detected species, are more robust. Nevertheless, the values derived from the maps are consistent with the core-averaged quantities.

Of the three species, MF is most spatially extended. Its emission closely follows the dust continuum emission and peaks near the continuum maximum (Fig. 9). The N(MF) map (Fig. 10) reveals a very similar morphology, with a clear column density gradient ranging from ∼10¹⁶ cm⁻² in the outer regions to ∼3.16× 10¹⁷ cm⁻² near the continuum peak, corresponding to the highest N among the three COMs. The gradient is clearly reflected in the radial profile (Fig. 11), which decreases steadily from the core centre towards larger radii. The T_ex(MF) map (Fig. 10) shows a robust temperature gradient, rising from ∼100 K in the outer parts of the core up to ∼250 K towards the centre. The radial profile presents a more rapid temperature increase inward of 0.″8. This combination of large spatial extent and wide T_ex dynamic range confirms MF as an excellent tracer of the gas temperature and its gradients in HMCs (e.g. Beltrán et al. 2018). These temperatures are consistent with the high rotational temperatures (120–300 K) expected in hot cores (Herbst & van Dishoeck 2009) and typical T_ex ∼ 150 K estimates Cesaroni 2005, Bisschop et al. 2007; Rivilla et al. 2017. Because MF has the most unblended transitions, its T_ex map has the lowest relative error $(\frac{\text{Δ}{T}_{\text{ex}}}{T_{\text{ex}}}\sim 10\%)$, which increases only in the outer regions with a lower S/N (Appendix B.1).

The distribution of EG is more compact. Its integrated-intensity and column density maps (Figs. 9 and 10) are more concentrated towards the inner core. The N(EG) increases from ∼10¹⁵ cm⁻² at the edges to ∼5 × 10¹⁶ cm⁻² towards the centre. The relative errors for N(EG) increase slightly in the centre of the emission (Fig. B.1). The T_ex(EG) map (Fig. 10) does not exhibit a well-defined peak as clearly as for MF. For clarity, the map is truncated at a radius of 0.″9 beyond which the uncertainties become too large for a reliable interpretation. Within this radius, the highest T_ex values are found within the core, coinciding with the region of strongest dust continuum emission. Overall, the temperature distribution appears to be relatively uniform throughout the inner region where EG emission is detected, with a relative error of ∼30% (Fig. B.2). This behaviour might partly reflect the greater difficulty in constraining the fits due to fewer unblended lines, as reflected in the higher relative error, which increases considerably in the outer colder region. Nevertheless, the radial profile clearly decreases in temperature from ∼170 K to ∼80 K with increasing distance from the core centre, although the profile appears to be less steep in the innermost regions.

The most compact species is GA, with emission confined to the innermost region of the core, close to the continuum peak. The N(GA) and T_ex(GA) maps peak centrally, with column densities ranging from ∼2 × 10¹⁶ cm⁻² to ∼4 × 10¹⁶ cm⁻² and an excitation temperature between ∼130 K to ∼160 K. The two map fits have a relative error of ∼10–15% (Appendix B.1). Over the radial range in which EG and GA are detected, the N gradient is comparable to that estimated in EG.

To evaluate whether the different spatial extents of MF, EG, and GA are driven by sensitivity limitations or reflect intrinsic chemical differentiation, we estimated the column density corresponding to a 5σ detection threshold (∼1.35 K) for each molecule with MADCUBA. For MF and EG, the 5σ limits (N ∼ 10¹⁶ cm⁻² and N ∼ 2.5 × 10¹⁵ cm⁻², respectively) closely match the minimum column densities reached by their radial profiles, indicating that the apparent truncation is largely driven by sensitivity. In contrast, the 5σ threshold for GA is N(GA) ≈ 2 × 10¹⁵ cm⁻², whereas the GA radial profile only extends to N(GA) ≈ 1.6 × 10¹⁶ cm⁻², almost one order of magnitude. This might suggest that the lack of detectable GA emission beyond the innermost regions of the core is not limited by sensitivity, but instead reflects its intrinsically compact distribution.

The compact emission of EG and GA suggests that they are excellent candidates for tracing the innermost high-density regions of the core. Due to their optically thin transitions, they might also be used to probe the kinematics within these regions. Furthermore, the fact that EG and GA are only detected in the warm inner region, unlike the more extended MF, might imply that a specific temperature threshold must be reached to sublimate them from the icy dust grains. This observation strongly supports a shared grain-surface formation pathway, consistent with models of dust-grain chemistry (e.g. Fedoseev et al. 2015; Simons et al. 2020; Garrod et al. 2022).

However, the innermost region must be interpreted with caution due to observational effects. At radii comparable to or smaller than the beam size (∼0.″21), the observed brightness temperature might be affected by beam dilution. In addition, dust opacity might become non-negligible towards the core centre, potentially attenuating line emission from the deepest layers. These effects might result in an underestimation of the true column densities and excitation temperatures in the central region, particularly for compact species such as EG and GA. Moreover, as shown in Fig. 11, the emission of the three species extends beyond the beam size, and therefore, it is expected to be hardly affected by beam dilution.

The observed spatial differentiation among the three species suggests that their distributions are not governed solely by thermal desorption, but also by their formation pathways and ice chemistry. MF, with a lower binding energy (∼4000 K; desorption at ∼100–120 K), is expected to be more extended and associated with CO-rich outer regions, whereas GA has a higher binding energy (∼5900 K; ∼140–190 K) and is therefore expected to trace warmer and denser regions, consistent with its more compact emission (Öberg et al. 2009). EG, with an even higher binding energy (∼7500 K; ∼150–180 K), does not show a correspondingly more compact distribution, suggesting that its spatial extent cannot be explained by desorption alone. Instead, its distribution might reflect efficient formation during the warm-up phase, such as via CH₂OH recombination and chemical stability once released into the gas phase (Öberg et al. 2009; Fedoseev et al. 2015). Overall, this chemical differentiation supports a scenario in which the temperature gradients and the initial ice composition can play a key role in shaping the spatial distribution of COMs (Garrod et al. 2008).

Fig. 8 Same as Fig. 6, but for MF, as described in Sect. 6.3.

Fig. 9 Integrated-intensity maps for EG (left), GA (middle), and MF (right). The white contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1, and the red cross shows its peak position. The red ellipses surrounded by white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown in the bottom left corner represents the synthesised beam of the interferometer.

Fig. 10 Column density (N) maps (top row) and excitation temperature (Tex) maps (bottom row) derived for EG (left), GA (middle), and MF (right). The corresponding relative errors are shown in Appendix B.1. The black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1, and the red cross shows its peak position. The red ellipses surrounded by white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown in the bottom left corner represents the synthesised beam of the interferometer.

Fig. 11 Radial profiles derived from the column density (N; top) and excitation temperature (Tex; bottom) maps. The blue, red, and green curves correspond to EG, GA, and MF, respectively. The shaded areas indicate the measurement uncertainties (see Appendix B.1). The bar in the lower left corner indicates the radius corresponding to the geometric mean of the synthetic beam.

7 Conclusions

We presented a comprehensive analysis of the chemical reservoir of the AG318.9477−0.1960 clump core 9 using data from the ALMAGAL project and a frequency range of ∼217–221 GHz. Throughout the ∼ 4 GHz setup, we identified 65 species, of which 44 are O-bearing species, 28 are N-bearing, 8 are S-bearing, and 2 are Si-bearing. Thirty-three of these species are COMs. The fact that most of the detections are O-bearing molecules suggests an efficient thermal sublimation of the ice mantles into the gas phase. For all the detected species, we derived their physical parameters and abundances with respect to H₂.

This spectral analysis allowed us to characterise weak and blended species, such as the three COMs of interest: EG, GA, and MF. For EG, we covered an E_UP range of 114–580 K and detected 8 unblended transitions. We found that the species that are more strongly blended with EG are MF, CH₃NCO, HNCO, HN¹³CO, C₂H₅OH, and CHD₂OH. For GA, we traced transitions with E_UP = 37–450 K including 3 unblended lines. The main blending species are CH₃COCH₃, ${\text{CH}}_3^{18}\text{OH}$, CH₃OCH₃, and CH₃NCO. Finally, MF exhibited transitions with E_UP = 100–457 K and 14 unblended lines, making it the species least affected by blending. All three species trace warm and dense gas, although they probe distinct thermal components. EG traces the warmest material, with a T_ex of ∼240 K. GA and MF trace slightly cooler gas, with T_ex ∼168 K and ∼167 K, respectively. Regarding N and abundance, MF has a column density of ∼2 × 10¹⁷ cm⁻² and is clearly the most abundant species, with X ∼ 3 × 10⁻⁸. This abundance is significantly higher than that of EG and GA, which have abundances of ∼7 × 10⁻⁹ and ∼2 × 10⁻⁹, and N of ∼5 × 10¹⁶ cm⁻² and ∼1.5 × 10¹⁶ cm⁻², respectively.

The chemical comparison with the chemically rich HMC G31.41+0.31 indicated that this core shows systematically higher column densities than AG318−c9 for O- and N-bearing COMs. This suggests that G31 is either intrinsically warmer or at a more advanced evolutionary stage than AG318−c9, having experienced more efficient thermal sublimation of its ice mantles. This comparison also revealed a chemical differentiation between these families, with a stronger enhancement of O-bearing species in G31, while N-bearing species showed a more moderate contrast between the two sources. Conversely, S-bearing species are more comparable between the two sources. Notably, AG318−c9 shows a significant overabundance of the prebiotic molecules NH₂CN and H₂CNH, which indicates specific physical conditions or precursor availability for its formation.

Our analysis of the T_ex and N maps for EG, GA, and MF within AG318−c9 revealed distinct morphologies. The emission of the three species extends beyond the beam size and is therefore not significantly affected by beam dilution. However, observational effects, in particular, dust opacity, might still lead to underestimated column densities and excitation temperatures in the innermost region. MF shows significantly higher column densities (N ≈ 3 × 10¹⁷ cm⁻²) and a spatially extended distribution, robustly tracing the thermal structure of the core over a broad temperature range (T_ex ∼ 100–250 K). In sharp contrast, EG and GA are compact and confined to the innermost core region, with GA being particularly compact. They both exhibit similar lower column densities (N ≈ 4 × 10¹⁶ cm⁻²). Their restricted distribution indicates that they only trace the warmest densest gas, implying a high sublimation temperature threshold. This might support a shared grain-surface formation pathway.

As a methodological output of this work, we developed two MADCUBA scripts that will be used for future molecular analysis studies within the ALMAGAL survey. Work is now underway on a complete study of EG, GA, and MF in all ALMAGAL cores to understand the correlations between these COMs, the physical conditions that lead to their formation, and their possible chemical pathways.

Data availability

The scripts used for the automated line identification and map generation, as well as the complete table of molecular transitions used for the MADCUBA fits of the molecular emission towards AG318−c9 (see Tables 2, 3, and D.1) are available at: https://github.com/JofreAllande/ALMAGAL/.

Acknowledgements

We thank the anonymous referee for the useful comments. The table of molecular transitions were made thanks to the Python library RichValues⁶ created by Andrés Megías. J.A. acknowledges the SKA PhD fellowship “Astrochemistry of prebiotic molecules in high-mass star-forming regions in preparation of the SKA observations” offered by INAF. C.Y.L. acknowledges financial support through the INAF Large Grant “The role of MAGnetic fields in MAssive star formation” (MAGMA). R.S.K acknowledges financial support from the ERC via Synergy Grant “ECOGAL” (project ID 855130) and from the German Excellence Strategy via the Heidelberg Cluster “STRUCTURES” (EXC 2181 – 390900948). In addition R.S.K. is grateful for funding from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206), and from DFG and ANR for project “STARCLUSTERS” (funding ID KL 1358/22-1). L.C. acknowledges support from the grant PID2022-136814NB-I00 by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICI-U/AEI/10.13039/501100011033 and by ERDF, UE. The project that gave rise to these results received the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PR25/12110012. L.B. gratefully acknowledges support by the ANID BASAL project FB210003. C.M. and the INAF-IAPS team acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 program through the ECOGAL Synergy grant (ID 855130). C.M. acknowledges funding from INAF Mini Grants RSN2 2024 "Zodyac" CUP C83C25000340005. R.K. acknowledges financial support via the Heisenberg Research Grant funded by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) under grant no. KU 2849/9, project no. 445783058. A.S.-M. acknowledges support from the RyC2021-032892-I grant funded by MCIN/AEI/10.13039/501100011033 and by the European Union ‘Next GenerationEU’/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and support from the PID2020-117710GB-I00 (MCI-AEI-FEDER, UE). A.L-G. acknowledges 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; 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). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

References

Appendix A Automatization with MADCUBA

We developed a MADCUBA macro to automate line identification across the ALMAGAL sample. Molecular species identified toward AG318−c9 using the detailed line identification described in Sect. 4.2 were added to the SLIM module using the CDMS⁷ (Müller et al. 2001, 2005; Endres et al. 2016) or the JPL⁸ (Pickett et al. 1998) catalogues. As previously mentioned, for each species only sufficiently strong (S/N ≥ 5) and partially unblended transitions were selected, ensuring that AUT-OFIT operates on high-quality lines and minimising potential errors from blended or contaminated transitions. In future studies where the entire ALMAGAL sample will be analysed and chemically poorer sources will be included, the presence of species detected in AG318−c9 will not pose a problem, since AUTOFIT will not converge for species absent in the observed spectra and will not affect the fitting of other species.

SLIM requires an initial input of N, T_ex, V, FWHM, and θ_s for each species to run AUTOFIT. If these initial values are absent or far from realistic estimates, the programme does not converge. For AG318−c9, we assumed as initial guesses V = −34.4 km s⁻¹, which is the systemic velocity of the clump, and FWHM = 5 km s⁻¹, which is the typical width of the lines for this source, as obtained from simple Gaussian fits to some species in our setup. The parameter θ_s was fixed to 0, which in MADCUBA means that the source emission size is equal to or larger than the synthesised beam and hence unaffected by beam dilution. This assumption is appropriate for AG318−c9. To reliably estimate T_ex, we require multiple and unblended transitions spanning a wide energy range. Because this is not possible for all the species, particularly if a species has a single transition, we decided to use a common T_ex as an initial input, estimated from a species with many transitions covering a broad range of energies. To achieve this automatically, we have estimated this common T_ex by calculating the rotational diagram of this species with MADCUBA. Although we initially wanted to use CH₃CN to estimate such a common initial T_ex, we realised that the lowest K transitions were saturated, indicating that they are optically thick, leading to an overestimate of the rotational temperature (T_rot) and an underestimate of the rotational column density (N_rot). Therefore, we used a total of 41 selected optically thin MF transitions that contribute significantly to the observed spectrum. Among these, 14 are non-blended and 27 are slightly blended. The selected transitions offer a broad coverage of energies (E_up99.72–429.4 K or 69.31 cm⁻¹–298.45 cm⁻¹) that allowed us to properly calculate the rotational diagram of MF and obtain T_rot, which has been used as initial seed for T_ex for all species. The rotational column density obtained from the diagram has also been used as initial guess for N to fit MF.

If AUTOFIT converges for MF with these initial parameters, the same T_ex, V, and FWHM values are applied, as initial guesses, for the remaining species, with the exception of N. Assuming starting common values of T_ex, V and FWHM for different species tracing the HMC is reasonable, as they are not expected to vary significantly. In contrast, N can change significantly and is the parameter for which the synthetic spectrum is most sensitive. To address this, SLIM provides the "SLIM_Estimate_LogN"⁹, that estimates log(N) for a given species by using the transition with the highest expected intensity based on the molecular parameters (T_ex, V, and FWHM).

Fig. A.1 Flowchart of the MADCUBA macro that automates the line identification process. The asterisk (*) indicates that physical parameters are kept free whenever possible, following the criteria described in the text.

This estimate is derived from the velocity-integrated intensity obtained by integrating the line emission in the velocity range defined by the linewidth.

Once appropriate inputs are obtained for all the species, a first AUTOFIT is performed in two rounds. In the first round, the FWHM remains fixed, and each species is fitted individually without taking into account the contribution from other species (disabling the ALL SPECIES flag in AUTOFIT). During this first step, brighter and unblended molecules converge easily and independently. The converged molecules provide constraints that will help in the second AUTOFIT iteration, which will try to fit all the species, including those blended. In this second round, the FWHM is left free and the contribution of all blended species is properly taken into account (enabling the ALL SPECIES flag in AUTOFIT).

We have tried to leave as many free physical parameters as possible in order to have the most unforced solution. However, for species with limited transitions, poor energy coverage, or blended lines, we adopted a hierarchical approach by fixing parameters when AUTOFIT does not converge: i) first FWHM, because usually this parameter can easily be obtained from the observed spectra; ii) V is fixed next, using the systemic velocity (V₀) as a reference to constrain velocity shifts, since large deviations from this value are not expected; and iii) T_ex is fixed when only one transition is available, as for example for${\text{H}}_2^{13}\text{CO}$ or C¹⁸O, or when the number of transitions with different enough energies is not sufficient to make the fit to converge. In our approach, we never fixed N. Although this procedure is largely automated, the resulting fit must be visually inspected to ensure consistency with the observations, and this last step should be done manually. For a flowchart of the process, see Fig. A.1. Note that if MF does not converge, the automated macro stops. This could be due to MF not being present in the spectra, being too weak to be properly fit, or due to incorrect initial guesses. In the latter case, we should investigate the reason and, if possible, adjust the initial guesses.

Appendix B Generation of the maps (N)

To generate the T_ex and N maps we used the SLIM module to fit the spectra on a pixel-by-pixel basis. MADCUBA provides a built-in tool to generate the integrated intensity map (moment 0), velocity field (moment 1), and velocity dispersion (moment 2) maps for selected transitions. To derive the moment maps, we selected an unblended transition for each of the three molecules of interest. For GA the transition 35_10,26 → 35_9,27 at 220.055 GHz, for EG the transition 24_1,24 → 23_1,23 at 217.45 GHz, and for MF the transition 17_3,14 → 16_3,13 at 218.298 GHz.

To generate the T_ex and N maps, we are using the spectral cube (spw0 and spw1) within a user-defined region of interest. The procedure performs the analysis only on pixels with a S/N above a given threshold in the moment 0 map, which in our case was S/N > 3. A preliminary rotational diagram analysis, for each pixel, is used to estimate the initial guesses for T_ex and log N, as detailed in Sect. A. The velocity (V) and linewidth (FWHM) parameters are extracted from the line velocity (moment 1) and velocity dispersion (moment 2) maps, respectively. Moment maps were generated with a built-in tool from MADCUBA. The AUTOFIT routine inside SLIM is then executed for each pixel to perform a full LTE fit. The MADCUBA procedure generates as output a SLIM product containing the individual pixel-by-pixel fits, and the log(N), T_ex, V, and FWHM maps. The only difference with the macro defined in Sect. A is that this procedure only takes into account the molecule that we are fitting, and not the contribution of all the other molecules. The provided maps have two channels: channel 1 contains the fitted values, while channel 2 stores the associated uncertainties.

Appendix B.1 Uncertainty N and T_ex maps

We present the relative uncertainty maps of the column density (N; Fig. B.1) and excitation temperature (T_ex; Fig. B.2) maps derived from the spectral fitting, in order to assess the reliability of the obtained physical parameters.

Fig. B.1 Relative error from the column density (N) maps derived for EG (top), GA (middle), and MF (bottom), respectively. Black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1. The red cross indicates the continuum peak position. The red ellipses outlined in white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown on the bottom left corner represents the synthesised beam of the interferometer.

Appendix B.2 Molecular transitions used for N and T_ex maps for ethylene glycol and glycolaldehyde

We present the molecular transitions used in the MADCUBA pixel-by-pixel fits for EG and GA.

Appendix C Zoomed-in line identification spectra

Figures C.1 and C.2 presents zoomed-in versions of the spw0 and spw1 spectra, respectively, shown in Sect. 5.1, aimed at better illustrating the weaker molecular lines used for line identification.

Appendix D Complete version of physical parameters table

Table D.1 presents the complete set of physical parameters derived from the SLIM analysis towards AG318−c9. Only a portion of this table is shown in the main text (Table 3).

Fig. B.2 Relative error from excitation temperature (Tex) maps derived for EG (top), GA (middle), and MF (bottom), respectively. Black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ and 260σ with σ = 0.28 mJy beam−1. The red cross indicates the continuum peak position. The red ellipses outlined in white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown on the bottom left corner represents the synthesised beam of the interferometer.

Fig. C.1 Same as Fig. 2 but with a zoomed-in y-axis to better highlight the weaker transitions and the quality of the spectral fit.

Fig. C.2 Same as Fig. 3 but with a zoomed-in y-axis to better highlight the weaker transitions and the quality of the spectral fit.

All Tables

All Figures

Fig. 1 ALMA map of the continuum emission at 1.38 mm of the AG318.9477−0.1960 clump (left) and a zoom-in on core 9 (right). The white contours correspond to 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ, with σ ∼ 0.28 mJy beam−1. The cyan ellipses surrounded by black show the position and size of each core. In particular, the size of core 9 is 0.″66 × 0.″57, with a position angle of PA=42°. The synthesised beam (0.″46 × 0.″37, PA=7°) is shown in the bottom left corner.

In the text

Fig. 2 Spectra of spw0 towards core 9 of the AG318.9477−00.1960 clump. The black histogram and its grey shadow are the observational spectrum, and the red curve indicates the LTE fit after taking into account the contribution of all the species (see Table 3 and its complete version in Appendix D.1). The horizontal dashed green line indicates the 5σ ≈ 1.35 K level, with σ being the noise (σ ∼ 0.28 K). Only transitions contributing at ≥ 5σ are labelled. The vertical dashed black lines indicate the position where ν was converged. U labels indicate unidentified transitions.

In the text

Fig. 3 Spectra of spw1 towards core 9 of the AG318.9477−00.1960 clump. The black histogram and its grey shadow are the observational spectrum, and the red curve indicates the LTE fit after taking into account the contribution of all the species (see Table 3 and its complete version in Appendix D.1). The horizontal dashed green line indicates the 5σ ≈ 1.35 K level, with σ being the noise (σ ∼ 0.28 K). Only transitions contributing at ≥5σ are labelled. The vertical dashed black lines indicate the position where ν was converged. U labels indicate unidentified transitions.

In the text

	Fig. 4 Column density comparison for all common species between AG318−c9 (x-axis) and G31 (y-axis). The solid black line indicates the 1:1 relation. The dark and light orange shaded regions represent deviations within 0.5 dex and 1.0 dex from equality, respectively.
In the text

Fig. 5 Top: comparison of molecular abundances (X = N/H2) towards G31 (blue bars) and AG318−c9 (orange bars), sorted by molecular groups and by the abundance of AG318−c9 inside each group. The species are grouped on the x-axis according to their detection status. Left: detected in both sources. Middle: detected in G31 alone. The orange triangles represent the corresponding upper limits for AG318−c9. Right: reported in AG318−c9. Bottom: logarithm of the ratio of the abundances of the two sources. HC3N* represents the group composed of HC3N, HC3N (v6 = 1), HC3N (v7 = 1), and HC3N (v7 = 2).

In the text

Fig. 6 Selection of GA transitions observed and fit (see Sect. 6.1). The black histogram and its grey shadow are the observational spectrum, and the red curve is the cumulative LTE fit considering all detected species. The blue curve is the fit of the individual transitions. The horizontal dashed green line indicates the 5σ level, corresponding to ∼1.35 K. Only transitions contributing at ≥5σ are labelled. The plots are sorted by decreasing intensity of the GA transitions (blue line). The upper left corner of each panel indicates the EUP of the corresponding GA transition. When multiple GA transitions are shown in the same panel, the EUP values are listed from top to bottom following the left-to-right order of the transitions.

In the text

	Fig. 7 Same as Fig. 6, but for a-EG, as described in Sect. 6.2.
In the text

	Fig. 8 Same as Fig. 6, but for MF, as described in Sect. 6.3.
In the text

Fig. 9 Integrated-intensity maps for EG (left), GA (middle), and MF (right). The white contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1, and the red cross shows its peak position. The red ellipses surrounded by white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown in the bottom left corner represents the synthesised beam of the interferometer.

In the text

Fig. 10 Column density (N) maps (top row) and excitation temperature (Tex) maps (bottom row) derived for EG (left), GA (middle), and MF (right). The corresponding relative errors are shown in Appendix B.1. The black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1, and the red cross shows its peak position. The red ellipses surrounded by white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown in the bottom left corner represents the synthesised beam of the interferometer.

In the text

	Fig. 11 Radial profiles derived from the column density (N; top) and excitation temperature (Tex; bottom) maps. The blue, red, and green curves correspond to EG, GA, and MF, respectively. The shaded areas indicate the measurement uncertainties (see Appendix B.1). The bar in the lower left corner indicates the radius corresponding to the geometric mean of the synthetic beam.
In the text

	Fig. A.1 Flowchart of the MADCUBA macro that automates the line identification process. The asterisk (*) indicates that physical parameters are kept free whenever possible, following the criteria described in the text.
In the text

Fig. B.1 Relative error from the column density (N) maps derived for EG (top), GA (middle), and MF (bottom), respectively. Black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ, and 260σ with σ = 0.28 mJy beam−1. The red cross indicates the continuum peak position. The red ellipses outlined in white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown on the bottom left corner represents the synthesised beam of the interferometer.

In the text

Fig. B.2 Relative error from excitation temperature (Tex) maps derived for EG (top), GA (middle), and MF (bottom), respectively. Black contours correspond to the continuum emission at 5σ, 10σ, 20σ, 30σ, 60σ, 90σ, 180σ, 220σ and 260σ with σ = 0.28 mJy beam−1. The red cross indicates the continuum peak position. The red ellipses outlined in white show the location and size of the dust continuum core AG318−c9 as estimated by Coletta et al. (2025). The hatched ellipse shown on the bottom left corner represents the synthesised beam of the interferometer.

In the text

	Fig. C.1 Same as Fig. 2 but with a zoomed-in y-axis to better highlight the weaker transitions and the quality of the spectral fit.
In the text

	Fig. C.2 Same as Fig. 3 but with a zoomed-in y-axis to better highlight the weaker transitions and the quality of the spectral fit.
In the text