© The Authors 2023

1 Introduction

Astronomical methanol (CH₃OH) masers were first reported in the Orion Kleinmann-Low (KL) nebula (Barrett et al. 1971), shortly after CH3OH emission had been discovered in the Galactic centre (Ball et al. 1970). Subsequent observations have resulted in the detection of numerous CH₃OH maser transitions, which have been proven to be important tracers of star formation regions (e.g. Menten 1991). These CH₃OH masers are divided into two categories, based on their different observational properties and pumping mechanisms (Batrla et al. 1987; Menten 1991). Class I CH₃OH masers are thought to be colli- sionally pumped (e.g. Leurini et al. 2016), and are often found to be offset from embedded infrared sources and ultracompact HII (UCHII) regions. In contrast, class IICH₃OH masers are thought to be radiatively pumped (e.g. Cragg et al. 2005) by the infrared radiation emitted by massive young stellar objects in whose environments they are found (Walsh et al. 1998). Only a small portion of these masers produce detectable radio continuum emission that can ionise UCHII regions (Nguyen et al. 2022).

Based on previous studies (see e.g. Table 1 in Leurini et al. 2016), the three methanol transitions at 84 GHz, 95 GHz, and 104.3 GHz are the brightest class I methanol masers in the 3 mm band. The methanol emission in the J_k = 5₋₁ − 4₀ E line near 84 GHz belongs to the same line series of strong class I masers as the 36 GHz 4₋₁ − 3₀ E transition (i.e. they connect the same k-ladders). Maser emission at 84 GHz was discovered by Batrla & Menten (1988) and Menten (1991) towards DR21(OH), NGC2264, and OMC-2. Since 2001, several extensive searches for this transition have been conducted towards sources associated with other class I masers, and associated with young bipolar outflows in regions where low- to intermediate-mass stars are formed (Kalenskii et al. 2001, Kalenskii et al. 2006; Rodríguez-Garza et al. 2018; Breen et al. 2019). The detection rate of the 84 GHz methanol line, which can show narrow maser-like features and broad quasi-thermal emission, in these targeted surveys is greater than 70%. Methanol emission in the 95 GHz J_K = 8₀ − 7₁ A⁺ transition belongs to the same line series (common K quantum numbers) as the strong and widespread 7₀ − 6₁ A⁺ maser at 44 GHz (Haschick et al. 1990). The number of known 95 GHz methanol masers has been significantly increased by systematic surveys, including towards (1) extended green objects (EGOs) identified in the 4.5 µm band of the Spitzer Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE), whose emission is dominated by rotationally excited shocked H₂ (Chen et al. 2011, 2013); (2) molecular outflow sources (Gan et al. 2013); (3) many sources associated with both GLIMPSE point sources and 1.1 mm dust continuum emission detected in the Bolocam Galactic Plane Survey (BGPS; Chen et al. 2012); (4) a large sample of BGPS sources (Yang et al. 2017, 2020); and (5) red Midcourse Space Experiment sources (Kim et al. 2018). Ladeyschikov et al. (2019) compiled an online database of class I methanol masers and reported a total of 129 and 534 maser detections at 84 and 95 GHz, respectively, which suggests that both masers are widespread across the Milky Way.

On the other hand, the 11₋₁ − 10₋₂ E methanol maser at 104.3 GHz is rarely seen. Maser emission in this line was first predicted by theoretical calculations (Voronkov 1999), and was successfully detected towards W33-Met (Voronkov et al. 2005). A follow-up Mopra 104.3 GHz survey led to only two detections out of 69 targets, suggesting that this maser line is rare (Voronkov et al. 2007). Theoretical calculations suggest that more energetic conditions (i.e. higher temperatures and densities) are required to produce these rare class I methanol masers (Sobolev et al. 2005; Voronkov et al. 2012) than the widespread class I methanol masers at 36, 44, 84, and 95 GHz. Prior to this work, only five sources (G019.61−0.23, G305.21+0.21, G357.97-0.16, IRAS16547−4247 and W33-Met) were known to harbour 104.3 GHz methanol masers (e.g. Voronkov et al. 2006). However, the overall incidence of the 104.3 GHz methanol masers in the Milky Way is still poorly constrained.

It has been established that class I methanol masers are closely related to shocked regions. Based on high angular resolution observations towards DR21/DR21(OH), Plambeck & Menten (1990) suggested that class I methanol masers trace the interface between outflows and ambient dense clouds. Subsequent observations supported this scenario and found a coincidence between class I masers and other molecular shock tracers (e.g. SiO, Kurtz et al. 2004; HNCO, Gorski et al. 2018; shocked H₂ emission at 4.5 µm and 2.12 µm, Cyganowski et al. 2009; Voronkov et al. 2006). However, the relationship between class I masers and shock tracers has not been well studied in a detailed statistical manner.

The APEX telescope large area survey of the Galaxy (ATLASGAL) is an unbiased 870 µm submillimetre continuum survey of the inner Galaxy (300° < l < 60°, |b| < 1.5°) (Schuller et al. 2009). This survey provides a large inventory of dense molecular clumps (over ~10 000, Contreras et al. 2013; Csengeri et al. 2014; Urquhart et al. 2014) that comprise a variety of early evolutionary stages related to high-mass star formation. The evolutionary stages of a large number (~5000) of ATLASGAL clumps have been determined and their physical properties (e.g. distance, dust temperature, bolometric luminosity, clump mass, and H2 column density) have been studied in detail (König et al. 2017; Giannetti et al. 2017; Urquhart et al. 2018, Urquhart et al. 2022). This well-characterised sample allowed us to carry out a targeted survey searching for new class I methanol masers and perform a statistical analysis of the properties of detected masers and their associated star-forming clumps at different evolutionary stages.

Silicon monoxide (SiO) emission is an excellent tracer of shock interactions in star formation regions. The abundance of SiO is enhanced by up to six orders of magnitude in shocked regions associated with molecular outflows with respect to that measured in quiescent gas (e.g. Jiménez-Serra et al. 2008) through grain-grain collisions and sputtering, after which the released atomic silicon forms SiO (e.g. Martin-Pintado et al. 1992; Caselli et al. 1997). The SiO emission towards our targeted ATLASGAL clumps has been well studied by Csengeri et al. (2016). This provides us a unique opportunity to statistically study the relationship between the class I methanol masers and the shocked gas traced by SiO emission.

In this work, we focus on the three class I methanol maser transitions by analysing 3 mm line survey data obtained using the Institut de Radioastronomie Millimétrique (IRAM) 30-meter telescope on 408 ATLASGAL clumps. The aims of this work are (i) to search for new class I maser sources; (ii) to statistically study the relationships between class I masers and a shock tracer; (iii) to investigate the properties between class I masers and their host clumps, as well as a function of their evolutionary stages; and (iv) to provide strong constraints on environments in which multiple class I masers occur. Also covered by our λ ~3 mm data are six additional methanol transitions that have detected maser features that belong to class II. These include the 85.5, 86.6, 86.9, 104.1, 107, and 108 GHz lines. We will focus on these class II methanol maser transitions, especially the maser emission and absorption features at 107 GHz, towards the same 408 ATLASGAL sources in a subsequent paper.

The structure of this paper is as follows. Our observations and data reduction are described in Sect. 2. The results are given in Sect. 3. In Sect. 4 we discuss the properties of the associated ATLASGAL clumps and SiO emission for our detections, and perform theoretical calculations to constrain the physical conditions of maser environments. A summary of this work and highlighted conclusions are presented in Sect. 5.

2 Observations and data reduction

The methanol transition data were extracted from unbiased spectral surveys with a frequency coverage of 83.8–115.7 GHz using the Eight MIxer Receiver (EMIR, Carter et al. 2012) at 3 mm (E090) of the IRAM 30-meter telescope (see Csengeri et al. 2016 for details, project IDs: 181-10, 049-11, and 037-12), from 2010 May to 2012 October. Flux-limited samples of ATLASGAL sources (Schuller et al. 2009; Csengeri et al. 2014) were selected, with additional infrared selection criteria that ensure coverage of a large range of evolutionary stages and luminosities. The sample includes (i) ~120 brightest ATLASGAL clumps that are infrared bright; (ii) ~50 bona fide massive young stellar objects that fulfil the Lumsden et al. (2002) infrared colour criteria for embedded massive (proto)stars and are associated with ATLASGAL clumps; (iii) ~120 brightest ATLASGAL clumps that are GLIMPSE 8 µm dark, but contain 24 µm sources; and (iv) ~ 120 brightest ATLASGAL clumps that are Spitzer 8 and 24 µm dark.

