© The Authors 2025

1 Introduction

The outer Galaxy (OG), namely the portion of the Galactic disk beyond the solar circle, is located at galactocentric distances (R_GC) approximately between 9 and 24 kpc, and it shows chemical properties significantly different from those of the inner Galaxy. Several chemical galactocentric gradient studies have been conducted, mapping the differences in elemental composition between the inner and outer Milky Way (e.g., C, N, O, Fe), using various tracers (i.e., H II regions, classical Cepheids, young OB stars, and planetary nebulae, among others). These works, focused on observational trends (e.g., Luck 2018; Bragança et al. 2019; Méndez-Delgado et al. 2022) and chemodynamical studies (Jacobson et al. 2016; Stanghellini & Haywood 2018; Palla et al. 2020; Santos-Peral et al. 2021), indicate that the overall metallicity in the OG is lower than the solar one (e.g., a factor of four lower at R_GC = 19 kpc; Shimonishi et al. 2021). The elemental abundances of carbon, oxygen, and nitrogen – the three most abundant elements in the Universe after hydrogen and helium – decrease as a function of R_GC (e.g., Esteban et al. 2017; Arellano-Córdova et al. 2020), similarly to other heavier elements, following the radial metallicity gradient. This observational result is linked to the traditional view of the OG as being poorer in molecules than the solar neighborhoods and the inner Galaxy. Moreover, based on the aforementioned scarcity of elements crucial for the formation of complex molecules and (rocky) planets, chemical evolution models excluded most of the OG from the so-called Galactic habitable zone (GHZ); that is, the region in the Milky Way where the probability for the development of habitable planets is highest: the optimal region is predicted to be an annulus of the Galactic disk from ~7 kpc extended up to R_GC ~ 9 kpc, with maximum habitability at ~8 kpc (Spitoni et al. 2014, 2017). This led to the consideration of the OG not being suitable for the formation of planetary systems in which Earth-like planets could be born and might be capable of sustaining life (Prantzos 2008; Ramírez et al. 2010).

Recent observational results challenged this prediction. First, the occurrence of Earth-like planets does not seem to depend on metallicity (e.g., Maliuk & Budaj 2020; Mulders et al. 2018). Second, the low rate of disruptive events in the OG suggests that it might be more favorable than the inner Galaxy for preserving life on potentially habitable planets (e.g., Piran & Jimenez 2014; Vukotić et al. 2016). Third, the detection of several complex organic molecules (COMs; e.g., CH₃OH, CH₃CCH, CH₃OCHO, CH₃OCH₃) of prebiotic interest in low-metallicity environments, such as star-forming regions in the OG (Shimonishi et al. 2021; Bernal et al. 2021; Fontani et al. 2022a, 2024; Koelemay & Ziurys 2025) and the Magellanic Clouds (Sewiło et al. 2018, 2022; Shimonishi et al. 2023; Golshan et al. 2024). Overall, these results suggest that environments prone to the presence of molecular precursors of life and Earth-like planets might be widespread in the Milky Way, including the OG.

The project CHEMical complexity in the outer Galaxy (CHEMOUT; Fontani et al. 2022a; 2022b, Colzi et al. 2022, Fontani et al. 2024), performed with the IRAM 30 m single-dish telescope, is aimed at deriving the molecular composition in 35 star-forming regions of the OG (R_GC ~ 9–24 kpc). The main results obtained in the previous works of the CHEMOUT project are (1) both the molecular detection rates and the molecular abundances of some organic species –namely HCO, H₂CO, and CH₃OH– do not vary significantly with R_GC (Fontani et al. 2022a; 2022b); (2) the ¹⁴N/¹⁵N ratio as a function of R_GC is consistent with chemical-evolution models in which ¹⁵N comes from novae (Colzi et al. 2022); (3) the high-angular-resolution maps of the target with the highest R_GC (WB89-670 at R_GC ~ 23.4 kpc) reveal chemical differentiation at core scales and, above all, a [C/H] abundance higher than expected at such a high galactocentric distance (Fontani et al. 2024).

The results obtained so far indicate a lack of decrease in the inspected molecular abundances, which, however, was limited to three species (HCO, H₂CO, and CH₃OH) and only 15 sources. In this paper, we analyze the molecules detected in more than 60% of the CHEMOUT sample according to Fontani et al. (2022a), namely HCN, HCO⁺, c-C₃H₂, H¹³CO⁺, HCO, and SO, and revise the molecular abundances of H₂CO and CH₃OH from Fontani et al. (2022b), using H₂ column densities from the new NIKA2 dust continuum maps. In Sect. 2, we describe the observational data used in this paper. In Sect. 3, we describe the methods implemented in the data analysis. The results are presented and discussed in Sect. 4 and 5. The conclusions and the outlooks are summarized in Sect. 6.

2 Observations

2.1 IRAM 30 m observations

This work is based on the observations done for the CHEMOUT project, described in Papers I and II (Fontani et al. 2022a; 2022b), performed with the Institut de RadioAstronomie Millimétrique (IRAM) 30 m telescope. The sample of sources analyzed in this paper is composed of the 35 star-forming regions, presented in Table 1 of Fontani et al. (2022a). In some spectra, two velocity features are identified, and the sources listed in Table A.1 in this paper include both velocity components. In the analysis conducted, we used spectral windows at 3 and 2 mm, collected with the EMIR receiver and described in Table 1. The spectral windows of the project 004-18 (see Table 1) have been optimized to observe, at 3 mm, the J = 1–0 transitions of ¹³C and ¹⁵N isotopologues of HCN and HNC and, at 2 mm, the J = 2–1 line of DNC. The spectral windows observed in the project 042-21, at 3 and 2 mm, have been designed to observe several methanol and formaldehyde lines. For further details (e.g., pointing and focus, and weather conditions), we refer the reader to Fontani et al. (2022a; 2022b).

2.2 NIKA2 bolometer observations

