© ESO, 2010

1. Introduction

Though they account for a tiny fraction of the galactic stellar population, massive stars (M ≥ 8 M_⊙) dominate the galaxies within which they form and evolve. The radiation and outflows from young massive stars act as feedback processes, altering their parent molecular cloud and thus influencing current and subsequent generations of nearby star formation. Eventually, the ionising fronts of H ii regions and flows of material driven by stellar winds from these stars can lead to the destruction of the molecular cloud itself, ending further star formation altogether. However, our understanding of these effects is still in its infancy, with further study limited by the lack of large well-selected samples from which to derive the general properties of young massive stars and to study how these evolve in terms of both stellar mass and age.

To this end, the Red MSX Source (RMS) survey (Hoare et al. 2005; Urquhart et al. 2008b) has carried out a series of follow-up observations of a sample of candidate massive young stellar objects (MYSOs) which was colour-selected from the MSX point source catalogue (Egan et al. 1999, 2003) by Lumsden et al. (2002). We define an MYSO as a radio-quiet mid-infrared (IR) point source which has reached a luminosity close to its final value, but where accretion is probably still ongoing and an ionised H ii region has yet to form. The mid-IR colours of MYSOs are similar to those of both young and evolved dust-enshrouded objects (such as ultra-compact (UC)H ii regions, low-mass YSOs, evolved stars, proto-planetary nebulae and planetary nebulae (PNe)). Therefore additional observations are required to identify contaminants and gain information about the candidates. The RMS survey has undertaken  ~1′′ radio continuum observations to identify radio loud UCH ii regions and PNe (Urquhart et al. 2007a, 2009),  ~1′′ mid-IR imaging outside the GLIMPSE survey region (Benjamin et al. 2003; Churchwell et al. 2009) to search for candidates which show multiple and/or extended sources at higher resolution (Mottram et al. 2007) and near-IR low-resolution spectroscopy in order to identify evolved stars and weak, unresolved UCH ii regions (e.g. Clarke et al. 2006). We have also undertaken pointed ¹³CO J = 1−0 or J = 2−1 observations to obtain kinematic velocities of our sources in order to calculate kinematic distances (Urquhart et al. 2007b, 2008a).

In a previous paper (Mottram et al. 2010) we presented new measurements of the fluxes of RMS sources in the far-IR (by which we mean 30 μm to 300 μm) using IRAS Galaxy Atlas (IGA, Cao et al. 1997) and Spitzer MIPSGAL (Carey et al. 2009) images. These are important if the luminosities of RMS sources are to be obtained, as the spectral energy distribution (SED) of young massive stars in both the MYSO and UCH ii region phases peak between 60 μm and 120 μm.

In this paper we will use the far-IR data from Mottram et al. (2010) along with other data from RMS observations and the literature to obtain the SEDs of these sources (Sect. 2). In Sect. 3 we discuss tests of the SED model fitter of Robitaille et al. (2007), which are used to inform the choice of input parameters for the fitter and evaluate what represents a good result. The results of the SED fitting and derived bolometric fluxes are presented in Sect. 4, along with the calculation of filter flux ratios and flux estimates for those sources without far-IR fluxes derived from these ratios. In Sect. 5 we present a comparison with other methods for obtaining the bolometric flux from SEDs, then discuss the luminosity distributions of RMS sources in terms of their source type, before we reach our conclusions in Sect. 6.

2. Input data and the model fitter

The SEDs of the candidate MYSOs presented in this paper are primarily based on fluxes from v2.3 of the MSX point source catalogue (PSC, Egan et al. 2003) used to define the RMS sample (see Lumsden et al. 2002) and far-IR fluxes from Mottram et al. (2010). Far-IR fluxes obtained from 18′′ resolution MIPSGAL (Carey et al. 2009) 70 μm data were used where available in preference to those from 1′  ×  1.7′ IRAS Galaxy Atlas (IGA, Cao et al. 1997) data. Sources without full detections in either far-IR data are not included, as there is therefore no data to constrain the far-IR peak of the SED for these objects (the result of fitting the SED without far-IR data will be discussed in Sect. 3.4). While the original coordinates for all sources were those from the MSX PSC, in some cases these have been updated to those from higher-resolution data.

The subsections below discuss additional data included, where possible, to better constrain the SEDs (Sects. 2.1−2.3), and provide distance information (Sect. 2.4).

2.1. Mid-IR imaging

In order to supplement these data, we also include fluxes from our  ~1′′ resolution 10.4 μm TIMMI2 observations (Mottram et al. 2007). The combined 3σ uncertainty in the offset between MSX and TIMMI2 positions for all sources observed with TIMMI2 was found to be 7.6′′. Therefore, where TIMMI2 observations show no source within a radius of 7.6′′ of the source coordinates, a non-detection upper limit flux is used based on the uncertainty in the 10.4 μm TIMMI2 flux.

From TIMMI2 and other high-resolution mid-IR imaging, it is clear that some of the flux detected within the MSX beam (18′′) in star formation regions is due to diffuse emission not directly relating to the source of interest. Therefore, the bolometric flux of larger-beam observations should be apportioned based on the higher-resolution data. In these cases, the apportion ratio was therefore calculated as the ratio of the 10.4 μm flux to the MSX 8 μm flux. Where the TIMMI2 flux is larger than the MSX 8 μm flux, the ratio is set to 1. The individual fluxes from large-beam observations (i.e. MSX, MIPSGAL, IGA, sub-mm and mm) were then multiplied by the apportion ratio before input into the model fitter and the uncertainties of these fluxes calculated as the addition in quadrature of the uncertainty in the apportion ratio and the measurement uncertainties.

