© The Authors 2026

1 Introduction

As the nearest benchmark for studying nuclear molecular zones in spiral galaxies, the Central Molecular Zone (CMZ) of the Milky Way contains a substantial reservoir of dense molecular gas, accumulated through bar-driven inflows. While strong shear and turbulence are thought to inhibit star formation in this region (e.g. Yusef-Zadeh et al. 2009; An et al. 2011), a growing body of evidence indicates that the CMZ continues to host significant star-forming activity. This is supported by the presence of compact H II regions (Gaume et al. 1995; Meng et al. 2019, 2022), abundant H₂O and methanol masers (Caswell et al. 2010; Longmore et al. 2017), and a population of young stellar objects (YSOs) (An et al. 2009, 2011, 2017; Jang et al. 2022), along with young massive star clusters containing O- and B-type stars (Mauerhan et al. 2010; Dong et al. 2012). More recently, ALMA observations have revealed a large number of deeply embedded protostars and associated bipolar SiO (5-4) outflows (Ginsburg et al. 2018; Walker et al. 2021), while JWST imaging has uncovered shocked molecular hydrogen structures and protostars, further confirming ongoing star formation in the region (Crowe et al. 2025).

Among the various tracers of star formation, YSOs and compact cloud cores (e.g. Hatchell & Fuller 2008; Heyer et al. 2016) not only mark ongoing star-forming activity but also provide insights into the chemical composition of the natal molecular clouds that will give rise to the next generation of stars. In particular, infrared (IR) spectroscopy reveals a rich set of absorption features arising from the vibrational and rotational transitions of ice species such as H₂O, CO, CO₂, and CH₃OH (e.g. Chiar et al. 1995; Gerakines et al. 1999; Whittet et al. 1998;

Pontoppidan et al. 2008; Boogert et al. 2011; Öberg et al. 2011; Boogert et al. 2015, 2022). These molecules form on icy dust grains under cold, dense conditions, beginning with the gravitational collapse of interstellar clouds, and are subsequently processed as the central protostar evolves. This processing includes thermal annealing, crystallisation, molecular segregation, and eventual sublimation, each reflecting changes in the thermal and chemical environment (Boogert et al. 2015). Therefore, observations of IR ice absorption bands provide a powerful diagnostic tool for probing the structural and evolutionary properties of YSOs and cloud cores.

In this context, both direct and indirect detections of methanol (CH₃OH) ice in point-like sources with extremely red colours in the CMZ (An et al. 2009, 2011, 2017; Jang et al. 2022) suggest that these objects are in an early stage of YSO evolution. The 15 μm CO₂ absorption band is commonly observed towards quiescent dark clouds, where it typically exhibits a relatively narrow profile, comprising contributions from pure CO₂ ice and mixtures with H₂O (peaked at 15.3 μm) and/or CO (15.1 μm). In contrast, significantly broader 15 μm features were found in the Spitzer/IRS mid-IR spectra of CMZ sources, due to the presence of a 15.4 μm shoulder component arising from CO₂ ice mixed with species like CH₃OH. It is attributed to Lewis acid-base interactions that shift the CO₂ asymmetric stretching mode to longer wavelengths (Ehrenfreund et al. 1999; Dartois et al. 1999). Since methanol ice forms under low temperatures (T < 15 K; Cuppen et al. 2009) and relatively high extinctions (A_V ≳ 9; Boogert et al. 2011; Chiar et al. 2011; Whittet et al. 2011), it serves as a sensitive tracer of cold, dense environments characteristic of the earliest stages of star formation. Accordingly, An et al. (2011) interpreted the observed features as indirect evidence of the presence of CH₃OH ice and identified 16 objects as YSOs and 19 as candidate YSOs. This interpretation was later confirmed by follow-up near-IR spectroscopy (An et al. 2017; Jang et al. 2022), which directly detected the 3.535 μm absorption band of solid CH₃OH in many of the same sources. These results demonstrated that methanol ice can comprise up to 17% ± 3% relative to H₂O, validating the earlier interpretation and confirming the presence of abundant CH₃OH ice in these red CMZ objects.

