Methanol maser associated outflows: detection statistics and properties

Abstract

We have selected the positions of 54 6.7 GHz methanol masers from the Methanol Multibeam Survey catalogue, covering a range of longitudes between 20° and 34° of the Galactic plane. These positions were mapped in the J = 3-2 transition of both the ¹³CO and C¹⁸O lines. A total of 58 ¹³CO emission peaks are found in the vicinity of these maser positions. We search for outflows around all ¹³CO peaks, and find evidence for high-velocity gas in all cases, spatially resolving the red and blue outflow lobes in 55 cases. Of these sources, 44 have resolved kinematic distances, and are closely associated with the 6.7 GHz masers, a subset referred to as Methanol Maser Associated Outflows (MMAOs). We calculate the masses of the clumps associated with each peak using 870 μm continuum emission from the ATLASGAL survey. A strong correlation is seen between the clump mass and both outflow mass and mechanical force, lending support to models in which accretion is strongly linked to outflow. We find that the scaling law between outflow activity and clump masses observed for low-mass objects, is also followed by the MMAOs in this study, indicating a commonality in the formation processes of low-mass and high-mass stars.

INTRODUCTION

Massive stars (>8 M_⊙) play a key role in the evolution of the Universe, as the principal sources of heavy elements and UV radiation. Their winds, massive outflows, expanding H ii regions and supernova explosions serve as an important source of enrichment, mixing and turbulence in the interstellar medium (ISM) of galaxies (Zinnecker & Yorke 2007). Our understanding of the formation and evolution of young massive stars is made difficult by their rarity, large average distances that demands observations at higher angular resolution, deeply embedded formation within dense clusters resulting in confusing dynamics and obscuration, and rapid evolution with short-lived evolutionary phases (Shepherd & Churchwell 1996b; Zinnecker & Yorke 2007).

The specific formation process of massive stars is not yet fully understood. These stars reach the zero-age main sequence while still accreting material from their parent molecular cloud. Due to their high mass, they radiate strongly. This radiation pressure exceeds the gravitational pressure, and should the formation process be similar to low-mass stars, the growing radiation pressure from the newborn stars will eventually become strong enough to stop the accretion, yielding an upper mass limit of ∼40 M_⊙ (Wolfire & Cassinelli 1987; Stahler & Palla 1993).

Previously, two solutions were proposed to overcome this problem: (i) a formation process involving multiple lower mass stars, either via coalescence of low- to intermediate-mass protostars (e.g. Bonnell, Bate & Zinnecker 1998; Bally & Zinnecker 2005), or competitive accretion in a clustered environment (e.g. Bonnell, Vine & Bate 2004), or (ii) a scaled-up version of the process found in low-mass star formation. The latter solution can be subdivided into the following main categories: (a) increased spherical accretion rates in turbulent cloud cores (order 10⁻⁴–10⁻³ M_⊙ yr⁻¹), high enough to overcome the star's radiation pressure (e.g. Norberg & Maeder 2000; McKee & Tan 2003) or (b) accretion via discs on to a single massive star (e.g. Jijina & Adams 1996; Yorke & Sonnhalter 2002).

A solution to overcome the radiation pressure barrier was proposed by Yorke & Sonnhalter (2002), that involved the generation of a strong anisotropic radiation field where an accretion disc reduces the effects of radiative pressure, by allowing photons to escape along the polar axis (the ‘flashlight effect’). However, these simulations showed an early end of the disc accretion phase, with final masses limited to ∼42 M_⊙. Krumholz et al. (2009) suggested that the early end of the accretion phase is because the disc loses its shielding property as it cannot be fed in an axially symmetric configuration. Contrary to the stable radiation pressure-driven outflows in Yorke & Sonnhalter (2002), they proposed a three-dimensional radiation hydrodynamic simulation with a Rayleigh–Taylor instability in the outflow region, allowing further accretion on to the disc.

Kuiper et al. (2010) took this further by introducing a dust sublimation front to their simulations. This preserves the shielding of the massive accretion disc and allows the protostar to grow to ∼140 M_⊙.

The easiest way to verify the disc accretion models, would be with the detection of accretion discs around massive protostars, but this is difficult without specialized techniques (e.g. Pestalozzi, Elitzur & Conway 2009), because they are small (at most several hundred au), short lived, and easily confused by envelopes (Kim & Kurtz 2006). Few clear examples of such discs exist (e.g. Cesaroni et al. 2007; Zapata, Tang & Leurini 2010).

However, we expect that if massive stars do form via accretion discs, they will generate massive and powerful outflows, similar to low-mass stars. These outflows are necessary to transport angular momentum away from a forming star (Shu et al. 1991, 2000; Konigl & Pudritz 2000; Chrysostomou et al. 2008). For massive stars, these outflows should be of much larger scale and easier to detect than the accretion discs (Kim & Kurtz 2006). Studying outflows offers an alternative approach to probe the embedded core.

There have been many studies that collectively suggest outflows are ubiquitously associated with massive star formation (e.g. Shepherd & Churchwell 1996a; Molinari et al. 1998; Beuther et al. 2002b; Xu et al. 2006).

Zhang et al. (2005) found outflow masses (∼tens to hundreds  M_⊙), momenta (10–100  M_⊙ km s⁻¹) and energies (∼10³⁹ J) towards their sample of luminous IRAS point sources about a factor 10 higher than the values of low-mass outflows (Bontemps et al. 1996). This suggests that outflows consist of accelerated gas that has been driven by a young accreting protostar, rather than swept-up ambient material (Churchwell 1999). It could also be material that originates from the accretion disc/young stellar object (YSO) and is funnelled out of the central system (e.g. Shepherd & Churchwell 1996a).

To date, CO observations of molecular outflows have been made using mainly two methods: (1) single-point CO line surveys towards samples of massive YSOs in search of high-velocity (HV) molecular gas (e.g. Shepherd & Churchwell 1996b; Sridharan et al. 2002) or (2) CO line mapping of carefully selected sources that exhibit HV wings (e.g. Shepherd & Churchwell 1996a; Beuther et al. 2002b). Unless outflows are mapped, it is difficult to determine their physical properties. Mapping outflows at sufficient sensitivity and high angular resolution is time-consuming, but the development of heterodyne focal plane arrays [e.g. HARP on James Clerk Maxwell Telescope (JCMT) or HERA on Institut de Radioastronomie Millimétrique (IRAM)] has made it possible to map statistically significant samples of massive star-forming regions to search for outflows l(e.g. López-Sepulcre et al. 2009; Gottschalk et al. 2012).

Outflows are one of the earliest observable signatures of star formation, and are believed to develop from the central objects during the infrared bright stage called the ‘hot core’ phase (Cesaroni, Walmsley & Churchwell 1992; Kurtz et al. 2000), just before the UCHii phase (Shepherd & Churchwell 1996a; Wu et al. 1999; Zhang et al. 2001; Beuther et al. 2002b; Molinari et al. 2002).

Another important signpost of the ‘hot core’ phase is the turn-on of radiatively pumped 6.7 GHz (class II) methanol masers, the second brightest masers in the Galaxy (Menten 1991; Sobolev, Cragg & Godfrey 1997; Minier et al. 2003). Observations indicate that these masers are rarely associated with H ii regions, but most of them are found to be associated with massive millimetre and submillimetre sources (e.g. Beuther et al. 2002b; Walsh et al. 2003; Urquhart et al. 2013a). It appears as if these masers occupy a brief phase in the pre-UCHii region, even as short as ∼10⁴ yr, and disappear as the UCHii region evolves (Hatchell et al. 1998; Codella & Moscadelli 2000; Codella et al. 2004; van der Walt 2005; Wu et al. 2010). They are also known to be mostly associated with massive star formation, making them important signposts of massive star formation (Minier et al. 2005; Ellingsen 2006; Breen et al. 2013; Caswell 2013).