In using a simple ratio of the 10.4 μm flux to the MSX 8 μm flux we assume that all the sources detected in the TIMMI2 observations have the same SED shape and that the 9.7 μm silicate feature affects both fluxes equally. While this is an obvious simplification, it is certainly an improvement over not performing this correction and the absence of higher resolution data at longer wavelengths for all our sources makes a more detailed determination impossible at this time. This method was required for 178 out of a total of 1183 sources.

Similarly, some RMS sources have detections at 8 μm in the GLIMPSE PSC, in which case this flux is used to calculate the apportion ratio as for TIMMI2 sources. 182 out of a total of 1183 sources were apportioned in this way. In addition, in some cases it is clear from higher-resolution mid-IR imaging that more than one YSO or H ii region contributes significantly to the flux within the MSX beam. For these sources, the individual components are given letter designations and apportioning carried out where possible. Where we do not have measured fluxes, such as cases which can be seen in the GLIMPSE 8 μm images but do not have 8 μm fluxes in the PSC, it is still clear that each individual source does not contain the total MSX flux. Therefore the remaining fraction of the total flux of the larger-beam observations is split equally between those sources which cannot be apportioned using TIMMI2 or GLIMPSE 8 μm fluxes. This is performed for 136 out of a total of 1183 sources.

2.2. 2MASS

In order to help constrain the near-infrared side of the SED, J, H and K_s band fluxes from the nearest 2MASS PSC (Skrutskie et al. 2006) entry within 7.6′′ (see above) are used. Though this radius may seem large, it is still less than half the full-width half maximum (FWHM) of MSX (18′′) and is motivated by the uncertainty in the source coordinates. In the case of sources with TIMMI2 observations, the 2MASS counterparts identified by Mottram et al. (2007) are used instead. Where there is no 2MASS detection within 7.6′′ of the MSX PSC coordinates, or where inspection of the 2MASS images during object classification indicated that the nearest source is unlikely to be the near-IR counterpart to the MSX target, the noise in each band is calculated in the following way in order to estimate an upper limit flux. All 2MASS PSC entries which are detected in other bands but not in the band in question within 1′ of the target coordinates are identified. An average of the upper limit fluxes in the band in question of all such sources is then calculated. While some 2MASS counterparts may be mis-identified by this method, the SEDs of young massive sources are dominated by emission in the far-IR, so this is unlikely to have an overriding influence on the calculated bolometric flux.

2.3. Sub-millimetre and millimetre continuum

SCUBA 450 μm and 850 μm data from the legacy archive (Di Francesco et al. 2008) are included, where available, for those sources with an offset between the SCUBA observations and the source coordinates  ≤ 14′′, the beam full-width at half maximum (FWHM) of SCUBA at 850 μm. The FWHM is used rather than the pointing error (usually of order ±5′′) as many of these SCUBA sources associated with RMS young massive stars are larger than the beam size. In cases where the offset was larger than 14′′ but smaller than the source size at 850 μm, the SCUBA fluxes were used as upper limits to the true flux at these wavelengths. The percentage errors on the SCUBA 450 μm and 850 μm fluxes are 50% and 20% respectively (Di Francesco et al. 2008).

Millimetre continuum observations at 1.2 mm with the SEST imaging Bolometer array (SIMBA) on the Swedish ESO submillimetre telescope (SEST) by Faúndez et al. (2004), Hill et al. (2005) and Beltrán et al. (2006), along with 1.2 mm observations with the Max-Planck Millimetre Bolometer (MAMBO) on the Institut de Radio Astronomie Millimetrique (IRAM) 30 m telescope by Beuther et al. (2002) were also included. Originally, all sources with offsets from the source coordinates less than or equal to the beam size (24′′ for SIMBA and 11′′ for MAMBO) were included as detections and those with offsets larger than the beam but smaller than the source size as upper limits based on the source flux. However, in some cases the source sizes as derived by Gaussian fits are large compared to the beam size and/or the source is irregular and extended at 1.2 mm, resulting in a much higher flux at this wavelength than is consistent with the rest of the data. All 1.2 mm images were therefore examined and compared to MSX images in order to determine which observations should be regarded as upper limits to the true source flux.

The uncertainties on detected sources are taken as the rms noise from the various data sets, i.e. 40 mJy Faúndez et al. (2004), 40 mJy for Beltrán et al. (2006) and 15 mJy for Beuther et al. (2002). The observations of Hill et al. (2005) were observed over multiple runs, so in order to be conservative the largest rms noise value from these runs (i.e. 150 mJy) was used for all data. For non-detections, or sources with offsets from the MSX position larger than the measured size of the detection at 1.2 mm, 3σ of the appropriate minimum background flux levels are used as an upper limit to the 1.2 mm source flux.

2.4. Kinematic distances

Kinematic distances have been obtained for the majority of RMS sources, calculated from kinematic velocities using the Brand (1993) rotation curve, assuming the galactocentric distance and velocity of the solar system to be 8.5 kpc and 220 km s^-1 respectively (Urquhart et al. 2007b, 2008a). All sources with a galactocentric radius  < 8.5 kpc have a kinematic distance ambiguity caused by the fact that two points along the line of sight, one near and one far, have the same line of sight velocity in the galactic rotation curve, though at the tangent point these are the same and so not an issue. We have solved this distance ambiguity in the northern hemisphere (Urquhart et al. 2010) by cross-matching our sources with a sample of clouds identified by Rathborne et al. (2009) in the ¹³CO Boston University Five College Radio Astronomy Observatory (BU-FCRAO) Galactic Ring Survey (GRS, 18° ≤ l ≤ 55.7°,  | b |  ≤ 1°; Jackson et al. 2006) for which the ambiguity has already been solved by Roman-Duval et al. (2009). In the southern hemisphere we are using data from the Southern Galactic Plane Survey (SGPS, McClure-Griffiths et al. 2005) and the H i self-absorption method (Liszt et al. 1981) to solve the ambiguity to our sources (Urquhart et al., in prep.).