The properties and observational parameters of these transitions are summarised in Table 1. The FFTS backend with a frequency resolution of ~200 kHz was used, resulting in a channel spacing of ~0.7 km s⁻¹ (multiplied by 1.16 to convert to velocity resolution; see Klein et al. 2012). Table 2 provides information on the observed sources. The half power beam width (HPBW) is ~30″ at 84 GHz. The observations were performed in position switching mode with an offset of 10’ in right ascension and declination for the reference position. The data were processed using the GILDAS/CLASS package (Pety 2005). To characterise the spectra, we employed a Gaussian fit method for each transition source by source, and a multi-Gaussian fit was used in the case of multiple velocity components.

3 Results

We consider emission above the 3σ level as a detection, which corresponds to a typical detection threshold of ≳0.6 Jy for our three methanol transitions. Table 3 gives an overview of detections for individual sources in the three class I maser lines. Overall, we detect 282, 224, and 29 sources with methanol emission at 84 GHz (282/404, 70%), 95 GHz (224/408, 55%), and 104.3 GHz (29/404, 7%)¹, respectively.

Figure 1 shows the number of sources detected in emission in the three class I methanol maser lines. A total of 212 sources (over 50% of the observing targets) are detected with methanol emission at both 84 and 95 GHz (see Fig. 1). Seventy sources show methanol detections only at 84 GHz and nine sources have detections only at 95 GHz. All 29 detected 104.3 GHz emission lines in our sample have associated methanol emission at both 84 and 95 GHz, corresponding to 14% of the 212 sources (see Fig. 1). The detection rate decreasing from the 84 GHz transition to the 104.3 GHz transition, is readily explained by their relative energy levels and critical densities (see Table 1).

Table 4 provides Gaussian fitting results of detected lines for each source. Figure A.1 shows example spectra of the three class I (and the other six class II methanol maser transitions at 3 mm) towards the first six ATLASGAL sources in our sample. The following subsections present the details of the detections of the widespread class I masers (at 84 and 95 GHz) and rare class I masers (at 104.3 GHz).

3.1 84 and 95 GHz methanol emission

Figure 2 shows diverse 84 and 95 GHz spectral profiles, from narrow maser-like features (e.g. G06.22−0.61; see Fig. 2a) to broad (quasi-) thermal-like components (e.g. G10.47+0.03; see Fig. 2b), and a combination of the two (e.g. G32.02+0.06; see Fig. 2c). The line profiles of the two transitions could be quite similar (Fig. 2a–c) or could be very different (see Figs. 2d–f). G24.79+0.08 (see Fig. 2d) is an example showing different numbers of narrow features in the two lines. Emission from both transitions is not always detected towards a given source (e.g. G36.90−0.41 and G37.48−0.10 in Figs. 2e–f).

Typical spectral maser profiles show strong and narrow emission features. However, with single-dish observations, we cannot definitively determine whether the observed sources show only maser emission or a combination of maser and thermal emission. In this work, an object showing narrow (≲2 km s⁻¹, i.e. three channels) and relatively strong (≳1 Jy, i.e. 5 σ) features (e.g. Fig. 2c for the 84 GHz line towards G32.02+0.06) is likely to be a maser source (flagged ‘Y’ in Table 3). Otherwise, objects showing other features are considered as maser candidates (flagged ‘y’). A lack of emission above the 3σ level is rated as a non-detection (flagged ‘N’). For example, G06.22−0.61 shows narrow and strong maser features in both transitions (Fig. 2a). On the other hand, the 84 GHz emission in G36.90−0.41 (Fig. 2e) is broader than 2 km s⁻¹. Thus, this source is considered a maser candidate.

We note that narrower emission (<1 km s⁻¹) could be better detected with higher spectral resolution. With our velocity resolution of ~0.7 km s⁻¹, such narrow emission could be smoothed out over channels, appearing broader and weaker. For example, the narrow and strong component of 95 GHz emission in G06.22−0.61 shows a line width of 0.8 km s⁻¹ and a peak flux density of 59 Jy at a velocity resolution of 0.11 km s⁻¹ significantly better than ours (Chen et al. 2013). Assuming that this maser component is not variable, our observations nevertheless yield a line width of 1.2 km s⁻¹ and a peak flux density of 38 Jy, indicating a constant integrated flux density, but blurred peak and line width parameters. Thus, it is not arbitrary, but potentially stringent, to use a line width of 2 km s⁻¹ as a criterion for determining narrow maser features in this work. Based on this criterion, a total of 54 and 100 sources show maser features at 84 and 95 GHz, respectively.

After cross-matching the maser database² (Ladeyschikov et al. 2019), we find that 260 and 101 sources with methanol emission at 84 and 95 GHz, respectively, are detected for the first time. Among them, 50 masers at 84 GHz and 29 masers at 95 GHz are new discoveries. For 84 GHz masers with known distances, 78% of them are located within 5 kpc, 16% of them are located between 5 and 10 kpc, and 6% of them are located farther than 10 kpc. The situation is similar for the 95 GHz masers. Since an angular resolution of 26″ corresponds to a linear scale of 0.6 pc at a distance of 5 kpc, our analysis primarily focuses on the scales of ~0.6 pc, which we refer to as the clump scale.

Due to the close proximity of some sources to each other (within 30″, nearly the beam size at 3 mm), the observing beams may overlap and the detected emission could be from the same region. Thirteen pairs of sources are in close proximity, flagged by a dagger (†) near their source names in Table 3. Among them, four pairs (G14.19−0.19 and G14.20−0.19; G14.63−0.57 and G14.63−0.58; G34.40+0.23 and G34.40+0.23A; G59.50−0.23 and G59.50−0.24) show narrow maser features at 95 GHz. However, their small number means that they make negligible contributions to the overall statistics.

Interferometric observations reveal that 95 GHz class I maser emission is stronger than 84 GHz maser emission detected in six spots towards IRAS16547−4247 (see Figs. 2 and 3 in Voronkov et al. 2006). Our observations allow us to determine whether the 84 GHz or 95 GHz masers are brighter in a statistical manner using single-dish data.

Figure 3 shows the histogram of the peak flux density ratios (S _pk,95/S _pk,84) of 95 and 84 GHz methanol emission. The median and mean values of peak flux density ratios for all 212 sources showing both emissions are 0.65 and 0.83, respectively. Half of the whole sample (107/212) have both maser candidates, and 92% of the sources have stronger 84 GHz emission than their 95 GHz counterparts. The S _pk,95/S_pk,84 ratios for these 107 sources range from 0.2 to 4.0 with a median value of 0.5 and a mean value of 0.6. If thermal emission dominates methanol detection in these sources, it is reasonable that 84 GHz emission is stronger than 95 GHz emission since the upper level energy and critical density of the former are much lower than for the latter.

On the other hand, 40 sources host both masers, and show a different distribution of peak flux density ratios compared with the sources hosting both maser candidates in Fig. 3. Among them, we find that 32 (80%) sources show stronger 95 GHz maser emission. The S _pk,95/S_pk,84 ratios range from 0.3 to 3.1 with a median value of 1.2 and a mean value of 1.3. This suggests that the 95 GHz methanol masers are usually stronger than their 84 GHz maser counterparts. The trend is different from what we find for the maser candidates, probably because the emission at 84 GHz may contain more emission of thermal origin for the candidates in Fig. 3.

Fig. 1 Sources detected in the three class I maser lines. The Venn diagram shows that all detected 104.3 GHz emission lines in our sample have associated methanol emission at both 84 and 95 GHz.

3.2 104.3 GHz methanol emission

We take the criteria of narrow (≲2 km s⁻¹) and relatively strong (≳1 Jy) features to determine maser features in the 104.3 GHz spectra, which are the same as those for the 84 and 95 GHz transitions. Among the 29 detected methanol sources at 104.3 GHz, four possess narrow and relatively strong maser characteristics (see Fig. 4), and the other 25 sources are regarded as maser candidates. The four 104.3 GHz masers are detected for the first time. Our work increases the number of known 104.3 GHz masers from five to nine.

As can be seen from Fig. 4, the 104.3 GHz maser emission generally aligns with the 84 and/or 95 GHz maser emission in velocity, indicating that for all three lines the emission arises from the same source. However, compared to the 84 and 95 GHz spectra, the 104.3 GHz maser spectra appear to show fewer components and narrower line widths, which is consistent with a previous study (Voronkov et al. 2007). An explanation may be that the spatial distribution of the 104.3 GHz masers is not as extensive as that of the 84 and 95 GHz masers.