However, there are limited simultaneous studies of methanol masers and outflow activity. Minier, Conway & Booth (2001) found that 10 out of 13 absolute positions for class II methanol maser sites coincided with typical tracers of massive star formation (e.g. UCHii regions, outflows and hot cores), while seven out of these ten were within less than 2000 au (∼10⁻² pc) from outflows. Their results supported the expected association between the occurrence of class II methanol masers and molecular outflows.

The Spitzer GLIMPSE survey (Churchwell et al. 2009) revealed a new signpost for outflows in high-mass star formation regions in the form of extended emission which is bright in the 4.5-μm band. These objects are generally referred to either as extended green objects (EGOs; Cyganowski et al. 2008) or green fuzzies (Chambers et al. 2009). The enhanced emission in this wavelength range is believed to be due to shock-excited H₂ and/or CO band-head emission (De Buizer & Vacca 2010). Cyganowski et al. (2008) found that many EGOs are associated with 6.7 GHz methanol masers, while Chen, Ellingsen & Shen (2009) showed a high rate of association with shock-excited class I methanol masers at 44 and 95 GHz. Sensitive, high-resolution searches for class II methanol masers towards a small sample of EGOs achieved a detection rates of 64 per cent (although this should be considered an upper limit since most targets had known 6.7 GHz methanol masers in their vicinity), with approximately 90 per cent of these sources also having associated 44 GHz class I methanol maser emission (Cyganowski et al. 2009). These results demonstrate a close association between methanol masers and young high-mass stars with active outflows.

Molecular outflows are more visible than the YSO or its disc, and because of the association of 6.7 GHz methanol masers with massive star formation, searching for outflows towards these masers and studying their physical properties can reveal information regarding the obscured massive cores they are associated with. Moreover, by selecting outflows that are only associated with methanol masers, deliberately biases the resulting sample towards a narrower, relatively well-defined evolutionary range which allows constraints to be placed on the ‘switch-on’ of the outflows and the study of their temporal development. In this paper, we focus on the study of the physical properties of the outflows and the relationship of these properties with those of their embedding clumps. In a following publication (de Villiers et al., in preparation), we will explore the effects of the maser selection bias in our sample and the resulting behaviour in the dynamical ages of our maser selected sample.

We present a survey of ¹³CO(J = 3-2) outflows towards a sample of 6.7 GHz Methanol Multibeam (MMB) masers (Green et al. 2009; Breen et al. in preparation) using the HARP (Heterodyne Array Receiver Programme) instrument on the JCMT. Observations and data reduction are described in Section 2. In Section 3, we describe the extraction and analysis of the spectra, as well as outflow mapping and outflow detection frequency. The results are presented in Section 4, where we demonstrate the calculation of the outflows’ physical properties and associated clump masses. The relation between the outflow and associated clump properties are examined, and compared with some low-mass relations found in the literature. We also inspect the correlation between outflow and 6.7 GHz maser luminosities, as well as between maser luminosity and clump masses, as a probe of the relationship between the physical properties of the driving force, outflow and associated maser. The main results are summarized in Section 6.

Although the study of the properties of massive molecular outflows and their relation with associated clump masses is not novel per se, the selection of the sources in this study is unique in terms of association with 6.7 GHz masers. This allows the selection of sources within a relatively well-defined evolutionary phase, which potentially could limit the scatter in parameter space compared to previous work. In this paper, we discuss and investigate the physical properties of the Methanol Maser Associated Outflows (MMAOs), and put them in context with other studies. In a second forthcoming paper, we discuss the effect and implications of the 6.7 GHz maser bias of our sample on our results.

OBSERVATIONS AND DATA REDUCTION

A sample of 6.7 GHz methanol masers were drawn from a preliminary catalogue of Northern hemisphere masers from the MMB Survey which has subarcsec positional accuracies (Green et al. 2009). The properties of these masers are described fully in Breen et al. (in preparation). The initial sample selection was chosen to have an even spread in maser luminosity, distance, association with UCHii regions and IR sources. A sample of 70 sources were observed between 20° < l < 34°.

The targets were observed with the JCMT, on the summit of Mauna Kea, Hawaii on seven nights between 2007 May 17 and 2008 July 22. Targets were mapped in the ¹³CO and C¹⁸O(J = 3-2) transitions (330.6 and 329.3 GHz), using the 16-receptor HARP. Only 14 of the 16 receptors were operational at the time of observation. The receptors are laid out in a 4 × 4 grid separated by 30 arcsec and the beam size of the individual receptors at 345 GHz is 14 arcsec. All the data were corrected for a main-beam efficiency of η_mb = 0.66 (Smith et al. 2008; Buckle et al. 2009). A HARP jiggle map (Buckle et al. 2009) produces a fully sampled, 16-point rectangular map with a pixel scale of 6 arcsec and a spectral resolution of 0.06 km s⁻¹. The field of view is approximately 2 arcmin × 2 arcmin. As the typical distance to the methanol maser target sources is >2kpc, and with an estimated maser lifetime of 2.5-4 × 10⁴ yr (van der Walt 2005), the expected outflows from the maser-associated YSOs should be sampled in a single JCMT HARP jiggle-map. The pointing accuracy of the JCMT is of order 2 arcsec or better. Pointing checks were carried out regularly during observation runs to ensure and maintain accuracy.

The weather during the observations was mostly in JCMT-defined band 3, which implies a sky zenith opacity τ₂₂₅ varying between 0.08 and 0.12 at 225 GHz as measured by the Caltech Submillimeter Observatory tipping radiometer.¹

Out of the 70 observed maser coordinates, 16 observations did not meet one or more of the quality thresholds due to (a) too low signal-to-noise (less than ∼2), (b) non-functioning receptors (report unreliable temperatures), or (c) target positioning too close to the field-of-view border or a dead receptor. The remaining 54 target coordinates are listed in Table 1 and occur in the shaded area in Fig. 1.

The ¹³CO and C¹⁸O maps were simultaneously obtained using the multiple subband mode of the back-end Auto-Correlation Spectral Imaging System (ACSIS; Dent et al. 2000). The raw ACSIS data are in a HARP time series cube, giving the response of the receptors (x-axis) as a function of time (y-axis). The third dimension is the velocity spectrum recorded at the time for that receptor. Data were reduced with the Starlink ORAC-DR pipeline (Cavanagh et al. 2008) using the |${\rm REDUCE\_SCIENCE\_NARROWLINE}$| recipe with minor modifications tailored for this data set.² The pipeline reduction process automatically fits and subtracts polynomial baselines. This was followed by truncation of the noisy spectral endpoints, removal of interference spikes and rebinning of the spectrum to a resolution of 0.5 km s⁻¹. Any receptors with high baseline variations compared to the bulk of the spectra, were flagged as bad in addition to those masked out by the pipeline. Lastly, the time series were then mapped on to a position–position–velocity cube. The reduced data antenna temperature (T_A) had an average rms noise level of 0.24 K (per 6 arcsec × 6 arcsec × 0.5 km s^{− 1} pixel), or a main-beam efficiency corrected average rms noise level of T_mb = 0.36 K.

DATA ANALYSIS

Finding the peak emission

¹³CO was used as an outflow tracer in this study. It is a useful probe of the cloud structure and kinematics, because it traces the higher velocity gas, but has a lower abundance than ¹²CO, and hence is less contaminated by other HV structures within the star-forming complex (Arce et al. 2010). Emission from the (J = 3-2) transition was observed (T_trans = 31.8 K; Curtis, Richer & Buckle 2010), which traces the warm, dense gas, close to the embedded YSO and also serves as a clearer tracer of warm outflow emission than lower J transitions. Targets were simultaneously observed in the optically thin C¹⁸O transition, which serves as a useful tracer of the column density. The C¹⁸O emission peak is most likely to coincide with the YSO core's position.