For sources where there are unsolved near-far distance ambiguities, fits are calculated for all distances. These 33 sources are included in the table of results (Table 1) presented in Sect. 4 as the average of these fits, but are not included in the analysis.

3. Fitting the spectral energy distributions of young massive stars

This section discusses how the bolometric flux is obtained from model fits to spectral energy distributions (Sect. 3.1), then describes a test designed to investigate how use of an incorrect distance in the fitting process affects the accuracy of the resulting bolometric flux (Sect. 3.2). The choice of the A_V range that the fitter is allowed is discussed (Sect. 3.3), followed by a study of the relative importance of the various data sets that make up the SED data (Sect. 3.4). The investigations in Sects. 3.2 and 3.4 were run using the data for G034.7569+00.0247 as the data set for this probable MYSO is one of the most complete.

3.1. The model SED fitter

The radiative transfer models of Robitaille et al. (2006) cover a range of different stellar, disk and envelope properties at a range of angles of inclination. The model fitter of Robitaille et al. (2007) then takes these models, convolves them with a user-specified set of filters, allows a range in distance and A_V and then minimises the χ² between the model filter fluxes and the input data. We use this model fitter in preference to other methods as it takes into account dust features, such as the 9.7 μm silicate feature, and consider a more realistic geometry for the circumstellar material than a simple sphere. It is important to mention that a good fit to the data does not prove that the geometry or conditions present in the model are a good description of the observed source, rather that the model is consistent with the observational data. In particular, it may be the case that the data are insufficient to fully constrain the 14 parameters used to construct the models, so it is not appropriate to infer knowledge of variables that the SED data alone do not specifically constrain (Robitaille 2008). de Wit et al. (2010) found that the best fitting models returned by the fitter for MYSOs usually have quite massive dust disks which are ruled out by mid-IR interferometric observations. In addition, due to the non-random nature in which the parameters were sampled care must be taken to establish that observed trends in the data are not a result of trends in the parameter space (see Robitaille 2008, for a full discussion).

The bolometric flux of the fit is, however, well constrained by the fits, as this relates directly to the SED (Robitaille 2008). The following analysis will therefore focus on this property, while comparisons between the bolometric flux obtained by SED fits with other methods will be addressed in Sect. 5.1.

It is important that values derived from the results of the model fitter include several models, in order to obtain a reasonable estimate of whether these results are constrained and unique with the data used. Therefore the number of fits within a Δχ² = 2 of the best fit was calculated (we use the same definition of χ² as Robitaille et al. 2007). If there were less than 10 fits included within ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}\mathrm{+}\mathrm{\Delta }{\mathit{\chi }}^{\mathrm{2}}$ then Δχ² was increased by integer steps until this was no longer the case. This method was used rather than simply finding the best 10 fits so that the spread in fits was not artificially limited if many models are good fits to the data with a range of values for distance, bolometric flux and A_V. A minimum of 10 fits were required in order to ensure that one or two fits did not strongly bias the results.

The model bolometric flux for the source was obtained by using the model luminosity and model distance derived for each individual fit to calculate the bolometric flux for each fit (f_i). The luminosity returned by the model fitter is the combination of that for the central source, accretion onto the central source and accretion in the disk, which is then corrected for foreground extinction consistent with the SED, and is therefore independent of viewing angle unlike the SED itself. Next the weighted mean and weighted standard deviation were calculated of all fits within ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}\mathrm{+}\mathrm{\Delta }{\mathit{\chi }}^{\mathrm{2}}$ to give the flux and flux uncertainty using: $\mathit{F̅}\mathrm{w}\mathrm{=}\frac{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{w}}_{\mathit{i}}{\mathit{f}}_{\mathit{i}}}{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{w}}_{\mathit{i}}}\mathit{,\delta }\mathit{F̅}\mathrm{w}\mathrm{=}\sqrt{\frac{{\mathit{N}}^{\mathrm{\prime }}}{\mathrm{(}{\mathit{N}}^{\mathrm{\prime }}\mathrm{-}\mathrm{1}\mathrm{)}}\frac{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{w}}_{\mathit{i}}\mathrm{(}{\mathit{f}}_{\mathit{i}}\mathrm{-}\mathit{F̅}\mathrm{w}{\mathrm{)}}^{\mathrm{2}}}{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{w}}_{\mathit{i}}}}$(1)where N′ is the number of non-zero weights, and the weights are calculated from the χ² as: ${\mathit{w}}_{\mathit{i}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\left({\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}\mathit{/}\mathrm{2}\right)}\mathit{.}$(2)The χ² values are used to weight the fits so that the best model fitted to the data has more impact on the results than poorer fits, but that the results are not heavily skewed away from the mean. This method is also used to calculate the mean model distance, the mean A_V, and mean filter flux ratios given the results for individual model fits to the data.

3.2. Model distance

A distance range is required by the model fitter within which to fit the observed data. Therefore the test described below was undertaken in order to explore the distance dependence of the bolometric flux obtained from the model fitter, and thus what effect an incorrect distance has on the model derived flux.

Fits were performed to the data for G034.7569+00.0247 for model distances d ± δd kpc, with d varying between 2 kpc and 13 kpc in 1 kpc integer steps. Example SED fits for 2.0 ± 1.0 kpc and 5.0 ± 1.0 kpc (the kinematic distance for G034.7569+00.0247 is 5.0  ±  0.7 kpc, Urquhart et al. 2008a) are shown in Fig. 1 (left and right respectively). For all runs of the model fitter, δd is set to 1.0 kpc as this is the order of the typical error in the kinematic distance for RMS sources.