This scenario is supported by interferometric observations towards IRAS16547−4247, where the rarer class I methanol maser emission in the 9.9, 25, and 104.3 GHz lines is confined to only one of the six 84 and 95 GHz maser spots, which is associated with the brightest shocked H₂ 2.12 µm emission (see Fig. 1 in Voronkov et al. 2006).

The sources G10.34−0.14 and G16.58−0.05 both show 104.3 GHz maser emission and both have associated EGOs (Cyganowski et al. 2008), which are believed to be caused by shocked molecular gas in active outflows (see above). Towner et al. (2017) detected class I methanol maser emission in the 5₂ − 5₁ E line at 25 GHz towards these two sources. We cross-matched our sample with the EGO sample associated with 25 GHz methanol maser observations (Towner et al. 2017). We find three sources (G12.68−0.18, G14.33−0.64, and G22.04+0.22) showing 25 GHz masers and 104.3 GHz maser candidates and five sources (G10.29−0.12, G11.92−0.61, G12.90-0.03, G35.03+0.35, and G45.47+0.05) showing 25 GHz masers, but without 104.3 GHz detections at the current sensitivity. This indicates that the 104.3 GHz masers require more energetic conditions than the 25 GHz masers.

The four detected 104.3 GHz maser sources are associated with nearby 6.7 GHz class II masers (e.g. Green et al. 2010; Article number, page 6 of 23 Breen et al. 2015; Hu et al. 2016), and three of them also have a 12.2 GHz class II maser association, with the exception of G28.20−0.05 (e.g. Breen et al. 2014, 2016; Song et al. 2022). Although G28.20−0.05 does not have an associated 12.2 GHz maser (Breen et al. 2010), it has been classified as a source in a later evolutionary stage due to the presence of an HII region (Urquhart et al. 2022). G16.58−0.05 and G27.37−0.17 are also classified as being in the HII region stage, while G10.34−0.14 is classified as a PDR+Embedded source because of the complex infrared background (Urquhart et al. 2022). Voronkov et al. (2006) detected a 104.3 GHz maser towards IRAS16547−4247, but no 6.7 GHz class II maser, which suggests that this source is either too young or too old to have a 6.7 GHz maser. Urquhart et al. (2022) classifies its host dust continuum source (AGAL343.128−0.062) to be in the HII region stage. Overall, these facts support that the 104.3 GHz masers are likely associated with dust clumps that may already host HII regions in Urquhart et al. (2022).

Fig. 2 Examples of diverse line profiles of observed class I methanol emission at 84 GHz (upper panel) and 95 GHz (lower panel), depicting narrow maser-like features (panel a), broad thermal-like features (panel b) and a combination of the two (panel c). The line profiles of the two transitions are sometimes quite similar (panels a to c) and sometimes very different (panels d to f). The frequency and the corresponding quantum numbers of each transition are shown in the upper right corner of each panel for six sources. The horizontal magenta dashed lines in G36.90−0.41 and G37.48−0.10 represent the 3σ noise level for the undetected transitions.

Fig. 3 Peak flux density ratio of 95 and 84 GHz emission. The grey shaded histogram represents the whole sample with both 95 and 84 GHz detections, the pink line shows the sources hosting both 84 and 95 GHz masers, and the black line shows the sources hosting maser candidates in both transitions. The red and blue dashed lines depict the median and mean values of the whole sample.

3.3 Systemic velocity

Several studies have revealed that class I methanol masers are a good tracer of systemic velocities (e.g. Jordan et al. 2017). In this work, we adopt V_LSR from Urquhart et al. (2022) as the ATLASGAL clump systemic velocities to study the relationship between the velocity of class I methanol masers and systemic velocity.

Table 5 gives the numbers of maser sources and all methanol detections showing single and multiple components in our three methanol transitions, as well as the fractions of masers (i.e. the number of masers divided by the number of all detection). It is clear that the majority of the methanol sources detected at 84, 95, and 104.3 GHz shows a single component, while the detection rate of masers in multi-component sources is higher than that in the single component sources. This is expected, because multiple components can be present and distributed in different spots within a telescope beam.

Figure 5 shows the histograms of peak velocity of methanol emission with respect to the systemic velocity ($V_{{\text{CH}}_3\text{OH}}-V_{\text{sys}}$) at 84 GHz, 95 GHz, and 104.3 GHz. Apart from some outliers (see description below), the distributions of the velocity difference show Gaussian profiles for both the 84 and the 95 GHz transitions. The value of $V_{{\text{CH}}_3\text{OH}}-V_{\text{sys}}$ ranges from −3.83 km s⁻¹ to 3.31 km s⁻¹ for the 84 GHz masers, and from −5.09 km s⁻¹ to 4.23 km s⁻¹ for the 95 GHz masers. The mean and median velocity difference of the 84 GHz (95 GHz) maser sample is 0.17 km s⁻¹ (0.12 km s⁻¹), and 0.03 km s⁻¹ (0.01 km s⁻¹), with a standard error of 0.19 km s⁻¹ (0.16 km s⁻¹). For all methanol detections (i.e. including masers and maser candidates) at 84 and 95 GHz, the mean and median values are similar to those of our maser samples. The standard deviations (1σ) of the velocity difference are 1.33 and 0.92 km s⁻¹ for methanol masers, and all detections at 84 GHz, and 1.54 and 1.40 km s⁻¹ for the masers and all detections at 95 GHz. All of the 84 GHz masers and 99% of the 95 GHz masers (with the only exception of G23.44−0.18 showing the greatest velocity difference) are consistent with the systemic velocity within 3e. Overall, the velocities of 84 and 95 GHz methanol masers are consistent with the systemic velocity, which confirms that the velocity of the 84 and 95 GHz methanol maser emission traces the systemic velocity well (e.g. Jordan et al. 2017).

There are some outliers (84 and 95 GHz methanol emission in G43.17+0.01, and 95 GHz maser emission in G31.02+0.26), showing large velocity differences (>10 km s⁻¹) relative to the systemic velocity listed by Urquhart et al. (2022). However, after examining our own 3 mm spectral line survey data, we found that the velocity of methanol emission is consistent with that of C¹⁸O (1−0) emission. G43.17+0.01 is located in the massive and luminous W49 region, so that complex velocity profiles are not surprising. Both the 84 and 95 GHz lines in this source show two blended broad components (ΔV > 4 km s⁻¹). In the case of the 84 GHz line, the stronger feature peaks at 13 km s⁻¹, and the weaker and broader feature peaks at 7.3 km s⁻¹ (4.6 km s⁻¹ in case of the 95 GHz transition), which is similar to what is found for the double-peak profile of C¹⁸O (1−0) emission, whose two components peak at 4.4 and 11.9 km s⁻¹. In G31.02+0.26, only a 95 GHz maser single component is detected. Its velocity coincides with that of the second strongest C¹⁸O (1−0) emission component (at ~78 km s⁻¹), but not with that of the strongest one (at ~96 km s⁻¹). Therefore, the consistency between the maser velocities and systemic velocities still holds. Previous observations of 84 and 95 GHz masers towards DR21(OH) suggest that they are formed in the interface between outflows and ambient dense gas (Batrla & Menten 1988; Plambeck & Menten 1990). The masers are thus expected at the cloud systemic velocities. Our statistical results indicate that this scenario should apply to most of class I methanol masers in the Milky Way.

Similar to class I methanol masers, water masers are also believed to be caused by collisional pumping in star formation regions (e.g. Elitzur 1992), but water masers show a different behaviour in that they can exhibit multiple maser features within a wide velocity range; in extreme cases such as W49N, the velocity of water masers deviates from the systemic velocity more than 250 km s⁻¹ (e.g. Morris 1976; Kramer et al. 2018), and they are thought to trace motions of gas bullets powered by stellar winds or jets emitted from massive YSOs (e.g. Sanna et al. 2010).

Due to the limited number of methanol masers detected at 104.3 GHz, the velocity differences shown in Fig. 5 do not follow a Gaussian profile. It is worth noting that the four detected 104.3 GHz methanol masers are all clearly red- shifted with respect to the systemic velocity (see also Fig. 4). G27.37−0.17 shows the maximum offset from the systemic velocity of 3.3 km s⁻¹, and G28.20−0.05 shows the minimum offset of 1.7 km s⁻¹. A faint and broad feature at 104.3GHz is also detected in G10.34−0.14 and G28.20−0.05 that could be of thermal origin and is aligned with the systemic velocity (see Figs. 4a and d). All four 104.3 GHz maser sources have nearby 22.2 GHz H₂O maser associations (e.g. Breen & Ellingsen 2011; Walsh et al. 2014).