Visual inspection indicated that the positions of peak emission in both ¹³CO and C¹⁸O did not always coincide with the maser coordinate. These offsets were larger than a beam size (14 arcsec) for seven maser coordinates, with a maximum offset of 1 arcsec. Although it is known that 6.7 GHz methanol masers are mostly associated with massive YSOs, two competing and unresolved formation hypotheses state that either methanol masers are embedded in circumstellar tori or accretion discs around the massive protostars (Pestalozzi et al. 2009), or that they generally trace outflows (De Buizer, Bartkiewicz & Szymczak 2012). It thus seems possible that although the masers are in the close vicinity of the YSO, some could be offset, as was found in this study.

Since the peak ¹³CO emission did not always coincide with the maser coordinates, and the exact coordinate of the peak emission was needed as the position from where the one-dimensional spectrum would be extracted as part of the outflow detection method, an alternative method was needed to pin-point this position. ClumpFind (Williams, de Geus & Blitz 1994), also used by Moore et al. (2007), Buckle et al. (2010) and Parsons et al. (2012), was used to carry out a consistent search for the position of peak emission in this study. The search was undertaken on two-dimensional images, intensity integrated over the emission peaks’ velocity ranges.

In a few cases, ClumpFind reported more than one clump coordinate per image, likely due to the irregular structure and crowded environment of massive star-forming regions. The purpose of using ClumpFind in this study was to find the position of the peak molecular emission in the vicinity of each methanol maser target. Multiple clumps were accepted if they were further than a beam width apart and not close to the edges of the image. Multiple spectra per field of view were extracted at these positions. In four cases (marked with asterisk in Table 1), ClumpFind did not detect any clumps, either due to low signal-to-noise or the physical area of the emission being too small to satisfy the ClumpFind criteria (minimum seven pixels). In these cases, we did detect some emission at the maser coordinate, hence we used the maser coordinates as the location for spectrum extraction.

Where clumps were detected close to a dead receptor or the edge of the map, they were rejected from further analysis, as any extracted spectra and derived results will be incomplete. This is the case for the maps of the targets associated with masers G 24.790+0.083 (clump 2), G 30.851+0.123 and G 31.182−0.148 (clump 2).

Of the original 70 targets observed, reliable clump detections were obtained in 54 maps (77 per cent), and because more than one clump was found in some images, a total of 58 positions were analysed. The positions of the observed clumps are summarized in Table 1. Intensity-integrated spatial maps were created for these targets, and are shown online in Appendix A as Supporting Information, where the integrated C¹⁸O emission is contoured over the background of ¹³CO emission, the latter integrated over velocity ranges v_low to v_high, listed in Table 2.

Contour intervals are shown in steps of the integrated noise rms, ϕ, where |$\phi = \sigma _{{\rm {{\rm rms}}}} \Delta v \sqrt{n_{{\rm chan}}}$|⁠, calculated for the same number of channels, n_chan, over which the C¹⁸O image is integrated, with |$\sigma _{{\rm {{\rm rms}}}}$| the average rms per channel and Δv = 0.5 km s⁻¹ the velocity range per channel. All values are listed in Table 1. Where contour intervals are larger than 2ϕ, this is indicated in the figure captions. The lowest level contour was manually selected for each image, because every image has a unique signal-to-noise and background. The lowest contour levels ranged between 3ϕ and 14ϕ. Both maser and clump coordinates are indicated on the maps, shown respectively as a star and circle symbol.

Sometimes, C¹⁸O noise levels were excessive due to poor atmospheric transmission, as this line is at the edge of the atmospheric window, which makes it susceptible to small changes in the water vapour column. Targets associated with masers G 30.370+0.482 (clump 2), G 31.182−0.148 (clump 1), G 33.317−0.360 and G 33.486+0.040 did not show sufficiently high signal-to-noise to isolate any clump emission above a 3 × ϕ threshold. For these targets, the C¹⁸O maps are not shown in Appendix A.

Spectrum extraction and wing selection

After locating the emission peak, the spectrum was extracted at this position from both the ¹³CO and C¹⁸O cubes. Table 2 lists the maser median velocity taken from the literature, or its associated IRDC velocity if the former was not available (Simon et al. 2006) and the literature references of these values, for each target. It also gives the measured peak main-beam efficiency-corrected temperatures (T_mb) and corresponding velocities for both ¹³CO and C¹⁸O. Sources marked with an asterisk exhibit self-absorption dips in their ¹³CO spectra. For these spectra, a Gaussian was fitted to the shoulders of the absorbed spectrum and the peak of this resultant profile was used as the estimate of peak temperature. The peak from the Gaussian fit showed on average a ∼13 per cent increase with respect to the peak T_mb of the original, absorbed spectrum, with three extreme cases of a 30–40 per cent increase. Two examples of these sources and their Gaussian fits are shown in Fig. 2. Plotted C¹⁸O spectra give an indication where the optically thin peak is expected. In the case of the double-peaked target G 23.010−0.411, the values marked with an asterisk in Table 2 represent peak values of fits to the individual peaks. The values of Δv_b and Δv_r in columns 9 and 10 are the blue and red velocities relative to the peak C¹⁸O velocity, measured respectively from each wing extreme, to be discussed in Section 4.4. The use of Int_b and Int_r in columns 11 and 12 will be discussed in Section 3.3.

Following Codella et al. (2004), the optically thin C¹⁸O profiles were used as tracers for the line cores of targets. The C¹⁸O spectra were scaled to the ¹³CO peak temperature. To avoid subtracting any emission from higher velocity features that may be present in the C¹⁸O if densities were sufficiently high, a Gaussian was fitted to the C¹⁸O peak to approximate the line core-only emission. This was done by gradually removing points from the outer (higher velocity) edges of the C¹⁸O spectrum until the peak could be fitted, following the same approach as van der Walt, Sobolev & Butner (2007); see Fig. 3(a). The scaled Gaussian fit was then subtracted from the ¹³CO spectra to show the velocity ranges in the line wings where there is excess emission in ¹³CO.

G 23.010−0.411 is a special case with a double-peaked profile. Assuming that this is caused by two separate but closely associated clumps, we used two Gaussians, each fitted just to the highest velocity shoulder of each C¹⁸O line peak. Whenever absorption dips occur in the ¹³CO profiles, no natural profile peak existed. Instead, the C¹⁸O spectra were scaled to the peak of the previous Gaussian fitted to the ¹³CO.

This Gaussian was then subtracted from the ¹³CO profile. The line wings are defined by the sections where the ¹³CO profile is broader than the scaled Gaussian representing the C¹⁸O line core emission, provided the ¹³CO corrected antenna temperature is higher than 3σ (σ is the noise per 0.5 km s⁻¹ channel, averaged over a 30 km s⁻¹ section of the emission-free spectrum). An example of this wing selection process is shown in Fig. 3(b), which shows the ¹³CO residual spectrum and discrete spectral points that satisfy the wing criteria (empty circles are blue, and solid circles are red).

There is a risk that some blue and red emission might be missed by analysing a single spectrum at the location of the clump peak. Therefore, when the position of peak intensity in both the blue and red integrated images was found (mapping of blue and red images is explained in Section 3.3), another two additional spectra, called the ‘red-wing spectrum’ and ‘blue-wing spectrum’, were extracted. Once again blue and red residual spectra were calculated. If broader wing emission was found, the initial wing ranges were expanded to incorporate the ranges covered by the red-wing and blue-wing spectra. The final velocity ranges for blue and red wings are listed in Table 2.