Fig. 1 Example SED fits with allowed distances of 2.0 ± 1.0 kpc (left) and 5.0 ± 1.0 kpc (right) respectively. The circular points indicate flux detections while the downward pointing triangles indicate upper limit fluxes. The solid line shows the best fit to the data, while the grey lines indicate other fits within a given Δχ2 of the best fit.

Fig. 2 Bolometric flux (left) and AV (right) vs. the mean model distance obtained from SED fits with specified model distances d ± 1.0 kpc with d varying between 2 kpc and 13 kpc for the data of G034.7569+00.0247. The black horizontal error bars indicate the uncertainty on the mean model distance while the dashed horizontal error bars indicate the range of allowed model distances.

The model fitter divides the user-specified distance range logarithmically in order to obtain model distances after which the model SED fluxes are convolved with common filter bands and interpolated to the apertures used to obtain the data. The logarithmic binning was set to 0.01 for this test, while 0.025 was used for the other tests and the results presented in Sect. 4. The results were checked for a few sources at the lower model distance binning step and the differences were found to be minimal. At small model distances the apertures may not encompass the whole model but only the inner region, the calculation of which can require large amounts of computer memory. Therefore, where the distance to the source was less than 2.4 kpc, the distance given to the fitter for the main results was reset to this minimum limit.

The results of the fits where d is varied for the data of G034.7569+00.0247 are shown in Fig. 2, in terms of bolometric flux (left) and A_V (right) vs. the mean model distance as calculated using the weighted mean and weighted standard deviation of the model distances for the individual fits within Δχ² of ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}$. The A_V is relatively constant with distance for these fits (see Fig. 2), suggesting it is constrained by the data. The increase in the bolometric flux and error at small distances is due to the source being so close that the A_V is insufficient to stop emission in the visual being observed (see Fig. 1). The model fitter allows this as the near-IR fluxes are upper limits, and so are treated as weaker constraints on the SED than the detections at longer wavelengths. Overall, though the bolometric flux appears to decrease with increasing distance, the change is of similar order or smaller than the uncertainties in the fluxes, and so the results can be considered distance independent.

3.3. Visual extinction

The visual extinction (A_V) between Earth and the edge of the circumstellar material considered by the model is an important free parameter in fitting model SEDs to observed data. Restricting this value too much could result in poor fits, particularly between the near/mid-IR and far-IR/sub-mm parts of the SED, while allowing too much freedom could lead to unphysical solutions. The extinction law model assumed was calculated by Robitaille et al. (2007) using the method of Kim et al. (1994) for a typical galactic ISM, taking into account the results of Indebetouw et al. (2005) with regards to mid-IR extinction.

We limit the A_V to be within the range 0 to 50, in line with typical values for UCH ii regions (e.g. Hanson et al. 2002), though as shown in Fig. 3 the majority of sources have fitted values below this.

Fig. 3 A histogram of AV derived by the model fitter for all RMS sources with and ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}\mathrm{\le }\mathrm{30}$ for Δχ2 ≤ 10.

Fig. 4 SED fits with different data sets included. Top: all data included. Middle: all data except MIPSGAL. Bottom: all data except sub-mm and mm. Bottom right: all data except 2MASS. The lines and points have the same meaning as in Fig. 1. The bolometric flux derived for each data set is indicated on the plot.

3.4. Data included

SED model fits for G034.7569+00.0247 including only subsets of the available data were undertaken in order to explore the effect of missing data on the bolometric flux results, examples of which, along with the fit using all data are shown in Fig. 4. Comparing the top plot, which includes all data, with the middle plot where the MIPSGAL data has been omitted, the effect of using the model fitter for sources without far-IR data is revealed. The lack of constraint at the peak of the SED leads to the fractional error in the bolometric flux being larger than 1, i.e. it is unconstrained, despite having an otherwise reasonably complete data set. As a result, as mentioned in Sect. 2 we do not fit SEDs for sources where we do not have far-IR flux data, though we will see in Sect. 4.1 that we can use the F₂₁ to F_Bol ratio as derived from fitted sources to estimate F_Bol for these sources. The lower plot of Fig. 4 illustrates that the lack of sub-mm and mm data for many of our sources does not change the results within the calculated errors.

4. Results

The SED model fitter of Robitaille et al. (2007) was run for all 1183 sources with far-IR data as discussed in Sect. 2 using kinematic distances with δd = 1 kpc obtained from molecular emission lines (Urquhart et al. 2007b, 2008a, 2010; Urquhart et al., in prep.) toward each source. There are 33 sources remaining with multiple solutions for the kinematic distance to some sources due to unresolved near/far distance ambiguity, for which the SED fitter was run for each possible distance and the flux and error are given as the mean and $\mathit{\sigma }\mathit{/}\sqrt{\mathit{n}}$ of the results from these runs. The maximum A_V allowed was 50 (see Sect. 3.3).