Millimeter continuum emission of 31 out of the 35 CHEMOUT targets was observed with the New IRAM KID Arrays 2 (NIKA2; Perotto et al. 2020) bolometer of the IRAM 30 m telescope in Pico Veleta, Spain. Observations were taken during the winter 2021 and 2023 pool seasons (February). The observations cover an almost circular region with a radius of ~9 arcmin, centered at the positions where the IRAM 30 m spectra have been observed (Table A.1). The dual-band capability of NIKA2 allowed us to simultaneously observe the dust’s thermal continuum emission at 2.0 and 1.2 mm (150 GHz and 260 GHz, respectively). The sensitivity achieved was ~1.5 mJy/beam at 1.2 mm and 0.5 mJy/beam 2.0 mm. A more complete description of the observations and a detailed analysis of the maps will be presented in a forthcoming paper of the CHEMOUT series (Fontani et al., in prep.,). Here, we have used the maps only to derive the H₂ column densities (Sect. 3.2) to compute the molecular fractional abundances with respect to H₂ (Sect. 4.2).

3 Analysis methods

3.1 Line fitting

The spectra were fit using the MAdrid Data CUBe Analysis (MADCUBA¹; Martín et al. 2019) software. The transitions were identified using the Spectral Line Identification and LTE Modelling (SLIM) tool of MADCUBA, which makes use of the Jet Propulsion Laboratory (JPL; Pickett et al. 1998) and Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2001) catalogs. The lines were fit using the AUTOFIT function of SLIM. This function produces the Gaussian synthetic spectrum that best matches the observed spectrum, assuming local thermodynamic equilibrium (LTE) conditions and the same fitting input parameters for all transitions of a given molecular species, which are total column density, N_tot; excitation temperature, T_ex; radial systemic velocity of the source, V_LSR; full width at half-maximum of the line, FWHM; and angular size of the emission, θ_S.

For the source size, we assumed that the molecular emission fills the telescope beam, except for methanol and formaldehyde, where an estimation of the emission size was made in Fontani et al. (2022b) and updated in this paper as discussed in Sect. 4.1.1. Such an assumption is justified by the spatial extent of the emission of the dust (see Figs. B.1, B.2, B.3), as well as by the updated H₂CO and CH₃OH emission sizes shown in Table C.2. The excitation temperature, whenever an estimate was not possible through the fit, has been fixed to the one shown in Table A.1, which was estimated from the methanol lines for the 15 targets observed in Fontani et al. (2022b). For the remaining sources, we adopted an average T_ex of 10.7 K, derived from those sources for which methanol transitions were available. This approximation was necessary because it was not possible to estimate the temperature individually for each source based on the observed transitions. To evaluate the impact of this assumption, we computed N_tot in the range of excitation temperatures derived from methanol (7 ≤ T_ex ≤ 15 K; see Table A.1). In this T_ex range, the densities vary from the values shown in Tables C.1 and C.2 by a factor of two at most for the species studied (lowest for HCN, HCO⁺, and H¹³CO⁺; highest for SO).

As a detection criterion, we considered a line as detected if the signal-to-noise ratio (S/N) is greater than three (an estimation of the rms noise is given in Fontani et al. 2022a), while the detection is considered as tentative for 2.5 ≤ S/N ≤ 3. If a line is undetected (e.g., S/N ≤ 2.5), MADCUBA provides a 3σ upper limit estimate to the line’s integrated intensity from Eq. (29) of Martín et al. (2019). The upper limit on the integrated intensity is converted into a column-density upper limit using the SLIM parameters T_ex and θ_S. For these estimates, the average FWHM for each species across the sample was adopted (Sect. 5.2). This approach assumes that different molecules trace gases with different physical, kinematic, and excitation conditions.

The total column densities derived from transitions observed at frequencies corresponding to beam sizes different from 28″ (i.e., the largest beam size among the IRAM 30 m telescope spectral windows) have been rescaled by a factor of $(\frac{{\theta }_{\text{beam }}}{{28}^{\mathrm{\prime }\mathrm{\prime }}}{)}^2$. For H₂CO and CH₃OH, the column densities were estimated assuming a specific angular emission size (see Table C.2, and Sect. 4.1.1). In cases where the emission size exceeds 28″, no correction was applied, under the assumption of uniform emission within the beam. Conversely, when the source size is smaller than 28″, the same scaling factor was applied to account for beam dilution effects. The estimation of the molecular fractional abundances has been performed using

$$X_{\mathrm{m}\mathrm{o}\mathrm{l}}=\frac{N_{\mathrm{t}\mathrm{o}\mathrm{t}}}{N_{{\mathrm{H}}_2}},$$(1)

where N_tot is the total molecular column density, and N_H2 is the H₂ column density (see Sect. 3.2).

3.2 Continuum analysis

For this paper, the analysis of these maps is limited to the extraction of the flux density from an angular region of 28 ″ centered on the pointing centers of the molecular spectroscopic data. The dimension of the extracted region has been decided to match the larger beam size of IRAM 30 m in the studied spectral windows (Table 1). A complete description and a deeper analysis of the NIKA2 observations will be given in a forthcoming paper (Fontani et al. in prep.).

The column densities of H₂, N_H2 were derived from the dust-continuum-emission maps. In the scales (and densities) traced by our observations, the thermal dust emission is expected to be optically thin. Under this assumption, we can use the following equation to derive N_H2 from the continuum flux density (e.g., Battersby et al. 2011):

$$N_{{\mathrm{H}}_2}=\frac{\gamma F_{\nu }}{{\mathrm{\Omega }}_{\mathrm{S}}{\kappa }_{\nu }{B}_{\nu }(T_{\mathrm{d}}){\mu }_{{\mathrm{H}}_2}{m}_{\mathrm{H}}},$$(2)