Mapping the outflows

The final blue- and redshifted velocity ranges are used to produce two-dimensional ¹³CO intensity integrated images corresponding to each wing. These are overlaid as solid blue- and dotted red contours on to the ¹³CO-integrated intensity image, representing the outflow lobes. Two examples are shown in Fig. 4, showing target G 24.493−0.039 with the maser and clump coordinates overlapping, and target G 29.956−0.016A with an offset between the maser and clump coordinates. The remainder of the maps are shown online in Appendix B as Supporting Information. Contours are plotted in 10 per cent intervals up to 90 per cent of the maximum intensity, Int_b or Int_r, for each integrated image (values listed in columns 10 and 12 in Table 2). The lowest contour is never lower than 30 per cent, but values differ for each image depending on the individual background brightness levels. The lowest contour is selected by eye as the level which encompassed the outflow lobe clearly.

As massive stars form in clusters, the observed targets often have contamination from similarly HV components as the outflow, but from different spatial structures in the field of view (Shepherd & Churchwell 1996a). This makes it difficult to isolate the outflow. Therefore, if identified as belonging to such structures, these pixels were flagged to be bad in any further analysis. Sometimes one or both of the outflows are partially cut off where they are situated close to the edge of the field of view or to a dead receptor. These sources are flagged as such in the second to last column in Table 5 and their calculated properties only serve as a lower limit because a fraction of the emission is not included in the analysis.

Three of the 58 analysed clumps have been too close to the edge of the field of view for any significant information to be derived and are excluded from further analysis. Out of the remaining 55 maps, 47 outflows are clearly bipolar (85), with the eight exceptions marked with a superscript in Table 2. For a sample of high-mass protostellar objects, Beuther et al. (2002b) had a bipolar outflow detection frequency of 81 in ¹²CO, comparable with what we find.

RESULTS

Detection frequency

All of the 58 spectra available for analysis (see Section 3.1) were found to have HV outflow signatures, either in the spectra or in the contour maps, resulting in a 100 per cent detection rate. Such a high detection rate of outflows towards massive YSOs is not uncommon. Shepherd & Churchwell (1996b) searched for ¹²CO(J = 1-0) HV line wings towards 122 high-mass star-forming regions and detected low-intensity line wings in 94 of them. Of these 94, 90 per cent were associated with HV gas in the beam. The argument has already been made at that stage, that if the HV gas is due to bipolar outflows, molecular outflows are a common property of newly formed massive stars. Sridharan et al. (2002) detected 84 per cent of sources with HV gas from a ¹²CO |${\rm (J = 2{\rm -}1)}$| survey of 69 protostellar candidates.

Zhang et al. (2001, 2005) observed a sample of 69 luminous IRAS point sources in CO (J = 2-1) and detected 39 molecular outflows towards them (57 per cent). They found the search for outflows hampered for Galactic longitudes <50° (due to confusion my multiple cloud components when observing in this transition). A total of 39 objects were outside of this region, towards which 35 outflows were detected, resulting in a 90 per cent outflow detection rate.

Kim & Kurtz (2006) observed 12 sources from the same Molinari et al. (1996) catalogue that Zhang et al. (2001) selected their sources from. They detected outflows in 10 sources and adding these sources to the detections from Zhang, results in a detection rate of 88 per cent ([35 + 10 − 3 = 42] out of [39 + 12 − 3 = 48]), taking into account that there are three sources in common between the two samples. More recently, López-Sepulcre et al. (2009) searched for molecular outflows towards a sample of eleven very luminous massive YSOs. They found HV wings, indicative of outflow motions, in 100 per cent of the sample.

Three further studies have dealt specifically with class II methanol masers. Codella et al. (2004) surveyed for molecular outflows towards 136 UCHii regions, out of which 56 positions showed either 6.7 GHz methanol or 22.2 GHz water maser emission. Their overall outflow detection rate from ¹³CO(J = 1-0) and (J = 2-1) transition lines was ∼39 per cent, but they found that in cases where observations were made towards 6.7 GHz methanol or 22.2 GHz water maser emission lines, the outflow detection rate increased to 50 per cent. As their observations were single pointings, they may have missed some outflows that were offset from the masers. Xu et al. (2006) studied molecular outflows using high-resolution CO (J = 1-0) mapping towards eight 6.7 GHz methanol masers closer than 1.5 kpc. They found outflows in seven of them, an 88 per cent detection rate. Wu et al. (2010) investigated the distinctions between low- and high-luminosity 6.7 GHz methanol masers via multiline mapping observations of various molecular lines, including ¹²CO(J = 1-0), towards a sample of these masers. They found outflows to be common among both sets of masers: of the low-luminosity masers, they found six outflows out of nine, and from the high-luminosity masers they found four outflows out of eight, an overall detection rate of 59 per cent. Note that the detection frequencies from both Xu et al. (2006) and Wu et al. (2010) are obtained from small number samples.

All these results suggest that the majority of massive YSOs have molecular outflows, and should 6.7 GHz methanol masers be present, they are closely associated with the outflow phase.

Maser distances

Green & McClure-Griffiths (2011) published the kinematic distances for about 50 per cent of the targets in this study, using the 6.7 GHz maser mid-velocity as an estimate for the systemic velocity. They used respectively the presence/absence of self-absorption in H i spectra in the proximity of the systemic velocity, to determine whether the source is at the near/far kinematic distance. However, methanol maser emission often consists of a number of strong peaks spread over several km s⁻¹. As differences of only a few km s⁻¹ in the velocity of the local standard of rest v_lsr can be enough to change the kinematic distance solution from near to far and vice versa in the H i absorption feature method of resolving the former, using estimated v_lsr values form the maser emission could lead to an incorrect distance solution. Therefore, molecular line observations provide more reliable measurements of the clump systemic velocity (Urquhart et al. 2014).

For this reason, and to prevent the additional uncertainties introduced by adopting distances calculated using different techniques by different authors, we decided to recalculate all distances of the methanol masers using their associated C¹⁸O peak velocities. This was done using the Galactic Rotation Curve (GRC) as fitted by Brand & Blitz (1993), with the Sun's Galactocentric distance, R₀, assumed as 8.5 kpc and its circular rotation Θ₀ as 220 km s⁻¹. Calculated values are listed in Table 3. The average difference between distances calculated using the GRC, versus the maser distances listed by Green & McClure-Griffiths (2011) for our targets, is 0.8 ± 0.6 kpc.

Alternative solutions to Brand & Blitz (1993) for calculating the kinematic distances are Reid et al. (e.g. 2009) and Clemens (1985). Urquhart et al. (2014) found that for sources in the inner Galaxy, the distances given, respectively, by these rotation curves, all agree within a few tenths of a kpc, which are smaller than their associated uncertainty of the order ±1 kpc due to streaming motions (Urquhart et al. 2011, 2012). Consequently, the statistical results are robust against the choice of model.

The galactocentric distances obtained using the GRC from Brand & Blitz (1993) were geometrically converted to heliocentric distances, of which two solutions exist within the solar circle, called the kinematic distance ambiguity (KDA). These distances are equally spaced on either side of the tangent position, and are generally referred to as the near and far distances. Sources with velocities within 10 km s⁻¹ of the tangent velocity are placed at the tangent distance (indicated by TAN in the reference column), since the error in the distance is comparable to the difference between the near/far distance and the tangent distance.