However, these follow-up studies consistently reported a lower abundance of CH₃OH with respect to H₂O in CMZ objects compared to values found in star-forming clouds within the Galactic disc. This discrepancy may reflect intrinsic differences in the chemical composition of CMZ clouds and/or variations in the chemical reaction networks operating under the extreme physical conditions of the CMZ interstellar medium (ISM). To better understand the origin of this systematic difference in the chemical properties of CMZ YSOs, it is useful to investigate the spatial distribution and structure of ice species within YSOs. For example, Pontoppidan et al. (2008) used observations of background stars in the Ophiuchus cloud to spatially resolve the distribution of ices within the extended envelope of a YSO. However, at a distance of ~8 kpc from the Sun (Reid et al. 2014), direct observations of the detailed structure of individual YSO envelopes in the CMZ are challenging; for example, most CMZ objects appear as point-like sources in the mid-IR bands of Spitzer/IRAC images (FWHM~2"; Ramírez et al. 2008).

On the other hand, near- and mid-IR spectroscopy can provide a unique opportunity to probe the internal structure of YSO envelopes through absorption features. From our near-IR spectroscopic observations (An et al. 2017; Jang et al. 2022), we found that many of the YSO candidates previously identified from Spitzer mid-IR spectra (An et al. 2011) exhibit strong 2.3 μm CO band-head absorption, a feature characteristic of (super-)giant stars. The extreme red colours of these objects further suggest that their near-IR fluxes are significantly attenuated by dust. These findings indicate that the observed near- and mid-IR spectra are, in fact, composites of two sources along the line of sight: the mid-IR emission is dominated by the YSO in the early evolutionary stage, while the near-IR absorption originates primarily from a background giant star situated behind the YSO’s extended envelope. Such a configuration, in which a background giant star backlights a foreground YSO, is plausible in the CMZ given its high source density. This backlighting scenario naturally could explain both the spectral energy distributions (SEDs) of Spitzer’s red point-like sources in the CMZ and their observed absorption features.

Building on the unique backlighting configuration of CMZ YSOs, we aim to expand the sample of such systems and use them as probes of how ice is distributed within YSO envelopes. To address this question, we carried out follow-up near-IR spectroscopic observations of 23 red point-like sources that exhibit a strong CO₂ shoulder feature, and measured the column densities (N) of H₂O and CH₃OH ice. By combining these measurements with our previous mid-IR spectral data, we present a comprehensive analysis of these CMZ objects, including extinction estimates from both foreground dust and the envelopes of YSOs. Assuming that our sample represents YSOs at a similar evolutionary stage, we investigated the internal distribution of methanol ice within their envelopes. This work presents the first effort to model an ensemble-averaged radial profile of CH₃OH ice abundance as a function of projected radial distance within YSO envelopes in the CMZ.

This paper is organised as follows. In Sect. 2, the Gem-ini/GNIRS observations of red point-like sources in the CMZ are described, along with the data reduction procedures. In Sect. 3, the GNIRS spectra are supplemented with existing near-IR data and combined with mid-IR spectra to produce the composite SEDs of the CMZ sources. In Sect. 4, the modelling of the overall SED from near- to mid-IR spectra is presented, from which the foreground dust extinction is derived. Column densities of H₂O and CH₃OH measured from the GNIRS spectra are also reported. The ice compositions of the CMZ objects are examined and compared with those in the Galactic disc. In Sect. 5, the derived ice abundances and extinction estimates are combined to construct an ensemble-averaged radial profile of methanol ice abundance. A summary of the results and a discussion of the observed trends are given in Sect. 6.

2 Gemini/GNIRS observations

2.1 Sample selection

We used the Gemini Near-IR Spectrograph (GNIRS; Elias et al. 2006a,b) on Gemini North to obtain K- and L-band spectra of 15 CMZ objects, as listed in Table 1. These targets were primarily selected from An et al. (2011), which reported mid-IR spectroscopic observations with Spitzer/IRS. Their near- and mid-IR colours are extremely red, with 1.6 ≤ K - [3.6] ≤ 4.8 and 2.4 ≤ [3.6] - [8.0] ≤ 5.3, which was the basis for their inclusion in the original Spitzer/IRS follow-up observations. These colours distinguish them clearly from other stars towards the CMZ, even though typical CMZ stars already appear significantly reddened due to large foreground dust extinction (A_V ~ 30 mag).

The GNIRS follow-up targets were primarily selected to maximise the likelihood of detecting a high CH₃OH ice abundance. The last column of Table 1 shows the classification of each object from An et al. (2011), which was determined based on the strength of the 15.4 μm shoulder absorption feature of CO₂ ice. Since the 15.4 μm feature is attributed to CO₂ ice mixed with CH₃OH (An et al. 2017; Jang et al. 2022), objects classified as YSOs are those for which An et al. (2011) assigned greater weight to the presence of a strong and statistically significant detection of this feature. In contrast, sources classified as possible YSOs show either weaker 15.4 μm absorption or less significant detections. In total, 12 of the brightest sources (L < 12.2 mag) from An et al. (2011) were observed, including 5 YSOs and 7 possible YSOs.

In addition, three objects (SSTGC 425399, 619964, and 696367) previously classified as known stars in An et al. (2011), based on a literature search, were also included in this study. Among them, SSTGC 619964 and SSTGC 696367 were identified as OH/IR stars in Lindqvist et al. (1992). However, our understanding of the red point-like objects in An et al. (2011) has since been significantly revised, following the realisation that many of these sources exhibit spectroscopic signatures of both young and evolved stars due to overlapping sightlines (An et al. 2017; Jang et al. 2022), as is further demonstrated in this work. Although these sources show relatively weak 15.4 μm shoulder features, their IR colours remain as red as those of the YSO candidates in An et al. (2011), providing an opportunity to explore the range of physical conditions and geometries of YSOs in the CMZ.

Most of our CMZ objects are likely Class 0 YSOs, i.e. in an early stage of evolution (Lada 1987; Andre et al. 1993), as is indicated by the large column of ices, including CO₂, along the line of sight. Moreover, As is discussed in Appendix D, neither H₂O nor CH₃ OH maser emission was detected from these targets, consistent with them being too young to excite Class I masers, which are collisionally pumped in shocks and outflows, or Class II masers, which are radiatively pumped by IR emission from warm dust close to the protostar.

2.2 Observations and data reduction

Queue observations were carried out at Gemini North in 2017 and 2022¹ to obtain L-band spectra for all selected targets. The GNIRS was used in long-slit mode with the 10.44 lines mm⁻¹ grating and the long red camera (0.05″pixel⁻¹). The slit width was matched to the seeing during natural seeing guided observations, providing a spectral resolving power of R = λ/∆λ ~ 400500 at ~3.5 μm, where the CH₃OH absorption band is located. In addition, three of the L-band targets were also observed in the K -band during the 2017 run (see Table 1). These observations employed the 31.7 lines mm⁻¹ grating with the short blue camera. For both L- and K -band observations, a standard ABBA dithering pattern was used to optimise sky subtraction and achieve a sufficient signal-to-noise ratio (S/N). The total on-source exposure time per exposure sequence was 16-32 min in K and 8-40 min in L, depending on the source brightness.

Standard data reduction procedures were performed using the IRAF/Gemini package (Tody 1986, 1993; Fitzpatrick et al. 2024)². Flat-fielding and wavelength calibration frames were acquired with GCAL, the facility calibration unit, as part of the observing sequences. Corrected individual frames obtained in ABBA dithering mode were combined. Wavelength calibration for the K-band spectra was performed using arc lamp lines, while the Maunakea atmospheric transmission model (for an airmass of 1.0 and a precipitable water vapour column of 1.0 mm) was employed for the L band. Telluric absorption features were corrected using observations of bright standard stars with spectral types B3-A5, taken with the same instrumental configuration as the science targets and observed immediately before or after the target to ensure airmass matching. Telluric corrections were applied using the Spextool package (version 4.0; Cushing et al. 2004)³.

Repeat L-band observations were obtained for four targets (SSTGC 563780, SSTGC 619964, SSTGC 696367, and SSTGC 799887) at different epochs to assess potential systematics in the observing conditions and data reduction. With the exception of SSTGC 563780, the L-band fluxes measured at two epochs differ systematically by up to ~50%, suggesting uncertainties arising from the extraction and calibration of the spectra. In particular, SSTGC 619964 and SSTGC 799887 lie in very crowded regions, where contamination from nearby sources likely affected the source extraction. Among the targets with only single L-band observations, SSTGC 524665 and SSTGC 770393 also have nearby sources that may have contributed to their measured fluxes.

The S/N of each spectrum was estimated by fitting a low-order polynomial around ~3.7 μm, where no strong absorption features are present. The typical S/N is ~40, which is sufficiently high to allow for the confident measurement of CH₃OH column densities above ~1 × 10¹⁷ cm⁻². From telluric absorption lines in the standard star spectra, the resolving power was measured to be in the range R ~ 400-500. For the repeat observations, we carried out independent modelling of spectra from each exposure sequence before combining the resulting parameter estimates, including column densities and foreground extinctions (see below).

Fig. 1 GNIRS L-band spectra of 12 red sources in the CMZ. The observed spectra (red) are vertically offset for clarity and approximately ordered by increasing strength of the broad H2O ice absorption band centred at 3 μm. To illustrate the varying depth of this feature, the best-fit IRTF/SpeX M5.5III template spectrum (Rayner et al. 2009) is shown in green for each source. These templates were reddened using the foreground extinction estimate derived from our SED fitting procedure (see Sect. 4.1). Together with the 3 μm ice absorption profile, they reproduce the overall continuum shape and slope of the observed spectra.

3 IR spectra of red CMZ objects

3.1 Composite spectra

Figure 1 presents the GNIRS spectra of the 12 red CMZ objects. The spectra are vertically offset with arbitrary zero points and are ordered by increasing strength of the broad H₂O ice absorption band centred at 3 μm, whose wing extends across much of the L band. Some objects, such as SSTGC 728480 and SSTGC 563780, exhibit nearly flat spectra with very weak or absent H₂O ice absorption. In contrast, the majority of our sample shows strong band absorption, suggesting significant extinction due to dense molecular clouds along the line of sight. Additionally, Figure 2 displays the spectra of the three objects for which both K- and L-band observations were obtained. It reveals CO bandhead absorption features at 2.3-2.4 μm, commonly associated with stellar photospheres. These features are indicative of (super-)giant stars, and their presence, together with the extremely red SED and strong molecular absorption bands, suggests the superposition of such stars behind a YSO envelope or a compact cloud core (see Jang et al. 2022).

In the following analysis, we combined the GNIRS spectra with mid-IR spectra obtained from Spitzer/IRS (An et al. 2011) to enable a comprehensive characterisation of the red pointlike sources in the CMZ. Figure 3 shows an example of such a composite spectrum (SSTGC 799887). The GNIRS spectra cover the wavelength range 2-4 μm, while the low-resolution Spitzer spectra, obtained with the short-low (SL) and long-low (LL) modules, span 5-35 μm. As is detailed below, this multiwavelength analysis allowed for the simultaneous estimation of key absorption features, including those of H₂O and CH₃OH, as well as the dust extinction along the line of sight.

However, the GNIRS spectra occasionally exhibited significant flux discrepancies with the IRS spectra near the 4-5 μm boundary. These offsets likely reflect a combination of observational limitations in the CMZ and variations in spatial resolution and aperture size between instruments. Spectral extraction in this environment is complicated by high source density, where crowding affects background subtraction and source isolation. The large apertures of the IRS modules (3.7" for SL and 10.7" for LL) are particularly prone to contamination from nearby sources, while GNIRS, with its higher spatial resolution, samples a more localised region. Morphological differences further contribute to these mismatches, as many CMZ sources appear compact in the near-IR but display extended emission at longer wavelengths. In contrast, our ground-based GNIRS observations benefit from higher spatial resolution but are subject to limitations in atmospheric correction, which inevitably lead to systematic uncertainties in the absolute flux calibration. The extended nature of the emission in the L band further increases slit-loss uncertainties, leading to flux mismatches between the K-and L-band spectra. Overall, both the GNIRS and IRS spectra are affected by distinct sources of systematic uncertainty, which can introduce offsets in the absolute flux calibration and continuum level across the full spectral range.

For these reasons, we anchored all spectra to broadband photometric measurements, both to mitigate source confusion in the large-aperture IRS data and to establish a robust zero-point flux calibration for the GNIRS spectra. To place the GNIRS and IRS spectra on a consistent flux scale, we matched them to broadband photometry from the Galactic plane survey conducted as part of the UKIRT IR Deep Sky Survey (UKIDSS; Lucas et al. 2008), and from the Spitzer/IRAC survey of the Galactic centre (GAL-CEN; Ramírez et al. 2008). We used K-band photometry to scale the K-band spectra, computing synthetic magnitudes using the corresponding filter response function available from the SVO filter profile service⁴. Similarly, the [3.6] photometry was used to set the flux scale for the L-band spectra. For the IRS spectra, we adopted the [8.0] photometry for scaling, as the other IRAC bands ([4.5] and [5.8]) fall in wavelength regions that are either only partially covered or entirely missing from the spectra.

Figure 3 shows the spectra calibrated in this way, with the broadband photometry over-plotted as filled circles at their respective effective wavelengths. We note that the [8.0] photometric point lies below the observed spectra, owing to the strong silicate absorption feature that affects the long-wavelength side of its passband. Because the K-, L-, and IRS spectral segments do not overlap, a direct assessment of their relative flux scaling is not possible. Nonetheless, the consistency of the flux calibration can be evaluated through SED modelling, introduced in Sect. 4.1. The solid blue lines in Fig. 3 show the corresponding best-fitting models, and any residual offsets with respect to the data are discussed therein.

As the flux uncertainties in the data products are often underestimated, we scaled them by a multiplicative factor derived from local pseudo-continuum regions free of strong absorption features. Specifically, we fitted a second-order polynomial to 2.10-2.25 μm for the K-band spectra and 3.50-3.70 μm for the L band. The scaling factor was chosen to match the median flux uncertainty to the standard deviation of the logarithmic flux in each range. This yielded more realistic random uncertainties of 1-4% for most L-band spectra, although the spectrum of SSTGC 524665 shows uncertainties as high as 17%. In addition, we applied a second-order Savitzky & Golay (1964) filter to the spectra to reduce noise while preserving the spectral resolution.

Fig. 2 Same as Figure 1, but for sources that have both K- and L-band spectra from GNIRS observations.

Fig. 3 Composite IR spectrum of a red CMZ object (SSTGC 799887) and modelling of the SED. The observed spectra are shown as a solid red line, constructed by combining near-IR data from Gemini/GNIRS with mid-IR data from Spitzer/IRS. Circles represent photometric measurements from UKIDSS and Spitzer/IRAC, which were used to set the flux scale of the spectra (see text). The solid blue line shows the best-fitting model to the observed spectra, which included the following components: (1) a template of an M5.5-type giant from the IRTF/SpeX spectral library (Rayner et al. 2009) (green line); (2) a two-component blackbody continuum (dashed blue and orange lines); (3) the 3 μm band of solid H2O ice, using the template from GC IRS 7 (Jang et al. 2022); and (4) the 10 μm and 18 μm silicate bands from the GC IRS 3 spectrum (Kemper et al. 2004). Other weaker absorption features were not included in the modelling. Spectra and their best-fitting models for the additional objects are shown in Appendix A.

3.2 NASA IRTF/SpeX

Jang et al. (2022) obtained K- and L-band spectra of red CMZ objects with a medium-resolution spectrograph SpeX (Rayner et al. 2003) on the 3.2m NASA IR Telescope Facility (IRTF). Their data are comparable to the GNIRS spectra presented here in both spectral resolution and wavelength coverage. Like our GNIRS sample, the IRTF data combined with IRS observations span the 2-35 μm range. Eleven of these objects are bright enough to be included in An et al. (2011), and they also exhibit extreme red colours ([3.6]-[8.0] ≳ 3.0). Two of these (SSTGC 653270 and SSTGC 696367) were also observed in this work; for these sources, we supplemented our GNIRS spectra with the K-band IRTF/SpeX data. The spectra of the remaining IRTF/SpeX objects, excluding one with very low S/N, are listed in Table 2, and were combined with the IRS spectra in the same manner as the GNIRS+IRS composite spectra.

The IRTF/SpeX observations were performed simultaneously in the K and L bands using a single instrumental set-up, thereby providing a consistency check on the reliability of the flux zero-point adjustments applied to match the photometric measurements. For the eight objects listed in Table 2, we found that the L-band spectra required a scaling factor roughly twice that of the K band to match the photometry. These offsets likely reflect the extended nature of the emission arising from YSO envelopes, which becomes more prominent at longer wavelengths. This interpretation is supported by the fact that the objects appear essentially point-like in the K band, but more extended in the L band (An et al. 2017). It is also possible that the lower spatial resolution of the IRAC [3.6] images used for L-band scaling may have included flux contributions from nearby sources, in contrast to the higher-resolution UKIDSS K-band images used for the K-band scaling; however, this effect is likely less significant than that of the extended emission.

In total, we have 23 objects with 3-35 μm composite spectra assembled from Gemini/GNIRS, IRTF/SpeX, and Spitzer/IRS observations. About half of these also have K-band spectra. For the remainder, we used UKIDSS K-band photometry to constrain their physical parameters, as is detailed below. The full set of composite spectra, corrected for zero-point flux offsets and including scaled flux uncertainties, is publicly available via Zenodo⁵.

4 Spectral analysis and derivation of physical parameters

Most of the observed CMZ objects in this study show deep silicate and H₂O ice absorption features, indicating the significant amount of molecular clouds along the line of sight to each object. While foreground dust extinction and column densities of ice species cannot be reliably constrained using broadband photometry alone, the SED traced by IR spectra provides a valuable diagnostic for quantifying such parameters. Below, we performed a global modelling of the SED to extract foreground extinctions, along with the column density of solid H₂O with a wide absorption band (Sect. 4.1). We combined this with additional analysis of more localised absorption features from CH₃OH ice (Sect. 4.2) to reveal complex structures of dust and ices along each line of sight towards the red CMZ objects.

4.1 Column densities of solid H₂ O and foreground extinction

As is shown in Fig. 3, the spectra of the red CMZ objects are characterised by several absorption bands. The near-IR absorption feature at 3 μm is primarily caused by the O-H stretching vibration in solid H₂O ice. The broad absorption features centred at 10 μm and 18 μm are characteristic of silicate dust and originate from the Si-O stretching and O-Si-O bending modes, respectively, in amorphous silicate grains. Superimposed on a continuum that would otherwise be nearly flat in λF_λ space, or slowly rising towards longer wavelengths, these absorption features are sufficiently broad to modify both the continuum level and the overall shape of the SED. In our previous SED modelling (An et al. 2017; Jang et al. 2022), we found that the overall SED can be well explained by the superposition of two sources of radiation: blackbody emission representing warm dust emission from the extended envelope of a YSO, and the spectrum of a red (super-)giant, both attenuated by significant columns of foreground dust. Given the above considerations, we modelled the composite IR spectra using the following components:

• Photospheric component: in the near-IR, many red CMZ objects exhibit CO bandhead absorption at 2.3 μm, indicating the presence of a red (super-)giant along the line of sight to the mid-IR source. Jang et al. (2022) performed detailed spectral fits using the NASA IRTF/SpeX library (Rayner et al. 2009), and showed that varying the template spectral type between late-K and late-M giants only mildly affects the parameter estimates, because, beyond 2 μm, the stellar spectra of late-type giants show little variation in slope. Based on this, we adopted the spectrum of an M5.5III giant (HD 94705) as a representative photospheric component. This choice was also motivated by the fact that K-band spectra are available for only 5 of the 15 GNIRS targets.

• H₂O ice bands at 3 μm and 13 μm: we utilised the observed IRTF/SpeX spectrum of GC IRS 7 from Jang et al. (2022) to model the 3 μm band, which arises from the O-H stretching mode. In addition, H₂O ice exhibits an even broader feature centred at 13 μm, associated with the librational mode. Since the GC IRS 7 spectrum is limited to the near-IR region, we complemented it with a laboratory spectrum of pure H₂O ice at 15 K from the Leiden Ice Database for Astrochemistry (LIDA; Rocha et al. 2022) for wavelengths beyond 3.8 μm. Although the laboratory data also include the 6 μm feature, we did not incorporate it, as it only partially accounts for the observed absorption in that region.

• Silicate bands at 9.7 μm and 18 μm: following An et al. (2009, 2011), we adopted the mid-IR silicate absorption profile from Kemper et al. (2004). Besides the continuum extinction, the strengths of the silicate features provide a strong constraint on the dust column density along the line of sight. The relative strengths of the two silicate modes were fixed during scaling to match the observed data. The bestfitting peak optical depth at 9.7 μm was converted to A_K by multiplying it by 0.99 (Roche & Aitken 1985), which is a factor of two lower than that towards the local ISM (Roche & Aitken 1984; Madden et al. 2022).

• Warm dust clouds: as demonstrated by An et al. (2017), the mid-IR continuum of red CMZ objects can be effectively reproduced by emission from two blackbody components, likely originating from the extended envelope of a YSO. According to our modelling (see below), the warmer component, with a temperature of ~300-400 K, dominates the spectral range between ~5 μm and 20 μm, where the photospheric contribution becomes negligible. In contrast, the cooler component, with a temperature of ~100 K, primarily accounts for the rising slope of the spectrum beyond ~20 μm.

We incorporated the components described above into our global SED modelling, assuming the continuum extinction curve from Boogert et al. (2011). In addition to these components, narrower but still significant absorption features are seen in our composite spectra, including the 6 μm and 6.8 μm complexes, and the 15 μm CO₂ ice band. Furthermore, the spectrum of SSTGC 799887 in Fig. 3 shows a narrower absorption feature at ~18.7 μm, which is in fact due to over-subtraction of background emission from [S III]. A similar effect is seen at ~10 μm, where over-subtraction at the bottom of the silicate band creates a spurious narrow absorption feature. On the other hand, the emission lines at 33.5 μm and 34.8 μm, from [S III] and [Si II], respectively, could trace the localised emission from the YSO itself. For our modelling, we selected spectral regions free of these features, using only clean continuum segments across the K-, L-, and mid-IR bands (see Appendix A).

The CH₃OH ice bands at 3.535 μm and 9.74 μm were not included in our global SED analysis and were instead masked, because both features are highly sensitive to the shape of the local continuum. The 3.535 μm band lies within the broader 3 μm H₂O and hydrocarbon absorption complex, while the 9.74 μm band is embedded in the deep 10 μm silicate absorption. Because a global SED fit provides only a coarse approximation to the continuum in these wavelength ranges, the continuum level in both regions is particularly uncertain. We therefore analysed the CH₃ OH bands separately, fitting a third-order polynomial to more accurately trace the local baseline (see Sects. 4.2 and C).

In modelling the water-ice absorption, we tied the strength of the 13 μm libration band to the peak optical depth of the 3 μm stretching mode. This choice was motivated by the strong and broad silicate absorption at 10 μm in our CMZ targets, which overlaps with the libration feature and prevents a reliable decomposition. In addition, while some studies isolated the 13 μm band using simplified polynomial continua (e.g. Boogert et al. 2011), we adopted a more physically motivated continuum model based on blackbody dust emission. As a result, the 13 μm feature could not be independently constrained in our fits and was instead anchored to the better-defined 3 μm absorption. We note, however, that the two bands may have probed slightly different lines of sight (see below).

The solid blue lines in Fig. 3 show the best-fitting model for SSTGC 799887, derived as is described above. The model included eight free parameters: a flux normalisation for the near-IR stellar template, temperatures and normalisations for the two blackbody components, foreground extinction affecting the near-IR source, foreground extinction affecting the mid-IR source (including both continuum and silicate absorption), and the column density of solid H₂O ice. To derive global parameters that reproduce the overall SED shape, we used curve_fit, the nonlinear least-squares optimisation routine from the SciPy library, applying weights based on the scaled uncertainties in logarithmic flux space. Results for the remaining objects are presented in Appendix A.

Overall, the models reproduce the observed SEDs of the red CMZ sources well, capturing the spectral curvature shaped by the major broad absorption bands. The 2.3 μm CO bandheads are also well matched. However, some local mismatches remain. For instance, the model slightly overestimates the depth of the 3 μm band (middle panel) and underestimates the flux on the shortwavelength side of the 9.7 μm silicate feature. The continuum at wavelengths shorter than 2.2 μm also appears systematically brighter in the model. Most notably, the model underpredicts the IRS flux between 5.5 and 7 μm. This mismatch may partly arise from source contamination in the IRAC photometry used to scale the IRS spectra, as the coarse spatial resolution of IRAC (~2") increases the risk of source contamination in crowded regions. Taken together, these deviations likely reflect a combination of factors, including uncertainties in spectral extraction in the dense CMZ environment, residual uncertainties in flux calibration and background subtraction, and possible limitations in the adopted templates for the 3 μm H₂O and 9.7 μm silicate bands in capturing the true continuum and band shapes. Nevertheless, as the model captures the global shape of the SED and the key spectral features, the derived parameters provide a reasonable approximation to the intrinsic spectral properties of the sources.

Results from the global SED fitting are summarised in Table 3 for both the Gemini/GNIRS and IRTF/SpeX samples. Only the three key parameters are reported. In the second and third columns, A^*_K denotes the foreground extinction in the K band associated with the near-IR photospheric component, while $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ represents the K-band extinction inferred from the mid-IR spectrum, including the contribution from the silicate absorption bands. The fourth column lists the peak optical depth of the 3 μm H₂O ice feature (τ_3.0). The corresponding column density is shown in the final column, computed by multiplying τ_3.0 by a fixed FWHM of 330 cm⁻¹ for amorphous H₂O ice and using a band strength of A = 2 × 10⁻¹⁶ cm molecule⁻¹ (Hagen et al. 1981).

The uncertainties in Table 3 represent the quadrature sum of fitting errors, derived from the covariance matrix, and various systematic uncertainties. The latter include the effect of varying the stellar template between M0IIIb (HD 213893) and M9III (IRAS 15060+0947), with the difference from the default case taken as an effective 2σ uncertainty. To evaluate the impact of the extinction law, we repeated the SED fitting with alternative curves from Indebetouw et al. (2005) and Fritz et al. (2011). The latter is notably steeper in the near-IR, leading to shifts in the derived parameters. The deviations from the default curve (Boogert et al. 2011) were treated as an effective 2σ uncertainty and included in the total error budget for both A_K and τ_3.0. For N(solid H₂O), we adopted a 10% uncertainty in the band strength and included it in the total error budget. For the targets with repeated observations, the differences in A^*_K or τ_3.0 between individual spectra exceed the statistical uncertainties but remain well below the combined systematic uncertainties. In such cases, half of the observed differences were added in quadrature to the total error budget. Upper limits are quoted at the 3σ level, calculated from the total error budget.

Several objects in this study overlap with Jang et al. (2022). Specifically, in addition to the eight sources listed as the IRTF/SpeX sample in Table 3, three objects (SSTGC 653270, SSTGC 696367, and SSTGC 719445) were also observed in the L band with Gemini/GNIRS. The unweighted mean and standard deviation of the differences are 2.5 ± 0.7 mag in A^*_K and -0.08 ± 0.44 in τ_3.0, indicating systematically higher extinction and lower optical depth in this study. These offsets primarily reflect our revised modelling approach, which incorporates the extended wavelength coverage up to 35 μm. Nevertheless, the dispersion remains well below the level of the total systematic uncertainties.

As is shown in the top panel of Fig. 4, $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ ranges from 2 to 5 mag, corresponding to A_V ~ 18-45 mag. In contrast, A^*_K spans a substantially larger range of 4-15 mag. The systematic excess of A^*_K over $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ indicates that the near-IR stellar template requires stronger attenuation than what is implied by the mid-IR silicate bands. Together with the presence of CO absorption bandheads in many of our CMZ spectra, this enhanced extinction supports the interpretation that our composite spectra combine two sources along nearly the same line of sight (Jang et al. 2022).

Figure 5 illustrates the light-ray geometry probed by our composite spectra. The photospheric near-IR emission from a background giant star passes through the extended dusty envelope of a protostar and therefore encounters a deeper column of material than the mid-IR emission emerging from the protostar’s core. At the distance of the CMZ, the two light rays are unresolved and their emissions blend after undergoing the same foreground extinction by dust in the Galactic disc and/or within the CMZ. Because the near-IR part of the spectrum is dominated by photospheric emission from the background star, while the mid-IR is dominated by thermal emission from the cooler protostellar core, this configuration naturally produces the distinct extinction estimates A^*_K and $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ obtained from the near-and mid-IR spectral regions, respectively. The absence of H₂O and Class I CH₃OH maser emission towards nearly all sources (Appendix D) is consistent with the early, envelope-dominated evolutionary stage required for this backlighting geometry.

The middle and bottom panels in Figure 4 show the H₂O column density from the 3 μm band plotted against the foreground extinction estimates. Of the 23 targets in our sample, eight exhibit strong H₂O ice absorption with N(solid H₂O) > 2 × 10¹⁸ cm⁻², while six show moderate absorption (1-2 × 10¹⁸ cm⁻²). The remaining sources display only weak features or yield upper limits. No clear correlation is found between N(solid H₂O) and the YSO classifications of An et al. (2011, see Table 1): for example, two OH/IR stars (SSTGC 619964 and SSTGC 696367) show significant detections, whereas another star known in the literature has only an upper limit. As is discussed further below, in the context of near-/mid-IR source superposition, the near-IR classifications of these sources as (super-)giants do not necessarily preclude them from being YSOs.

For sources with significant detections of H₂O ice, the measured optical depths are τ_3.0 ≈ 0.5-4.0. A minor contribution could arise from intervening Galactic disc clouds, but this effect is expected to be negligible, with Moultaka et al. (2015) finding τ_3.0 ≈ 0.03 towards stars in the vicinity of Sgr A^*. Even if this full amount were present along our sightlines, it would account for only ~6% of the weakest detection and <1% of the strongest. We therefore consider the foreground H₂O ice absorption negligible and apply no correction.

For comparison, the solid lines in Fig. 4 show the empirical relation for Taurus derived by Whittet et al. (2001), who found that H₂O ice appears in dense clouds once the extinction exceeds A_K ≈ 0.35 mag (see also Whittet et al. 2011). Since the CMZ sightlines include foreground diffuse ISM clouds contributing A_K ~ 2.2 mag (A_V ~ 20 mag), we shifted the Taurus relation horizontally by this amount. The shifted curve thus reflects the extinction-ice column density relation expected for CMZ sources. However, the shifted Taurus relation provides a poor match to either of the observed distributions. In the middle panel of Fig. 4, most objects with detected H₂O ice (open circles) lie above the N(solid H₂O) values predicted from their $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$. In the bottom panel, A^*_K extends to much higher values, emphasising the exceptionally large extinctions along the near-IR path. The Whittet et al. relation lies between the two distributions defined by $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A^*_K, indicating that a single Taurus-like relation does not adequately describe our CMZ sample. This is not surprising, since CMZ sightlines are considerably more complex than those in nearby molecular clouds, as is illustrated by the following effects:

• The threshold extinction for the onset of H₂O ice formation is not consistently reproduced by either $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ or A^*_K. In the middle panel, sources with upper limits in N(solid H₂O) broadly define a CMZ threshold near $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ ≈ 2.55 at the minimum detected N(solid H₂O), yet one object (SSTGC 653270) already shows substantial ice at this value. In the bottom panel, the apparent threshold in A^*_K is nearly twice as large. This behaviour is naturally explained by the highly patchy foreground extinction towards the CMZ, which varies on ~ 1′ scales and shifts the effective extinction at which ice appears along a given line of sight.

• Even where the threshold extinction appears roughly consistent in the middle panel, the shifted Taurus relation still fails to reproduce the observed $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$-N(solid H₂O) distribution. Nearly half of the sources lie above the Taurus prediction in N(solid H₂O) at a given $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$, which would imply either that $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ underestimates the relevant dust column or that H₂O ices form more efficiently in the CMZ than in Taurus. The latter scenario seems unlikely in the strongly UV-irradiated, warmer CMZ environment, suggesting instead that $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ does not trace the same column of cold dust that dominates the ice absorption.

• Despite the measurement uncertainties, an overall increase in N(solid H₂O) with A^*_K is evident in the bottom panel. Some of the scatter may again be driven by patchy foreground extinction. However, because our N(solid H₂O) estimates are mainly constrained by the near-IR 3 μm ice band, while the librational band provides only a weak additional constraint, N (solid H₂O) is expected to correlate more tightly with A^*_K than with $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$, in agreement with the observed trend.

• This behaviour reflects the different optical paths that define $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A^*_K through the same YSO envelope (Fig. 5): the mid-IR path ($A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$, middle panel) samples a shorter route through warmer, more deeply embedded regions, whereas the near-IR path (A^*_K, bottom panel) traverses cooler, more extended parts of the envelope. If A^*_K and $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ therefore probe different dust columns, then thermal processing of ices and radial gradients in ice abundance within the envelope naturally lead to departures from a simple, single-column Taurus-like dust-N(solid H₂O) relation.

Fig. 4 Distribution of extinction and H2O ice column densities derived from the global fits to the composite spectra. The two extinction estimates, $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A*K, come from our composite spectra but rely on different wavelength regimes: the mid-IR spectral shape for $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and near-IR template matching for A*K. Objects with large solid CH3OH column densities (N > 3 × 1017 cm−2; see below) are marked with open blue circles. The dashed line in the top panel indicates unity. For comparison, the solid lines in the middle and bottom panels show the empirical relation derived for the Taurus molecular cloud by Whittet et al. (2001), shifted horizontally by AK = 2.2 mag to account for the average foreground extinction towards the CMZ.

Fig. 5 Schematic illustration of the light-ray geometry probed by the composite spectra. A*K reflects extinction constrained mainly by the near-IR, while $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ is determined primarily from the mid-IR. The quantity $A_K^{\mathrm{f}\mathrm{g}}$ represents the foreground extinction estimated from nearby source-free regions, which accounts for dust located along the line of sight in the Galactic disc and within the CMZ. The envelope-only extinction, $A_{K,0}^{\ast }$, is defined as $A_K^{\ast }-A_K^{\mathrm{f}\mathrm{g}}$, isolating the extinction produced by material within the protostellar envelope along the sightline to the background star.

4.2 Column densities of solid CH₃ OH

In addition to the SED fitting, we analysed the 3.535 μm absorption feature from solid CH₃OH (Hudgins et al. 1993; Schutte et al. 1996). The band is sufficiently narrow that the local continuum was fitted with a third-order polynomial over the 3.40-3.50 μm and 3.60-3.70 μm intervals. The feature was modelled in optical-depth space with the laboratory spectrum of pure CH₃OH ice at 10 K (Hudson et al. 2024), fitted using the curve_fit routine from the SciPy library. As it lies on the long-wavelength wing of a stronger, broader band centred at 3.42 μm, it was extracted using the same polynomial fit and wavelength ranges adopted for the optical-depth spectra. The model spectrum was then smoothed to match the spectral resolution of our observations.

Table 4 summarises the peak optical depth of the 3.535 μm absorption in the best-fit model (τ_3.535) for both the Gem-ini/GNIRS and IRTF/SpeX samples. The column density N(solid CH₃OH) was computed using a band strength of 5.19 × 10⁻¹⁸ cm molecule⁻¹ (Hudson et al. 2024). The value is close to 5.3 × 10⁻¹⁸ cmmolecule⁻¹ in Hudgins et al. (1993); however, a conservative 12% uncertainty in the band strength was adopted to account for the impact of the 3.42 μm feature. We further assessed systematic effects by re-fitting the spectra with a second-order polynomial continuum and by modelling the feature with a Gaussian profile to evaluate model dependence in τ_3.535. The uncertainties listed in Table 4 are the quadrature sum of these contributions and the statistical errors derived from the covariance matrix of the fit. Upper limits are quoted at the 3σ level, calculated from the total error budget.

Of the targets with repeated observations, only SSTGC 696367 exhibits a significant detection of the 3.535 μm band, and its individual τ_3.535 measurements are consistent within the statistical uncertainty. The three objects (SSTGC 653270, SSTGC 696367, and SSTGC 726327) with significant detections were also included in Jang et al. (2022); the differences in τ_3.535 are well within the measurement uncertainties (< 0.01). Since the comparison for SSTGC 653270 and SSTGC 696367 was made based on independent observations, it demonstrates the robustness of our τ_3.535 estimates. In case of SSTGC 726327, the IRTF/SpeX spectra have been re-analysed in this work, which also yields a result consistent with the previously reported value of 0.10 ± 0.01 in An et al. (2017).

Figure 6 shows the optical-depth spectra of objects with detections at >2.5σ significance, including four with >3σ confidence. Spectra of the remaining objects are presented in Appendix B. The red shading denotes the ±1σ range of the total uncertainty. Among these six objects, SSTGC 726327 and SSTGC 817031 have previously been reported to show CH₃OH (Jang et al. 2022), while the 3.535 μm detections for the remaining four objects are new from this work. As is noted above for H₂O ice, the detection of CH₃OH ice is not correlated with the mid-IR classifications in An et al. (2011): for example, the previously known stars (SSTGC 696367 and SSTGC 817031) exhibit moderate to significant detections.

Figure 7 shows N(solid CH₃OH) plotted against our foreground extinction measurements: $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ in the top panel and A^*_K in the middle panel. As for water ice in Fig. 4, N(solid CH₃ OH) measurements show a large upward scatter above the local relation with respect to $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ (top panel). The scatter of N(solid CH₃OH) is more pronounced against A^*_K (middle panel). For reference, the local relation from Whittet et al. (2011) is shown as a solid line, shifted by A_K = 2.2 mag to account for the extinction of the foreground diffuse clouds in front of the CMZ. Methanol ice tends to form in denser regions and deeper layers of clouds, where the extinction is generally higher than in the environments where H₂ O ice forms. Thus, the extinction threshold for detecting CH₃OH ice (A_K ≈ 1.87 ± 0.50) is higher than that for water ice (see Boogert et al. 2011). However, the threshold extinction of N(solid CH₃ OH) in the top panel is less apparent than that for water ice in Fig. 4.

In the bottom panel of Fig. 7, N(solid CH₃OH) is shown against N(solid H₂O). The absence of a clear correlation reflects the significant scatter observed in the top and middle panels. Nevertheless, in local clouds, both N(solid H₂O) and N (solid CH₃OH) typically increase with extinction (Boogert et al. 2011), suggesting higher N (solid CH₃OH) at higher N (solid H₂O). Therefore, the lack of a positive trend in our sample, particularly considering upper limits in these measurements, may indicate that the ice abundances have been modified by environmental effects, such as thermal processing within YSO envelopes. These possibilities and their implications for ice chemistry in the CMZ are explored in the following section.

4.3 Mixing ratio of CH₃ OH and CO₂

The primary feature distinguishing our sample objects from the other CMZ sources is the strong, broad 15 μm absorption from CO₂ ice. As is described in Sect. 1, CH₃ OH-rich CO₂ ice grains are responsible for the long-wavelength shoulder of this absorption band. The greater shoulder absorption compared to the polar (H₂ O-rich) and apolar (CO-rich) components of CO₂ ice is uniquely seen towards massive YSOs in the CMZ, suggesting that distinct chemical networks operate in CMZ star-forming clouds, producing unique chemical mixtures.

The upper left panel of Figure 8 shows N(solid CH₃ OH) plotted against the CO₂ shoulder column densities from An et al. (2011) for all 23 objects in our CMZ sample, with red circles indicating detections and arrows marking 3σ upper limits. The ratio of these quantities represents the CH₃OH:CO₂ mixing ratio in ice grains, spanning 0.7-5.9 for the six objects with useful 3.535 μm detections (>2.5σ). A linear fit to the CMZ detections (excluding upper limits) yields a slope of N(solidCH₃OH)/N(CO₂ shoulder) = 3.3 ± 0.7 (dashed line), consistent with 2.9 ± 0.6 reported by Jang et al. (2022). For comparison, the upper right panel shows measurements for YSOs in the Galactic disc (Boogert et al. 2008; Pontoppidan et al. 2008), which span a broader range (1.3-52.2) and have a mean ratio of 12.8 ± 2.4 (dashed line), significantly higher than the CMZ value. This contrast is also evident in the lower panels, where individual probability density functions (thin lines) are modelled as bivariate normal distributions using the measured column densities and their 1σ uncertainties; their sum (thick line) represents the overall distribution.

According to the LIDA laboratory experiment of CO₂ ice absorptions (Ehrenfreund et al. 1999), even a relatively low abundance of CH₃OH in CH₃OH:CO₂ mixtures can produce the shoulder absorption of solid CO₂ centred at 15.4 μm at temperatures below 130 K. The band strength extracted from the laboratory data is nearly the same as in more methanol-rich mixtures, such as those appropriate for the Galactic disc population. On the other hand, the shoulder component could also arise in mixtures with species other than CH₃ OH. In that case, the presence of additional species would enhance the strength of the 15.4 μm shoulder, while N(solid CH₃OH) would remain unchanged because the 3.535 μm band originates primarily from CH₃OH. Thus, if other species that act as bases in Lewis acid-base interactions, such as ethanol or related molecules, contribute to the CO₂ shoulder (e.g. Dartois et al. 1999), the inferred CH₃OH:CO₂ ratio would decrease, since part of the 15.4 μm absorption would not be due to CH₃OH. Such a scenario would make the CMZ objects distinct from those in the Galactic disc, but only if the CMZ uniquely hosts a large reservoir of these more complex molecules.

If the 15.4 μm shoulder arises entirely from CO₂ ice mixed with CH₃OH, the relatively low abundance ratios shown in Fig. 8 may point to a distinctive chemical environment in the CMZ. In particular, they could suggest that CH₃ OH formation is less efficient relative to CO₂ in CMZ clouds, even though both species are produced through CO reactions involving hydrogenation or oxidation (see Jang et al. 2022). However, the CH₃OH:CO₂ mixing ratio in the CMZ displays substantial scatter. Among the six CMZ sources with significant detections, SSTGC 817031 shows the highest ratio (~27.7), whereas SSTGC 719445 has the lowest (~2.1). Including upper limits for the other 17 objects broadens the range even further, implying large variations across the sample. Such a wide spread in abundance ratios would imply a high level of chemical diversity among clouds in the CMZ.

Instead, we argue in the next section that the observed scatter primarily reflects variations in the lines of sight, with each spectrum probing different regions of the YSO envelope. This geometry offers a rare opportunity to investigate the internal chemical structure of CMZ YSOs: CO₂ absorption is measured from mid-IR spectra, whereas N(solid CH₃OH) is derived from near-IR data. By contrast, the consistently higher abundance ratios found for YSOs in the Galactic disc may stem from tracing the same mid-IR sightlines. In that case, N(solid CH₃ OH) estimates were based on the 9.7 μm C-O stretching mode (Skinner et al. 1992) in the same mid-IR spectra (Pontoppidan et al. 2008). We revisited this feature in our own data in Appendix C, but the CMZ sources are considerably fainter, leading to much lower S/N at the bottom of the deep 9.7 μm silicate bands.

Fig. 6 Optical depth spectra of CMZ sources exhibiting 3.535 μm solid CH3OH absorption at >2.5σ significance. The red line shows the best-fit model for pure solid CH3OH ice, and the shaded area denotes the ±1σ confidence interval. Spectra for the remaining sources are provided in Appendix B.

5 Ensemble average of CMZ YSOs

5.1 Reduced foreground extinction

Our CMZ objects appear as point-like sources in Spitzer/IRAC mid-IR images (FWHM ~2"). At a distance of 8.3 kpc to the Galactic centre (Reid et al. 2014), this corresponds to ≲8300 au in projection, implying that the near-IR sightlines towards the background giant through the YSO envelope are unresolved within this scale and can be regarded as the near-IR counterparts of the mid-IR sources. Consequently, the projected separation between each near-IR sightline and the YSO centre cannot be directly measured. As an alternative, foreground extinction can serve as a proxy for projected separation, as the integrated dust column density within the YSO envelope generally increases towards its centre. However, the complex nature of foreground extinction fields of the CMZ - arising from patchy extinction in intervening Galactic disc clouds, additional obscuration by CMZ clouds, and dust within the YSO envelope - means that such an approach requires a careful understanding of the foreground dust properties.

In order to isolate the extinction contributed by a YSO envelope from the rest of the foreground extinction, we defined a reduced extinction as $$A_{K,0}^{\ast }\equiv A_K^{\ast }-⟨A_K^{\mathrm{f}\mathrm{g}}⟩$$(1)

where $⟨A_K^{\mathrm{f}\mathrm{g}}⟩$ is the mean foreground extinction from dust lying between the Sun and the YSO’s outer surface. Since A^*_K measures the total A_K from the Sun to the background giant, subtracting $⟨A_K^{\mathrm{f}\mathrm{g}}⟩$ yields a reduced A_K that more closely represents the extinction through the YSO envelope alone. For $A_K^{\mathrm{f}\mathrm{g}}$, we adopted the foreground extinction estimates from An et al. (2013), who measured the depth of the 9.7 μm silicate feature in spectra of the diffuse ISM, deliberately selected to exclude point sources. For most of our sample, such measurements are available at four orthogonal positions, each located approximately 1′ (equivalent to 2.4 pc in projection) from the target. We averaged the A_K values converted from the silicate depths and incorporated their standard deviation into the total $A_{K,0}^{\ast }$ uncertainty budget. These results are listed in the second and third columns of Table 5.

Alternatively, the fourth column of Table 5 shows the differential extinction, obtained by subtracting the foreground extinction (second column) from the mid-IR extinction $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$. Although $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ does not trace the column density towards a background giant but rather probes the line of sight into the YSO itself, one might expect that subtracting the foreground ISM extinction from the total mid-IR extinction would yield a positive residual, corresponding to the local envelope material. The presence of prominent mid-IR absorption features, such as the 15 μm CO₂ band, supports this expectation. However, even in sources that exhibit significant CH₃OH ice absorption, the residual extinction listed in the fourth column often appears consistent with zero or even slightly negative values. This result does not contradict the observed mid-IR absorption but rather reflects the limitations imposed by the measurement uncertainties.

Figure 9 shows the column densities of H₂O and CH₃OH ices as a function of $A_{K,0}^{\ast }$. As in the case of A^*_K (Fig. 4), N(solid H₂O) increases with extinction (Pearson’s correlation coefficient of r = 0.565), but the zero-point shift in $A_{K,0}^{\ast }$ lowers the minimum extinction to ~1 mag. The decline of N(solid CH₃OH) with $A_{K,0}^{\ast }$ closely parallels that with A^*_K (Fig. 7), although the trend is less significant than in the H₂O case (r = −0.367). Interpreting $A_{K,0}^{\ast }$ as a proxy for projected distance from the source centre, these trends suggest that N(solid H₂O) rises and N(solid CH₃OH) falls towards the source centre. The linear fits yield dispersions of ∆N = 1.2cm⁻² and 1.8 cm⁻² in the top and bottom panels, respectively, close to the expected uncertainties of 1.0 cm⁻² and 1.4 cm⁻². The small excess may reflect intrinsic variations in ice abundances due to differing temperature structures among YSOs of varying masses or evolutionary stages.

Nonetheless, $A_{K,0}^{\ast }$ extends from ~1 mag to nearly 12 mag, indicating a wide range of dust column densities along the various lines of sight in our sample. To assess whether dust extinction can be quantitatively linked to the projected distance from the source centre, we constructed a spherically symmetric gas model representing the envelope in the early phase of YSO evolution. Assuming a density power-law index of 1.5 (e.g. Jørgensen et al. 2002) and a total mass of 10-20 M_⊙ for our CMZ YSOs, as is inferred from multi-wavelength SED fitting (An et al. 2011), and adopting the gas-to-dust ratio from Bohlin et al. (1978), we estimated that the dust extinction through the envelope is 4 ≲ A_K ≲ 8 mag at 1" from the centre. At 0.5″, the extinction increases to 6 ≲ A_K ≲ 11 mag. This range of A_K is broadly consistent with $A_{K,0}^{\ast }$ in Fig. 9, supporting the view that our observed near-IR sightlines probe different regions of the extended envelopes of the CMZ YSOs.

Fig. 7 Methanol-ice column densities against foreground extinction estimates ($A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A*K in the top and middle panels) and water-ice column densities (bottom panel). The solid line in the top and middle panels shows the empirical relation for local clouds from Whittet et al. (2011), shifted horizontally by the average CMZ foreground extinction (∆AK = 2.2 mag).

Fig. 8 Mixing ratios of CH3OH ice to the shoulder component of CO2 ice. Circles and arrows in the top panels indicate measurements and 3σ upper limits, respectively, for the red CMZ objects in this study (left) and for YSOs in the Galactic disc (Pontoppidan et al. 2008, right). In the bottom two panels, probability density functions of the mixing ratio N(solidCH3OH)/N(CO2 shoulder) are shown for the CMZ and disc YSOs, respectively. Thin lines correspond to individual measurements, assuming a bivariate normal distribution for the measured column densities and their 1sgr; uncertainties; the thick line shows the sum of all distributions.

5.2 Average chemical structure of YSOs

The top panel of Figure 10 shows the fractional abundance of CH₃OH ice, defined as N(solid CH₃OH) divided by N (solid H₂O), for sources in our sample. The ratios are shown as a function of reduced extinction $A_{K,0}^{\ast }$, which serves here as a proxy for the projected distance from the source centre. The ratio peaks at ~0.3 near $A_{K,0}^{\ast }$ ≈ 2 mag, but generally stays below 0.1 for $A_{K,0}^{\ast }$ > 3 mag. To highlight the overall trend, the thick red line represents the weighted mean abundance ratio in four adaptive extinction intervals, including the single upper limit at $A_{K,0}^{\ast }$ ~ 1 mag. When viewed as an ensemble, the CMZ objects show a trend in which a typical CMZ YSO exhibits a sharp rise in methanol ice abundance in the outer envelope, while the abundance remains low (< 0.1) in the inner region.

For comparison, the thick blue line in the top panel of Fig. 10 shows the methanol abundance profile towards SMM 4, a nearby, well-characterised low-mass Class 0 YSO projected against the young cluster SVS 4 (Pontoppidan et al. 2004). To construct this profile, we used their ice abundance measurements, which include H₂O column densities along ten sightlines towards background stars spanning projected distances of ~4500-19000 au from the centre of SMM 4. The H₂O abundance relative to H₂ was found to be nearly constant across most of this range, except for a single data point at ~4700 AU that showed a significantly higher abundance. CH₃OH ice was detected in four sightlines within ~ 12 000 au, and the inferred abundance remained nearly constant across these positions, with a strong upper limit placed at the outermost radius (~ 19 000 au). In constructing the CH₃OH/H₂O ratio profile, we adopted their abundance values and assumed a constant CH₃OH abundance inwards of the innermost detection, based on their best-fitting envelope model. This results in a broad peak in the CH₃OH/H₂O ratio at intermediate radii, caused by the combination of flat CH₃OH abundance and a centrally increasing H₂O abundance.

As in SMM 4, the CH₃OH abundance of the CMZ ensemble exhibits a peak at intermediate radii within the envelope. In the cold and less dense outer envelope, methanol formation is suppressed due to the reduced efficiency of CO hydrogenation, likely caused by the limited mobility and availability of atomic hydrogen at low temperatures (e.g. Watanabe et al. 2003). While water ice can still accumulate in these conditions, methanol production becomes less efficient, leading to a declining N(solid CH₃OH)/N(solid H₂O) ratio with increasing distance from the source centre.

The peak CH₃OH abundance in the CMZ ensemble (≈0.3) is comparable to that observed towards SMM 4, suggesting that methanol formation via CO hydrogenation on grain surfaces may impose a common upper limit on CH₃OH ice production in starforming cores. However, the abundance peak in the CMZ sources spans a narrower A_K range than in SMM 4, potentially reflecting their higher total masses (~10-20 M_⊙; An et al. 2011) compared to ~5.8 M_⊙ for SMM 4 (Pontoppidan et al. 2004). The shorter evolutionary timescales of high-mass YSOs likely result in hotter internal structures, which in turn enhance CH₃OH ice sublimation in the inner envelope. As a consequence, methanol ice is efficiently thermally desorbed near the protostar, producing a characteristic abundance peak at intermediate depths.

This interpretation can further explain the observed low N(solidCH₃OH)/N(CO₂ shoulder) values of CMZ YSOs (Fig. 8). The bottom panel of Fig. 10 shows the trend of N(solid CH₃OH) relative to the CO₂ shoulder component, with the ratio peaking at lower $A_{K,0}^{\ast }$ (i.e. in the outer envelope regions) and gradually decreasing towards higher $A_{K,0}^{\ast }$ (i.e. the warmer, inner regions of the YSOs). The declining ratio with increasing $A_{K,0}^{\ast }$ suggests that CH₃OH is more efficiently depleted than CO₂ from the ice mantles in the inner envelopes. This trend is somewhat counter-intuitive, given the higher sublimation temperature of CH₃ OH compared to CO₂. However, CH₃OH may be efficiently destroyed or chemically transformed into more complex organic species as the ices are subjected to UV irradiation followed by thermal processing. UV photons dissociate CH₃OH into radicals, which then become mobile and recombine during warm-up to form complex molecules. This combined processing reduces the CH₃OH abundance relative to CO₂, which is more photo-stable and less reactive under similar conditions (e.g. Gerakines et al. 1996; Öberg et al. 2009). In addition, desorption driven by increasing gas temperature may remove both CO₂ and CH₃ OH; however, CH₃OH may appear preferentially depleted, as it is a minor component in the ice mixture. In contrast, CO₂ is more abundant and remains detectable through its prominent 15.4 μm shoulder feature.

Previously, Jang et al. (2022) proposed that the reduced CH₃ OH abundance could result from elevated oxygen abundance in the CMZ, which favours oxidation of CO over hydrogenation, thereby limiting CH₃OH production. However, the present result suggests that the lower methanol ice abundance may not directly arise from suppressed CH₃OH formation pathways. Instead, it may reflect a sample bias towards higher-mass YSOs, in which methanol ice has been more significantly depleted.

Fig. 9 Column densities of H2O ice (top) and CH3OH ice (bottom) as a function of reduced extinction, $A_{K,0}^{\ast }$. Since $A_{K,0}^{\ast }$ traces the dust column integrated along the line of sight within a YSO envelope, it serves as a proxy for projected distance from the centre of the system, with higher $A_{K,0}^{\ast }$ corresponding to smaller impact parameters.

Fig. 10 Relative abundance of CH3 OH ice in the CMZ sources, shown with respect to H2O ice (top panel) and to the CO2 shoulder component (bottom panel), as a function of reduced extinction $A_{K,0}^{\ast }$. The fractional abundance is defined as the column density of CH3OH ice divided by the sum of N(solid CH3OH) and either H2O or CO2 shoulder column densities, respectively. Circles show our measurements, and arrows mark 3σ upper limits. Only sources with significant detections (>2.5σ) of 3 μm H2O ice are included. In the top panel, the thick red line shows the weighted mean abundance in four adaptive intervals, while the thick blue line shows the abundance measured in the envelope of SMM 4, a low-mass YSO in the Serpens cloud (Pontoppidan et al. 2004). The approximate length scale shown at the top corresponds to the projected distance from the centre of SMM 4.

6 Summary

We have presented a spectroscopic study of 23 extremely red, point-like CMZ sources, combining new L-band data from Gem-ini/GNIRS (15 objects) with previously published IRTF/SpeX L-band data (8 objects), and supplementing these with archival Spitzer/IRS spectra to trace ice chemistry and foreground attenuation from 2 to 35 μm. Many targets show K-band CO band heads consistent with (super-)giant photospheres, while their mid-IR continua and strong ice bands are characteristic of embedded YSOs (An et al. 2017; Jang et al. 2022). This supports a backlighting geometry in which a background giant star illuminates a foreground YSO envelope (or dense core), enabling us to probe the envelope material through absorption features. We modelled each composite SED and derived both near- and mid-IR extinction estimates, along with robust ice column densities.

Water ice absorption is strong and ubiquitous among our sources, with N(solid H₂O) ~ (0.1-6.5) × 10¹⁹ cm⁻². The mid-IR extinction ($A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$) is typically 2-5 mag, whereas the extinction towards the background star (A^*_K) spans 4-15 mag, with A^*_K systematically exceeding $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ at fixed N(solid H₂O). This offset is consistent with the background star being behind additional envelope material. Solid CH₃OH is detected at > 2.5σ towards six objects; in the remainder we report meaningful upper limits.

To compare N(solid CH₃OH) measured in the near-IR with the CO₂:CH₃OH mixtures traced in the mid-IR (An et al. 2009, 2011), we revisited the 15.4 μm CO₂ ‘shoulder’ component. In the CMZ, the bulk ratio N(CH₃OH)/N(CO₂ shoulder) is ~2-5%, with a mean of 3.3 ± 0.7%, confirming earlier findings (Jang et al. 2022) that it is systematically lower than typical Galactic disc values (5-15%). This indicates either intrinsically reduced CH₃OH production or preferential depletion (or processing) of CH₃OH in CMZ envelopes. Exploiting the backlighting geometry, we introduced a reduced extinction ($A_{K,0}^{\ast }$) as a proxy for projected distance through an envelope. In this ensemble framework, N(CO₂ shoulder) increases with $A_{K,0}^{\ast }$, while N (solid CH₃OH) shows a declining trend; the ice abundance ratio CH₃OH/CO₂ is ~10% in inner lines of sight (large $A_{K,0}^{\ast }$) and rises steeply to ~30% in the outer envelope (small $A_{K,0}^{\ast }$). As our CMZ sample is likely biased towards more massive YSOs, the trends suggest that internal heating preferentially sublimates and/or chemically alters CH₃ OH relative to CO₂ in the inner envelope, lowering the line-of-sight CH₃OH/CO₂ while leaving high methanol abundances in the outer regions. This thermal/UV processing scenario offers an alternative to purely chemicalnetwork explanations (e.g. distinct elemental ratios in the ISM) for the CMZ’s low integrated CH₃OH/CO₂ ratios.

Our results demonstrate that (i) red CMZ point sources often comprise superposed near- and mid-IR components; (ii) joint modelling of near- and mid-IR spectra enables an estimate of the foreground extinction towards the backlit giant and thus the line-of-sight integrated column through the YSO envelope; and (iii) backlit spectroscopy can reveal ensemble-averaged radial abundance trends that would otherwise be inaccessible at the distance of the CMZ. Our conclusions are tempered by systematic uncertainties in cross-instrument flux scaling and by crowding and/or background-subtraction effects; moreover, a limited S/N at the base of the 9.7 μm silicate band precludes robust use of the CH₃OH C-O stretching mode. High-S/N, higher-resolution mid-IR spectroscopy and spatially resolved mapping (e.g. with JWST) of additional backlit sightlines will enable tests of the inferred abundance gradients and quantify the roles of sublimation, photochemistry, and mixing in setting the distinct ice chemistry of CMZ YSOs.

Acknowledgements

We thank the anonymous referee for constructive comments that helped improve the clarity and presentation of this manuscript. D.A. thanks Kris Sellgren and Solange Ramírez for their advice and support in the preparation of the Gemini observations. D.A. acknowledges the hospitality and support provided by the staff at the Gemini Observatory’s Hilo office during a research visit in 2017, which contributed to the development of this work. D.A. also thanks Minjin Kim for his assistance with the data acquisition during that visit. Y.K. and D.A. are grateful to Young-Jun Kim for discussions on YSO modelling. We acknowledge support provided by the National Research Foundation (NRF) of Korea grant funded by the Ministry of Science and ICT (RS-2021-NR058093). Based on observations obtained at the international Gemini Observatory, a programme of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comu-nicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea. We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance. NOIRLab IRAF is distributed by the Community Science and Data Center at NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation.

References

Appendix A Additional spectra of the sample

Figure A.1 presents the full set of observed composite spectra for all sources in our sample. Each panel follows the same format as the top panel of Fig. 3, which illustrates the representative case of SSTGC 799887. The observed data are shown in red, and the coloured curves represent the best-fitting model components as is described in Sect. 4.1, including the stellar continuum, the foreground extinction, and thermal dust emission, as well as major molecular ice absorption features. To avoid strong molecular features and ensure reliable continuum estimation, the model fitting was performed over the following wavelength intervals: 1.90-3.30 μm, 3.65-6.5 μm, 7.8-9.3 μm, 10.2-12.4 μm, 14.0-14.9 μm, 15.8-16.7 μm, 17.8-18.4 μm, 19.0-27.5 μm, and 28.8-30.0 μm. These plots are provided to illustrate the consistency and diversity of spectral shapes across the sample and to complement the parameter values listed in Table 5.

Fig. A.1 Same as Fig. 3, but showing the composite spectra of the remaining red CMZ sources.

Fig. A.1 Spectra of the sample sources (continued).

Fig. A.1 Spectra of the sample sources (continued).

Appendix B Additional optical depth spectra centred at 3.53 μm

Figure B.1 presents the remaining set of observed optical depth spectra centred on the 3.535 μm CH₃OH band. Each panel follows the same format as Fig. 6, which shows only the sources with significant detections. The observed spectra are shown in grey, and the red curves indicate the best-fitting models described in Sect. 4.2. These plots complement the fitting results listed in Table 4.

Fig. B.1 Same as Fig. 6, but showing the optical-depth spectra of the remaining CMZ sources with non-significant 3.535 μm absorption.

Appendix C Methanol ice column densities from spectra around 9.7 μm

In addition to the CH₃OH ice derived from the 3.535 μm absorption band, we used Spitzer/IRS spectra (An et al. 2011) to measure the peak optical depth and column density of CH₃OH ice from the 9.74 μm band. Because the band depth is highly sensitive to the precise shape of the local continuum and the feature lies within the core of the stronger silicate absorption, we adopted a third-order polynomial to define the continuum across the band rather than relying on the global SED silicate profile (Sect. 4.1). The results are listed in Table C.1, with uncertainties that include both statistical errors and systematics, such as those due to continuum placement. Owing to the very strong 10 μm silicate absorption towards our extremely red CMZ sources, the S/N at the band centre is low, and meaningful constraints could be derived for only a few objects.

Figure C.1 presents the optical-depth spectra of the three CMZ objects with significant (> 2.5σ) detections of the 9.74 μm CH₃OH absorption. The 3.535 μm measurements (Table 4) show a mild detection for SSTGC 653270, but SSTGC 425399 and SSTGC 770393 exhibit non-detections. The discrepancy between the near-IR (Table 4) and mid-IR (Table C.1) abundances may reflect the low S/N at the bottom of the silicate band affecting the 9.7 μm measurements, or could arise from dual sightlines, with the near- and mid-IR bands tracing different columns along the line of sight.

Fig. C.1 Optical depth spectra of CMZ objects showing the 9.74 μm absorption feature of CH3OH ice. Only the three sources detected at > 2.5σ are displayed.

Appendix D Maser observations of the full An et al. (2011) survey sample

We conducted a simultaneous survey of H₂O 623 - 5₂₃ (22.23508 GHz) and class I CH₃OH 7₀ - 6₁ A⁺ (44.06943 GHz) masers towards all 107 sources in June 2011 and April 2012. Observations were carried out with the Korean VLBI Network (KVN) Tamna telescope in single-dish mode (Kim et al. 2011; Lee et al. 2011). The main beam sizes were 130" at 22 GHz and 65" at 44 GHz, with corresponding aperture efficiencies of 0.70 and 0.69. Pointing and focus checks were performed every 2-3 h, with pointing accuracy better than 5". A digital spectrometer backend provided 4096 channels across a 64 MHz bandwidth. All spectra were obtained using position switching. Typical system temperatures were 100-200 K at 22 GHz and 150-350 K at 44 GHz. The on-source integration time was typically 10 min, yielding an rms noise level of ~0.5 Jy at a velocity resolution of 0.2 km s⁻¹. Data were calibrated with the standard chopper-wheel method and placed on the T^*_A scale. The conversion factors from T^*_A to flux density were 11.1 Jy K⁻¹ at 22 GHz and 11.6 Jy K⁻¹ at 44 GHz.

We detected H₂O maser emission in 7 sources (7%) - SSTGC 304239, 354683, 358370, 360559, 511261, 524665, and 381931 - and CH₃OH maser emission in 6 sources (6%) — SSTGC 511261, 524665, 543691, 564417, 536969, and 381931. Three sources (SSTGC 381931, 511261, and 524665) exhibit both water and methanol maser emission. The detection rates are much lower than those reported in previous surveys of high-mass protostellar objects and ultra-compact H II regions; for example, 42% for H₂O masers (Sridharan et al. 2002) and 31% for 44 GHz CH₃OH masers (Fontani et al. 2010). The detection rates of these two maser species tend to increase as the central objects evolve (Kim et al. 2018; Kim et al. 2019).

No H₂O or Class I CH₃OH maser emission was detected towards any of the YSOs in our sample except SSTGC 524665, despite their initial classification as massive YSO candidates. H₂O masers were detected both in YSOs and evolved stars, while CH₃OH masers were predominantly detected towards high-mass YSOs, with Class I masers tracing shocked gas in outflows (Menten 1991; Ellingsen 2006). The almost complete absence of maser emission in our sample therefore suggests that these objects may be less evolved than previously thought - possibly in the deeply embedded Class 0/I stage - or that their masses are overestimated, placing them in a regime where the necessary shock conditions for maser pumping are absent.

All Tables

All Figures

Fig. 1 GNIRS L-band spectra of 12 red sources in the CMZ. The observed spectra (red) are vertically offset for clarity and approximately ordered by increasing strength of the broad H2O ice absorption band centred at 3 μm. To illustrate the varying depth of this feature, the best-fit IRTF/SpeX M5.5III template spectrum (Rayner et al. 2009) is shown in green for each source. These templates were reddened using the foreground extinction estimate derived from our SED fitting procedure (see Sect. 4.1). Together with the 3 μm ice absorption profile, they reproduce the overall continuum shape and slope of the observed spectra.

In the text

	Fig. 2 Same as Figure 1, but for sources that have both K- and L-band spectra from GNIRS observations.
In the text

Fig. 3 Composite IR spectrum of a red CMZ object (SSTGC 799887) and modelling of the SED. The observed spectra are shown as a solid red line, constructed by combining near-IR data from Gemini/GNIRS with mid-IR data from Spitzer/IRS. Circles represent photometric measurements from UKIDSS and Spitzer/IRAC, which were used to set the flux scale of the spectra (see text). The solid blue line shows the best-fitting model to the observed spectra, which included the following components: (1) a template of an M5.5-type giant from the IRTF/SpeX spectral library (Rayner et al. 2009) (green line); (2) a two-component blackbody continuum (dashed blue and orange lines); (3) the 3 μm band of solid H2O ice, using the template from GC IRS 7 (Jang et al. 2022); and (4) the 10 μm and 18 μm silicate bands from the GC IRS 3 spectrum (Kemper et al. 2004). Other weaker absorption features were not included in the modelling. Spectra and their best-fitting models for the additional objects are shown in Appendix A.

In the text

Fig. 4 Distribution of extinction and H2O ice column densities derived from the global fits to the composite spectra. The two extinction estimates, $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A*K, come from our composite spectra but rely on different wavelength regimes: the mid-IR spectral shape for $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and near-IR template matching for A*K. Objects with large solid CH3OH column densities (N > 3 × 1017 cm−2; see below) are marked with open blue circles. The dashed line in the top panel indicates unity. For comparison, the solid lines in the middle and bottom panels show the empirical relation derived for the Taurus molecular cloud by Whittet et al. (2001), shifted horizontally by AK = 2.2 mag to account for the average foreground extinction towards the CMZ.

In the text

Fig. 5 Schematic illustration of the light-ray geometry probed by the composite spectra. A*K reflects extinction constrained mainly by the near-IR, while $A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ is determined primarily from the mid-IR. The quantity $A_K^{\mathrm{f}\mathrm{g}}$ represents the foreground extinction estimated from nearby source-free regions, which accounts for dust located along the line of sight in the Galactic disc and within the CMZ. The envelope-only extinction, $A_{K,0}^{\ast }$, is defined as $A_K^{\ast }-A_K^{\mathrm{f}\mathrm{g}}$, isolating the extinction produced by material within the protostellar envelope along the sightline to the background star.

In the text

	Fig. 6 Optical depth spectra of CMZ sources exhibiting 3.535 μm solid CH3OH absorption at >2.5σ significance. The red line shows the best-fit model for pure solid CH3OH ice, and the shaded area denotes the ±1σ confidence interval. Spectra for the remaining sources are provided in Appendix B.
In the text

	Fig. 7 Methanol-ice column densities against foreground extinction estimates ($A_K^{\mathrm{M}\mathrm{I}\mathrm{R}}$ and A*K in the top and middle panels) and water-ice column densities (bottom panel). The solid line in the top and middle panels shows the empirical relation for local clouds from Whittet et al. (2011), shifted horizontally by the average CMZ foreground extinction (∆AK = 2.2 mag).
In the text

Fig. 8 Mixing ratios of CH3OH ice to the shoulder component of CO2 ice. Circles and arrows in the top panels indicate measurements and 3σ upper limits, respectively, for the red CMZ objects in this study (left) and for YSOs in the Galactic disc (Pontoppidan et al. 2008, right). In the bottom two panels, probability density functions of the mixing ratio N(solidCH3OH)/N(CO2 shoulder) are shown for the CMZ and disc YSOs, respectively. Thin lines correspond to individual measurements, assuming a bivariate normal distribution for the measured column densities and their 1sgr; uncertainties; the thick line shows the sum of all distributions.

In the text

	Fig. 9 Column densities of H2O ice (top) and CH3OH ice (bottom) as a function of reduced extinction, $A_{K,0}^{\ast }$. Since $A_{K,0}^{\ast }$ traces the dust column integrated along the line of sight within a YSO envelope, it serves as a proxy for projected distance from the centre of the system, with higher $A_{K,0}^{\ast }$ corresponding to smaller impact parameters.
In the text

Fig. 10 Relative abundance of CH3 OH ice in the CMZ sources, shown with respect to H2O ice (top panel) and to the CO2 shoulder component (bottom panel), as a function of reduced extinction $A_{K,0}^{\ast }$. The fractional abundance is defined as the column density of CH3OH ice divided by the sum of N(solid CH3OH) and either H2O or CO2 shoulder column densities, respectively. Circles show our measurements, and arrows mark 3σ upper limits. Only sources with significant detections (>2.5σ) of 3 μm H2O ice are included. In the top panel, the thick red line shows the weighted mean abundance in four adaptive intervals, while the thick blue line shows the abundance measured in the envelope of SMM 4, a low-mass YSO in the Serpens cloud (Pontoppidan et al. 2004). The approximate length scale shown at the top corresponds to the projected distance from the centre of SMM 4.

In the text

	Fig. A.1 Same as Fig. 3, but showing the composite spectra of the remaining red CMZ sources.
In the text

	Fig. A.1 Spectra of the sample sources (continued).
In the text

	Fig. A.1 Spectra of the sample sources (continued).
In the text

	Fig. B.1 Same as Fig. 6, but showing the optical-depth spectra of the remaining CMZ sources with non-significant 3.535 μm absorption.
In the text

	Fig. C.1 Optical depth spectra of CMZ objects showing the 9.74 μm absorption feature of CH3OH ice. Only the three sources detected at > 2.5σ are displayed.
In the text