where γ is the gas-to-dust ratio, F_ν is the continuum flux density at frequency ν, Ω_S is the source solid angle, κ_ν is the dust-mass opacity, T_d is the dust temperature, and μ_H2 is the mean molecular weight for molecular hydrogen that is assumed to be 2.8 (Kauffmann et al. 2008). The H₂ column density was estimated with both a constant gas-to-dust ratio (γ = 100) and a gas-to-dust ratio increasing with distance from the Galactic center from Giannetti et al. (2017) (γ ~ 3000 at R_GC = 23.4 kpc). T_d was derived under the assumption of thermal coupling between gas and dust, adopting the excitation temperatures obtained from the molecular line analysis of CH₃OH (Table A.1) or, for those sources where methanol is not detected, adopting 10.7 K (the average T_ex from CH₃OH). The dust-mass opacity was derived at the observed frequencies from the relation κ_ν/κ_ν0 = (ν/ν₀)^β, assuming κ_ν0 = 0.899 at ν₀ = 230 GHz (Ossenkopf & Henning 1994). The spectral index, β, was estimated, for each source, from the ratio between the flux densities at the two frequencies observed (150 and 260 GHz) with NIKA2 (Appendix B).

4 Results

4.1 Detection summary

A comprehensive molecular identification was performed over the whole sample to identify the species detected in the greatest number of targets. This allowed us to analyze abundance gradients as a function of galactocentric distance with the highest possible statistics. Over 30 molecular species have been detected, and the most frequently detected ones are HCN (39/39), HCO⁺ (39/39), c-C₃H₂ (38/39), H¹³CO⁺ (30/39), HCO (29/39), and SO (25/39). The detection statistics over these species are already listed in Fontani et al. (2022a). We also included H₂CO (18/18) and CH₃OH (17/18) among the species considered. The analysis of the HCO, H₂CO, and CH₃OH gradients was conducted by Fontani et al. (2022b); however, here, we present updated trends based both on the revised fits performed with MADCUBA (Sect. 4.1.1) and on our new H₂ column-density estimates from NIKA2 observations. The spectroscopic parameters of the transitions detected in these molecules are displayed in Table 3. In Appendix C, the best-fit parameters obtained from the analysis with MADCUBA, and the spectra of the fit transitions for HCN, H¹³CO⁺, and SO are shown in detail. The spectra of c-C₃H₂ and HCO⁺ are shown in Fontani et al. (2022a), and the ones of HCO, H₂CO, and CH₃OH in Fontani et al. (2022b). The results regarding the molecular column densities (e.g., their galactocentric gradients) are provided in Appendix D.

In Table A.1, the N_H2 values calculated assuming a gas-to-dust ratio from Giannetti et al. (2017) are presented. These values are used hereafter for the analysis. All the data, results, and discussions concerning fractional abundances with respect to H₂ are based on these estimates, rather than on those obtained assuming a fixed γ = 100. The NIKA2 dust continuum maps and the galactocentric gradients for H₂ column densities are shown in Appendix B.

Six sources (19383+2711, 19571+3113, WB89-006, WB89-060, WB89-080, and WB89-380) show multiple velocity features also reported in Fontani et al. (2022a). Four of them (19383+2711, 19571+3113, WB89-006, and WB89-380) exhibit double-peak emission in most of the observed lines. Because the peak velocities of these two features are consistent in the different transitions in which they are detected, they probably arise from two clumps within the telescope beam that have an angular separation smaller than the beam itself. This is confirmed by the NIKA2 maps (Figs. B.1–B.3), for 19383+2711 and 19571+3113, showing two different emission components in the beam of the IRAM 30 m telescope. WB89-006 and WB89-380 do not show a clear double peak in the continuum maps, but the two clumps could be non-resolved with NIKA2 or aligned along the line of sight. In principle, each of these clumps may have distinct physical and chemical properties and will therefore be analyzed separately in the following sections. For these, the peak velocities were estimated using the V_LSR of the c-C₃H₂ J_Ka,Kc = 2_1,2−1_0,1 line, and we assumed the same Galactic longitude, ℓ, for both clumps, to estimate their galactocentric radius, R_GC, and heliocentric distance, d. The remaining two sources (WB890-60 and WB89-080) are characterized by a double peak, but only in some more intense and potentially optically thick lines (e.g., HCN and HCO⁺) and do not show coherence in the peak velocities. Consequently, in these cases, the double-peaked line shapes are more likely due to high-optical-depth effects, such as (self-)absorption, rather than multiple velocity features. These have thus been analyzed as single sources, and the absorbed lines discarded from the analysis. Moreover, some transitions showed non-Gaussian, high-velocity wings. We decided to limit our analysis to the narrower component of the molecular emissions related to the line emission around the intensity peak, and thus arising from the bulk emission of the cores.

4.1.1 Formaldehyde and methanol updated fits

The fits for H₂CO and CH₃OH have been revised due to an incorrect assumption in Fontani et al. (2022b), which originated from a conversion between $T_{\mathrm{A}}^{\ast }$ and T_MB. Since the beam efficiency η_eff of the telescope was not explicitly specified in Bernal et al. (2021), the values were taken from Appendix C.3 of the ARO 12m user manual; however, outdated efficiencies were inadvertently used. Updated efficiencies² have now been adopted to estimate the T_MB of these two species. For H₂CO at 2 mm we divided $T_{\mathrm{A}}^{\ast }$ by 0.81, and for CH₃OH at 3 mm we divided by 0.91. The source sizes for H₂CO and CH₃OH were re-estimated using Eq. (3) in Fontani et al. (2002). Consequently, the MADCUBA fits were recalculated, and the resulting parameters are reported in Appendix C. The excitation temperatures of CH₃OH are listed in Table A.1.

4.2 Fractional abundances

The molecular fractional abundances, X_mol, with respect to H₂ were computed through Eq. (1), with the molecular column densities from Tables C.1 and C.2, and the N_H2 from Table A.1. The values of X_mol are in ranges presented in Table 2. All the values are shown in Table 4. These values are estimated for 31 star-forming regions of the sample, i.e., the ones observed in the dust emission, for which we have the H₂ column densities. For sources exhibiting a double clump structure (see Sect. 4.1), the fractional abundances were obtained by summing the molecular column densities of the two components, and the galactocentric distances were determined as the average R_GC of the individual clumps.