Fig. 4 Spectra of the 84, 95, 104.3 GHz lines (from top to bottom) for 4 sources showing 104.3 GHz methanol maser features. The vertical magenta dashed lines represent the systemic velocity adopted from Urquhart et al. (2022). The frequency and the corresponding quantum numbers of each transition are shown in each panel.

4 Discussion

4.1 Relationship between class I methanol masers and SiO emission

SiO is one of the best tracers of shocks in star formation regions. The abundance of SiO is enhanced in shocked regions (e.g. Jiménez-Serra et al. 2008) and the emission profile varies with shock velocity. In low-velocity shocks, the SiO emission shows a narrow profile (with a FWHM line width of ≤1–2 km s⁻¹, Jiménez-Serra et al. 2009), while the SiO profiles are broader in high-velocity shocks.

Cross-matching with the SiO data of Csengeri et al. (2016), we find that nine sources showing narrow and strong methanol maser characteristics are not associated with SiO emission. Among these SiO non-detections, four sources (G10.75+0.02, G13.87+0.28, G14.18-0.53, and G17.64+0.15) harbour 84GHz masers and seven sources (G08.05−0.24, G10.75+0.02, G13.87+0.28, G27.55−0.94, G31.02+0.26, G31.10+0.26, and G37.48−0.10) harbour 95 GHz masers. The abundance of gas phase methanol may be enhanced by low-velocity shocks in which the high-velocity shock tracer SiO is harder to detect. G10.75+0.02, G14.18−0.53, and G31.02+0.26 are classified as sources in a quiescent stage of evolution (Urquhart et al. 2022), suggesting that class I methanol masers could be a unique signpost of protostellar activity in extremely embedded objects at the earliest evolutionary stage.

Figure 6 shows 84, 95, and 104.3 GHz integrated intensities of methanol masers and maser candidates plotted against integrated SiO (2−1) and SiO (5−4) intensities. Their linear fitting results in log-log form for three methanol transitions are summarised in Table 6. Given that the correlation coefficients and slopes of the fitting lines are almost identical within the errors, we infer that there is no significant difference between masers and maser candidates in the 84 and 95 GHz lines. On the other hand, the positive slopes and correlation coefficients for the two lines indicate that stronger class I methanol maser emission is typically associated with higher SiO integrated intensity. Despite the small number of sources showing the 104.3 GHz emission, the integrated intensities of the 104.3 GHz emission appear to also be linearly correlated with the integrated SiO (2−1) and SiO (5−4) intensities. This is consistent with previous studies reporting similarly positive correlations between class I methanol masers at 36, 44, and 84 GHz and SiO emission (Jordan et al. 2017; Breen et al. 2019). In order to minimise the contribution of thermal emission of a maser source, we choose the 95 GHz maser sources and maser candidates with a single component, and do least-squares fitting separately. Figure 7 shows no significant difference between the least- squares fitting results for all masers and those for masers with a single component, suggesting that the potential contribution of thermal emission is negligible. Hence, the scaling relationships support the close link between the three class I methanol masers and SiO emission.

SiO emission with high-velocity wings indicates shocked gas due to fast material ejection. Csengeri et al. (2016) detected high- velocity SiO (2−1) wings towards 167 sources in our sample, while 120 sources do not show such wings. Figure 8 reveals that there are more detections and higher detection rates of the three class I methanol transitions (and for methanol maser emission) towards sources showing wings than sources without wings. The uncertainties of the detection rates in this work are calculated under the assumption that the detection rates satisfy the binomial distribution, that is $\sqrt{p\times \left(1-p\right)/n}$, where p is the detection rate and n is the sample size. Figure 9 shows that methanol masers have higher integrated intensities in the sources showing SiO (2−1) wings than in the sources without SiO (2−1) wings. This indicates that brighter class I methanol masers are produced by higher shock velocities. In order to further investigate the relationship, we adopt the full width at zero power (FWZP) for this analysis in the following.

Beuther & Sridharan (2007) and Li et al. (2019) divided their observed SiO spectra by the FWZP into several regimes, and considered that sources with FWZP >20 km s⁻¹ indicate gas at relatively high velocities with respect to the ambient gas. We use the FWZP of SiO (2−1), which is taken from Csengeri et al. (2016), to roughly characterise shock speed, and consider that the broader FWZP traces the faster shock speed. We perform a Kolmogorov-Smirnov (K-S) test to study the FWZP distribution for the sources with and without methanol detection. Figure 10 shows that the FWZP of SiO emission for sources with and without methanol detection are significantly different, indicated by a very small p-value (≪0.0013) from the K-S test. This difference is shown in all three of our methanol transitions. Moreover, the sources with methanol detection show broader FWZPs of SiO emission than sources without methanol detection. The median value of the FWZP of SiO emission for sources showing 84 or 95 GHz emission is about 30 km s⁻¹, while the SiO emission generally exhibits a larger FWZP with a median of about 40 km s⁻¹ in sources detected in the 104.3 GHz transition. This suggests that the excitation of 104.3 GHz methanol masers requires faster shock velocities than the 84 and 95 GHz masers. Figure 11 shows that the differences of SiO (2−1) FWZP between masers and maser candidates are not statistically significant at both 84 and 95 GHz. This insignificant difference indicates that the shock speeds at which maser emission and maser candidate emission are created are similar on clump scales. It may also be due to our conservative classification that a significant portion of maser candidates host masers.

Figure 12 shows the detection numbers and detection rates of the three methanol transitions in different FWZP ranges of SiO (2−1) emission. For the 84 and 95 GHz transitions, the detection numbers first increase with increasing FWZP of SiO emission and peak at a FWZP range of 20–30 km s⁻¹, then decrease as the FWZP increases. When the FWZP < 10 km s⁻¹, few sources show methanol emission, and only two out of four 95 GHz detections show maser features. With FWZP > 30 km s⁻¹, the numbers of detected 84 and 95 GHz emission gradually decrease, the detection rates for both lines keep rising (due to a smaller number of sources), and a detection rate of 100% for FWZP > 50 km s⁻¹ is reached. For the 104.3 GHz transition, the detection rate also shows a tendency to increase with increasing FWZP, when FWZP < 50 km s⁻¹. With FWZP > 50 km s⁻¹ the detection rate drops slightly. This indicates that the faster shock speeds lead to the higher detection rates of methanol emission when shocks have not yet completely destroyed SiO and methanol molecules. This is consistent with the astrochemical modelling results that the methanol abundance is correlated with shock velocities (e.g. see Fig. 3 in Holdship et al. 2017).

We further quantify the relationship between methanol emission and the SiO (2−1) FWZP in Fig. 13, which shows upward trends of 84, 95, and 104.3 GHz methanol integrated intensity with increasing FWZP values of the SiO (2−1) emission. This further supports the perspective that the integrated intensities of class I methanol emission increase with increasing shock speeds. The linear fitting results in log-log form are summarised in Table 6.

Fig. 5 Peak velocity of methanol emission with respect to systemic velocity (from Urquhart et al. 2022) for the 84 GHz (top), 95 GHz (middle), and 104.3 GHz (bottom) lines. The orange filled and purple open histograms represent methanol masers and all methanol detections (including masers and maser candidate), respectively. The black solid and dashed lines depict the Gaussian fitting results for masers and all methanol detection, respectively.

Fig. 6 Methanol integrated intensity as a function of SiO integrated intensity. The top panels (a, b, c) respectively show the integrated 84, 95, and 104.3 GHz methanol emission vs. SiO (2−1) integrated intensity. The bottom panels (d, e, f) respectively show the integrated 84, 95, and 104.3 GHz methanol emission vs. SiO (5−4) integrated intensity. The orange and purple triangles represent the clumps hosting masers and maser candidates, respectively. The orange solid and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

Fig. 7 95 GHz methanol integrated intensity vs. SiO (2−1) integrated intensity. The orange and purple dots represent the clumps hosting masers and maser candidates, respectively. The orange and purple circle are the sources with a single component. The orange and purple dashed lines depict the least-squares fitting results for methanol masers and maser candidates with a single component, respectively. The orange solid line represents the least-squares fitting results for methanol masers, as in Fig. 6b.

Fig. 8 Distribution of detection number and detection rate of three methanol transitions for SiO (2−1) emission with or without wings. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the filled and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

4.2 Physical properties on a clump scale