As Green & McClure-Griffiths (2011) could resolve the KDA for many masers using corresponding H i self-absorption profiles, when available, we used their values to resolve the KDAs for our targets. Where maser distances were not published in Green & McClure-Griffiths (2011), alternative sources were used where available, being other publications of 6.7 GHz masers (Purcell et al. 2006; Caswell & Green 2011), EGOs (Cyganowski et al. 2009), OH-masers (Fish et al. 2003), associated IRDCs (Simon et al. 2006), or molecular clouds (Roman-Duval et al. 2009) – if the maser position fell within the cloud as well as within ∼ 5 km s^{− 1} of the cloud's v_lsr – were used to resolve our KDAs. The literature reference used to resolve the distance for each source, should it exist, is also listed in Table 3.

For four out of the 58 targets, no published distance could be found to resolve the KDA. In these cases, both the near and far values are listed. Three targets are rejected from this list in Table 3: G 24.790+0.083B, because of too high noise, and G 30.851+0.123 and G 31.182−0.148B, because the clumps are mostly cut off by the edge of the field of view. The columns after the distance columns, list outflow lobe surface areas and lengths, values used for calculations discussed in Section 4.4.

Dealing with uncertainties

Calculation of the physical properties of molecular outflows can provide useful information on the obscured driving source. These calculations are subject to a number of uncertainties, most prominent of which is the outflow orientation. However, this is not easily determined (e.g. Shepherd & Churchwell 1996b; Curtis et al. 2010) and, as such, no correction is applied in this study. Although we will not apply any corrections, we discuss in the following paragraphs the corrections usually applied in the literature for outflow orientation. We also report the effects of such corrections on the calculated outflow properties, as well as additional contributors to uncertainties.

Our observations were some of the first to be carried out with the HARP instrument and some of the receptors exhibited poor performance and did not yield useful data. Often more than two receptors had to be switched off. At times this resulted in some of the clump/outflow emission being missed. Potentially, this could also result in outflow lobes not being detected at all. Blue and red contour levels were determined by eye, since each image is uniquely characterized by the background noise, emission brightness and available receptors. The 14 arcsec beam of the telescope places a limit on the size of outflows that can be resolved, especially for the more distant targets.

The most significant of the above uncertainties, is θ, the angle of the outflow axis with respect to the line of sight. As only a projection of the outflow is observed, any inclination with respect to the plane of sky will reduce the length of the outflow (not the width) by sin(θ), and increase the observed Doppler broadening by cos(θ). Cabrit & Bertout (1990) give a detailed discussion of the effect of inclination angle. Due to the lack of a specific orientation for each outflow, many authors assume a mean inclination angle for their sample to correct the calculated outflow parameters. The most commonly used angle is 57| $_{.}^{\circ}$|3, determined using the assumption that outflows are distributed uniformly and with random inclinations to the line of sight (Bontemps et al. 1996; Beuther et al. 2002b; Hatchell, Fuller & Richer 2007; Curtis, Richer & Buckle 2010). Table 4 summarizes the corrections due to inclination for the outflow parameters calculated (see Section 4.4). Unknown inclinations mostly cause the outflow parameters to be underestimated. Time-scales t_d are thus likely to represent a lower limit to the true age of the outflows (Parker, Padman & Scott 1991) and hence also to the time over which the embedded protostars responsible for the outflows have been accreting from their surroundings (Beuther et al. 2002b).

Other contributors to uncertainties are possible difficulty separating the outflowing gas from the ambient gas; higher interstellar extinction towards the molecular ring in the inner Galaxy (0° < l < 50°) in addition to their internal extinction (Zhang et al. 2005); different CO/H₂ abundance ratios used by different authors (e.g. Rodriguez et al. 1982; Cabrit & Bertout 1992; Herbst & van Dishoeck 2009).

Some authors use the mean atomic weight of the mixture of hydrogen and helium gas (Garden et al. 1991), while others only consider pure hydrogen molecular gas (Snell et al. 1984), resulting in a difference of 0.36 amu in the mean atomic weight (Wu et al. 2004). The excitation temperature, T_ex, is assumed to range from 30 to 50 K for high-mass sources (Shepherd & Churchwell 1996b; Beuther et al. 2002b; Wu et al. 2004). However, a constant temperature assumption will underestimate the kinetic energy for an outflow with high jet/ambient density contrast (Downes & Cabrit 2007). Contamination from additional unrelated velocity components within the telescope beam could make it difficult to isolate the outflow, unless the components have a different spatial distribution from the outflow gas (Shepherd & Churchwell 1996a). Finally, even though we used ¹³CO as a tracer, we note that for their ¹²CO observations, Cabrit & Bertout (1990) estimated typical errors in the outflow parameters that reflect uncertainties in ¹²CO/H₂, distance determinations, T_ex, inclination angles, optical depth effects and possible low-level contamination of ¹²CO emission in the reference position. These error values are a factor ∼3 on outflow mass M_out, a factor ∼10 on mechanical force F_m, and a factor ∼30 on mechanical luminosity L_m.

Calculation of outflow physical properties

¹³CO is a less abundant molecule than ¹²CO, thus exhibiting a narrower spectral profile. A sample of 56 sources for which both ¹²CO and ¹³CO spectra were published has been investigated and the average ¹²CO/¹³CO full width zero intensity ratio is found to be ∼2 with a standard deviation of 1.3 (Cabrit, Goldsmith & Snell 1988; Shepherd & Churchwell 1996a; Su, Zhang & Lim 2004; Bronfman et al. 2008; Narayanan, Snell & Bemis 2012; Ortega et al. 2012; Xu & Wang 2013). All calculations containing wing velocities relative to the systemic velocity are scaled by this factor, implying a factor 2 increase in p and factor 4 increase in E.

For sources with unresolved distances, both values are shown and distinguished by the numbers next to the target names (1 = far and 2 = near). Exclusions are then made in further analyses for targets which have (i) kinematic distance ambiguities, hence uncertainties in calculated physical parameters, or, (ii) offsets of more than 3 pixels (18 arcsec, of the order of a beam size) between the maser coordinate and peak CO emission. These targets are marked as such in the notes of Table 5. Following the exclusions, we are left with 44 outflows in our sample that are positionally associated with methanol masers and for which we can calculate physical properties that are unaffected by distance ambiguities. We refer to this sample as MMAOs, indicated as such in Table 5, and base all further discussion on these outflows.

Clump masses

The evolutionary sequence of massive stars begins with prestellar clumps and cores, which are gravitationally bound overdensities inside a molecular cloud that show signs of inward motion before a protostar starts forming (Zinnecker & Yorke 2007; Dunham et al. 2011). Most massive stars form in star clusters (e.g. Clarke, Bonnell & Hillenbrand 2000; Lada & Lada 2003), which are part of a hierarchical structure, defined by Williams, Blitz & McKee (2000) and summarized by Bergin & Tafalla (2007). In this approach, the largest structure is a molecular cloud, with masses of the order 10³–10⁴ M_⊙ and diameters ranging from 2 to 15 pc. Clouds contain subunits of enhanced density gas and dust, called clumps, wherein the earliest stages of massive star formation take place. Clumps will typically form stellar clusters (Williams et al. 2000).

Studies of massive star formation regions showed that clumps generally have sizes of the order of ∼1 pc, and masses ranging from order 10 M_⊙ to ∼10³–10⁴ M_⊙ (Kurtz et al. 2000; Smith, Longmore & Bonnell 2009). They are defined to be coherent in position–velocity space. Smith et al. showed that the gravitational potential of these clumps causes global collapse, which channels mass from large radii towards the centre of the cluster, where protostars with the greatest gravitational radius accrete it, causing them to become massive. Stars (or multiple systems such as binaries) eventually form from gravitationally bound subunits in the clumps, called cores (Williams et al. 2000). Cores have sizes typically ≤0.1 pc and masses ranging from 0.5 M_⊙ up to ∼10²–10³ M_⊙ (Kurtz et al. 2000; Smith et al. 2009).