The plots showing X_mol against R_GC for each molecule are shown in Fig. 1. For comparison, the fractional abundances calculated using a gas-to-dust ratio γ = 100 are also shown. A linear fit is performed for each plot, and the Pearson correlation coefficient³, ρ, is calculated to inspect possible (anti-)correlations. The linear fit applied to the data of fractional abundances, calculated using a gas-to-dust ratio from Giannetti et al. (2017), gives:

$$[\mathrm{H}\mathrm{C}\mathrm{N}]=(-0.105\pm 0.013)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(8.79\pm 0.20)\mathrm{d}\mathrm{e}\mathrm{x},$$(3)

$$\left[{\mathrm{H}\mathrm{C}\mathrm{O}}^{+}\right]=(-0.078\pm 0.018)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(9.50\pm 0.28)\text{ dex},$$(4)

$$[c-{\mathrm{C}}_3{\mathrm{H}}_2]=(-0.058\pm 0.021)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(9.95\pm 0.32)\text{dex},$$(5)

$$\left[{\mathrm{H}}^{13}{\mathrm{C}\mathrm{O}}^{+}\right]=(-0.039\pm 0.025)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(11.32\pm 0.38)\mathrm{d}\mathrm{e}\mathrm{x},$$(6)

$$[\mathrm{H}\mathrm{C}\mathrm{O}]=(-0.057\pm 0.026)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(9.83\pm 0.41)\mathrm{d}\mathrm{e}\mathrm{x},$$(7)

$$[\mathrm{S}\mathrm{O}]=(-0.080\pm 0.023)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(8.98\pm 0.36)\mathrm{d}\mathrm{e}\mathrm{x},$$(8)

$$\left[{\mathrm{H}}_2\mathrm{C}\mathrm{O}\right]=(-0.103\pm 0.050)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(8.42\pm 0.75)\mathrm{d}\mathrm{e}\mathrm{x},$$(9)