A sample of the results from these fits are shown in Table 1, while the full version of this table is available online at the CDS. The RMS sample sources have been given one of a number of classifications, shown in Col. 4 of Table 1, which are discussed in Sect. 2.1 of Mottram et al. (2010). “Young/Old Star” is abbreviated to “Y/O Str” and an additional classification for more evolved Compact or diffuse H ii regions, where the ionised region fills much of the MSX beam, is abbreviated to “DH ii”. For details on individual sources see the RMS database (http://www.ast.leeds.ac.uk/RMS). As distance solutions are subject to change as new information becomes available, these are not included in the table. However, up to date distance information is available on the database and once distances are available for all sources, these will be the subject of a future publication. In addition, as discussed in Sect. 3.2 the flux results are distance independent. The data included in the SED of a source is given by a series of flags in Col. 6 of Table 1. These are from left to right 2MASS (Y or N), TIMMI2 (Y or N), whether the far-IR was from the IGA (I) or MIPSGAL (M), submillimetre (Y or N) and millimetre data (N, 1 = Beuther et al. 2002, 2 = Faúndez et al. 2004, 3 = Hill et al. 2005, 4 = Beltrán et al. 2006). Table 1 contains results for 1173 sources, while for a further 10 sources the fitter was run but returned a flux smaller than the error in the flux (i.e. unconstrained).

We do not discuss in this paper the results of the 14 free parameters in the models of Robitaille et al. (2007) as these are not uniquely constrained by the data, and the models do not cover the parameter space uniformly. Indeed, the distributions of most of the parameters follow trends in the models (see Mottram 2008, for more details).

A total of 1069 sources have an unique assigned kinematic distance and have fits with $\mathit{\delta }\overline{{\mathit{F}}_{\mathrm{w}}}\mathit{/}\overline{{\mathit{F}}_{\mathrm{w}}}\mathrm{\le }\mathrm{0.5}$ and ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}\mathrm{\le }\mathrm{30}$ for Δχ² ≤ 10, of which 488 are identified as YSOs, 20 as H ii/YSOs, 45 as Young/Old Star, 76 as CH ii and 440 as H ii regions. Histograms of the number of sources as a function of their luminosity for YSOs and H ii regions are shown in Fig. 5. These sources have typical uncertainties in the bolometric flux of  ~17%, have a median distance of 5.1 kpc, a typical distance uncertainty of  ~19% (assuming δd = 1.0 kpc), and have typical uncertainties in the luminosity of 34%. The sources with MIPSGAL far-IR fluxes do not show any significant difference in terms of distribution to those with IGA data, implying that the improvement in resolution provided by the MIPSGAL data does not result in a lack of luminous sources. Figure 6 shows the plots for YSOs and H ii regions from Fig. 5 for both the assigned distances (black) and with all near/far determinations for all assigned kinematic distances inverted. This is undertaken to examine the effect of incorrect near/far distance ambiguity assignments on the luminosity distributions. As can be see, this does not change the overall shape of the luminosity distributions appreciably. Therefore the remaining individual uncertainties in the distances to the RMS sample do not significantly affect the global results on a statistical level. These results will be discussed in detail in Sect. 5.

Fig. 5 Plots of the number of sources vs. luminosity for RMS candidates with classifications of YSO (488 sources, black) and H ii (440 sources, red). The blue arrows indicate the luminosity of a B0 type star using the surface temperature of Boehm-Vitense (1981) and the luminosity/temperature/mass relationship of Meynet & Maeder (2003).

Fig. 6 Plots of the number of sources vs. luminosity for YSOs (left) and H ii regions (right) with the assigned distances (colours as in Fig. 5) and with all near/far distance determinations swapped (blue and green respectively). The blue arrows have the same meaning as in Fig. 5.

Fig. 7 Left: the distribution of FBol/FMSX   21 as a function of FBol for all YSOs and H ii regions (red and blue respectively in the online version of this figure). Error bars are shown in black in order to give an idea of the total variation in values. Right: the histogram of FBol/FMSX   21 for all sources with MIPS data. A Gaussian is fitted to the data is shown by the dashed line.

Fig. 8 Plots of the number of sources vs. luminosity calculated from MSX 21 μm fluxes for RMS candidates without far-IR data with classifications of YSO (110 sources, left) and H ii region (103 sources, right). The blue arrows have the same meaning as in Fig. 5. The dot-dashed histograms show the normalised luminosity distributions from Fig. 5.

4.1. Sources without far-IR data

Though far-IR fluxes were obtained by Mottram et al. (2010) for the majority of young RMS sources, this was not possible for some sources due to confusion within the available data, which in turn prevented SED fitting of these sources (see Sect. 3.4). However, the properties of those sources where fits were obtained can be used to calculate the mean ratio between the bolometric flux and the MSX 21 μm flux for all sources with good SED fits. This mean ratio can then be used to obtain estimates of the bolometric flux of RMS sources which do not have far-IR data, since all have measured MSX 21 μm fluxes.

In addition to the total flux, filter-band fluxes were output for the SED model fits to the data. The MSX 21 μm band filter fluxes in W m^-2 can be calculated by multiplying the SED MSX 21 μm flux in Janskys by the bandwidth of the filter (4.041  ×  10^-14 Hz, Cohen et al. 2000). The weighted mean ratio of the total bolometric flux to the MSX 21 μm band filter flux can therefore be obtained for each source, in a similar manner to that used for the bolometric flux (see Sect. 3), with the uncertainty given by the weighted standard deviation. The results of these calculations do not appear to be dependent on source flux or type (see left-hand plot of Fig. 7), though the scatter in sources is relatively large and so may mask such relationships. The ratio itself has a roughly log-normal distribution, as shown in the right-hand plot of Fig. 7. Sources with MIPSGAL data generally have a slightly smaller spread than the general distribution, so the properties of the Gaussian fits to the mean ratio for these 613 sources are used.

This mean ratio was therefore used to obtain estimates of the bolometric flux of RMS sources which do not have far-IR data, since all have measured MSX 21 μm fluxes. The bolometric fluxes obtained for sources with TIMMI2 (69) or GLIMPSE PSC (20) fluxes, as well as split sources (19), were apportioned using the method discussed in Sect. 2.1. An example of these results is presented in Table 2, while the full version of this table is available online at the CDS.

Fig. 9 Comparison of bolometric fluxes derived by the model SED fitter with those obtained by simple trapezium rule integration (left) and two component greybody fits (right). The green dashed line indicates y = x. The blue dashed lines indicate equality  ±50%.