We have corresponding C¹⁸O maps for 51 out of the 55 ¹³CO maps. The optically thin C¹⁸O serves as a useful tracer of the central clump (e.g. López-Sepulcre et al. 2009). With a median source distance of 7.2 kpc and telescope beam of 14 arcsec, our resolving power is of the order of 0.5 pc, which, given the above definitions, implies the traced structures are more likely clumps than cores.

However, López-Sepulcre et al. (2009) noted that their dust clump masses are a factor ∼5 larger than the corresponding C¹⁸O masses, and stated that this difference might be explained by the fact that C¹⁸O and the submm continuum are tracing different parts of the clump. They found the angular full width at half-maximum (FWHM) measured in the submm continuum surveys to be a factor ∼2.5 larger than what is mapped by C¹⁸O (J = 2–1), and speculated that this played a main role in the difference between the two mass estimates. Hofner et al. (2000) also concluded from their survey that masses derived from submm dust emission, tend to be systematically higher than masses derived from C¹⁸O by a factor ∼2. They pointed out that contributing sources of uncertainty to this discrepancy could be C¹⁸O abundance, optical depth estimates, and the dust grain emissivity adopted.

Therefore, we also use continuum measurements to calculate the clump masses associated with MMAOs, and use the latter in all further discussions. The 870 μm flux measurements were obtained from the the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL) survey (Schuller et al. 2009; Contreras et al. 2013), using offsets within a FWHM beam (beam size 19 arcsec) as matching criteria. Using the matching fluxes from Csengeri et al. (2014) for the targets in this study, clump masses were calculated following Urquhart et al. (2013a), with a gas-to-dust mass ratio assumed to be 100, dust absorption coefficient κ_ν of 1.85 cm² g⁻¹ and dust temperature of 20 K. All values are listed in Table 6. Two clump masses are listed for sources with distance ambiguities, marked with (1) next to the name for the far distance, and (2) for the near distance. Excluding all targets with distance ambiguities, these clump masses range from ∼30 to 1.4 × 10⁴ M_⊙, have a mean value of ∼(2.5 ± 0.5) × 10³ M_⊙ and median of ∼1.3 × 10³ M_⊙. The clump masses associated with the MMAO subset have a mean of ∼2.8 × 10³ M_⊙. For the targets with resolved distances, 96 per cent have masses >10² M_⊙ and 49 per cent have masses of the order of 10³–10⁴ M_⊙. This confirms that the majority of these targets are likely classified as clumps, as per definitions given above.

We find that the clump masses derived from dust measurements for our targets, are on average a factor 8 higher than masses derived using their C¹⁸O emission, in agreement with López-Sepulcre et al. (2009) and Hofner et al. (2000).

DISCUSSION

Clump and outflow mass relations

Here, we investigate the relationships between outflow properties and the mass of the clumps that they are associated with. While it is not possible to resolve the contribution from individual stars or protostellar cores in our data we can at least infer a relationship between clump mass and the most massive star present in the clump (e.g. Urquhart et al. 2013b).

McKee & Tan (2003) derived a relation for the accretion rate of a free-falling envelope as a function of time, wherein it is proportional to the clump surface density Σ^0.75. This implies higher mass accretion in the most massive and dense clumps, compared to those forming in lower mass clumps. The former will arrive at the main-sequence sooner, and hence form an H ii region much more quickly than less massive stars forming in lower mass clumps. Urquhart et al. (2014) found their results to confirm this hypothesis, being consistent with the decreasing time-scale with increasing massive YSO luminosity found by Mottram et al. (2011).

In addition to this, Urquhart et al. also found that the most massive stars form predominantly towards the centres of their spherical, centrally condensed host clumps, where the highest densities exists and the gravitational potential is the deepest. These stars evolve much quicker than those outside the central region, and reach the main sequence in ∼10⁵ yr, ahead of the lower mass stars which can take ∼10 times longer. Thus, while the most massive star, traced by the 6.7 GHz methanol maser in this study, is near to joining the main sequence, it is likely that the lower mass members of the protocluster still need to evolve to a stage where they make a significant contribution to the observed luminosity.

Urquhart et al. (2013b) found that for the H ii regions they investigated, the bolometric luminosities (which are effectively a measure of the whole protocluster luminosity), were very similar to values estimated from the radio continuum flux (which only trace the most massive stars). This suggests that it is likely that the bolometric luminosity is actually dominated by the most massive stars. Urquhart et al. (2014) also found that there is a strong correlation between the clump masses and the bolometric luminosities of the most massive stars, but that the total clump luminosities are much lower than would be expected from the fully formed cluster. These findings agree with the stated hypothesis, that the most massive stars have very rapid evolution times, and are consequently likely to dominate the observed clump properties. Thus, if the luminosity of a clump is dominated by the most massive stars, then it is reasonable to assume that so too, is the luminosity, and consequently the energetics, of the outflow.

In the following, we compensate for the loss of one lobe's mass in monopolar targets by doubling the values of the detected lobe for M_out, p, E, |$\dot{M}{_{\rm out}}$|⁠, F_m and L_m, as they all depend on the outflow mass. Figs 5(a) and (b) shows the relation between, respectively, outflow masses and mass-loss rates, and clump masses for the MMAO sample (solid circles). In order to compare these relations with low-mass YSOs, we obtained associated outflow masses from Wu et al. (2004), Narayanan et al. (2012) and Davis et al. (2010) for 13 of the core envelope masses listed by Bontemps et al. (1996, empty circles). The pink squares represent the clump masses from Beuther et al. (2002a), derived from 1.2 mm dust continuum data, together with their associated outflow masses and mass-loss rates, as calculated in Beuther et al. (2002b). These values are corrected according to the erratum that was later published by Beuther et al. (2005), where the authors explain that their grain emissivity approximation should be a factor 2 higher than initially calculated, which would cause their derived clump masses to be a factor 2 lower than reported.

In Fig. 5(a), a best-fitting power law to the MMAO sample (circles) is given by log(M_out/M_⊙) = (−0.8 ± 0.3) + (0.8 ± 0.1)log(M_clump/M_⊙), holding over three orders of magnitude for massive outflows, and extending towards the low-mass regime to cover six orders of magnitude in total.

López-Sepulcre et al. (2009) converted their clump masses derived from C¹⁸O to dust masses, and found a tight correlation between M_clump and M_out, described by M_out = 0.3M^0.8_clump. Sánchez-Monge et al. (2013) did a similar fit, but they calculated their outflow masses from SiO observations and clump masses from SED fits to the Hershel infrared Galactic Plane Survey (Hi-GAL) data. They found the same relation as López-Sepulcre et al. (2009). The power law found for MMAOs agrees with these authors within uncertainties. It is interesting to note that the outflow masses and mass-loss rates for MMAOs are generally higher than those estimated by Beuther et al. (2002b). We will discuss this property more detail in an upcoming second paper (de Villiers et al., in preparation).

Fig. 5(b) shows the best-fitting linear relation between the logarithmic values for the mass-loss rates and clump masses for MMAOs, given by |${\rm log}(\dot{M}_{{\rm out}}/{\mathrm{M}_{{\odot }} \,\mathrm{yr}^{-1}}) = (-5.0 \pm 0.4)+(0.6 \pm 0.1){\rm log}(M{_{\rm clump}}/{\mathrm{M}_{{\odot }}})$|⁠. It again holds over three orders of magnitude for massive outflows, and extend to the low-mass regime to cover six orders of magnitude in total.