We made use of the latest catalogue of ATLASGAL physical clump properties (Urquhart et al. 2022) to study the relationship between the physical parameters of clumps and associated class I methanol (maser) emission detections. The clump parameters used for analysis include the bolometric luminosity, L_bol; the FWHM clump mass, M_fwhm; the luminosity-to-mass ratio, L_bol/M_fwhm; the dust temperature, T_dust; the mean H₂ FWHM volume density, n_fwhm(H₂); and the peak H₂ column density, N(H₂) (Urquhart et al. 2018). The FWHM clump mass is calculated using the integrated 870 µm flux density within the 50 per cent contour of the peak of the ATLASGAL continuum emission, and the mean H₂ FWHM volume density is calculated using the M_fwhm and the FWHM clump size, which is determined by the number of pixels within the 50 per cent contour (see details in Eqs. (3) and (5) in Urquhart et al. 2022). These FWHM parameters eliminate observational bias due to source evolution, as demonstrated by Billington et al. (2019).

Figure 14 presents the cumulative distributions of the clump physical properties for sources with and without 84 GHz methanol detection. The K-S tests reveal that except for the FWHM clump mass, the properties of clumps with and without 84 GHz methanol detection are significantly different. The situation is similar in our sample of 95 and 104.3 GHz sources; their cumulative distributions are shown in Figs. B.1 and B.2. Previous research towards 6.7 GHz class II methanol masers (Billington et al. 2019) and 95 GHz class I methanol masers (Ladeyschikov et al. 2020) also shows a very similar trend in that only the FWHM clump mass shows no difference between maser detection and non-maser detection or for the entire sample.

Table 7 provides a statistical summary of the physical properties for the maser subsample, the detection subsample (i.e. masers and maser candidates), and the non-detections for the 84, 95, and 104.3 GHz transitions. It is clear that the L_bol, L_bol/M_fwhm, T_dust, n_fwhm(H₂), and N(H₂) of clumps with methanol detection typically show higher values (i.e. the minimum, maximum, median and mean) than those without a detection for all three methanol transitions, as also suggested in Figs. 14, B.1, and B.2; instead, for the clumps with methanol maser and methanol detections the values of these properties are comparable. Figures B.3 and B.4 show the results of K-S tests of the cumulative distributions of the clump physical properties for sources with masers and maser candidates for 84 and 95 GHz transitions. The p-values of the K-S tests for each property are not smaller than 0.0013, indicating that the two sub- samples belong to the same distribution. We conclude that the warmer, brighter and denser clumps are more easily detectable in methanol and are more likely to host the class I methanol masers.

Figure 15 shows the cumulative distributions of the clump physical properties for sources with 84, 95, and 104.3 GHz methanol detection, which are similar to Fig. 14. The K-S tests reveal that except for the FWHM clump mass, the properties of clumps with 104.3 GHz methanol detection are significantly different from the clumps with 84 or 95 GHz (p-values ≪0.0013). In contrast, the properties of clumps with 84 and 95 GHz methanol detection are similar (p-values ≫0.0013), due to the fact that a large portion of sources have both maser species. This supports a significant overlap in the physical conditions that excite 84 and 95 GHz masers.

As shown in Fig. 15, the clumps with 104.3 GHz masers generally show brighter luminosities, warmer dust temperatures, higher luminosity-to-mass ratios, and denser environments than the clumps with 84 and 95 GHz masers (see also Table 7). The narrower coverage (i.e. maximum minus minimum in Table 7) of these parameters for the clumps hosting 104.3 GHz masers also indicates that this maser species occurs only under more demanding physical conditions. The minimum values of these parameters of the clumps hosting 104.3 GHz masers suggest that this maser species would not arise in environments (on clump scales) where L_bol ≲ 10⁴ L_⊙, M_fwhm ⊙ 200 M_⊙, T_dust ⊙ 22 K, N(H₂) ⊙ 10²³ cm⁻², and n_fwhm(H₂) ⊙ 10⁵ cm⁻³.