Estimate total fluxes were obtained using this method for an additional 280 young RMS sources. The uncertainty in the total flux for each source is a combination in quadrature of the uncertainty in F_Bol/F_MSX   21 and the uncertainty in the MSX 21 μm flux. The fact that the Gaussian fit to F_Bol/F_MSX   21 is performed in log space results in asymmetric uncertainties, which are propagated through to the total flux.

The distributions of luminosities for sources without far-IR data identified as YSO or H ii region with uniquely assigned distances is shown in Fig. 8, with the distributions of SED derived luminosities from Fig. 5 also shown for comparison. Due to the lower number of sources (110 YSOs and 103 H ii regions), the bin size for the distributions for luminosities derived using F_Bol/F_MSX   21 is twice that used for the SED luminosities. The spread of sources is larger than for the SED luminosities, probably caused by the larger uncertainties in the bolometric fluxes derived using the MSX 21 μm flux and the fact that these sources are in regions which are confused at far-IR wavelengths. However the general shape of the distributions is similar, and the peaks in the distributions are near those for the SED luminosities. This method of estimating the far-IR flux therefore seems to provide a reasonable estimate when no far-IR information is available.

4.2. Filter effects

As with the ratio of F_Bol/F_MSX   21 discussed above, various other ratios both of filter fluxes to the bolometric flux and filter fluxes to other filter fluxes can be calculated from the model SED fits. However the SED models of Robitaille et al. (2006) do not include PAH emission features, so the flux ratios for bands where these are important may be different, e.g. the 8.0 μm PAH band may affect ratios involving the MSX 8 μm band. The results of these fits are presented in Table 3.

Fig. 10 Example SED (left) and greybody (right) fits to the same data. For the SEDs, the lines and points have the same meaning as in Fig. 1. The solid line in the greybody images indicates the best two-component greybody fit to the data.

5. Discussion

5.1. Comparison with other methods

In order to explore the relative accuracy of other methods commonly used to obtain bolometric fluxes compared to the SED fitter results, we performed calculations using both simple trapezium rule integration and a combination of two greybody fits (e.g. Minier et al. 2005; Hill et al. 2009) with model rather than observed fluxes. The SED model “observed” flux for each fit was first corrected for the extinction derived by that fit. Next the weighted mean and weighted standard deviation were calculated in order to obtain a corrected model flux with error for each filter where data were input to the fitter for that source. The greybody fit was obtained using two temperature components given by: $\begin{array}{c}\\ {\mathit{F}}_{\mathit{\nu }}\mathrm{=}{\mathit{B}}_{\mathit{\nu }}\mathrm{(}{\mathit{T}}_{\mathrm{1}}\mathrm{)}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }\mathrm{1}}}\right){\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathit{\nu }}\mathrm{(}{\mathit{T}}_{\mathrm{2}}\mathrm{)}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }\mathrm{2}}}\right){\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}\\ \mathrm{where} {\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{h}{\mathit{\nu }}^{\mathrm{3}}}{{\mathit{c}}^{\mathrm{2}}\left({\mathrm{e}}^{\frac{\mathit{h\nu }}{\mathit{kT}}}\mathrm{-}\mathrm{1}\right)} \mathrm{and} {\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{\tau }}_{\mathrm{ref}}{\mathrm{(}\mathit{\nu }\mathit{/}{\mathit{\nu }}_{\mathrm{ref}}\mathrm{)}}^{\mathit{\beta }}\end{array}$(3)where θ is the source size, T is the dust temperature, τ_ν is the optical depth at frequency ν and β is the dust emissivity index such that τ can be evaluated at any frequency relative to a reference frequency ν_ref. A total of five free parameters were used for the fits: the temperatures corresponding to the peaks of the two greybody components at lower and higher temperature (T₁ and T₂), the reference optical depth (τ_ref) at 850 μm, the source size in arcseconds of the high-temperature component (θ₂) and β₁. β₂ was kept constant at 1, while θ₁, the source size of the low-temperature component, was set to the beam size of 1.2 mm SEST data (i.e. 24′′). β₁ was restricted to only be between 1 and 2. Only sources with at least five data points were fitted using this method, due to the number of free parameters, but this resulted in the exclusion of only a few sources which tend to have only upper limit 2MASS fluxes, an upper limit in the MSX 12 μm band and MIPSGAL data. Though these sources have less than 5 detected data points, the SED fitter is also constrained by the upper limits.

Figure 9 shows a comparison between bolometric fluxes derived by the model SED fitter with those obtained by simple trapezium rule integration (left) and two component greybody fits (right). There is quite a large scatter, with a significant fraction more than 50% away from the SED fitter derived flux. Both alternative methods are more likely to over-estimate the flux, with mean values of F_grey/F_fitter = 1.35 and F_trap/F_fitter = 1.16. If the Emerson flux (Emerson 1988) is calculated with IRAS fluxes from the models instead of a full trapezium integration, the mean value of F_Emerson/F_fitter = 1.09 though this does not include the effect of contamination which is often a major factor for IRAS fluxes of sources in the galactic plane. If IRAS PSC fluxes are used, F_Emerson/F_fitter = 1.82 with considerably more scatter. Many of the greybody fits do not fit both the near and mid-IR well, as shown by the example in Fig. 10. The group of sources with considerably larger greybody fluxes in Fig. 10 are all sources with MIPS far-IR data but no submm or mm points where the fit has the far-IR peak at λ > 70 μm, and so overestimates the flux from cold dust. Neither method will deal with the 9.7 μm silicate feature accurately. Figure 11 shows a histogram of the fractional errors of fluxes derived by the SED fitter, greybody fits and trapezium integration. These indicate that errors in fluxes derived using trapezium integration are systematically larger than those from obtained using the other methods, and that the errors using greybody fits are on average larger than those using the SED fitter. The median fractional error is 0.17 for the SED fitter fluxes, 0.24 for the two-component greybody fluxes and 0.59 for trapezium integration. Overall, therefore, the SED fitter produces more consistent and accurate results than either greybody fitting of simple integration.