As both M_out and |$\dot{M}_{{\rm out}}$| depend on mass and distance, a nonparametric measure of the statistical dependence between these parameters and M_clump is needed. The Spearman-rank test is used, similarly to Ridge & Moore (2001), where a perfect positive or negative correlation exists between the ranks when r_s is ±1. No correlation exists when r_s = 0. The Spearman parameter r_s is related to the two-tailed Student's t-test via, |$t_{\rm s} = r_{\rm s} \sqrt{(n-2)/(1-r^{2}_{\rm s})}$| for a sample size n. For no significant correlation, |$|t_{\rm s}| < |t_{s{_{\rm crit}}|}$|⁠. For a sample size of 43 (MMAO sample of 44 had 43 matches from ATLAGAL for clump masses), |$t_{s{_{\rm crit}}}$| is ±2.02. The relation between M_out and M_clump for MMAOs is statistically significant with r_s = 0.59 (t_s = 4.66). The same is true for the relation between |$\dot{M}_{{\rm out}}$| and M_clump, where r_s = 0.64 (t_s = 5.28).

The result that both relations in Fig. 5 are found to hold over six orders of magnitude when extrapolated to the low-mass regime, suggests that a similar process causes outflows in both low- and high-mass star formation.

Fig. 6(a) shows the relation between the outflow mass and dynamical time-scales for MMAOs. A power law fitted to the data resulted in a poor correlation given by log(M_out/M_⊙) = (0.3 ± 1.2) + (0.3 ± 0.3)log(t/ yr^{− 1}), with r_s = 0.18 (t_s = 1.18, for a complete MMAO sample of 44) implying no significant correlation, as also found by Wu et al. (2004). Our range of dynamical times is possibly too small for any significant evolutionary trends, such as the increase in outflow mass, to be observed.

Shepherd & Churchwell (1996a) plotted the same relation for their 10 massive star-forming regions. They divided their plot into three regions based on the bolometric luminosities of the outflow sources in Cabrit & Bertout (1992). We show the same luminosity boundaries in Fig. 6(a). Using data from Cabrit & Bertout (1992), Shepherd & Churchwell (1996a) found that sources with L_* < 10² L_⊙ were located below the bottom line and sources with L_* > 2 × 10⁴ L_⊙ above the top. This sectioning contributes a useful estimate of the expected bolometric luminosity of the YSO associated with each outflow studied – a property that could not be derived directly from the available data. Consistent with what is expected for massive YSO's, the majority (95 per cent) of the sources occur above the L_* = 10² L_⊙ boundary, with ∼79 per cent of them between L_* = 10² L_⊙ and L_* = 2 × 10⁴ L_⊙, and ∼16 per cent above L_* = 2 × 10⁴ L_⊙. Also, we find that higher luminosity (and hence more massive) YSOs will drive the higher mass outflows we study in this paper. A more massive accreting YSO is likely to have higher angular momentum, hence power larger outflows. In addition to this, star formation regions that form O and B stars tend to have larger reservoirs of material available to be entrained into a molecular outflow.

Fig. 6(b) shows the relation between the CO momentum flux/mechanical force F_CO and the clump masses for MMAOs (circles), together with the values from Beuther et al. (2002b) (pink squares), as well as with 33 low-mass sources from Bontemps et al. (1996, triangles). F_CO is the inclination corrected mechanical force, derived from F_m by applying a correction factor of 2.9, corresponding to a mean inclination angle of 57| $_{.}^{\circ}$|3 (i.e. assuming random outflow orientations, Beuther et al. 2002b and Table 4). Although we do not apply inclination corrections to the data in this study, in this specific case we correct F_m values in order to compare the CO momentum flux values from MMAOs with those from Beuther et al. (2002b) and Bontemps et al. (1996). Once again, as in Fig. 5, MMAOs’ momentum flux values are generally higher than those estimated by Beuther et al. (2002b, also to be discussed by de Villiers et al., in preparation).

The outflow's mechanical force is a very important parameter in studying the early phases of star formation. It is a measure of the rate at which momentum is injected from the underlying driving agent, via interactions with the molecular gas in the core, into the envelope. In other words, it is a measure of the outflow's strength and used to understand the driving mechanism of outflows (Bachiller & Tafalla 1999; Downes & Cabrit 2007). Hatchell et al. (2007) suggests that, although it is difficult to separate this effect from contamination due to the initial conditions, the correlation between F_CO and M_clump suggests that the outflow activity declines during the later stages of the accretion phase (when the clump mass decreases).

Many authors have studied this relation before, both for low-mass outflows (e.g. Bontemps et al. 1996; Hatchell et al. 2007; Takahashi et al. 2008; Curtis et al. 2010; van der Marel et al. 2013) and high-mass outflows (e.g. Henning et al. 2000; Beuther et al. 2002b). For our massive MMAOs, we find the relation to be log(F_CO/M_⊙ km s⁻¹ yr⁻¹) = (−4.5 ± 0.4) + (0.8 ± 0.1)log(M_clump/M_⊙) (black solid line) with a Spearman-rank coefficient of r_s = 0.66(t_s = 5.62). This relation found for MMAOs corresponds within uncertainties to what Bontemps et al. (1996) found for low-mass sources, log(F_CO/M_⊙ km s⁻¹ yr⁻¹) = −(4.2 ± 0.1) + (1.1 ± 0.2)log(M_clump/M_⊙) (blue dotted line).

If, as stated in the beginning of this section, the outflow energetics are dominated by the most massive cores in a clump harbouring massive star formation, the presence of the same relationship between the outflow parameters and the clump masses for both high and low masses, is consistent with the theory that massive star formation is an scaled-up version of the low-mass scenario.

Mass-loss rate and accretion rates

If massive outflows are produced by the same mechanism as low-mass outflows, it implies that the outflows are momentum driven by the jet coming from the central YSO which entrains the surrounding molecular gas and forms the outflow. One of the proposed solutions, enabling disc accretion to be a possible formation process for massive stars, is when accretion proceeds through a disc with a high enough mass accretion rate to overcome the radiation pressure of the central massive star (Jijina & Adams 1996; Yorke & Sonnhalter 2002). Even more involved disc accretion models including asymmetric configurations and Rayleigh–Taylor instabilities (Krumholz, Klein & McKee 2007; Krumholz et al. 2009), assume accretion rates of the order of 10^{− 4} M_⊙ yr^{− 1}. High accretion rates (∼10⁻⁵–10⁻³ M_⊙ yr⁻¹) are also used in more recent disc accretion models that incorporate a dust sublimation front and consequently yielded the growth of the highest mass stars ever formed in multidimensional radiation hydrodynamic simulations (Kuiper et al. 2010).

The mean |$\dot{M}{_{\rm out}}$| for MMAOs is ∼1.7 × 10⁻³ M_⊙ yr⁻¹, which would lead to accretion rates of ∼3 × 10⁻⁴ M_⊙ yr⁻¹, being of the same order of magnitude as the ∼10⁻⁴ M_⊙ yr⁻¹ found by both Beuther et al. (2002b) and Kim & Kurtz (2006). Our approximate accretion rate for MMAOs is also much higher than typical accretion rates of |${\rm 10^{-7}\hbox{--}10^{-5} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}}$| expected for low-mass YSOs (Shu 1977), and agrees with the theoretical rates used in disc accretion models for massive stars (Krumholz et al. 2007; Kuiper et al. 2010).

Comparison with methanol masers

Integrated 6.7 GHz flux densities (S_i) are obtained from the MMB catalogue (Breen et al., in preparation) for 38 of the MMAO masers. Using the distance of each target, S_i is converted to an integrated spectral luminosity L_6.7 GHz.

Fig. 7 shows the relation between the outflow mechanical luminosities (or energy supply rates) and maser luminosities for 38 masers from the MMAO subset, given as log(L_m/ L_⊙) = (3.2 ± 1.3) + (0.2 ± 0.2)log(L_6.7  GHz/ L_⊙) and a Spearman-rank coefficient of r_s = 0.28(t_s = 1.78), which, although not quite, is at the margin of being statistically significant as the 5 per cent acceptance interval for t is ±2.02.