The isotropic luminosity of CH₃OH masers can be calculated via the following equations: $$L_{84}=8.81\times {10}^{-8}{L}_{\odot }{\left(\frac{d}{1\text{\hspace{0.17em}kpc}}\right)}^2\left(\frac{{\displaystyle \int S_{84}\text{d}V}}{1\text{\hspace{0.17em}Jy\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right),$$(1) $$L_{95}=9.92\times {10}^{-8}{L}_{\odot }{\left(\frac{d}{1\text{\hspace{0.17em}kpc}}\right)}^2\left(\frac{{\displaystyle \int S_{95}\text{d}V}}{1\text{\hspace{0.17em}Jy\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right),$$(2) $$\begin{array}{l}{L}_{104.3}=1.09\times {10}^{-7}{L}_{\odot }{\left(\frac{d}{1\text{\hspace{0.17em}kpc}}\right)}^2\left(\frac{{\displaystyle \int S_{104.3}\text{d}V}}{1\text{\hspace{0.17em}Jy\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right).\\ \end{array}$$(3)

Here d is the distance in units of kpc and ∫ S_vdV is the total integrated flux density in units of Jy km s⁻¹ at the corresponding frequency v. Figure 16 shows the correlations between 84 GHz methanol (maser) luminosity or integrated intensity and ATLASGAL properties. Similarly, the correlations for the 95 and 104.3 GHz lines are shown in Figs. C.1 and C.2. Their linear fitting results in log-log form for three methanol transitions are summarised in Table 8.

From these figures and Table 8 it can be seen that (1) for all three methanol transitions the class I maser luminosity shows significant positive correlations with the bolometric luminosity of the embedded protostellar objects and the FWHM clump mass, which is consistent with previous studies of the 95 GHz masers towards BGPS sources (Chen et al. 2011, 2012; Yang et al. 2020); (2) for the 84 and 95 GHz methanol transitions the total integrated intensity of class I methanol masers has a significant positive correlation with peak H₂ column density; (3) for all three methanol transitions the methanol luminosity of class I methanol masers has a very weak, or a statistically insignificant, correlation with luminosity-to-mass ratio, dust temperature or the mean H₂ FWHM volume density. The methanol luminosity independence on the mean H₂ FWHM volume density agrees with the results in Chen et al. (2012).

Fig. 9 Integrated intensity for 84 GHz (left) and 95 GHz (right) masers, with and without SiO (2−1) wings. In each box plot the horizontal coloured line represents the median value, the box represents the interquartile range between the 25th and 75th percentiles, and the vertical lines (the ‘whiskers’) show the ranges from the minimum value to the 25th percentile and from the 75th percentile to the maximum value.

Fig. 10 Cumulative distribution functions of the FWZP of SiO (2−1) emission for sources with and without methanol detections. From panel a to c, the distributions are shown for 84, 95, and 104.3 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S test are presented in the panel.

Fig. 11 Cumulative distribution functions of the FWZP of SiO (2−1) emission for masers and maser candidates. Panels a and b show the distributions for 84 and 95 GHz methanol masers (orange lines) and maser candidates (purple lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S test are presented in the panel.

Fig. 12 Distribution of detection number and detection rate of three methanol transitions under different FWZP of SiO (2−1) emission. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the solid and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

Fig. 13 Methanol integrated intensity as a function of the FWZP of SiO (2−1) for the 84, 95, and 104.3 GHz detection. The orange and purple triangles represent the clumps hosting masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least- squares fitting results for methanol masers and all methanol detections, respectively.

Fig. 14 Cumulative distribution functions of the physical properties for the clumps with and without 84 GHz methanol detections. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mas s ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 84 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

4.3 Class I methanol masers in different evolutionary stages

As massive proto-stars evolve, the surrounding material can be dramatically changed due to powerful stellar radiation and jet or outflow activity, and the changing physical and chemical properties may successively meet the excitation conditions for different maser species. Ellingsen (2007) suggested using the presence and absence of interstellar masers to infer an evolutionary timeline for high-mass star formation regions, and proposed a possible evolutionary sequence for common maser species, including class I and II methanol, water, and OH masers. Breen et al. (2010) then refined and quantified the maser-based evolutionary timeline (see their Fig. 6). Because the infrared colours are redder in sources solely associated with class I methanol masers than in sources associated with both class I and II masers, class I methanol masers are generally believed to trace an earlier stage than class II masers (Ellingsen 2006). However, some studies have revealed that in some cases class I methanol masers can also be associated with sources in more evolved stages (e.g. Voronkov et al. 2010, 2014; Gómez-Ruiz et al. 2016), for example with expanding HII regions (Voronkov et al. 2010, 2014). Here we use the properties of the ATLASGAL clumps to study the evolutionary stages in which class I methanol masers may arise.

ATLASGAL clumps have been classified into four evolutionary stages, namely quiescent, protostellar, young stellar objects (YSOs), and HII regions (König et al. 2017; Giannetti et al. 2017; Urquhart et al. 2018, 2022). In this work we adopt the latest evolutionary stage classification given by Urquhart et al. (2022) for each source. Table 2 lists the evolutionary stage of each target. Out of 408 ATLASGAL sources, a total of 309 (~76%) are assigned to one of the four evolutionary stages: 58 sources are starless (i.e. in the quiescent stage), 69 sources belong to a protostellar stage, 90 sources are classified as YSOs, and 92 sources host HII regions. The other 99 sources are classified as ’ambiguous’, ’complicated’, or ‘PDR+embedded source’, or remain undefined in Urquhart et al. (2022), due to the complexity of the submillimetre/infrared emission from their environments, which precludes a reliable determination of the current level of star formation in them. Consequently, these 99 sources are excluded from the following statistical analysis. Figure 17 shows the proportions of the clumps in the four evolutionary stages.

Figure 18 shows that the detection numbers and detection rates increase with evolution for all three of our studied class I methanol detections (including maser sources and candidates). In particular, we find that class I methanol masers are favoured to be detected at the HII stage on a clump scale. Previous VLA observations result in a detection rate (54%) towards UCHII regions (Gómez-Ruiz et al. 2016), supporting our findings on an even smaller scale. Different from the 84 GHz and 95 GHz masers that are seen in four evolutionary stages, methanol emission at 104.3 GHz never occurs before the YSO stage, and the rare 104.3 GHz maser emission only arises at the HII stage. This is expected because the 104.3 GHz line has an upper energy level of 158.6 K and is not easily excited in cold environments.

In addition, the bolometric luminosity-to-mass (L_bol/M_fwhm) ratio can be used as an evolutionary indicator in the star formation process (e.g. Molinari et al. 2016; Elia et al. 2021; Urquhart et al. 2022). High L_bol/M_fwhm ratios are related to evolved sources where the luminosity increases and the envelope mass is consumed, while low ratios are related to young sources. Our observing targets cover a wide range of L_bol/M_fwhm values from 0.07 to 1037 L_⊙/M_⊙ spanning more than four orders of magnitude.

Billington et al. (2020) presented a box plot of the central 95 per cent of L_bol/M_fwhm ratios for 6.7, 12.2 GHz class II methanol masers, OH masers, and H₂O masers (see their Fig. 14). Ladeyschikov et al. (2022) updated this analysis with larger samples for these masers. Figure 19 shows a similar box plot of the L_bol/M_fwhm ratio distributions for ATLASGAL clumps associated with class I masers, comparing them to the four maser species studies by Ladeyschikov et al. (2022). From the figure, we can see that class I masers at 84 and 95 GHz masers have a great degree of overlap in the L_bol/M_fwhm ratios with H₂O, OH, and class II methanol masers. These two maser species can trace a similar evolutionary stage as H₂O maser, and appear prior to 6.7 and 12.2 GHz methanol and OH masers. Despite the small number, the 104.3 GHz class I masers appear to trace a short and more evolved stage compared to the other class I maser species.

Fig. 15 Cumulative distribution functions of the physical properties for the clumps with 84, 95, and 104.3 GHz methanol detections. Panels (a) to (f) are similar to the corresponding panels in Fig. 14. The vertical dashed lines in the corresponding colours depict the median values of each sample. The p-values from the K-S tests between each two methanol emission are presented in each panel.

4.4 Physical conditions of class I masers

As indicated in Leurini et al. (2016), the intensity ratio of two maser transitions usually gives much stronger constraints on physical conditions than one single transition. If there are three coincident methanol maser transitions detected in one source, theoretical calculations of the maser intensity ratios would put stronger constraints on the physical conditions (such as kinetic temperature and H₂ volume density) where all masers are excited.

We use RADEX (van der Tak et al. 2007) and myRadex³ (Du 2022) to create model grids for the methanol lines. MyRadex solves the same problem as RADEX, except that a different approach is used to solve the statistical equilibrium problem (i.e. using an ordinary differential equation solver to evolve the system towards equilibrium under an initial distribution). In our work we show the modelling results from the myRadex code instead of RADEX, because RADEX calculations do not reach convergence when the specific methanol column densities (N/AV) are higher than 10¹⁴ cm⁻² km⁻¹ s, and such high specific methanol column densities are required to generate bright class I methanol masers (Leurini et al. 2016).

We make use of the molecular data of CH₃OH from the Leiden Atomic and Molecular Database (LAMDA⁴; Schöier et al. 2005), where the energy levels, transition frequencies and Einstein A coefficients are adopted from the CDMS (CDMS; Müller et al. 2005; Endres et al. 2016) and the latest collisional rates are obtained from Rabli & Flower (2010). The collisional rates for the torsional ground states of A- and E-type methanol are only known for collisions with para-H₂ in a temperature range from 10 to 200 K and including the rotational states up to J = 15 for collisions (Rabli & Flower 2010). We note that myRadex is able to account for a variable ortho-to-para ratio, which can be calculated under a local thermodynamic equilibrium assumption and only depends on local temperature (see e.g. Eq. (1) in Takahashi 2001). However, significant conversion from para- to ortho-H only starts at 700 K in C-shocks and at even higher temperatures in J-shocks (Wilgenbus et al. 2000). That means that the abundance of ortho-H₂ is very low at a kinetic temperature of <200 K. Thus, it is reasonable to consider that collisions with ortho-H₂ are negligible at kinetic temperatures of <200 K.

Methanol molecules can exist in one of two different symmetry species, A- and E- type CH₃OH. Since the very long timescale of proton exchange reactions to interconvert the two types of methanol, we can consider A-type and E-type methanol as two different molecules. The A/E methanol abundance ratio is related to the temperature where methanol was formed. The value of the ratio is about unity at a spin temperature (characterising the relative population of non-interacting spin types of a molecule) of 30−10 K and increasing at lower temperatures (see Wirström et al. 2011). Based on previous studies (Sobolev et al. 2005; Voronkov et al. 2006; Leurini & Menten 2018), class I methanol masers are usually excited at kinetic temperatures higher than 40 K. Therefore, it is reasonable to assume a methanol A/E ratio of unity in our calculations.

Figure 20 shows myRadex calculations for 84, 95, and 104.3 GHz methanol lines over an equally spaced grid with 39 temperatures from 10 to 200 K and 41 para-H₂ densities from 10⁴ to 10⁸ cm⁻³. The specific column densities of methanol are fixed at 10¹⁵ cm⁻² km⁻¹ s (Fig. 20a) and 10¹⁶ cm⁻² km⁻¹ s (Fig. 20b). A plane-parallel slab geometry is adopted, which is appropriate for shocks, and the corresponding escape probability β = (1−e^−3τ)/3τ, where τ is the optical depth. No external radiation field is taken into account except for the cosmic microwave background, because class I methanol masers are usually offset from infrared sources and UCHII regions, and interferometric observations reveal that they reside in the interface regions between outflows and ambient dense material (Plambeck & Menten 1990; Voronkov et al. 2006). The black and white contours represent the peak intensity ratios of T_r,104.3/T_r,84 and T_r,95/T_r,84, respectively. We assume that the maser components in these three transitions with the same velocity are co-spatial, and have the same small masing area (i.e. emission size ≪ beam size). In the calculations of peak intensity ratios, different beam dilution factors of each transition have been taken into account. For example, T_r,104.3/T_r,84 is equal to (T_mb,104.3/T_mb,84) · $\left({\theta }_{\text{beam,104}.3}^2/{\theta }_{\text{beam,84}}^2\right)$, where T_mb and θ_beam are the observed brightness temperature and the beam size at the corresponding frequency.

From our four detected 104.3 GHz masers, only G10.34−0.14 shows distinguishable and velocity-aligned (at ~14.5 km s⁻¹) maser features in the three transitions. We adopt this feature as an example to determine under what physical conditions the three masers can be excited simultaneously. The red star in Fig. 20 is the intersection of the observed two pairs of ratios, and the red dashed lines represent the uncertainties of the observed ratios. Under an assumed specific column density of methanol of 10¹⁵ cm⁻² km⁻¹ s, a kinetic temperature of 62±3 K and a para- H₂ number density of 4.4(±1.5)×10⁶ cm⁻³ are needed to excited the three class I methonal masers. When the assumed specific column density of methanol increases to 10¹⁶ cm⁻² km⁻¹ s, the masering region in G10.34−0.14 should have a kinetic temperature of 57 ±3 K and a para-H₂ number density of 7.9(±2.5)×10⁵ cm⁻³. Comparing the two panels in Fig. 20, it can be seen that (1) the 104.3 GHz methanol masers are usually produced in a dense (≳10⁵ cm⁻³) environment and (2) with increasing specific column densities of methanol the absolute values of optical depth become higher, which means stronger maser emission, and the strongest maser occurs in a less dense environment.

Fig. 16 Distribution of 84 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. From panel (a) to (e) the methanol isotropic luminosity is plotted against the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, and the mean H2 FWHM volume density. In panel (f) the total integrated intensity of the 84 GHz methanol transition is plotted against the peak H2 column density. The orange and purple triangles represent the ATLASGAL clumps that host masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

Fig. 17 Proportions of 309 ATLASGAL clumps classified into the four evolutionary stages.

Fig. 18 Distribution of methanol detection number and detection rate of three methanol transitions under four evolutionary stages. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the solid and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

Fig. 19 Box plot of distributions of bolometric luminosity-to-mass ratios for clumps associated with different masers. Except for class I maser data, the plotting data for H2O, OH, and class II methanol masers at 6.7 and 12.2 GHz are from Fig. 7(A) in Ladeyschikov et al. (2022). The boxes filled in red, blue, light yellow, and green indicate the class I and II methanol masers, H2O, and OH masers, respectively. Each box represents the interquartile range between the 25th and 75th percentiles, with a thick black line denoting the median value. The whiskers show the full range of the data.

Fig. 20 Statistical equilibrium calculations of optical depths and line ratios. The left and right panels show that the calculations are performed under the specific column densities of 1015 and 1016 cm−2 km−1 s, respectively. The coloured background shows the regions with negative optical depths of the 104.3 GHz methanol transition where masers can be generated. The black and white contours represent the peak intensity ratios of Tr,104.3/Tr,84 and Tr,95/Tr,84, respectively. The red star is the intersection of the observed two pairs of ratios in G10.34−0.14, and the red dashed lines represent the 1σ uncertainties of the observed ratios.

5 Summary

We performed a survey of three class I methanol masers in the 84, 95, and 104.3 GHz transitions using the IRAM 30-meter telescope towards 408 ATLASGAL clumps. The main results are summarised as follows:

1. In 282 (70%), 224 (55%), and 29 (7%) sources we detect methanol emission at 84, 95, and 104.3 GHz, respectively. The 104.3 GHz emission is only found in sources with both 84 and 95 GHz detections. A total of 54, 100, and 4 sources show maser-like features at 84, 95, and 104.3 GHz, respectively. Among them, 50 masers at 84 GHz, 29 masers at 95 GHz, and 4 masers at 104.3 GHz are new discoveries. Our work increases the number of known 104.3 GHz masers from five to nine.

2. In our sample the 95 GHz class I methanol maser is generally stronger than its 84 GHz maser counterpart. For both the 84 and 95 GHz class I masers, the relative velocities between maser velocity and systemic velocity are within ~5 km s⁻¹, confirming that the class I maser velocity can trace the systemic velocity well.

3. All detected 104.3 GHz methanol masers are redshifted with respect to the systemic velocity. Based on the properties of their associated ATLASGAL clumps, this maser species would not arise in environments where L_bol < 10⁴ L_⊙, M_fwhm ≲ 200 M_⊙, T_dust ≲ 22 K, N(H₂) ≲ 10²³ cm⁻², and n_fwhm(H₂) ≲ 10⁵ cm⁻³.

4. We find nine sources that show class I methanol masers but no SiO emission, indicating that class I methanol masers might be the only signpost of protostellar activity in extremely embedded objects at the earliest evolutionary stage. More and stronger class I masers were detected towards sources showing SiO line wings than towards sources without SiO wings. The total integrated intensity of class I masers is positively correlated with SiO integrated intensity and FWZP of SiO (2−1) emission. These facts strongly suggest that the properties of class I masers are regulated by shock properties also traced by SiO.

5. The properties of class I methanol masers show positive correlations with the following properties of associated ATLASGAL clumps: bolometric luminosity, clump mass, and peak H₂ column density. There is no statistically significant correlation between the luminosity of class I methanol masers and the luminosity-to-mass ratio, dust temperature, or mean H₂ volume density.

6. Our results show that the 84 GHz and 95 GHz methanol masers exist in the quiescent, protostellar, YSO, and HII region stages. In contrast, the 104.3 GHz methanol masers are only detected in the HII region stage. Based on the distribution of bolometric luminosity-to-mass ratios associated with different masers, we suggest that class I methanol masers at 84 and 95 GHz can trace an evolutionary stage similar to that of H₂O masers, and appear prior to 6.7 and 12.2 GHz methanol and OH masers. The 104.3 GHz class I masers appear to trace a shorter and more evolved stage compared to the other class I maser species.

7. Our calculations show that the physical conditions can be better constrained in a region where multiple class I methanol masers arise. With an assumed specific column density of methanol of 10¹⁵ cm⁻² km⁻¹ s (10¹⁶ cm⁻² km⁻¹ s), a kinetic temperature of 62±3 K (57±3 K), and a para-H₂ number density of 4.4(±1.5)×10⁶ cm⁻³ (7.9(±2.5)×10⁵ cm⁻³) are needed to simultaneously excite the three class I methanol masers in G10.34−0.14.

We present a systemic study of class I methanol masers in the inner Galaxy and their relation with SiO emission on clump scales. In future work, more in-depth statistical analysis such as mathematical classification models can be performed to understand which properties of the ATLASGAL clump are the best predictors of an associated class I methanol maser. This analysis can then be used to predict which ATLASGAL clumps that have not been observed should be targeted for class I maser searches. In addition, higher angular follow-up observations will be able to pinpoint the location of class I methanol masers on smaller scales, which will provide more constraints on the relation between the masers and their driving stars.

Acknowledgements

T. Cs. has received financial support from the French State in the framework of the IdEx Université de Bordeaux Investments for the future Program. This work made use of Python libraries including Astropy (https://www.astropy.org/) (Astropy Collaboration 2013), NumPy (https://www.numpy.org/) (van der Walt et al. 2011), SciPy (https://www.scipy.org/) (Jones et al. 2001), Matplotlib (https://matplotlib.org/) (Hunter 2007). This research has made use of the VizieR catalogue, operated at CDS, Strasbourg, France.

Appendix A Spectra of nine methanol maser transitions for 408 observed ATLASGAL sources

Fig. A.1 Observed spectra of nine methanol maser transitions at 3 mm towards the first six ATLASGAL sources in our sample. The frequency and the quantum numbers of class I and II CH3OH maser transitions are labelled in red and blue, respectively, and their Gaussian fitting results are plotted in corresponding colours. The horizontal magenta dashed lines represent the 3σ noise level for the non-detection transitions. For marginal detections with S/N values of ~3, both Gaussian fitting results and 3σ levels are presented. The spectra of the whole sample are available in electronic form at the Zenodo via: https://doi.org/10.5281/zenodo.7442831.

Appendix B K-S test results of ATLASGAL properties

Fig. B.1 Cumulative distribution functions of the physical properties for the clumps with and without 95 GHz methanol detections. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 95 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

Fig. B.2 Cumulative distribution functions of physical properties for clumps with and without 104.3 GHz methanol detections. Panels (a) to (f) are similar to the corresponding panels in Fig. B.1.

Fig. B.3 Cumulative distribution functions of the physical properties for the clumps with 84 GHz methanol masers and maser candidates. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mas s ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 84 GHz methanol masers (orange lines) and maser candidates (purple lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

Fig. B.4 Cumulative distribution functions of the physical properties for the clumps with 95 GHz methanol masers and maser candidates. Panels (a) to (f) are similar to the corresponding panels in Fig. B.3.

Appendix C Distribution of 95 and 104.3 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps

Fig. C.1 Distribution of 95 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. In panels (a) to (e) total methanol isotropic luminosity is plotted against the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, and the mean H2 FWHM volume density. In panel (f) the total integrated intensity of the 95 GHz methanol transition is plotted against the peak H2 column density. The orange and purple triangles represent the ATLASGAL clumps that host masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

Fig. C.2 Distribution of 104.3 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. Panels (a) to (f) are similar to the corresponding panels in Fig. C.1.

References

All Tables

All Figures

	Fig. 1 Sources detected in the three class I maser lines. The Venn diagram shows that all detected 104.3 GHz emission lines in our sample have associated methanol emission at both 84 and 95 GHz.
In the text

Fig. 2 Examples of diverse line profiles of observed class I methanol emission at 84 GHz (upper panel) and 95 GHz (lower panel), depicting narrow maser-like features (panel a), broad thermal-like features (panel b) and a combination of the two (panel c). The line profiles of the two transitions are sometimes quite similar (panels a to c) and sometimes very different (panels d to f). The frequency and the corresponding quantum numbers of each transition are shown in the upper right corner of each panel for six sources. The horizontal magenta dashed lines in G36.90−0.41 and G37.48−0.10 represent the 3σ noise level for the undetected transitions.

In the text

	Fig. 3 Peak flux density ratio of 95 and 84 GHz emission. The grey shaded histogram represents the whole sample with both 95 and 84 GHz detections, the pink line shows the sources hosting both 84 and 95 GHz masers, and the black line shows the sources hosting maser candidates in both transitions. The red and blue dashed lines depict the median and mean values of the whole sample.
In the text

	Fig. 4 Spectra of the 84, 95, 104.3 GHz lines (from top to bottom) for 4 sources showing 104.3 GHz methanol maser features. The vertical magenta dashed lines represent the systemic velocity adopted from Urquhart et al. (2022). The frequency and the corresponding quantum numbers of each transition are shown in each panel.
In the text

Fig. 5 Peak velocity of methanol emission with respect to systemic velocity (from Urquhart et al. 2022) for the 84 GHz (top), 95 GHz (middle), and 104.3 GHz (bottom) lines. The orange filled and purple open histograms represent methanol masers and all methanol detections (including masers and maser candidate), respectively. The black solid and dashed lines depict the Gaussian fitting results for masers and all methanol detection, respectively.

In the text

Fig. 6 Methanol integrated intensity as a function of SiO integrated intensity. The top panels (a, b, c) respectively show the integrated 84, 95, and 104.3 GHz methanol emission vs. SiO (2−1) integrated intensity. The bottom panels (d, e, f) respectively show the integrated 84, 95, and 104.3 GHz methanol emission vs. SiO (5−4) integrated intensity. The orange and purple triangles represent the clumps hosting masers and maser candidates, respectively. The orange solid and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

In the text

Fig. 7 95 GHz methanol integrated intensity vs. SiO (2−1) integrated intensity. The orange and purple dots represent the clumps hosting masers and maser candidates, respectively. The orange and purple circle are the sources with a single component. The orange and purple dashed lines depict the least-squares fitting results for methanol masers and maser candidates with a single component, respectively. The orange solid line represents the least-squares fitting results for methanol masers, as in Fig. 6b.

In the text

Fig. 8 Distribution of detection number and detection rate of three methanol transitions for SiO (2−1) emission with or without wings. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the filled and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

In the text

Fig. 9 Integrated intensity for 84 GHz (left) and 95 GHz (right) masers, with and without SiO (2−1) wings. In each box plot the horizontal coloured line represents the median value, the box represents the interquartile range between the 25th and 75th percentiles, and the vertical lines (the ‘whiskers’) show the ranges from the minimum value to the 25th percentile and from the 75th percentile to the maximum value.

In the text

Fig. 10 Cumulative distribution functions of the FWZP of SiO (2−1) emission for sources with and without methanol detections. From panel a to c, the distributions are shown for 84, 95, and 104.3 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S test are presented in the panel.

In the text

Fig. 11 Cumulative distribution functions of the FWZP of SiO (2−1) emission for masers and maser candidates. Panels a and b show the distributions for 84 and 95 GHz methanol masers (orange lines) and maser candidates (purple lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S test are presented in the panel.

In the text

Fig. 12 Distribution of detection number and detection rate of three methanol transitions under different FWZP of SiO (2−1) emission. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the solid and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

In the text

	Fig. 13 Methanol integrated intensity as a function of the FWZP of SiO (2−1) for the 84, 95, and 104.3 GHz detection. The orange and purple triangles represent the clumps hosting masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least- squares fitting results for methanol masers and all methanol detections, respectively.
In the text

Fig. 14 Cumulative distribution functions of the physical properties for the clumps with and without 84 GHz methanol detections. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mas s ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 84 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

In the text

	Fig. 15 Cumulative distribution functions of the physical properties for the clumps with 84, 95, and 104.3 GHz methanol detections. Panels (a) to (f) are similar to the corresponding panels in Fig. 14. The vertical dashed lines in the corresponding colours depict the median values of each sample. The p-values from the K-S tests between each two methanol emission are presented in each panel.
In the text

Fig. 16 Distribution of 84 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. From panel (a) to (e) the methanol isotropic luminosity is plotted against the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, and the mean H2 FWHM volume density. In panel (f) the total integrated intensity of the 84 GHz methanol transition is plotted against the peak H2 column density. The orange and purple triangles represent the ATLASGAL clumps that host masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

In the text

	Fig. 17 Proportions of 309 ATLASGAL clumps classified into the four evolutionary stages.
In the text

Fig. 18 Distribution of methanol detection number and detection rate of three methanol transitions under four evolutionary stages. The labels for the detection numbers are on the left side, while those for the detection rates are on the right side. The red, blue, and green bins present the detection numbers of methanol emission at 84, 95, and 104.3 GHz, respectively, where the solid and hatched parts represent masers and maser candidates. The red, blue, and green triangles indicate the detection rates of methanol emission at 84, 95, and 104.3 GHz.

In the text

Fig. 19 Box plot of distributions of bolometric luminosity-to-mass ratios for clumps associated with different masers. Except for class I maser data, the plotting data for H2O, OH, and class II methanol masers at 6.7 and 12.2 GHz are from Fig. 7(A) in Ladeyschikov et al. (2022). The boxes filled in red, blue, light yellow, and green indicate the class I and II methanol masers, H2O, and OH masers, respectively. Each box represents the interquartile range between the 25th and 75th percentiles, with a thick black line denoting the median value. The whiskers show the full range of the data.

In the text

Fig. 20 Statistical equilibrium calculations of optical depths and line ratios. The left and right panels show that the calculations are performed under the specific column densities of 1015 and 1016 cm−2 km−1 s, respectively. The coloured background shows the regions with negative optical depths of the 104.3 GHz methanol transition where masers can be generated. The black and white contours represent the peak intensity ratios of Tr,104.3/Tr,84 and Tr,95/Tr,84, respectively. The red star is the intersection of the observed two pairs of ratios in G10.34−0.14, and the red dashed lines represent the 1σ uncertainties of the observed ratios.

In the text

Fig. A.1 Observed spectra of nine methanol maser transitions at 3 mm towards the first six ATLASGAL sources in our sample. The frequency and the quantum numbers of class I and II CH3OH maser transitions are labelled in red and blue, respectively, and their Gaussian fitting results are plotted in corresponding colours. The horizontal magenta dashed lines represent the 3σ noise level for the non-detection transitions. For marginal detections with S/N values of ~3, both Gaussian fitting results and 3σ levels are presented. The spectra of the whole sample are available in electronic form at the Zenodo via: https://doi.org/10.5281/zenodo.7442831.

In the text

Fig. B.1 Cumulative distribution functions of the physical properties for the clumps with and without 95 GHz methanol detections. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 95 GHz methanol detections (red lines) and non-detections (blue lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

In the text

	Fig. B.2 Cumulative distribution functions of physical properties for clumps with and without 104.3 GHz methanol detections. Panels (a) to (f) are similar to the corresponding panels in Fig. B.1.
In the text

Fig. B.3 Cumulative distribution functions of the physical properties for the clumps with 84 GHz methanol masers and maser candidates. In panels (a) to (f) are plotted the cumulative distributions of the bolometric luminosity, the FWHM clump mass, the luminosity-to-mas s ratio, the dust temperature, the peak H2 column density, and the mean H2 FWHM volume density for 84 GHz methanol masers (orange lines) and maser candidates (purple lines), respectively. The vertical dashed lines in the corresponding colours depict the median values of the two samples. The p-values from the K-S tests are presented in each panel.

In the text

	Fig. B.4 Cumulative distribution functions of the physical properties for the clumps with 95 GHz methanol masers and maser candidates. Panels (a) to (f) are similar to the corresponding panels in Fig. B.3.
In the text

Fig. C.1 Distribution of 95 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. In panels (a) to (e) total methanol isotropic luminosity is plotted against the bolometric luminosity, the FWHM clump mass, the luminosity-to-mass ratio, the dust temperature, and the mean H2 FWHM volume density. In panel (f) the total integrated intensity of the 95 GHz methanol transition is plotted against the peak H2 column density. The orange and purple triangles represent the ATLASGAL clumps that host masers and maser candidates, respectively. The orange solid lines and black dashed lines depict the least-squares fitting results for methanol masers and all methanol detections, respectively.

In the text

	Fig. C.2 Distribution of 104.3 GHz methanol luminosity or integrated intensity against properties of ATLASGAL clumps. Panels (a) to (f) are similar to the corresponding panels in Fig. C.1.
In the text