Fig. 11 Histogram of the fractional errors (i.e. δFBol/FBol) for fluxes obtained using the SED fitter (black), trapezium integration (red) and greybody fitting (blue).

Fig. 12 SED (left) and MIPSGAL 70 μm image (right) for the H ii region G045.4543+00.0600, for which we derive a luminosity of (3.7 ± 1.0) × 105 L⊙ while Wood & Churchwell (1989) derive a luminosity, after correction to our distance of 7.3 kpc, of 8.2 × 105 L⊙ using IRAS far-IR data. The RMS source is merged with the source to the east (IRAS 19120+1103) in IRAS data and has a PSC 60 μm flux of 5340 Jy compared to a MIPS 70 μm flux of 2402 Jy, accounting for much of the remaining difference in luminosity measurements. Saturated pixels are indicated in yellow in the online version of this figure. For the SED, the lines and points have the same meaning as in Fig. 1.

5.2. YSOs

The upper left plot of Fig. 5 shows that the number of YSO sources peaks around 5 × 10³ L_⊙, and decreases sharply above  ~2 × 10⁴ L_⊙. The gradual drop in the number of YSOs at luminosities below the peak is caused by incompleteness in the RMS sample due to the 21 μm flux-limit of the MSX survey. Though there are some sources which have L ≳ 10⁵ L_⊙ (which corresponds to an  ~30 M_⊙ O7 type star using the temperature/logQ scale of Martins et al. 2005, corrected for the luminosity/temperature/mass relationship of Meynet & Maeder 2003), inspection reveals that many of these sources should be apportioned based on GLIMPSE images, and thus have lower luminosities, but lack PSC detections at 8 μm. As already discussed in Sect. 4, the shape of the overall distribution is robust against these issues.

5.3. H ii regions

The left-hand plot of Fig. 5 shows that the number of H ii regions in the RMS sample peaks around 3 × 10⁴ L_⊙. There are few sources which are identified as H ii regions but have a luminosity derived from the SED lower than the limit for formation of an H ii region (L ~ 10³ L_⊙, spectral type  ~B3, Boehm-Vitense 1981; Meynet & Maeder 2003). Inspection of these sources suggest that they are correctly identified as H ii regions, but suffer from individual distance determination issues.

While the models used by the SED fitter do not include additional sources of heating of the dust inside the ionised zone, in particular Lyman alpha (e.g. Hoare et al. 1991), the SEDs of UCH ii regions are similar to those of YSOs except at radio wavelengths, which the models do not probe. As the primary concern of this work is the bolometric flux, the fits are acceptable for this purpose.

The most luminous UCH ii region observed by Wood & Churchwell (1989) has a luminosity of  ~3 × 10⁶ L_⊙, while no RMS source has a luminosities above 10⁶ L_⊙. However, some of the Wood & Churchwell (1989) sources are now known to contain multiple objects and in some cases the RMS kinematic distance measurements are significantly less than the ones used by those authors. For example Wood & Churchwell (1989) use d = 9.7 kpc for G045.4543+00.0600 but the RMS kinematic distance to this source is 7.3 kpc, which would change their 1.44 × 10⁶ L_⊙ to 8.2 × 10⁵ L_⊙. The SED fit for this source gives L = (3.7 ± 1.0) × 10⁵ L_⊙ but MIPSGAL rather than IRAS far-IR fluxes are used and the MIPSGAL 70 μm image shows multiple bright regions nearby (see Fig. 12). The SED does not fit the 850 μm point well, but this source exhibits significant surrounding low-level emission in the MIPSGAL image and strong radio emission. Both of these factors could increase the observed SCUBA flux above that consistent with the SED models based on the lower wavelength emission, which is well fit.

5.4. Comparison of YSOs and H ii regions

In comparing the luminosity distributions of YSOs and H ii regions (Fig. 5), the difference in the peak of these distributions is apparent. Indeed, a Kolmogorov-Smirnov (KS) test of the cumulative distribution functions, indicates that the probability of these distributions being the same is of the order 10^-26. It is unlikely that this is due to an under-detection of B spectral type H ii regions in the sample as the RMS radio observations (Urquhart et al. 2007a, 2009) were deep enough to detect H ii regions down to spectral type B1 to distances of 20 kpc assuming an electron temperature of 10⁴ K. In addition, many of these MYSOs are visible at near-IR wavelengths, so are unlikely to be heavily embedded hyper-compact (HC) H ii regions.

Main-sequence stars with a spectral type of B3 (L ≥ 10³ L_⊙) or earlier are capable of producing H ii regions, and evolutionary tracks indicate that core hydrogen burning should have started in YSOs of spectral type B3 or earlier despite the fact that accretion is still ongoing (e.g. Yorke & Sonnhalter 2002). It has therefore been unclear until recently why MYSOs, where the central source is massive enough to form an H ii region but no H ii region is observed, are detected at all, rather than all young massive stars transitioning directly from intermediate mass YSOs to UCH ii regions as they gain mass (Hoare & Franco 2007). Some mechanism must either quench or mask the existence of an H ii region (e.g. Walmsley 1995; Keto 2003), or increase the radius of the star such that the stellar T_eff, and thus UV flux, are decreased (e.g. Yorke & Bodenheimer 2008; Hosokawa & Omukai 2009; Hosokawa et al. 2010).