Furthermore, following their positions in Fig. 6(a), we divide all the outflows plotted in Fig. 7 into three categories; outflows associated with clumps whose luminosities above 2 × 10⁴ L_⊙, between the L_* = 10² L_⊙ and L_* = 2 × 10⁴ L_⊙, and below L_* = 10² L_⊙. The data points in Fig. 7 are marked accordingly. Although bright outflows are mostly associated with bright clumps, and low-luminosity outflows associated with lower luminosity clumps, there is no such preference for maser brightness.

Finally, we plot the relation between the 6.7 GHz maser luminosities against clump masses in Fig. 8 and found the best-fitting linear relation between the logarithmic values of these parameters to be log(L_6.7 GHz/ L_⊙) = ( − 10.6 ± 0.7) + (0.8 ± 0.2)log(M_clump/M_⊙). The associated value for r_s is 0.59 (t_s = 4.36), which implies a statistically significant positive correlation.

Breen et al. (2010) studied the relation between 101 1.2 mm dust clumps and their associated 6.7 GHz methanol masers. Their study showed that more luminous 6.7 GHz methanol masers are likely to be associated with dust clump sources that have bigger radii, higher mass and lower clump densities (the latter possibly due to the assumption of a constant dust temperature). These results are represented graphically in their fig. A2. Wu et al. (2010) also found molecular clumps (derived from NH₃ lines) associated with high-luminosity 6.7 GHz masers to be larger and ∼10 times more massive than those associated with low-luminosity masers. They also found that outflows associated with high-luminosity masers have wider line wings and larger sizes than those associated with low-luminosity masers. This lead to their interpretation that masers with higher luminosities are associated with YSOs with larger masses.

Urquhart et al. (2013a) found a weak linear correlation between the maser luminosity and clump masses for 442 ATLASGAL–MMB associations. They speculated that this correlation may be related to the fact that the most massive clumps are likely to form more massive stars, and that a higher isotropic maser luminosity is somehow related to a higher stellar luminosity.

The results from Figs 8 and 7 agrees with above authors and suggest that there is some correlation between the mass (hence also brightness) of a massive YSO and both the luminosity of the outflow it generates as well as the total 6.7 GHz maser luminosity it pumps. Note that in this case, we know that a higher maser luminosity could not be due to a larger number of masers in more massive regions, as all our MMB-sample selected masers have supporting interferometric observations that reveal they are principally single-maser spots.

SUMMARY AND CONCLUSIONS

We analysed the ¹³CO and C¹⁸O spectra extracted at 58 ¹³CO clump coordinates towards 6.7 GHz methanol maser coordinates between 20° < l < 34°. All spectra showed HV outflow signatures. Of these, the HV emission was mapped for 55 with 47 showing bipolar structures. A subset containing 44 targets is referred to as MMAOs, with the criteria that all members have to have resolved kinematic distances and be closely associated with a methanol maser. Only MMAOs were used for further analysis.

The wing spectra and spatial maps were used to calculate the physical properties of all 55 mapped outflows, generally following Beuther et al. (2002b), with a few adjustments in their approach for this study. The associated clump masses for MMAOs were calculated from the 870 μm flux measurements towards these targets, of the ATLASGAL survey (Csengeri et al. 2014).

The main results for this study can be summarized as follows. (i) A statistically significant relation over three orders of magnitude was found between the outflow and clump masses for the MMAOs, given by log(M_out/M_⊙) = (−0.8 ± 0.3) + (0.8 ± 0.1)log(M_clump/M_⊙). This relation agreed, within uncertainties, with Beuther et al. (2002b), López-Sepulcre et al. (2009) and Sánchez-Monge et al. (2013). Low-mass sources (Bontemps et al. 1996; Wu et al. 2004; Davis et al. 2010; Narayanan et al. 2012) were found to follow a similar trend as the massive sources, expanding the power-law relationship to over six orders of magnitude.

(ii) The relation between the outflow mass-loss rate and clump masses was described by |${\rm log}(\dot{M}_{{\rm out}}/{\mathrm{M}_{{\odot }} \,\mathrm{yr}^{-1}}) = (-5.0 \pm 0.4)+(0.6 \pm 0.1){\rm log}(M{_{\rm clump}}/{\mathrm{M}_{{\odot }}})$|⁠, which held over six magnitudes down to the low-mass regime.

(iii) When plotting outflow masses against dynamical time-scales, 95 per cent of the MMAOs occurred above the L_* = 10² L_⊙ boundary for YSO luminosities, with ∼79 per cent of them between L_* = 10² L_⊙ and L_* = 2 × 10⁴ L_⊙, and ∼16 per cent above L_* = 2 × 10⁴ L_⊙. This suggests that higher luminosity (and hence more massive) YSOs will generate higher mass outflows.

(iv) The relationship between the mechanical force (CO momentum flux) F_CO and the clump masses was described by log(F_CO/M_⊙ km s⁻¹ yr⁻¹) = (−4.5 ± 0.4) + (0.8 ± 0.1)log(M_clump/M_⊙), a trend which held over six orders of magnitude if we included low-mass outflows.

(v) We derived an approximate accretion rate of ∼3 × 10⁻⁴ M_⊙ yr⁻¹ from our mean mass-loss rate. This is of the same order of magnitude than both the mean accretion rates found by Beuther et al. (2002b) and Kim & Kurtz (2006) and agrees with the theoretical rates used in disc accretion models for massive stars (Krumholz et al. 2007; Kuiper et al. 2010).

(vi) An investigation of the relation between the outflow mechanical luminosities and maser luminosities for 38 MMAO targets showed a weak positive correlation wherein bright outflows were associated with brighter masers and clumps. However, low-luminosity outflows and clumps showed no preference for maser brightness. The relation between the 6.7 GHz maser luminosities and clump masses was given by log(L_6.7 GHz/ L_⊙) = ( − 10.6 ± 0.7) + (0.8 ± 0.2)log(M_clump/M_⊙). Although weakly correlated, these relations suggest that there is a correlation between the mass (hence also brightness) of a massive YSO and both the luminosity of the outflow it generates as well as the 6.7 GHz maser it pumps. It agrees with the speculation from Urquhart et al. (2013a) that such a correlation may be related to the fact that the most massive clumps are likely to form more massive stars, and that a higher isotropic maser luminosity is somehow related to a higher stellar luminosity.

Results (i), (ii) and (vi) indicate that both the outflow mass relations and energetics follow a common relation for low- and high-mass YSOs. This lends evidence to the hypothesis that a similar process is responsible for outflows in both mass regimes, i.e. high-mass star formation is a scaled-up version of the low-mass scenario. In addition to this, result (v) shows that the approximate mass accretion rate for MMAOs is sufficiently high to overcome the radiation pressure of a massive central star.

Although this evidence suggests that massive stars form in a similar fashion to low-mass stars, to determine whether each massive YSO indeed drives its own outflow, high-resolution imaging is required.

In a following paper (de Villiers et al., in preparation), the outflow property distributions are compared to those from other surveys, focusing on the evolutionary sequence of 6.7 GHz methanol masers and molecular outflows during the hot core phase.

We want to thank James Caswell for useful comments and contributions to the paper. Timea Csengeri acknowledges financial support for the ERC Advanced Grant GLOSTAR under contract no. 247078. The JCMT is operated by the Joint Astronomy Centre on behalf of the STFC of the United Kingdom, the National Research Council of Canada, and (until 2013 March 31) the Netherlands Organization for Scientific Research. The programme ID for the observations is M07AU20. We also thank the anonymous referee for constructive comments on the paper.

REFERENCES

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

Appendix A

Appendix B (Supplementary Data).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.