$$\left[{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\right]=(-0.146\pm 0.054)\mathrm{d}\mathrm{e}\mathrm{x}{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\cdot R_{\mathrm{G}\mathrm{C}}-(7.67\pm 0.82)\text{ dex.}$$(10)

4.3 Line widths

From the fit of the lines, we also obtained the full widths at half maximum (FWHMs). The measured FWHMs for the molecules studied are in the ranges shown in Table 2. The mean values of line width estimated for each molecule have been used in the estimation of the upper limits of the column densities, as an educated guess (2.1 ± 0.6 km s⁻¹ for HCN; 1.8 ± 0.9 km s⁻¹ for HCO⁺; 2.1 ± 1.0 km s⁻¹ for c-C₃H₂; 2.1 ± 1.0 km s⁻¹ for H¹³CO⁺; 2.2 ± 1.0 km s⁻¹ for HCO; 2.0 ± 0.8 km s⁻¹ for SO; 2.3 ± 0.7 km s⁻¹ for H₂CO; 1.9 ± 0.8 km s⁻¹ for CH₃OH).

The plots of FWHM of HCN, HCO⁺, c-C₃H₂, H¹³CO⁺, HCO, SO, H₂CO, and CH₃OH as a function of the galactocentric distance are shown in Fig. 2. A linear fit is performed, and the Pearson correlation coefficient, ρ, is estimated. We find that, for all species, the line width shows an overall negative trend with R_GC, along with a weak negative correlation (ρ ≲ −0.2.

5 Discussion

5.1 Galactocentric gradients of fractional abundances

The trends obtained for the molecular abundance (with respect to H₂) were compared with the galactocentric abundance gradients of carbon (both ¹²C and ¹³C), oxygen, nitrogen, and sulfur. These are constituting elements of the studied species, and the gradients were determined in Méndez-Delgado et al. (2022) through the analysis of the spectra of Galactic H II regions. The galactocentric gradient of ¹³C was estimated by coupling the carbon gradient with the galactocentric trend of the ¹³C/¹²C ratio from Yan et al. (2023). Among 8.34 kpc, the galactocentric distance of the Sun, and 24 kpc, which is the distance of the most distant source in our sample, the abundance ratios [C/H], [¹³C/H], [N/H], [O/H], and [S/H] decrease by factors of approximately ~16, ~36, ~10, ~5, and ~4, respectively. Extrapolating the gradients derived in the OG (Eqs. (3)–(8)) to the Solar circle, over the same range of R_GC, we find that the abundance of the studied species, i.e., HCN, HCO⁺, c-C₃H₂, H¹³CO⁺, HCO, SO, H₂CO, and CH₃OH decrease by a factor of ~44, ~17, ~8, ~4, ~8, ~18, ~41, and ~194.

The decreasing rates of HCN, HCO⁺, c-C₃H₂, HCO, H₂CO, and CH₃OH are consistent, within the error bars, with that of [C/H]. This implies that, considering the dispersion and uncertainties in estimating these quantities, the fractional abundances of these carbon-bearing species scale with the elemental abundance of carbon ([C/H]) between 8.34 kpc and 24 kpc. This suggests that the efficiency in forming these molecules, which scales with the availability of the parent element, C, is at least as high as in the local Galaxy.

In contrast, the decreasing rate of SO appears steeper than that of [S/H]. For SO, the higher decreasing rate with respect to [S/H] indicates a lower efficiency of the formation of these molecules in the OG with respect to the local Galaxy. We speculate that a steeper decreasing trend for SO, which is the best shock tracer among those that we analyzed in this work, could indicate that, in the OG, outflows have a reduced influence on the gas composition compared to the local Galaxy. This result would be consistent with the decreasing line width with R_GC (Sect. 4.3). Another shock tracer is CH₃OH, which, although its decreasing trend is marginally consistent with that of carbon, exhibits the most pronounced reduction among the species investigated in the OG. The steep decreasing trend of CH₃OH could also arise from a less efficient grain-surface chemistry in the OG. The results concerning HCO, H₂CO, and CH₃OH, obtained with the updated abundances (and distances), confirm the ones of Fontani et al. (2022b). In Fontani et al. (2024), a correlation between the abundances of CH₃OH and SO was found, which is confirmed over the entire CHEMOUT sample, although no correlation between the line widths of these two species is observed, hinting that they are not tracing the same gas components.

The decreasing rate of H¹³CO⁺, although lower than that of [¹³C/H], aligns well with the decreasing rate observed for oxygen. This correlation is not observed for HCO⁺, whose decreasing rate is consistent with carbon rather than oxygen. This discrepancy could suggest a possible revision of the galactocentric gradient of the ¹³C/¹²C ratio, with a shallower trend at high R_GC, in agreement with the Galactic chemical evolution discussed in Colzi et al. (2022). Moreover, the detailed analysis of the ¹³C/¹²C ratios from the HCO⁺and H¹³CO⁺ lines will be the subject of a forthcoming paper (Colzi et al. in prep.). Therefore, the results regarding c-C₃H₂, HCO, and especially H¹³CO⁺ lead to the puzzling conclusion of observing more carbon-bearing molecules than the available carbon atoms in the OG, their decreasing rates being lower than that of [C/H] (or [¹³C/H]).

Fig. 1 Galactocentric gradients of fractional abundances with respect to H2, Xmol. The plot shows the trends for HCN (a), HCO+(b), c-C3H2 (c), H13CO+ (d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of the galactocentric radius (RGC). The light blue data represent the abundances calculated using a constant gas-to-dust ratio (γ = 100), while the dark blue data illustrate those estimated using a non-constant gas-to-dust ratio (γ = γ(RGC) from Giannetti et al. 2017). For both datasets, the linear regression results are shown as the light blue and dark blue lines, respectively. The 1σ error bars over the slope of the gradients are displayed for each fit. The upper-limit values are represented with triangles. The gradients of the elemental abundances of carbon (C and 13C), nitrogen (N), oxygen (O), and sulfur (S), as reported by Méndez-Delgado et al. (2022), are plotted as dashed lines. The product of the parent elements of the species is represented, in each subplot, by a dash–dotted black line. All elemental trends are plotted with a vertical offset to align with the starting point of the molecular gradients derived using γ(RGC), to facilitate comparison of their slopes with the linear fit of the molecular abundances. In the upper-right side of each subplot, the Pearson correlation coefficient, ρ, is shown (estimated only for the abundances estimated using the non-constant gas-to-dust ratio).

5.2 Galactocentric gradients of line widths

Considering the low temperature (Table A.1) of the sources, the thermal broadening of the emission lines can be neglected (~0.1–0.2 km s⁻¹); therefore, the line width is a good tracer of the gas turbulence of the sources. The decreasing galactocentric trends observed for the FWHM of the studied molecules suggest that the sources tend to have lower turbulence as the distance from the Galactic center increases. Such a result would be in line with an environment becoming more and more quiescent at increasing R_GC. The nearest source to the Sun, WB89-002, is the only outlier that deviates from these trends.

This result is also consistent with the study of Lin et al. (2025), where the velocity dispersion of low-J ¹³CO lines has been studied toward molecular clouds in the low-metallicity part of the Galactic disk. A decreasing trend of the line width is observed with R_GC, and the sub-virial dynamics observed in the metal-poor clouds hint at how the gas turbulence is not sufficient to counterbalance the cloud’s self-gravity. As stated in Lin et al. (2025), in low-metallicity environments, the magnetic field may play a more dominant role, replacing gas turbulence in sustaining the clouds.

Fig. 2 Galactocentric gradients of line widths, FWHM. The plot shows the trends for HCN (a), HCO+ (b), c-C3H2 (c), H13CO+(d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of the galactocentric radius (RGC). The linear fit results are shown as dark blue lines. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

5.3 Galactocentric gradients of molecular ratios

Studies on the formation pathways of HCO⁺, H¹³CO⁺ (Klemperer 1970; Millar et al. 1991), HCO, H₂CO, and CH₃OH (Watanabe & Kouchi 2002; Fuchs et al. 2009; Minissale et al. 2016; Rivilla et al. 2019; Jiménez-Serra et al. 2025) have shown that their formation is linked to the presence of CO. In contrast, c-C₃H₂ is formed from “free” carbon (not bound in CO) (Thaddeus et al. 1985; Sipilä et al. 2016). Moreover, it is known that [C/O] decreases with distance from the Galactic center (Méndez-Delgado et al. 2022). Therefore, we expect to have less free carbon in the OG to produce c-C₃H₂, because most of it would be locked in CO. This would lead to a prediction of an increasing trend for the CO/c-C₃H₂ ratio.

Considering the lack of direct CO measurements for our sample, we can use species containing C and O that are formed from CO, such as HCO⁺, H¹³CO⁺, HCO, H₂CO, and CH₃OH. In Fig. 3, the galactocentric gradients of these ratios are presented. No clear trends arise from the data, which are also characterized by |ρ| < 0.2 (except for methanol, for which we have ρ ~ −0.3). The ratio between H¹³CO⁺ and c-C₃H₂ is the only one presenting a marginally increasing trend (in agreement with the formation pathways and the C/O abundance); however, it is too shallow to be considered reliable.

Fig. 3 Galactocentric gradients of molecular ratios. The plot shows the trend of the ratios between CO molecule products (i.e., HCO+, H13CO+, HCO, CH3OH, and H2CO) and c-C3H2, as a function of galactocentric radius (RGC). The 1σ error bars over the slope of the gradients are displayed for each molecular fit. The upper limit values are represented with triangles. The linear fit results are shown as the dark blue lines. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

5.4 Caveats and uncertainties

The comparison between observations and chemical models by Fontani et al. (2024) for WB89-670 suggests that the methanol could be sub-thermally excited: the kinetic temperatures predicted by chemical modeling are in the ~20–40 K range, while the excitation temperatures are on the order of ~10 K. Therefore, an important caveat arises from the assumed T_ex for the species for which we cannot estimate it directly from the observations.

Furthermore, the sources are located in a range of heliocentric distances, d, between 3 and 15.1 kpc, which implies that the linear scale probed is not the same in all targets. Fontani et al. (2024) demonstrated the presence of a high differentiation at smaller angular scales in WB89-670, and hence the average column densities derived in this work refer to linear scales with (eventually) different average conditions. However, this caveat is typical of source surveys, and in any case, we compare emissions arising from core envelopes in all targets.

Some results may also be influenced by low statistics, particularly within certain R_GC ranges, where the data characterizing a few sources have a bigger impact on the gradients. This concern is more significant for formaldehyde and methanol, given their even lower sampling rates (only 18 molecular clumps of 39 are observed in H₂CO and CH₃OH spectral windows), over a smaller R_GC range (Fontani et al. 2022b).

Column densities of species such as HCN or HCO⁺might be affected by large optical-depth effects in the observed lines, which could impact the derived Galactocentric gradients. However, for HCN, the column densities were corrected by the line opacity (Table C.1) obtained from the hyperfine structure of the observed transition, and both HCN and HCO⁺ were observed also in their less abundant ¹³C-bearing isotopologues (see Appendix E for the gradient of H¹³CN). Therefore, a gradient can also be inferred for these species using these isotopologues, expected to be optically thin.

Moreover, in the comparison between molecular and elemental trends, a further caveat comes from the fact that the elemental abundance gradients are extrapolated in the extreme OG (R_GC ≳ 16 kpc), where observations are sparse and thus the gradients are poorly constrained, and could be different from those at R_GC ≤ 16 kpc. Indeed, recent studies have reported evidence of a flattening of the metallicity gradient at large galactocentric radii, which could influence the observed molecular abundances and potentially explain the unexpectedly high chemical complexity in the outskirts of the Galaxy (Donor et al. 2020; Myers et al. 2022; Magrini et al. 2023; Da Silva et al. 2023).

6 Conclusions and outlooks

In this work, we studied the chemical composition of 35 star-forming regions in the OG (9 ≤ R_GC ≤ 24 kpc), using IRAM 30 m telescope data from the CHEMOUT project. We identified 39 molecular clumps and analyzed the molecular species detected in the greatest number of targets across our sample; i.e., HCN (39/39), HCO⁺ (39/39), c-C₃H₂ (38/39), H¹³CO⁺ (30/39), HCO (29/39), and SO (25/39). We estimated their column densities and calculated the fractional abundances of these species with respect to H₂. We also updated the abundances of H₂CO (18/18) and CH₃OH (17/18) from Fontani et al. (2022b), using the H₂ column densities obtained from the new NIKA2 continuum maps.

The fractional abundances of the studied species with respect to H₂ estimated assuming a gas-to-dust ratio dependent on R_GC (Giannetti et al. 2017) show a decreasing trend with the galactocentric radius. For HCN, HCO⁺, c-C₃H₂, HCO, SO, H₂CO, and CH₃OH the abundances scale as the elemental abundance of carbon ([C/H]) between the solar circle and ~24 kpc. This indicates that the efficiency in the formation of these molecules, which scales with the availability of the parent element (i.e., C), is at least as high as that found in the local Galaxy. Instead, the decreasing rate of SO is higher than that of [S/H], indicating a lower efficiency in the formation of this species in the OG. Conversely, the decreasing rate of H¹³CO⁺presents the puzzling result that it is lower than that of [¹³C/H], resulting in a higher formation efficiency in the OG compared to the inner Galaxy.

The galactocentric gradients of the abundances indicate how the metallicity of the environment, which scales as effectively with R_GC, does not affect the formation efficiency of most of the molecules studied, even for so-called prebiotic molecules such as formaldehyde. Metallicity appears to act mostly as a scaling factor for the molecular abundances relative to the elemental ones. The role of metallicity in the efficiency of molecule formation may become evident when studying molecules more complex than methanol, and/or it could be hidden by other environmental factors such as the radiation field, gas temperature, or cosmic ionization rate, which may counterbalance its effect. These results confirm (after the previous papers of the CHEMOUT series) that the presence of organic molecules and tracers of protostellar activity is ubiquitous in the low-metallicity environment of the OG.

The next step for the CHEMOUT project is to infer a more precise kinetic temperature of the sources to better estimate the column densities from the line analysis. Determining this parameter will also help to reduce the space parameters on a chemical modeling analysis (e.g., using UCLCHEM; Holdship et al. 2017), which is crucial to inferring the environmental differences of the star-forming regions in the different areas of the Galaxy and understanding the parameters driving the chemistry (Vermariën et al. 2025). The comparison between high-resolution observations and chemical models is crucial to disentangling the several contributions to these formation processes.

It is furthermore crucial to enhance the statistics acquired as much as we can. In fact, given the scarcity of star-forming regions in the extreme OG, the data we have for R_GC ≳ 20 kpc weigh significantly on the overall trends. Increasing the observed sample is essential to constraining the gradients obtained and securing more precise results over the ongoing chemistry in this region of the Milky Way.

Data availability

Figures C.1, C.2, C.3, and F.1 are available on Zenodo (DOI: https://doi.org/10.5281/zenodo.17187165).

Acknowledgements

We thank the anonymous referee for providing valuable comments and suggestions that improved the scientific content of this paper. F.F. is grateful to the IRAM 30 m staff for their precious help during the observations. L.C. and V.M.R acknowledge support from grant no. 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 Cientificas (CSIC) and the Centro de Astrobiología (CAB) through the project 20225 AT015 (Proyectos intramurales especiales del CSIC); and from the grant CNS2023-144464 funded by MICIU/AEI/10.13039/501100011033 and by “European Union NextGenerationEU/PRTR”. The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement nos. 730562 and 101004719 (ORP: https://www.orp-h2020.eu). SV and GV acknowledge support from the European Research Council (ERC) Advanced grant MOPPEX 833460. 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 PID2023-146675NB-I00 (MCI-AEI-FEDER, UE).

References

Appendix A The CHEMOUT sample

Table A.1 presents the observed sample together with its main properties.

Appendix B Millimeter continuum maps and H₂ column densities

In this appendix, we show the millimeter continuum maps obtained with the IRAM 30 m bolometer NIKA2 (Figs. B.1-B.3). We observed 31 out of the 35 CHEMOUT targets. A detailed analysis of the maps will be presented in a forthcoming paper (Fontani et al. in prep.).

In Table B.1, all the parameters derived from the continuum maps, and used to estimate the H₂ column densities, are presented. The fluxes have been extracted from the NIKA2 maps using a 28 ″ beam, matching the widest beam size of the spectral data centered at the same position from which we observed the IRAM 30 m spectra. The spectral index β is estimated using the following equation:

$$\beta =\frac{\log \left(\frac{F_{{\nu }_2}}{F_{{\nu }_1}}\frac{B_{{\nu }_1}\left(T_{\mathrm{d}}\right)}{B_{{\nu }_2}\left(T_{\mathrm{d}}\right)}\right)}{\log \frac{{\nu }_2}{{\nu }_1}}$$(B.1)

where ν₂ and v₁, are 150 and 260 GHz (i.e., the frequencies corresponding to 2.0 and 1.2 mm), and T_d is the assumed dust temperature (see Sect. 3.2). In Fig. B.4, the Galactocentric trends of H₂ are shown. The dataset estimated with a non-constant gas-to-dust ratio γ = γ(R_GC) (Giannetti et al. 2017) exhibits a positive trend with increasing values of R_GC, and a weak, positive correlation (ρ = 0.19), which is due to the increasing γ (Sect. 3.2).

Fig. B.1 Preliminary reduced NIKA2 images (color-scale in mJy/beam) at 2.0 mm (150 GHz) and 1.2 mm (260 GHz) of 12 out of the 31 CHEMOUT sources observed (Fontani et al. in prep.). The contours start from the 5σ rms level in each map. The NIKA2 angular resolution is shown in the bottom-left corner of each frame. The black circles on the maps represent the 28 ″ beam, centered on the position targeted for the molecular line observations analyzed in this paper.

Fig. B.2 Same as Fig. B.1 for 12 additional targets.

Fig. B.3 Same as Fig. B.1 for 7 additional targets.

Fig. B.4 Gradients of column densities of H2, NH2, as a function of Galactocentric radius (RGC). The light blue dataset represents the column densities estimated using a constant gas-to-dust ratio (γ = 100), while the dark blue dataset shows the column densities calculated with a gas-to-dust ratio function of RGC from Giannetti et al. (2017). Both datasets have been fitted with a linear regression. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper-right side of the plot, the Pearson correlation coefficient, ρ, is shown (only for the dark blue dataset).

Appendix C Line-fitting results

In this appendix, we present the results of the line fit from MADCUBA. Tables C.1 and C.2 present the best-fit parameters from MADCUBA for HCN, HCO⁺, c-C₃H₂, H¹³CO⁺, HCO, and SO. The fit results for H₂CO and CH₃OH have already been presented in Fontani et al. (2022b). The upper-limits values are also estimated for the undetected transitions. For the (self-) absorbed lines, no best-fit parameters are given.

In Figs. C.1, C.2, and C.3 (available on Zenodo), the observed spectra of the whole sample are shown with the Gaussian fits performed with MADCUBA for HCN, H¹³CO⁺, and SO. The plots showing the spectra of the remaining species are shown in Fontani et al. (2022a) (c-C₃H₂, HCO⁺) and in Fontani et al. (2022b) (HCO, H₂ CO, CH₃ OH).

Appendix D Molecular column densities

The derived column densities of the 39 molecular clumps for the studied molecules are within the ranges shown in Table 2. The trends of N_tot as a function of the Galactocentric radius for all molecules are shown in Fig. D.1. A linear fit is performed for each plot, and the Pearson correlation coefficient, ρ, is calculated. As illustrated, almost constant trends and no- or weak correlations (i.e. |ρ| ≲ 0.2) are found for HCN, HCO⁺, c-C₃H₂, H¹³CO⁺, HCO, SO, and H₂CO. A higher anti-correlation is found between CH₃OH and the Galactocentric distance.

Fig. D.1 Galactocentric gradients of column densities, Ntot. The plots show the trends for HCN (a), HCO+ (b), c-C3H2 (c), H13CO+ (d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of Galactocentric radius (RGC). The upper limit values are represented with triangles. The dark blue lines represent the linear regressions computed over the datasets. In the upper-right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

Appendix E Galactocentric gradients of H¹³CN

The Galactocentric gradients of H¹³CN fractional abundances are presented in Fig. E.1. This isotologue of HCN is characterized by a decreasing rate of ~114, between the Solar circle and 24 kpc, higher than that of ¹³C (and N). A steeper decrease of this species suggests a lower formation efficiency in the OG, with respect to the inner Galaxy. It is important to highlight the contrasting behavior between HCN and HCO⁺, compared to their ¹³C isotopologues: while HCO⁺ shows a steeper decrease than H¹³CO⁺, HCN displays a shallower gradient compared to its ¹³C-bearing isotopologue.

Fig. E.1 Galactocentric gradients of fractional abundances w.r.t. H2, Xmol, of H13CN. The light blue data represent the abundances calculated using a constant gas-to-dust ratio (γ = 100), while the dark blue data illustrate those estimated using a non-constant gas-to-dust ratio (γ = γ(RGC) from Giannetti et al. 2017). For both datasets, the linear regression results are shown as the light blue and dark blue lines, respectively. The 1σ error bars over the slope of the gradients are displayed for each molecular fit. The upper limit values are represented with triangles. The gradients of the elemental abundances of carbon (13C) and nitrogen (N), as reported by Méndez-Delgado et al. (2022), are plotted as dashed lines. The product of the parent elements of the species is represented by a black dash–dot line. All the elemental trends are plotted shifted to the start of the molecular gradients estimated with γ(RGC) for reference to allow comparison of their slopes with the linear fit of the molecular abundances. In the upper-right side of each plot, the Pearson correlation coefficient, ρ, is shown (estimated only for the abundances estimated using the non-constant gas-to-dust ratio).

Appendix F Correlation between fractional abundances and between line widths

Figure F.1 (available on Zenodo) shows plots to highlight possible correlations between fractional abundances and between line widths of different species. In particular, we analyzed the correlation between HCO⁺, H¹³CO⁺, HCO, H₂CO, CH₃OH (i.e., CO-product molecules), and c-C₃H₂, produced from the free carbon.

We found a high correlation between the abundances of c-C₃H₂ with HCO and H¹³CO⁺, with ρ ~ 0.8–0.9, and a good correlation between their line width, indicating how these two couples of molecules trace the same gas. The correlations with CH₃OH and H₂CO are less constrained, for both abundances and line widths. As said in Sect. 5.3, we would expect a low correlation coefficient between the species chemically linked to the CO and c-C₃H₂ because of their different formation pathways, and their weak chemical link. The correlation between these species may not result from a direct chemical link. Still, it could arise from their simultaneous formation in the star-forming region or from the similar gas they trace.

All Tables

All Figures

Fig. 1 Galactocentric gradients of fractional abundances with respect to H2, Xmol. The plot shows the trends for HCN (a), HCO+(b), c-C3H2 (c), H13CO+ (d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of the galactocentric radius (RGC). The light blue data represent the abundances calculated using a constant gas-to-dust ratio (γ = 100), while the dark blue data illustrate those estimated using a non-constant gas-to-dust ratio (γ = γ(RGC) from Giannetti et al. 2017). For both datasets, the linear regression results are shown as the light blue and dark blue lines, respectively. The 1σ error bars over the slope of the gradients are displayed for each fit. The upper-limit values are represented with triangles. The gradients of the elemental abundances of carbon (C and 13C), nitrogen (N), oxygen (O), and sulfur (S), as reported by Méndez-Delgado et al. (2022), are plotted as dashed lines. The product of the parent elements of the species is represented, in each subplot, by a dash–dotted black line. All elemental trends are plotted with a vertical offset to align with the starting point of the molecular gradients derived using γ(RGC), to facilitate comparison of their slopes with the linear fit of the molecular abundances. In the upper-right side of each subplot, the Pearson correlation coefficient, ρ, is shown (estimated only for the abundances estimated using the non-constant gas-to-dust ratio).

In the text

Fig. 2 Galactocentric gradients of line widths, FWHM. The plot shows the trends for HCN (a), HCO+ (b), c-C3H2 (c), H13CO+(d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of the galactocentric radius (RGC). The linear fit results are shown as dark blue lines. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

In the text

Fig. 3 Galactocentric gradients of molecular ratios. The plot shows the trend of the ratios between CO molecule products (i.e., HCO+, H13CO+, HCO, CH3OH, and H2CO) and c-C3H2, as a function of galactocentric radius (RGC). The 1σ error bars over the slope of the gradients are displayed for each molecular fit. The upper limit values are represented with triangles. The linear fit results are shown as the dark blue lines. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

In the text

Fig. B.1 Preliminary reduced NIKA2 images (color-scale in mJy/beam) at 2.0 mm (150 GHz) and 1.2 mm (260 GHz) of 12 out of the 31 CHEMOUT sources observed (Fontani et al. in prep.). The contours start from the 5σ rms level in each map. The NIKA2 angular resolution is shown in the bottom-left corner of each frame. The black circles on the maps represent the 28 ″ beam, centered on the position targeted for the molecular line observations analyzed in this paper.

In the text

	Fig. B.2 Same as Fig. B.1 for 12 additional targets.
In the text

	Fig. B.3 Same as Fig. B.1 for 7 additional targets.
In the text

Fig. B.4 Gradients of column densities of H2, NH2, as a function of Galactocentric radius (RGC). The light blue dataset represents the column densities estimated using a constant gas-to-dust ratio (γ = 100), while the dark blue dataset shows the column densities calculated with a gas-to-dust ratio function of RGC from Giannetti et al. (2017). Both datasets have been fitted with a linear regression. The 1σ error bars over the slope of the gradients are displayed for each fit. In the upper-right side of the plot, the Pearson correlation coefficient, ρ, is shown (only for the dark blue dataset).

In the text

Fig. D.1 Galactocentric gradients of column densities, Ntot. The plots show the trends for HCN (a), HCO+ (b), c-C3H2 (c), H13CO+ (d), HCO (e), SO (f), H2CO (g), and CH3OH (h) as a function of Galactocentric radius (RGC). The upper limit values are represented with triangles. The dark blue lines represent the linear regressions computed over the datasets. In the upper-right side of each subplot, the Pearson correlation coefficient, ρ, is shown.

In the text

Fig. E.1 Galactocentric gradients of fractional abundances w.r.t. H2, Xmol, of H13CN. The light blue data represent the abundances calculated using a constant gas-to-dust ratio (γ = 100), while the dark blue data illustrate those estimated using a non-constant gas-to-dust ratio (γ = γ(RGC) from Giannetti et al. 2017). For both datasets, the linear regression results are shown as the light blue and dark blue lines, respectively. The 1σ error bars over the slope of the gradients are displayed for each molecular fit. The upper limit values are represented with triangles. The gradients of the elemental abundances of carbon (13C) and nitrogen (N), as reported by Méndez-Delgado et al. (2022), are plotted as dashed lines. The product of the parent elements of the species is represented by a black dash–dot line. All the elemental trends are plotted shifted to the start of the molecular gradients estimated with γ(RGC) for reference to allow comparison of their slopes with the linear fit of the molecular abundances. In the upper-right side of each plot, the Pearson correlation coefficient, ρ, is shown (estimated only for the abundances estimated using the non-constant gas-to-dust ratio).

In the text