Ascertaining further how selection effects influence the difference in the distribution of RMS YSOs and UCH ii regions will require the luminosity functions of these sources, which will be the subject of a future paper (Mottram et al., in prep.).

6. Conclusions

After preforming tests of the influence of the input parameters to the model SED fitter of Robitaille et al. (2007), we have obtained the luminosities of 1173 young RMS sources, of which 1069 have unique kinematic distances and good fits.

The luminosity distributions of these sources in terms of their RMS classification have been presented, and show that there are few MYSOs with L ≥ 10⁵ L_⊙, though we detect UCH ii regions up to  ~7 × 10⁵ L_⊙, a difference which we find to be statistically significant. We also find that using the SED fitter to obtain bolometric flux measurements is more consistent and has less scatter, than either using simple trapezium rule integration or greybody fits for mid-IR bright young massive stars. Finally, we obtained the flux ratios consistent with our SED fits, which allow us to estimate the bolometric flux, of 280 sources for which far-IR fluxes could not be obtained.

Having obtained the luminosity distributions of YSOs and H ii regions in the RMS sample, we can go on to calculate the luminosity functions for these sources. This will be discussed in a forthcoming paper.

Acknowledgments

The authors would like to thank the referee, Barbara Whitney, for her helpful comments and suggestions which improved the quality and clarity of this paper. J.C.M. is partially funded by a Postgraduate Studentship and by a Postdoctoral Research Associate grant from the Science and Technologies Research Council of the United Kingdom (STFC). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This paper made use of information from the Red MSX Source survey database at www.ast.leeds.ac.uk/RMS which was constructed with support from the Science and Technology Facilities Council of the UK.

References

All Tables

All Figures

	Fig. 1 Example SED fits with allowed distances of 2.0 ± 1.0 kpc (left) and 5.0 ± 1.0 kpc (right) respectively. The circular points indicate flux detections while the downward pointing triangles indicate upper limit fluxes. The solid line shows the best fit to the data, while the grey lines indicate other fits within a given Δχ2 of the best fit.
In the text

	Fig. 2 Bolometric flux (left) and AV (right) vs. the mean model distance obtained from SED fits with specified model distances d ± 1.0 kpc with d varying between 2 kpc and 13 kpc for the data of G034.7569+00.0247. The black horizontal error bars indicate the uncertainty on the mean model distance while the dashed horizontal error bars indicate the range of allowed model distances.
In the text

	Fig. 3 A histogram of AV derived by the model fitter for all RMS sources with and ${\mathit{\chi }}_{\mathrm{best}}^{\mathrm{2}}\mathrm{\le }\mathrm{30}$ for Δχ2 ≤ 10.
In the text

	Fig. 4 SED fits with different data sets included. Top: all data included. Middle: all data except MIPSGAL. Bottom: all data except sub-mm and mm. Bottom right: all data except 2MASS. The lines and points have the same meaning as in Fig. 1. The bolometric flux derived for each data set is indicated on the plot.
In the text

	Fig. 5 Plots of the number of sources vs. luminosity for RMS candidates with classifications of YSO (488 sources, black) and H ii (440 sources, red). The blue arrows indicate the luminosity of a B0 type star using the surface temperature of Boehm-Vitense (1981) and the luminosity/temperature/mass relationship of Meynet & Maeder (2003).
In the text

	Fig. 6 Plots of the number of sources vs. luminosity for YSOs (left) and H ii regions (right) with the assigned distances (colours as in Fig. 5) and with all near/far distance determinations swapped (blue and green respectively). The blue arrows have the same meaning as in Fig. 5.
In the text

	Fig. 7 Left: the distribution of FBol/FMSX   21 as a function of FBol for all YSOs and H ii regions (red and blue respectively in the online version of this figure). Error bars are shown in black in order to give an idea of the total variation in values. Right: the histogram of FBol/FMSX   21 for all sources with MIPS data. A Gaussian is fitted to the data is shown by the dashed line.
In the text

	Fig. 8 Plots of the number of sources vs. luminosity calculated from MSX 21 μm fluxes for RMS candidates without far-IR data with classifications of YSO (110 sources, left) and H ii region (103 sources, right). The blue arrows have the same meaning as in Fig. 5. The dot-dashed histograms show the normalised luminosity distributions from Fig. 5.
In the text

	Fig. 9 Comparison of bolometric fluxes derived by the model SED fitter with those obtained by simple trapezium rule integration (left) and two component greybody fits (right). The green dashed line indicates y = x. The blue dashed lines indicate equality  ±50%.
In the text

	Fig. 10 Example SED (left) and greybody (right) fits to the same data. For the SEDs, the lines and points have the same meaning as in Fig. 1. The solid line in the greybody images indicates the best two-component greybody fit to the data.
In the text

	Fig. 11 Histogram of the fractional errors (i.e. δFBol/FBol) for fluxes obtained using the SED fitter (black), trapezium integration (red) and greybody fitting (blue).
In the text

Fig. 12 SED (left) and MIPSGAL 70 μm image (right) for the H ii region G045.4543+00.0600, for which we derive a luminosity of (3.7 ± 1.0) × 105 L⊙ while Wood & Churchwell (1989) derive a luminosity, after correction to our distance of 7.3 kpc, of 8.2 × 105 L⊙ using IRAS far-IR data. The RMS source is merged with the source to the east (IRAS 19120+1103) in IRAS data and has a PSC 60 μm flux of 5340 Jy compared to a MIPS 70 μm flux of 2402 Jy, accounting for much of the remaining difference in luminosity measurements. Saturated pixels are indicated in yellow in the online version of this figure. For the SED, the lines and points have the same meaning as in Fig. 1.

In the text