Contents

• The evolution of the spectral energy distribution in massive young stellar objects
• 1 Introduction
• 2 The source sample
• 3 The approach
• 4 Results
• 5 An evolutionary sequence for massive YSOs
• 6 Discussion
• 7 Summary
• References
• 8 Online Material
• Appendix A: Notes and data analysis log for individual objects

A&A 481, 345-365 (2008)
DOI: 10.1051/0004-6361:20078661

The evolution of the spectral energy distribution in massive young stellar objects

S. Molinari¹ - S. Pezzuto¹ - R. Cesaroni² - J. Brand³ - F. Faustini¹ - L. Testi^2,4

1 - INAF-Istituto Fisica Spazio Interplanetario, Via Fosso del Cavaliere 100, 00133 Roma, Italy
2 - INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
3 - INAF - Istituto di Radioastronomia, Via Gobetti 101, 40129 Bologna, Italy
4 - ESO Headquarters, Karl-Schwarzschild-Strasse 2, 85748, Germany

Received 12 September 2007 / Accepted 26 December 2007

Abstract
Context. The mechanism of formation of massive stars is still a matter of debate. It is not yet clear if it can be considered to be a scaled-up analogue of the low-mass star regime, or if there are additional agents like merging of lower-mass forming objects or accretion from initially unbound material. Most of the uncertainties come from the lack of diagnostic tools to evolutionarily classify large samples of candidate massive protostellar objects that can then be studied in more detail.
Aims. We want to verify whether diagnostic tools like the SED shape and the relationship between envelope mass and bolometric luminosity can be extended to the study of high-mass star formation.
Methods. The 8-1200 $\mu$m SED of YSOs in 42 regions of massive star formation has been reconstructed using MSX, IRAS, and submm data partly available from previous works. Apart from IRAS catalogue fluxes, the fluxes in the Mid-IR and sub-mm/mm were derived directly from the images. The SEDs were fitted to an extensive grid of envelope models with embedded ZAMS stars, available from the literature. Sources that could not be fitted with a single model were then fitted with a two-component model composed of an embedded ZAMS for the mid-IR part and a single-temperature optically thin greybody for the longer wavelength emitting component. Sources were classified as ``IR'' if they were fitted with an embedded ZAMS envelope, and ``MM'' if they could only be fitted with a greybody with a peak at high $\lambda$; further subclassification was based on being the most massive object in the field (``P'', for primary) or not (``S'', for secondary).
Results. The different classes of sources identified in our analysis have very different SEDs and occupy distinct areas in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram; by analogy with the low-mass regime, we see that MM-P, IR-P and IR-S objects could be interpreted as the high-mass analogue of Class 0-I-II. Evolutionary tracks obtained from a simple model based on the turbulent core prescriptions show that the three classes of sources possibly mark different periods in the formation of a massive YSO. The IR-P objects are consistent with being at the end of the main accretion phase, while MM-P sources are probably in an earlier evolutionary stage. The timescales for the formation decrease from $\sim$ $4\times 10^5$ to $\sim$ $1\times 10^5$ years with stellar mass increasing from $\sim$6 to $\sim$40 $M_{\odot }$; these timescales, and the association with young clusters with median stellar age of a few 10⁶ years suggest that the most massive objects are among the last ones to form.
Conclusions. Our results are consistent with the high-mass star formation being a scaled-up analogue of the traditional accretion-dominated paradigm valid for the low-mass regime.

Key words: stars: formation - stars: pre-main sequence

   
1 Introduction

The widely accepted paradigm for the formation of solar-type stars via spherical accretion (Shu et al. 1987) predicts an evolution from cores to protostars, and finally pre-main sequence stars that is well matched with distinctive characteristics of their IR-mm SED (Lada & Wilking 1984). The empirical classification of the SED of low mass YSOs has thus been used as a powerful tool to constrain theoretical models. In particular, is was shown that the peak of the SED shifted from wavelengths longward of 100 $\mu$m for Class 0 objects toward shorter wavelengths the more the YSO approaches the Main Sequence (MS). An additional tool that proved to have good diagnostic power is the relationship between the bolometric luminosity of the YSO and the envelope mass out of which it forms. Objects known to be in different stages according to the SED shape classification also occupy distinct regions in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram, and the entire Class 0-II path can easily be followed even with very simple models (Saraceno et al. 1996).

Figure 1: A view of two regions prototypical for distinct evolutionary phases of massive star formation: IRAS 23385+6053 in panel  a) ( left) for massive Class 0s, and IRAS 20126+4104 in panel  b) ( right) for Hot Cores. In both panels the color image (grey scale in the printed version) dipicts the 21 $\mu$m continuum observed with MSX in band E, the green contours (white in the pr. ver.) represent the cold dust emission as traced in the sub-millimeter continuum, and the blue contours (black in the pr. ver.) are the centimetric continuum as observed at the VLA.

Establishing a similar scenario for intermediate and high-mass YSOs is much more difficult due to the clustered environments in which they are born and to the large ($d\geq 1$ kpc) characteristic distances. Qualitative schemes have been proposed (Evans et al. 2002); a more quantitative assessment especially in terms of lifetimes of the duration of the different phases is still elusive. IRAS data have been mainly used to select samples of massive YSOs using color criteria, but such investigations can only go as far as providing statistical indications about the evolutionary stage that different groups of sources exhibiting specific global SEDs properties may be in. The complexity of the situation clearly emerges when these bright IRAS sources are observed at higher resolution in the mid-IR, and are often resolved in clusters of several members with a diffuse emission component local to the star forming region. The same region observed in the sub-millimeter even with the limited spatial resolution attainable with single-dish facilities may or may not offer a scenario consistent with the mid-IR picture. Two textbook examples in this respect, IRAS 20126+4104 and IRAS 23385+6053, have been the subject of extensive studies which clarified their nature as very young massive forming objects yet in distinct evolutionary stages.

IRAS 23385+6053 is a luminous ($L\sim10^4$ $L_{\odot}$) source originally proposed as a template for massive Class 0 protostars (Molinari et al. 1998b), where molecular line properties suggest an evolutionary stage prior to HC (Thompson & Macdonald 2003; Fontani et al. 2004). A strong dust continuum peak traces the presence of a compact massive object, driving a molecular outflow. This compact object is completely absent in the mid-IR shortward of 15 $\mu$m. Figure 1a clearly shows that the structured mid-IR emission is completely anti-correlated with the peak sub-mm emission and is partially resolved in several compact structures likely marking additional signposts for YSOs. The presence of thermal free-free emission (Molinari et al. 2002) and of a $K_{\rm s}$-detected pre-MS stellar population (Faustini et al., in prep.) only in association with the mid-IR structures and not with the central mm core, suggests for the latter a very early stage of formation, possibly prior to the ignition of nuclear H burning.

IRAS 20126+4104 (Cesaroni et al. 1997, 1999) is a similarly luminous Hot Core (HC) exhibiting a strong molecular and dust continuum peak, where temperatures of $\sim$200 K are revealed in its inner envelope using high-density tracers like CH₃CN. It exhibits thermal free-free radio emission, probably arising from a jet. A disk-envelope structure which collimates a well-developed bipolar molecular outflow is clearly detected in mm-wavelength interferometric observations. The HC is also revealed in the mid-IR as a strong compact peak (see Fig. 1b), and is coincident with a conspicuous cluster of young low mass YSOs revealed in $K_{\rm s}$-band observations.

It is clearly possible that given the large average distance of massive YSOs in general (few kpc on average) the dominant millimeter peak may itself be a compact cluster of sources which is unresolved even with interferometric facilities; the mean separation between the trapezium stars in the inner ONC is $\sim$0.01 pc (or 5 $^{\prime \prime}$ at a distance of 500 pc), at the limit of the currently available sub-mm arrays for distances in excess of 2 kpc. Multi-wavelength comparative analysis carried out with present state-of-the art facilities is therefore not able to address the evolutionary stage of single massive YSOs (unless for very nearby objects), but it would be sensitive to evolutionary differences of small groups of massive YSOs within arcmin-sized star forming regions. Our first aim is therefore to verify in a large sample of sources if this approach has the potential to improve our understanding of massive star formation.

A detailed SED characterization of massive YSOs within individual star forming regions is also mandatory for a careful estimate of the bolometric luminosity, which for low-mass YSOs is the principal discriminant between Class 0 and Class I sources. The presently available data at the peak of the SEDs are the IRAS 100 $\mu$m fluxes, but the resolution of 5$^{\prime}$ is insufficient for the purpose; the structure present in the field of IRAS 23385+6053 as reported in Fig. 1a, for example, is entirely encompassed by one 100 $\mu$m IRAS beam. Our second aim is therefore to devise an approach to split the bolometric luminosity of a group of bright YSOs between its main members in order to provide a more firm placement of massive YSOs in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram.

   
2 The source sample

We selected 42 sources from the original sample of Molinari et al. (1996) and from the complementary sample collected (Fontani et al. 2005) for the portion of the southern sky ( $\delta \leq30$$^{\circ}$) that was left out in the Molinari et al. (1996) work. The samples contain sources grouped into High and Low sources mainly according to their [25-12 $]\;\lower.6ex\hbox{$<$ }\kern-7.75pt\raise.65ex\hbox {$>$ }\;0$.57 IRAS color value. The adopted threshold was recommended by Wood & Churchwell (1989) to select sources with a high probability of association with UCHII regions, and these are the High sources in our terminology. A robust case has been built in the past to show that the two groups of sources exhibit different properties; the Low sources appear to be less strongly characterised than High  but they show convincing evidence of being on average relatively younger (Molinari et al. 1996, 1998a, 2000, 2002; Brand et al. 2001). The present source list contains 16 High and 26 Low sources and is reported in Table 1.

 

 

Table 1: List of sources.

IRAS field	$\alpha$(J2000)	$\delta$(J2000)	[25-12]	$D_{\rm kin}$ (kpc)
05137+3919	05:17:13.3	+39:22:14.0	0.63	11.5
05345+3157	05:37:47.8	+31:59:24.0	0.19	2.1
05373+2349	05:40:24.4	+23:50:54.0	0.47	1.6
08211-4158	08:22:52.1	-42:07:53.7	0.46	1.7
08470-4243	08:48:47.7	-42:54:24.0	0.82	2.2
08477-4359	08:49:32.9	-44:10:47.0	0.47	1.8
08563-4225	08:58:12.5	-42:37:34.0	0.53	1.6
08589-4714	09:00:40.4	-47:25:54.0	0.81	1.5
09131-4723	09:14:55.5	-47:36:13.0	0.55	1.7
10019-5712	10:03:40.4	-57:26:39.0	0.86	1.8
10095-5843	10:11:15.8	-58:58:15.0	0.47	1.1-2.9
10184-5748	10:20:14.9	-58:03:39.0	1.04	5.4
10501-5556	10:52:10.8	-56:12:28.0	0.69	2.5
10555-6242	10:57:34.3	-62:59:03.2	0.92	3.0
11380-6311	11:40:27.6	-63:27:56.0	0.56	1.3-6.0
12063-6259	12:09:01.2	-63:15:56.0	0.91	10.1
12127-6244	12:15:24.0	-63:01:20.0	0.97	10.8
14000-6104	14:03:38.0	-61:17:56.0	0.54	5.6
14057-6032	14:09:23.7	-60:46:56.0	1.03	0.8-10.6
15068-5733	15:10:43.3	-57:44:47.0	1.30	4.2-8.9
15454-5335	15:49:16.3	-53:44:59.0	1.15	5.2-9.1
15519-5430	15:55:50.4	-54:38:58.0	0.35	2.6-11.7
15557-5337	15:59:38.2	-53:45:32.0	0.52	3.3-11.2
16148-5011	16:18:35.2	-50:18:53.0	0.53	3.3-11.9
16232-4917	16:27:03.3	-49:23:59.0	0.37	3.4-11.9
16501-4314	16:53:41.1	-43:19:23.0	0.55	6.7-9.5
16535-4300	16:57:05.5	-43:05:20.0	0.35	6.8-9.4
17082-4114	17:11:46.2	-41:18:03.0	0.43	2.5-14.0
17425-3017	17:45:45.1	-30:18:51.0	0.24	8
18024-2119	18:05:25.4	-21:19:40.0	0.24	0.12-16.6
18144-1723	18:17:24.2	-17:22:13.0	0.56	4.3-12.1
18278-1009	18:30:35.2	-10:07:12.0	0.52	5.7
18511+0146	18:53:38.1	+01:50:27.0	0.45	3.9
19092+0841	19:11:37.4	+08:46:30.0	0.51	4.5-7.8
20062+3550	20:08:09.8	+35:59:20.0	1.42	4.9
20106+3545	20:12:31.3	+35:54:46.0	0.47	1.6-3.2
21307+5049	21:32:31.5	+51:02:22.0	0.51	6.2
21519+5613	21:53:38.8	+56:27:53.0	1.04	7.3
22172+5549	22:19:09.0	+56:04:45.0	0.32	5.0
22305+5803	22:32:24.3	+58:18:58.2	1.05	5.6
22506+5944	22:52:38.6	+60:00:56.0	0.73	5.4
23385+6053	23:40:53.3	+61:10:19.1	0.54	5.0

The table reports the IRAS name, the J2000 coordinates, the [25-12] color and the distance in kpc. The value of the color is reported in normal typeface for Low sources (where this value is below 0.57), and in boldface in High (where it is higher that this threshold). In several inner Galaxy sources the distance ambiguity could be tentatively solved using a variety of methods, but we chose here to leave this open because in the analysis we perform the distance will be a fitting parameter (see the rest of the paper, and in particular Sect. 3.3).

   
3 The approach

We outline below the general method, and in detail the various steps in the following subsections. We want to identify YSOs in the target fields and estimate their fluxes both in the mid-IR and in the millimeter, where the spatial resolution of the data is relatively higher, trying to isolate as much as possible the contribution of the YSOs out of the total that may be due to other nearby objects or extended components. The bulk of the SED is however emitted in the far-IR, where the resolution is worse. Far-IR emission is mostly due to relatively cold dust which is also traced by the millimeter and submillimeter emission, rather than in the mid-IR; so the information gathered from the millimeter maps is used to formulate assumptions on the amount of IRAS 60 and 100 $\mu$m flux that is to be assigned to the various YSOs.

Compact peaks identified in the mid-IR and in the millimeter often do not coincide, and the question then arises whether they are counterparts of the same object, as in the case of IRAS 20126+4104 (Fig. 1 right), or different objects, as in IRAS 23385+6053 (Fig. 1 left). We made generous allowances for the inaccuracies of the astrometry for both the mid-IR and millimeter maps, and we considered as non-associated those mid-IR/mm counterparts more than 30 $^{\prime \prime}$ away from each other; where the separation is less than that value, they were assumed to be associated at first. At this point a first tentative SED is build for each identified source.

SEDs were then compared to models of envelopes surrounding ZAMS stars spanning a wide parameter space in terms of spectral classes of the embedded star, the density and geometry of the circumstellar envelopes. Sources were then divided into two groups, those where models were found providing a good fit to the data across the entire wavelength range from the mid-IR to the millimeter as in IRAS 20126+4104, and those where this did not happen, as in IRAS 23385+6053. In the latter case we assumed that the mid-IR and millimeter counterparts were in fact not associated, and the SED was then fitted with a two-component model consisting of the sum of an embedded ZAMS model (as for the first group) for the mid-IR part, and of a modified greybody for the far-IR/mm part of the SED.

Sections 3.1, 3.2 and 3.3 below explain in detail the various steps of the analysis procedure.

   
3.1 Building the source SEDs

   
3.1.1 The millimeter

Sub-mm/mm continuum maps are available for all fields in Table 1. Large maps covering several square arcminutes at 1.2 mm ( ${\it HPBW}=24$ $^{\prime \prime}$) were obtained at the SEST telescope with SIMBA by Beltrán et al. (2006) for all fields with $08~{\rm h}\leq\alpha({\rm J}2000)\leq18$ h. For all other fields small ( ${\it FOV}\sim2$$^{\prime}$²) 850 $\mu$m (and occasionally 450 $\mu$m) jiggle-maps were acquired in service observing at the JCMT telescope using the SCUBA instrument ( ${\it HPBW}\sim 15$ $^{\prime \prime}$ @850 $\mu$m). In each field the peaks were visually identified and labelled with numbers; the brightness distribution in the observed fields was analyzed using cuts in specific directions to verify the presence of changes in the functional form of the brightness profile that would be caused by a core+halo structure. This allowed us to identify and quantify the level, if any, of underlying diffuse emission, and to provide reliable estimates of the FWHM of the identified peaks. This information was then used with the AIPS JMFIT routine to fit the brightness profiles of the compact components to 2D Gaussian functions (simultaneously, if identified peaks were blended) and estimate the integrated flux. In all cases we verified a posteriori that the sum of the integrated fluxes of all sources plus the level of diffuse emission was compatible with the total integrated flux over the entire area encompassing the identified sources. We report in Table 2 for each field the computed fluxes for all identified components. The location of the mm peaks is indicated as offsets in ( $\alpha,\delta$) with respect to the IRAS positions reported in Table 1.

This method to estimate millimeter fluxes differs from that used in Beltrán et al. (2006); in that context the aim was to quantify the flux (and hence the mass) of the entire clumps associated with IRAS sources and dense gas, so they chose to adopt the CLUMPFIND algorithm in its 2D version. This algorithm performs an integration of the contiguous emission area above a certain flux threshold and calls that a ``clump''. The aim of the present work, however, is to determine the flux quantity that can be assigned to the massive forming object only, discarding what is more likely from the diffuse clump medium. This is the reason why our flux estimates are always between 30% and 100% of the Beltrán et al. ones; if we include also the contribution of the extended halos wherever they are identified, we recover the same fluxes.

Some of the cores in our source list were also observed, and most of them detected, at 3.4 mm with the OVRO interferometer (Molinari et al. 2002); the integrated fluxes are considered here in the millimeter SED to be model-fitted.

   

Table 2: Sub-mm and mm fluxes.

IRAS field	MM	$\Delta\alpha$	$\Delta\delta$	F450	F850	F1200
		( $^{\prime \prime}$)	( $^{\prime \prime}$)	(Jy)
05137+3919	1	+8	+7		1.15
	2	-7	-22		0.33
05345+3157	1	+69	+40		1.95
	2	+74	+12		1.00
	3	+37	-11		0.5
05373+2349	1	+7	+1		3.27
08211-4158	1	+7	-11			1.45
08470-4243	1	+1	-11			4.2
08477-4359	1	+29	-71			1.4
08563-4225	1	-6	-10			1.6
08589-4714	1	-3	-17			0.93
09131-4723	1	-9	+5			1.5
	2	-44	+16			1.1
10019-5712	1	+3	-6			1.15
	2	-23	+11			0.47
10095-5843	1	0	-10			0.94
10184-5748	1	+17	-23			3.2
	2	+13	-19			0.2
10501-5556	1	+12	-14			1.9
10555-6242	1	-8	-5			1.7
11380-6311	1	+24	-3			1.7
12063-6259	1	-3	-11			3.3
12127-6244	1	-4	-12			5.5
	2	+101	-36			0.7
14000-6104	1	-18	-26			3.27
	2	-80	+1			1.2
	3	-30	+191			1.25
	4	-11	+136			1.15
	5	+215	-111			0.65
	6	+275	-60			1.9
	7	-216	+150			0.15
	8	-233	-55			1.0
14057-6032	1	+13	-13			0.93
15068-5733	1	+12	-13			2.35
15454-5335	1	+30	-20			6.5
15519-5430	1	-20	-20			3.9
15557-5337	1	+3	+7			6.5
	2	-13	-26			8.7
16148-5011	1	+36	-10			1.7
	2	+72	-300			0.6
16232-4917	1	-2	-1			0.8
16501-4314	1	+16	-8			0.83
	2	-15	0			0.76
16535-4300	1	+16	+3			0.91
17082-4114	1	-27	-22			2.8
	2	-10	+30			1.6
	3	-61	+6			1.8
17425-3017	1	-3	+12			1.0
18024-2119	1	+7	+15	31.2	3.9
18144-1723	1	0	+1	28.3	7.6
18278-1009	1	+10	+4	12.3	1.6
18511+0146	1	0	+3	53.4	3.6
	2	0	-15	7.7	1.5
19092+0841	1	+30	+1	44.6	5.4
20062+3550	1	+16	-1		3.1
20106+3545	1	+20	+2	5.55	0.84
	2	+11	+10	6.5	0.89
21307+5049	1	-2	-3		1.5
21519+5613	1	+9	-2		1.2
22172+5549	1	-1	+14		1.9
22305+5803	1	-5	+4			0.44a
22506+5944	1	+1	-10	17.5	3.5
	2	+1	-11	12.9	b
	3	+31	-23	4.0	0.35
23385+6053	1	+12	+10		1.8

$\parbox{7.2cm}{ {$^a$\space Obtained at IRAM-30~m (${\it HPBW}\sim10$$^{\prime \prime}$ ).} \\ {$^b$\space Blended with MM1 @850~$\mu$ m.}}$

   
3.1.2 The mid-infrared

For the mid-IR wavelength range we used the MSX Galactic Plane imaging survey (Price et al. 2001) carried out in 4 bands centered at 8.3 (band A), 12.1 (band C), 14.7 (band D) and 21.4 $\mu$m (band E). Images at all wavelengths are publicly available and have been reconstructed with a beam of 20 $^{\prime \prime}$. In all examined fields the 8.3 $\mu$m image is richer in sources than it appears at the other wavelengths; this is not surprising since only YSOs exhibit a flat or rising spectrum in the mid-IR, while a wealth of other source types with generally decreasing spectra is found in the Galactic Plane. As a general procedure we have then restricted our analysis for each field only to those sources that were detected at least at 21.4 $\mu$m. The integrated fluxes for each isolated selected source were computed using the IRAM-Grenoble GRAPHIC package and using an aperture hand-tailored to the level where the source brightness value reaches the local background, and the value of the local background was then subtracted. The sources have been labeled with letters to distinguish them from submm counterparts that may or may not be physically associated; this labelling should not generate confusion with the naming of the MSX photometric bands since the context in which they are used in the text should be clear. The fluxes computed in this way are reported in Table 3. We chose not to apply here the same method used for the millimeter and submillimeter maps because here, and especially in the A and C MSX bands, the morphology of the emission is often very complex and cannot be fitted by a simple plateau+Gaussian combination (see the examples below).

The component identification and consequent flux estimate was done directly on the images, and not using the point source catalogue. The reason is that the PSC photometry was compiled using a fixed aperture equal to the reconstruction beam centered on the peak of each detected source. While this may be satisfactory for true point sources, it may lead to significant flux underestimates in the case of compact and yet resolved sources like dusty envelopes around massive YSOs.

As an example we can consider the field of IRAS 10501-5556, a High source, shown in Fig. A.49 in Sect. A.13. This is a relatively ``simple'' field where one component (A) is detected in all four bands of MSX. The analysis of the 8.3 $\mu$m brightness profiles reveals that source A is a resolved and yet compact source with a ${\it FWHM}\sim 4$0 $^{\prime \prime}$, superimposed on an arc-like plateau with the tip of the convexity centered on the point-source, and the concavity oriented S-E. This diffuse structure is also visible in band C (12.1 $\mu$m) but not in bands D and E, suggesting that most of it may be due to PAH emission. To estimate the flux of the source we then manually define a circular aperture centered on source and extending until the signal level of the diffuse arc plateau is reached; this corresponds to a $\sim$2 FWHM diameter aperture, definitely larger than the fixed 20'' aperture used in the automatic aperture photometry generating the MSX Point Source Catalogue and which appears to be inadequate for this type of sources. As concerns the background, it is clear from the morphology of the underlying diffuse plateau that adopting a circular annulus just outside the aperture would lead to severe underestimates of the background; instead, we choose to integrate the emission from the east tip of the arc-like structure and subtract it, after normalizing to the different areas, from the aperture integral carried out on source. This procedure is followed for bands A and C, while for bands D and E (where the diffuse plateau vanishes) we adopt a large aperture freely encompassing the entire compact source down to the uniform background of the field, and a neighboring area is considered for the background estimate. The resulting mid-IR photometry (see Table 3) considerably differs from the one provided by the MSX Point Source Catalogue. The values we compute for bands A and C are a factor 2$\div$3 smaller than PSC numbers, most probably because the background was severely underestimated in the MSX PSC photometry; on the contrary our band D and E values are 30$\div$40% higher than PSC values, most likely due to the larger aperture used in the present work.

The situation illustrated above for IRAS 10501-5556 is representative of most of the sources in the other fields in the present sample. A straightforward adoption of the MSX PSC fluxes would have led to a completely wrong SED characterization for source A; in particular in the $8\rightarrow15$ $\mu$m portion the SED slope would have been decreasing, while a careful analysis of the images and determination of the photometry (admittedly difficult in an automatic photometry pipeline) reveals that the SED is instead steeply rising with wavelength. The implications for the characterization of the source nature are clearly critical. Extreme care should be exercised when studying the mid-IR portion of the YSOs SED using the MSX database, and any conclusion based on the straightforward adoption of the MSX PSC photometry should be regarded with caution.

 

 

Table 3: Mid-IR MSX fluxes.

IRAS field	IR	$\Delta\alpha$	$\Delta\delta$	F8.3	F12.1	F14.7	F21.4
		( $^{\prime \prime}$)	( $^{\prime \prime}$)	(Jy)
05137+3919	A	+12	+8	3.1	4.9	6.3	15
	B	-4	-25	1.0	1.7	1.1	5.3
05345+3157	A	+18	-7	8	11.1	4	8.2
	B	+62	+36	-	-	-	2.7
05373+2349	A	+1	-4	2.2	3.5	6.6	14.4
08211-4158	A	+2	+6	7.6	10.8	8.2	17
	B	-88	-125	1.0	1.4	2.1	6.9
08470-4243	A	+1	+12	13.6	18.8	12.6	56.7
	B	-33	-9	2.0	3.0	1.4	1.2
	E	-86	29	-	-	1.0	4.5
08477-4359	A	+3	+3	12.0	12.3	6.3	14
08563-4225	A	+3	+1	12.8	15.3	14.1	29.2
08589-4714	A	+12	-9	0.6	0.3	0.7	5.2
	B	+52	-72	0.3	1.0	-	1
09131-4723	A	-23	+21	10.9	11.9	6.9	16.1
	B	-51	+98	4.5	6.3	5.8	7.1
10019-5712	A	+12	-2	4.2	6.2	4.6	31.1
	B	-11	+11	3.6	4.1	2.5	11.5
10095-5843	A	+2	+6	6.7	8.8	5.5	12.5
	B	+162	-89	1.1	2.1	0.94	0.26
10184-5748	A	+14	-19	3.6	10.1	44	240
	B	+9	+26	1.3	6.5	15	36
	C	+10	-98	4.5	4.0	5.9	17
	G	-67	-155	1.8	3.1	1.6	8.9
10501-5556	A	+3	-6	2.6	2.1	5.6	34
10555-6242	A	-2	+13	3.1	3.7	3.9	17.8
11380-6311	A	-5	+4	3.3	2.4	2.8	7.8
12063-6259	A	+3	-1	36	74	91	331
12127-6244	A	+2	-4	39	94	162	590
14000-6104	A	-19	-33	14	13	14	35
	B	+212	-126	7.7	8,.4	4.6	11
	C	+91	-30	3.1	4.8	8	10.6
	D	-229	+146	8.4	13.7	18.8	22.7
	E	+313	-92	0.2	0.3	0.9	1
	F	-87	+101	1.6	3.6	3.3	10.7
	G	-51	+68	-	0.7	1.7	4.1
	H	-28	+123	0.1	-	0.5	2.7
	I	-31	+172	0.5	0.75	1	2.6
14057-6032	A	+12	-2	0.6	1.1	2.3	11.2
15068-5733	A	+6	0	2.8	3.9	5.5	52
	B	-91	+59	1.4	2.2	1.7	0.7
	C	+168	-7	1.1	1	1.1	4.2
	E	+83	-125	-	-	1.3	3.4
15454-5335	A	+28	-13	2	1.8	4.3	22.2
	B	-93	+47	0.6	0.9	3.8	20.2
15519-5430	A	0	+10	3.5	5.1	2.1	4
	B	-253	+144	10	13	10.4	12.3
	G	-87	-154	-	1	1.5	5
15557-5337	A	-1	+11	64.5	187	208	713
	B	-19	-23	36.9	108	129	440
	C	+41	-51	-	32	33	145
	D	+23	+39	8.1	13	10.6	79
16148-5011	A	+20	-2	32	46	26	72
	B	+382	-330	48	87	71	57.3
	C	+215	-301	7	9	7.5	36.7
	D	+126	-254	2.6	2.2	2.6	2.4
	E	+33	+322	3	3.8	6	5.5
16232-4917	A	-7	-5	10.6	14	12.1	26.7
16501-4314	A	+3	-4	5.8	7.6	6	13
16535-4300	A	+12	-1	1.2	2.3	1.7	2.2
	C	+72	-256	4.1	6	11.7	12.7
17082-4114	A	-8	-6	12.4	14.5	8.7	17
	B	-53	-60	7.7	8.7	4.9	10.8
17425-3017	A	-3	-4	2.9	2.4	1.6	3.1
	F	-228	-52	0.9	0.9	1	1.8
18024-2119	G	+3	+9	-	-	-	1.2

IRAS field	IR	$\Delta\alpha$	$\Delta\delta$	F8.3	F12.1	F14.7	F21.4
		(s)	( $^{\prime \prime}$)	(Jy)
18144-1723	A	+6	-3	1.4	3.4	5.4	18.2
	B	-178	-12	8.2	4	2.1	1.7
	C	+155	+87	1.2	1.9 6 1	0.9
	D	-63	-159	0.8	0.5	0.4	0.4
	E	-32	-86	-	-	0.7	1.2
18278-1009	A	+15	+3	5	9.7	15.6	15.5
	B	-68	-38	1.3	2.7	1.8	1.9
	C	-12	-169	2.5	3.6	3.6	2.8
	D	+140	-45	4.5	4.8	3.9	2.2
18511+0146	Aa	+2	+5	18.3	26.5	36.1	47.1
19092+0841	Aa	+15	+16	3	3.1	1.7	4.1
	C	-172	+64	2.7	4	2.8	2.5
	D	+224	-108	1.7	1.5	1.1	4.4
20062+3550	A	+17	+3	1	1.3	2.7	17.2
20106+3545	A	+9	+4	2.6	2.5	2.5	10.4
	B	+34	+65	8	3.5	3.1	3.1
21307+5049	Aa	+3	+2	1	1.1	2.7	6.5
21519+5613	A	+10	-5	0.75	1.4	2.2	10.8
22172+5549	A	+2	+18	2	5.8	3.2	7
22305+5803	A	+3	+5	1.2	1.1	2.4	12.1
22506+5944	A	+2	+5	4.7	6.7	5.7	18.4
23385+6053	Aa	+25	+27	2.5	6.4	0.5	0.7
	Ba	-16	+23	0.7	1.5	0.2	0.1
	C	-87	+12	0.9	0.9	-	0.8

^a Source also observed with ISOCAM in LW2 and LW3 filters (see Table 4).

A word of caution should also be given for spectroscopic contamination of the mid-IR photometric bands due to PAHs. Emission from PAHs is generated as non-equilbrium thermal irradiation after absorption of an energetic photon; as such, its presence is related to the intensity of the incident UV field which is quite abundant in regions of massive star formation where several young objects are commonly found which are already on the ZAMS and are able to generate copious amounts of UV photons. PAH broad features are present over the entire wavelength range from 6 to 13 $\mu$m and therefore contaminate both MSX bands A and C. This leads to an overestimate of the dust continuum emission which is difficult to parametrize for the sources of our sample. A quick analysis of unpublished Spitzer/IRS spectra of the IRAS 23385+6053 region reveals that the continuum-subtracted integral of PAHs emission contributes between 50% and 70% of the photometric fluxes in bands A and C. For the SED model fitting (see below) we will then assume that the true dust continuum in the MSX band A and C may be up to a factor 4 lower than the photometric fluxes.

It may be objected that the MSX data may not the optimal choice to characterize the mid-IR spectrum of YSOs in crowded areas given the coarse spatial resolution. The use of GLIMPSE survey data (Benjamin et al. 2003) might be more appropriate in principle, given Spitzer's better resolution in the mid-IR, but there are two limitations. First GLIMPSE does not cover the entire plane, being limited to $\vert l\vert \leq 65$$^{\circ}$. Second, its 3.6-8 $\mu$m wavelength range is not optimally suited to catalogue objects according to traditional SED-based schemes. Detailed modelling (Whitney et al. 2004) shows that for luminous objects ( $T_{\star} \ge 10^4$ K) Class I and II spectra can be easily distinguished beyond 8$\div$10 $\mu$m; the comparison of their Figs. 3b and 5 shows that the 1-10 $\mu$m SED portion is rising for both Class I and II. The 10-20 $\mu$m portion shows a clear turn-down for Class II, while it keeps rising for Class I. MSX data are therefore more suited to pinpoint the nature of the mid-IR emitting object; the MSX fluxes will most likely contain the flux contribution from unresolved close-by objects, but they will certainly be dominated by the brightest and reddest of the objects.

Four fields were also observed with the ISOCAM instrument, on-board the ISO satellite, in the two broad-band filters LW2 (6.75 $\mu$m) ad LW3 (15 $\mu$m) with a ${\it PFOV}=3$ $^{\prime \prime}$). Data for field IRAS 23385+6053 were published in Molinari et al. (1998b); the same calibration and reduction procedures were adopted for the other fields. Background-subtracted fluxes were estimated just as for the MSX fluxes, and are reported in Table 4.

 

 

Table 4: Mid-IR ISOCAM fluxes.

IRAS field	IR	$\Delta\alpha$	$\Delta\delta$	F6.75	F15
		( $^{\prime \prime}$)	( $^{\prime \prime}$)	(Jy)
18511+0146	A1	-2	-2	19	27.6
	A2	9	-20	2.5	2.1
19092+0841	A	13	3	1	0.8
21307+5049	A	-14	-4	0.7	1.4
23385+6053	A1	21	30	0.4	0.2
	A2	24	24	0.5	0.2
	A3	24	5	0.5	0.12
	B1	-21	17	0.25	0.1
	B2	-21	9	0.1	0.4

   
3.1.3 The far-infrared

IRAS PSC fluxes are computed over an area ( ${\it FWHM}=5$$^{\prime}$ @100 $\mu$m) that encompasses the entire source plus plateau structure, so it is necessary to quantify and subtract the contribution of the diffuse emission before using the IRAS fluxes in the fit to SED models of YSOs. As noted above, the MSX maps show that the mid-IR diffuse emission, due either to warm dust or more probably to well-known PAH features is only visible in bands A and C, but readily disappears in bands D and E and hence does not contribute also at longer wavelengths. On the other hand the plateau visible in the 1.2 mm maps suggests the presence of a cold and diffuse ISM component in which the YSO is located; the properties of dust within star forming regions have been investigated in nearby star forming regions, where mean temperatures of $\sim$18 K have been determined (Schnee et al. 2005). Given the relatively higher luminosities of the sources in the present sample we choose to adopt a mean dust temperature of 20 K for the diffuse dust, as also suggested by estimates derived using IRAS colors in the outer layers of high-mass star forming regions (e.g., Tej et al. 2006). After verifying that the mean background level in the 1.2 mm map is 0 within the rms, we can obtain an estimate of the emission in the diffuse plateau by integrating the flux of the whole source and subtracting the integrated flux of the compact peak obtained previously in AIPS/JMFIT. We then use a greybody function with $\beta =2$ and T=20 K to estimate the global contribution of this diffuse component to the 60 and 100 $\mu$m IRAS fluxes. These values are then subtracted from the IRAS PSC2 fluxes; if there is more than one millimeter core detected in an area of a few arcminutes in diameter, these residual 60 and 100 $\mu$m fluxes will be split among the various millimeter cores in proportion to their millimeter fluxes. We acknowledge the arbitrariness of the procedure but we believe it is the most conservative approach.

   
3.2 The SED models

It has been shown (Molinari et al. 1998a, 2000; Brand et al. 2001) that the sample used to extract the sources for the present work selects massive YSOs in an early evolutionary stage. Masses of circumstellar envelopes traced in the sub-millimeter are usually in the range from a few to hundreds of solar masses (e.g. Molinari et al. 2000). The ubiquitous association of massive YSOs with strong molecular outflows exhibiting mass loss rates of the order of 10 $^{-4}{\div}10^{-3}$ $M_{\odot }$ yr^-1 (Zhang et al. 2001; Beuther et al. 2002b) probably requires that seemingly high-rate accretion must be ongoing. Besides, the presence of molecular outflows also suggests that the envelopes are not spherical and that equatorial disks and polar cavities should be present. The direct evidence for the presence of equatorial disks also in high-mass YSOs was established with detailed interferometric studies toward specific targets like IRAS 20126+4104 (Kurtz et al. 2000). As for the nature of the embedded accreting object it is commonly believed that we are in the presence of ZAMS stars not only in case of UCHII regions, but probably also in the case of Hot Cores (HCs). The question of whether this is also true for the so-called UCHII precursors is still under debate. The absence of thermal free-free radio emission is not unequivocal evidence in this respect (Molinari et al. 1998a); extreme compactness of a newborn HII region, strong UV absorption by dust close to the ZAMS star and/or quenching by residual accretion can all account for the lack of radio emission. For the present analysis we use ZAMS stars as embedded objects in model envelopes. Clearly, this becomes a critical assumption to be questioned if it is impossible to obtain a proper fit to the data for any reasonable combination of envelope parameters.

Two codes are publicly available to generate synthetic SEDs suited for star formation: DUSTY (Ivezic & Elitzur 1997) and the 3D Radiation Transfer Code by Whitney et al. (2003). We adopted the latter due to its 3D nature which allows us to model the disc around the central star, to take into account the polar cavities due to outflows, and to derive, for the same input parameters, how the SED changes when seen through different lines of sight.

When we completed the first data analysis on the target fields, the impressive model grid produced by Robitaille et al. (2006) was not yet publicly available. We then computed an extensive grid of SED models for massive YSOs spanning a wide range of envelope radii, accretion rates, polar cavity angles and luminosities. Almost 400 parameter combinations have been used, and for each one of them the model computes the SED for a set of ten inclination angles of the system with respect to the observer; our library of models will then consist of about 4000 SEDs. Table 5 summarizes the extent of parameters space mapped with our model grid.

 

 

Table 5: Parameter space for the SED model grid.

ZAMS	$M_{\rm D}$	$R_{\rm D}$	$R_{\rm env}$	$\theta _{\rm p}$	$\dot{M}_{\rm\rm acc}$
Sp. Cl.	($M_{\odot }$)	(AU)	(104 AU)	($^{\circ}$)	(10-4 $M_{\odot }$/yr)
B5	0.3	100	[1,2,4,6,8]	0, 20, 40	[0.1,0.5,1]
B3	0.6	200	[1,2,4,6,8]	0, 20, 40	[0.1,0.5,1]
B1	0.9	300	[1,2,4,6,8,10]	0, 20, 40	[0.5,1,5]
B0	2.0	500	[1,2,4,6,8,10]	0, 20, 40	[0.5,1,5]
O8	2.0	500	[2,4,6,8,10 12.5,25]	0, 20, 40	[0.5,1,5,10a]
O5	5.0	1000	[4,6,8,10,12.5,15,30]	0, 20, 40	[1,5,10,50b]
O3	5.0	1000	[4,6,8,10,12.5,15,30]	0, 20, 40	[1,5,10,50b]

^a $\dot{M}$ = 10^-3 $M_{\odot }$ yr^-1 only for models with $R_{\rm env}$ =  $2.5\times 10^5$ AU.
^b $\dot{M}$ =  $5\times 10^{-3}$ $M_{\odot }$ yr^-1 only for models with $R_{\rm env}$ =  $3\times 10^5$ AU.

Here Col. 1 is the spectral class of the embedded YSO; Cols. 2-3 are the mass and radius of the circumstellar disk; Col. 4 is the envelope radius; Col. 5 is half of the opening angle of the envelope polar cavity; Col. 6 is mass accretion rate which eventually defines the mass of the envelope. Disk parameters were not changed since the SED is insensitive to it when a massive envelope is present (Robitaille et al. 2006). The a posteriori comparison with Robitaille's model grid shows that we extended the parameter space to higher stellar masses, and to larger and denser envelopes. We limited the parameter exploration to embedded ZAMS stars only; the quantity of free parameters in the model is by far larger than the number of data points at hand, and we prefer to have a clear-cut choice of the star's nature rather than an unreliably testable spread of ages. The words ``ZAMS models'' will be used in the rest of the paper to signify our model grid.

Figure 2: Variation in the shape of the SED produced by the Whitney et al. (2003) models for different combinations of input parameters. For each panel the title shows the parameter held fixed, and the legend (or the labels in the figure) indicate the parameter being varied.

Given the Monte Carlo nature of the models, the accuracy of the output model SED depends on the number of photons adopted in each run and to keep computing time to a manageable level, we limited the models to 10⁷ photons. This results in a very ``noisy'' model SED longward of $\sim$500 $\mu$m, and this required a linear extrapolation (in log-log space) of the SED from shorter wavelengths down to the millimeter range to enable a useable comparison to the submm and mm data. This procedure was carried out selecting manually for each model run the SED portion to be used for extrapolation; we experimented selecting different portions to be used for extrapolation between 200 and 500 $\mu$m  and we obtained 1.2 mm fluxes (the wavelength where we compare models and data for most of the southern sources) spanning a factor of 4 at most. In the subsequent SED fitting (see below) we will then consider a threshold of 50% as the maximum tolerable amount for model departures from the data (in excess or in deficit).

To evaluate the diagnostic power of the SED model grid, Fig. 2 illustrates the effect of varying various parameters on the appearance of the SED. An increase in the accretion rate (panel a) is mainly reflected in a more massive envelope and results in a SED shifting toward longer wavelengths due to the higher content of cold dust which also increases the extinction toward the warmer regions of the envelope. Increasing the envelope radius (panel b) has a similar but less pronounced effect, because the density reaches its minimum at the outer envelope and so the increase in total mass of the envelope is not as important compared to an increase in volume density as implied by a rise in $\dot{M}_{\rm\rm acc}$. The polar cavity aperture angle (panel c) is particularly important in providing relatively clear lines of sight toward the inner and warmer regions of the envelope; the main change in SED for increasing cavity aperture angles is to greatly enhance emission shortward of 20 $\mu$m. Finally, panel d shows the relation between the emerging SED and the spectral class of the embedded star; as it could be easily anticipated, the change is both in luminosity and peak $\lambda$, since an envelope of fixed radius is consistently warmer for a harder input radiation field.

   
3.3 SED fitting and model selection

The fitting of the object SEDs to our nearly 4000 models was done by changing the model distance using a weighted least-squares algorithm in which the weights are not related to the uncertainties of the data points, but are assigned so that a certain balance is achieved between different wavelength regions of the SED. Our SED data points are split into three groups: i) ``MIR'' with 8 $\mu$m  $\leq\lambda\leq$ 25 $\mu$m, ii) ``FIR'' with 60 $\mu$m  $\leq\lambda\leq$ 100 $\mu$m, and iii) ``mm'' with $\lambda \geq300$ $\mu$m. The quantity that we minimize in the fit is

$$\xi=\sum _{i=1} ^{N_{\rm MIR}} {{(F_i-M_i)^2}\over{\zeta _i ^... ... + \sum _{k=1} ^{N_{\rm mm}} {{(F_k-M_k)^2}\over{\zeta _k ^2}}$$ (1)

where F_i,j,k are the fluxes of each data point (excluding upper limits) and M_i,j,k are the model values at the same wavelengths. The choice we make is that the weights of the MIR and mm portions of the SED should be equal, and that of the FIR part should be 1/3 of that. The reason for this choice is that the FIR portion of the SED only relies on IRAS 60 and 100 $\mu$m data which have much worse spatial resolution compared to the MSX and submm maps. The assignment of FIR fluxes to compact components is mostly based on assumptions and extrapolations from other wavelengths, and we transfer this uncertainty into the fitting procedure by assigning a lower weight to the FIR with respect to adjacent wavelengths ranges. If we make 1 the total of the weights for the whole SED, the individual $\zeta$ values in each wavelength group will therefore be assigned so that

$${{1}\over{\zeta_{[i,j,k]}^2}} = \left[ \frac{3}{7},\frac{1}{7},\frac{3}{7} \right].$$ (2)

In order to ensure that the values of $\xi$ for different sources can be compared to one another, no further renormalization of the $\zeta$ is done if one or two sums in Eq. (1) are missing. In some sources we may miss the mm portion, or the FIR data are only upper limits; in this way we are sure that each wavelength portion weighs always the same independently of the other portions.

For each model, the distance $D_{\rm mod}$ for which the minimum of $\xi$ is reached is then compared to the known kinematic source distance $D_{\rm sou}$, and the model is accepted if the difference relative to $D_{\rm sou}$ is below a minimum threshold of 30% adopted to allow for adjustments due to the finite steps in the parameters' values in the models grid. For each source we also computed the kinematic distance spread corresponding to a departure from pure Galactic rotation of $\pm$5 km s^-1 (the cloud-cloud velocity dispersion; see, e.g., Brand & Blitz 1993); if this relative spread was higher than 30% then the model acceptance threshold was raised accordingly. In the case of inner Galaxy sources with near-far distance ambiguity we tried a fit assuming both distances to avoid that a prior, possibly wrong, assumption made to resolve the distance ambiguity may prevent a fit being found.

For each source, all the accepted models were $\xi$-ranked and the overall minimum $\xi$ model was selected. Figure 3 reports the distribution of the minimum $\xi$ obtained for each fitted source SED.

Figure 3: Distribution of minimum $\xi$ for all fitted source SEDs.

The distribution shows a first peak with $0\leq\xi\leq0.1$, a second peak at $0.1<\xi\leq0.25$, and a more or less uniformly spread distribution for $\xi>0.25$. All minimum $\xi$ models were visually examined to verify the correspondence of the $\xi$ values to the effective ability or impossibility to find a good quality model fit. The presentation of all best model fits to all sources is given for each source in the appendix; here we give a brief overview of the results to illustrate the general guidelines used to accept or reject a model fit.

Figure A.45 in the appendix shows a model fit to source IRAS 10184-5748 A; visually, the quality of the fit is very high, and corresponds to $\xi=0.09$. Since $\xi$ repesents the cumulative weighted difference between the model and the data, a low value of $\xi$ could also result from an exceptional model-data agreement on most of the data points with the exception of one or two points where the disagreement is clearly apparent. The noticeable case is represented in Fig. A.29, where a $\xi=0.1$ value is obtained, but where the fit appears to depart by a factor of $\sim$4 from the data at 8 $\mu$m, and it is also $\sim$50% off at 100 $\mu$m. The potential impact of PAHs emission on the MSX A band (see Sect. 3.1.2) might explain the model deficiency at 8 $\mu$m; on the other hand the uncertainties in the estimate of the contamination of IRAS 100 $\mu$m fluxes by extended emission might have led in this case to an incorrect source flux estimate, and the relatively lower weight assigned to the FIR portion of the SED explains that a low value of $\xi$ is still obtained. We conclude that all $\xi\leq0.1$ model fits are acceptable.

On the opposite side is the best model fit obtained for source IRAS 14000-6104 A, and presented in Fig. A.67 as a dashed-triple-dotted line. The fit has $\xi=0.31$, and visual inspection shows that the model fails by more than a factor of 3 in fitting the data longward of 60 $\mu$m. The departure at 100 $\mu$m, although severe, could still be justified by an incorrect source flux assignment given the complexity of the field and the other sources that may contribute in the IRAS beam (see Fig. A.66); the 1.2 mm flux estimate, however, is much more reliable and it is then not possible to reconcile it with the model deficiency in the submm. All models with $\xi>0.3$ are of similar, or lower, quality and have been rejected; the only exceptions are sources IRAS 14000-6104 H and I (Figs. A.71 and A.72), which are tentatively accepted given the impossibility to fit them otherwise (see below) due to lack of reliable data in the FIR. These sources are also of very limited importance for the results of this paper so the acceptance or rejection of their model fit does not require further attention.

The situations where the correspondence between the $\xi$ value and quality of the fit is not immediately apparent are found in the intermediate range $0.1<\xi\leq0.3$. The typical fit with these $\xi$ values is one where the departures of the model from the data are present in both the MIR and mm portions of the SED, which have equally high weight in the fitting procedure, like source IRAS 15454-5335 A (Fig. A.81) with $\xi=0.27$. The deficiency of the model at 1.2 mm is within the threshold (50%) which we consider acceptable for the fit (see Sect. 3.1.1); likewise, the data excess of a factor $\sim$3 at 8 $\mu$m is still considered acceptable given the level of the potential PAH contamination (Sect. 3.1.2); this is a typical situation for the majority of the fits with $\xi$ in this range, which are all accepted. Only in a handful cases does the fit severely differs in one point only, and the decision for acceptance or rejection of the model required a careful scrutiny of the data quality (see the discussion for each source in the appendix).

 

 

Table 6: Derived parameters from model fits.

IRAS			$\xi$		type	Sp.	$L_{\rm bol}$	$M_{\rm env} ^c$	$R_{\rm env}$	$\theta _{\rm p}$	id
						Cl.		($M_{\odot }$)		($^{\circ}$)	($^{\circ}$)
05137+3919	L	A/1	0.12		IR-P	O8	255	2890	25	20	66
		B/2	0.11		IR-S1	B1	5.5	120	10	20	72
05345+3157	L	A/3	0.09		IR-S1	B5	0.5	7	6	40	60
		B/1	0.11f		MM-P		1.7	35
		2			MM-S		0.8	20
05373+2349	L	A/1	0.04		IR-P	B5	0.4	26	6	0	60
08211-4158	L	A	0.15g		IR-S2	B3	1.0	3	8	20	53
		1			MM-P		1.1	80
		B	0.01		IR-S1	B5	0.4	8	6	20	53
08470-4243	H	A/1	0.22		IR-P	B1	4.8	85	8	20	66
08477-4359	L	A	0.06		IR-S1	B5	0.45	0.1	1	40	60
		1			MM-P		0.7	85
08589-4714	H	A/1	0.11		IR-P	B5	0.3	20	8	40	0
08563-4225	L	A/1	0.06		IR-P	B3	1.0	23	8	20	60
09131-4723	L	A/1	0.15		IR-P	B5	0.45	26	8	20	66
		B	0.02		IR-S1	B5	0.3	4.3	8	40	60
		2			MM-S		0.47	25
		3			MM-S		0.3	21.5
10019-5712	H	A/1	0.13		IR-P	B3	0.65	11	6	40	37
		B/2	0.09		IR-S1	B5	0.3	7	6	40	37
10095-5843b	L	A	0.17l		IR-P	B1	4.8	85	8	20	66
10184-5748	H	A/1	0.09		IR-P	O5	255	2890	30	20	66
		B	0.02		IR-S1	B0	15	54	8	40	37
		C	0.02		IR-S1	B1	4.5	21.6	4	40	60
		G/2	0.14		IR-S1	B1	3.5	62	8	40	45
10501-5556	H	A/1	0.07		IR-P	B1	3.3	88	8	40	37
10555-6242	H	A/1	0.18		IR-P	B1	3.3	88	8	40	37
11380-6311a	L	A/1	0.10		IR-P	B5	0.3	20	8	40	37
12063-6259	H	A/1	0.04		IR-P	O5	373	2890	30	20	78
12127-6244	H	A/1	0.13		IR-P	O3	888	2050	30	20	72
		2			MM-S		105	1352
14000-6104	L	A	0.31		IR-S2	B0	19	54	8	40	53
		1			MM-P		21	1965
		B	0.13		IR-S1	B1	3.6	88	10	40	45
		C	0.01		IR-S1	B3	1.0	2.5	8	0	53
		D	0.01		IR-S1	B1	8.6	21.5	4	40	78
		H	0.31h		IR-S1	B1	4.9	132	10	0	0
		I	0.34h		IR-S1	B1	5.0	132	10	0	45
14057-6032a	H	A/1	0.02		IR-P	B5	0.25	6.5	4	40	25
15068-5733b	H	A	0.21		IR-P	O5	176	2145	30	40	37
15454-5335b	H	A	0.27		IR-P	O5	225	2890	30	20	37
		B	0.02		IR-S1	O8	50	720	25	20	53
15519-5430a	L	A	0.56		IR-S2	B3	1.1	4	8	40	53
		1			MM-P		9.8	485
		G	0.01		IR-S1	B5	0.5	3	8	20	60
15557-5337a	L	A/1	0.05		IR-P	O8	65.4	721	25	20	73
		B/2	0.06		IR-S1	O8	57	721	25	20	66
		C	0.01		IR-S1	B0	13	4.4	6	40	25
		D	0.05		IR-S1	B1	4.1	1	2	40	53
16148-5011a	L	A	0.03		IR-S1	O8	71.5	20	10	0	72
		C	0.16		IR-S1	O5	255	2890	30	40	66
		1			MM-Pe		240	2650
		2			MM-S		43.2	1160
16232-4917a	L	A/1	0.06		IR-P	B1	4.25	88	10	40	53
16501-4314a	L	A	0.26i		IR-S2	B1	6.1	26	10	0	78
		1			MM-P		24.7	610
		2			MM-S		24.7	610
16535-4300a	L	A	0.54		IR-S2	B3	2	0.4	1	40	84
		1			MM-P		13.4	750
		C	0.01		IR-S1	B1	7.5	1.2	2	0	84
17082-4114a	L	A	0.36		IR-S2	B3	2.1	0.1	1	20	84
		1			MM-P		2.5	282
		2			MM-S		1.4	166
		3			MM-S		1.5	185
17425-3017	L	A	0.53		IR-S2	B1	5.1	13	10	40	53
		1			MM-P		23.1	1250
18024-2119b	L	G/1	0.53		MM-P		84.2	5850
18144-1723b	L	A	0.20		IR-P	O5	176	2145	30	40	37
18278-1009	L	A	0.13		IR-P	B1	9.1	120	10	20	84
18511+0146	L	A1	0.17m		IR-S2	B1	6.4	1.2	2	0	72
		1			MM-P		9.2	330
		A2	0.21		IR-S1	B1	4.5	120	10	20	60
19092+0841a	L	A	0.43		IR-S2	B3	1.1	8	6	0	66
		1			MM-P		7.8	738
20062+3550	H	A	0.15		IR-P	B1	4.8	85	8	20	66
20106+3545b	L	A	0.11n		IR-P	B1	3.6	88	10	40	45
		1			MM-S		1.8	45
21307+5049	L	A/1	0.16		IR-P	B1	3.6	88	10	40	45
21519+5613	H	A/1	0.18		IR-P	B1	5.4	85	8	20	75
22172+5549	L	A/1	0.16		IR-P	B1	3.6	88	10	40	45
22305+5803	H	A/1	0.05		IR-P	B1	4.8	86	8	20	66
22506+5944	H	A	0.31		IR-S2	B1	7	25	10	20	78
		1			MM-P		8.9	671
23385+6053	L	1	0.36		MM-P		17.5	222

^a Near distance adopted.
^b Far distance adopted.
^c A correction factor of 0.4 has been applied to all envelope masses of IR sources to account for different mass opacity assumptions in Whitney et al.'s models and in our greybody mass estimates.
^d Between polar axis and plane of the sky.
^e Source MM1 is assumed as primary instead of C because source C is not part of the A/MM1/MM2 complex (it's a different IRAS source), although we assume it connected to the same star forming region.
^f Model rejected, see Sect. A.2.
^g Model rejected, see Sect. A.4.
^h Model accepted, see Sect. A.18.
ⁱ Model rejected, see Sect. A.26.
^l This is not the minimum-$\xi$ model (see Sect. A.11), but it is chosen among the 10 best-$\xi$ models. The minimum $\xi$ model ($\xi=0.14$) provides a better fit to the mid-IR portion, but misses the 100 $\mu$m by slightly more than 50%.
^m Model rejected, see Sect. A.33.
ⁿ Same as footnote l; see Sect. A.36.

Thus, we conclude that the $\xi$ parameter offers a sufficiently accurate representation of the model fit quality to decide on acceptance, $\xi\leq0.3$, or rejection, $\xi>0.3$ of a model fit. Only in a few cases with $0.1<\xi\leq0.3$ a careful analysis of the data quality led us to reject the model fit (noted in Table 6). The top 10 models were also visually inspected to make sure there is no higher $\xi$ model which ``to the eye'' gives a better overall fit; this happened only in two cases (noted in Table 6). In all cases where a ZAMS model fitting the entire wavelength range could not be found, we assumed that the mid-IR and the mm counterparts are not physically related, in spite of the positional coincidence within the uncertainties. The procedure followed was then to try a two-component fit.

We start by fitting ZAMS models to the mid-IR portion of the SED ( $\lambda \leq25$ $\mu$m) following the same procedure outlined above. We then subtract the FIR fluxes of this model SED from the FIR and mm observed fluxes to obtain the SED of the component supposedly responsible for most of the FIR-mm fluxes. Fitting ZAMS models to these ``cold'' SEDs does not provide any result; the observed levels of mm flux always require models that produce too much flux shortward of 100 $\mu$m even for the most massive and opaque envelopes in the parameter range spanned by our model grids. We would then need to produce less bolometric luminosity for a given envelope mass, which inevitably translates into lowering the luminosity of the central object; it is indeed well known from detailed models of pre-MS evolution (e.g. Palla & Stahler 1993) that for $M_{\star}$ $\geq$ 3 $M_{\odot }$ the luminosity generated by the forming objects (either by radiative contraction or H-burning if already on the ZAMS) starts to dominate over accretion (this mass threshold applies for $\dot{M}$ = 10^-5 $M_{\odot }$ yr^-1, and increases with $\dot{M}$). It is then clear that a lower luminosity object, for a given mass, is incompatible with its being on the ZAMS.

Releasing the ZAMS assumption for the central object complicates the analysis because one has to model the input radiation field from an accreting and contracting object; not only does this involve more free parameters, but the physical structure of the object is strongly influenced by the accretion flow. Since our present purpose is to characterize the nature of these FIR-mm condensations from a statistical point of view, we feel it is more appropriate to apply a simple analysis where these condensations are treated as single-temperature optically-thin envelopes which can then be modeled via a standard optically thin greybody, where

$$F\propto M\cdot B_{\nu}(T)\cdot \nu ^{\beta}.$$ (3)

The mass is derived from the dust column density assuming a gas/dust ratio =100 and the opacities from Priebisch et al. (1993); $\beta$ is the usual exponent for the power-law frequency dependence of dust emissivity. Mass, temperature and $\beta$ are the three free parameters in our greybody fits and provide a decent representation of the global properties of the FIR-mm clumps not associated with mid-IR emission and which cannot be fitted by the Whitney et al. (2003) models. There are a handful of caseswhere the SED of a source is only known longward of 20 $\mu$m. In those cases no two-component fit is obviously possible and the modelling is directly done with the greybody fit. In the case of a distance ambiguity we always first try to find a good ZAMS fit assuming the ``near'' distance. If a fit is not found we then try with the ``far'' distance, and if a fit is not found we will revert to the ``near'' distance for the ZAMS+GB composite fit.

   
3.3.1 General considerations

The paragraphs above give a detailed description of the method followed to characterize the sources in each field. Clearly each field has its own peculiarities which often require dedicated consideration; we report in the appendix, for each examined field, a brief data analysis log complete with images and SEDs, documenting the particular assumptions and problems encountered for each source.

Clearly the quality of the data, mainly in terms of spatial resolution, and for the assumptions we necessarily had to make in many steps of the science analysis are such that the results should not be regarded as valid on a source-by-source basis, but rather in a statistical sense. In other words the number of examined fields is such that the global trends emerging from the analysis can be trusted, although the conclusion for a specific field may be incorrect.

An indication that we find is that the positional association of a mid-IR and a submm counterpart does not seem to guarantee per se that the two sources are physically the same. There are cases (see the appendix for a complete report) where the counterparts have a separation of less than 20 $^{\prime \prime}$ and yet the SED fitting using a unique embedded ZAMS model is unsuccessful. There may be clearly peculiarities, effects or physical properties which the model does not take into account and that might explain the inability to find a plausible fit; yet the possibility that the non-association is real is certainly plausible. For example, how would a source like IRAS 23385+6053 in Fig. 1 left appear from a line of sight looking from the east in the plane of the sky? The mid-IR emission patch and the millimeter core would be positionally associated and yet it would not be possible to model them with a single ZAMS model.

   
4 Results

The methodology described in the previous section has been applied to analyse all the fields listed in Table 1. The results concerning the nature of the examined sources are reported in Table 6, where for each component we report in Cols. 1 and 2 the IRAS name and its High/Low type. Column 3 is the component identifier built from the source IDs in MSX (letters) and mm images (numbers); both types of codes are present if the SED fitting procedure yields a satisfactory overall match, and the counterparts are hence to be considered associated. Column 4 identifies if the source could be fitted with a Whitney et al. model (``IR'') or with a simple greybody (``MM''); in the former case Col. 5 also report the obtained ZAMS spectral class. Columns 6-7 are the bolometric luminosity and envelope mass as derived from the fitted models/greybodies; only in the case of ``IR'' sources, Cols. 8-9 report the model envelope radius (in AU) and the semi-aperture (in degrees) of the envelope polar cavity.

It is also apparent from Tables 2 and 3, that there are several cases where more than one mid-IR and/or mm source are considered in a 0.2$^{\circ}$ $\times$ 0.2$^{\circ}$ field. Sources that were neglected in the present analysis are i) mid-IR sources clearly separated from any sign of mm emission and with rapidly decreasing MSX mid-IR fluxes (and that can therefore be reasonably assumed to be field stars not associated with the SFR pinpointed by the IRAS source), and ii) isolated millimeter cores with no sign of mid-IR or far-IR emission whatsoever and which are impossible to fit in a reliable way (see appendix for details).

The procedure adopted for the fitting of models to the data points is such that a good fit is often obtained for a model that is scaled to a distance that sometimes can be up to 50% different (in either direction) from the source kinematic distance (see Sect. 3.3 for a discussion of the acceptance threshold). It is then clear that the physical parameters of the best fitting model will not really represent the fitted source. This may be due to the fact that the grid of models does not span the parameter space continuously, so that a different distance compensates for parameter variations finer than the resolution of the model grids. To verify this we ran the fitting procedure for one source (we chose IRAS 05137+3919 B, with the best $\xi=0.11$) without imposing any distance threshold, and looked at the parameters of all the models that provide a $\xi \leq 0.15$ (and so, acceptable in principle). Good models were found with spectral class varying from B1 to B5, $R_{\rm env}$ from from $4\times 10^4$ to 10⁵ AU, polar cavity aperture $\theta _{\rm p}$ between 20 and 40$^{\circ}$, and accretion rates between 10^-5 and 10^-4 $M_{\odot }$ yr^-1. On one side this shows that there are several combinations of model parameters that can conspire to produce a good SED model, but on the other it also says that if a fit is not found this is a strong indication that the models and the observed SED are not compatible. To mitigate the former aspect we can say that the model parameters for the best $\xi$ fit should not be taken at face value for each individual source; the global parameter range spanned by all fitted sources is however significant since the best fit model distances will be both lower and higher than the kinematic distances.

Concerning the macroscopic parameters used later in the discussion, the envelope mass and bolometric luminosity, Fig. 4 shows the relationship between the two parameters for the models selected in the above simulation.

Figure 4: $L_{\rm bol}$- $M_{\rm env}$ plot for all models fitted to source IRAS 05137+3919 B with $\xi \leq 0.15$; the log-log fit has a slope of 1.14. The two arrows starting at the plus sign illustrate the extent of luminosity and mass variations for a difference of 30% in distance from the nominal (kinematic) value.

Both mass and luminosity vary by large amounts, but on average their variations are correlated (as could be easily anticipated) with a slope of 1.14. Both parameters are found to vary with the square power of the distance (as expected). The two arrows show the extent of the variations both in mass and luminosity for a distance difference of 30% with respect to the nominal (kinematic) value.

4.1 Source classification

Modeled sources were classified depending on whether they can be fitted with an embedded ZAMS model, called ``IR'' sources, or with a greybody model only, called ``MM'' sources. A further subclassification is made to identify in each field the source with the most massive envelope; such a source will have a flag ``P'' (for primary) while all the others will have a flag ``S'' (for secondary). Clearly, there will be only one ``P'' source per field.

A finer subdivision is useful within IR-S sources to distinguish between those that can be consistently fit with an embedded ZAMS model from the mid-IR throughout the millimeter just like the IR-P (but they are not the ones with the most massive envelope in the field), called ``IR-S1'', and those fitted with a ZAMS model limited to the mid-IR portion of the SED in case of a composite two-components fit, called ``IR-S2''.

Table 6 reports the results of the fitting procedure for the various components modeled in each field. The clearest entries are those where the entire SED could be fit by a single model; in these cases both the mid-IR and millimeter component codes are listed for the same entry under the ``Comp'' column. In most cases the fitted model is an embedded ZAMS (i.e. an ``IR'' source), so that the source is classified either as an IR-P or as an IR-S1, and the corresponding value of $\xi$ is reported in boldface (which means that the model was accepted); typical examples are Mol008/IRAS 05137+3919 A/1 and B/2. In only two cases, IRAS 18024-2119 G/1 and Mol011/IRAS 05345+3157 B/1, which are only detected at 21.4 $\mu$m in their mid-IR SED portion, the entire SED could not be fit by a ZAMS model (and so the $\chi$ value is not in boldface), and a greybody fit was used.

MM-P objects are mostly found in ZAMS+greybody composite fits and so in Table 6 they appear together with IR-S2 (emphasized by chain brackets). The attempt to fit the entire SED with an embedded ZAMS model was clearly not successful (the $\chi$ value is not in boldface), and the proximity of the mid-IR and millimeter components is a projection effect. An example is IRAS 08211-4158. As already mentioned, the composite fit was always attempted when the millimeter and mid-IR counterparts were less than $\sim$30 $^{\prime \prime}$ away (to account for potential astrometric problems). In the case of IRAS 08477-4359 and IRAS 16148-5011 the separation is larger so that the two counterparts are fitted separately.

4.1.1 IR sources

Figure 5 shows the distribution of the values for the relevant model parameters for IR sources.

Figure 5: Distribution of model parameters for IR-P (full lines), IR-S1 (dotted line) and IR-S2 (dashed line) sources. Top-left: model inclination angle of the polar axis i with respect to the plane of the sky. Top-right: semi-amplitude of the polar cavity angle $\theta _{\rm p}$. Bottom-left: envelope radius. Bottom-right: mass accretion rate.

From the fitting method it is clear that the SED of IR-P and IR-S2 sources are different, with the latter having the SED peak at lower wavelengths (as can be seen from the SED plots in the appendix). There are several parameters that can induce this difference (see Fig. 2); the distributions in Fig. 5 confirm that IR-P and IR-S1 cover a common range of model parameters, and that the reason for the different mean SED shapes between IR-P and IR-S2 is not caused by a systematically different inclination angle of the models, or a wider aperture of the outflow-excavated polar cavity. However, IR-S2 envelopes are less dense (i.e. with lower accretion rates); together with the average similarity in fitted spectral classes between IR-P and IR-S2 sources, this would seem to suggest that IR-S2 might represent a more evolved stage of an IR source, where the circumstellar envelope has been experiencing a dispersion process for some time.

The envelope masses provided by the Whitney et al. (2003) models were multiplied by a factor of $\sim$0.4 to enable a consistent comparison with the masses derived from our greybody fits due to the different dust opacities used; the values reported in Table 6 for IR sources include this correction factor. This factor was determined by extrapolating the frequency dependence of the dust opacities used in the Whitney et al. models to $\lambda=1.3$ mm and normalising it to the value we adopted at the same wavelength in our greybody fits (Preibisch et al. 1993). A correction was also applied to the model luminosities to account for the fraction of the flux that is scattered out of the line-of-sight in the presence of a polar cavity; the correction is recommended by Whitney et al. (2003) and depends both on the polar cavity aperture and the inclination of the model envelope with respect to the line of sight.

4.1.2 MM sources

The parameters of MM sources are derived from simple greybody fits to three points, on average: the 60 and 100 $\mu$m IRAS fluxes, tentatively corrected for the contribution of extended emission and the submm/mm data. We report in Fig. 6 the distribution of fitted temperatures and dust emissivity index for the MM-P (full lines) and MM-S (dashed lines) sources.

Figure 6: Distribution of greybody parameters for MM-P (full line) and MM-S (dashed line) sources. Left: dust temperature. Right: dust emissivity index $\beta$.

Dust temperatures are found in the range between 20 and 60 K for all MM sources, with a median value of $\sim$30 K. For the values of $\beta$, the dust emissivity index, we find that the entire range between 0 and $\sim$2.5 is spanned by our MM sources, with a median value of 1.5. No significant systematic differences are found between MM-P and MM-S sources.

4.1.3 IR-P vs. MM-P sources

While the SED differences between IR-S2 and MM sources do not require explanation, being inherent to the SED analysis method, the comparison of MM with IR-P is more instructive. We restrict the comparison to MM-P sources to outline and characterise the SED appearance of the most massive source in each of our target fields.

Figure 7: Color-color diagrams for all-P sources in the target fields; full and empty circles represent IR-P and MM-P sources, respectively. The thick line marks the locus of greybody emission for various indicated temperatures and dust $\beta =2$. The colors are determined from the models (Whitney's models for the IR-P and greybody models for the MM-P sources).

The SED dichotomy is apparent in Fig. 7 where color-color diagrams are presented for MM-P and IR-P sources. For the present (MIPSGAL) and future (with Herschel) Far-IR surveys, the colors are computed at the wavelength characteristics of the MIPS (Spitzer), PACS and SPIRE (Herschel) photometric cameras. The figure shows that longward of $\sim$100 $\mu$m there is no clear difference between IR-P and MM-P sources. The dynamic range in the color distribution increases the shorter we go in wavelength below $\sim$100 $\mu$m (left panel). While the Far-IR portion of the SED is where the YSO emits the bulk of its energy, it is in the 20$\div$70 $\mu$m that the existence of two distinct classes of SED clearly stands out. This difference cannot be easily interpreted based on the SED shape alone, but it requires consideration of the absolute value of the SED as well, i.e. the luminosity.

In principle, the ability to use color-color diagrams as a rapid tool to classify IR-P and MM-P objects would require that MM-P and IR-S2 objects are spatially separated so that the contribution of the IR-S2 object does not contaminate the spectrum of the MM-P. In other words, in order to have an MM-P object with [24-70] color in the range ``predicted'' for MM-P objects we have to make sure that no contribution to the 24 $\mu$m flux comes form the nearby IR-S2 source. While the real spatial separation of millimeter and mid-IR counterparts in MM-P/IR-S2 combinations could be resolved by Spitzer/MIPS at 24 $\mu$m and Herschel/PACS at 70 $\mu$m, additional indications can be derived from a careful comparative analysis of the global SED from IR-P objects against a combined MM-P/IR-S2 SED (as it is seen at relatively poor spatial resolution, comparable to those of the data used in the present work).

Figure 8: Median average IRAS (asterisks) and MSX (diamonds) fluxes, normalised to the IRAS 100 $\mu$m flux. Data points for IR-P sources are connected with a solid black line; a black dashed line is used to connect points for the MM-P/IR-S2 combination. Dotted lines are used to connect the IRAS 25 $\mu$m and 12 $\mu$m fluxes.

Figure 8 reports the 100 $\mu$m-normalised median averaged IRAS (asterisks) and MSX (diamonds) fluxes. Data points for IR-P sources are connected by a full black line, while a dashed black line is used for the MM-P plus related IR-S2 sources. The dashed line is the SED we would assign to the same physical object, had we not performed the detailed model fitting analysis described in the sections above which instead suggests the presence of two distinct components; the comparison between the full and dashed black lines in Fig. 8 is suggestive of a two-component SED. However, the 8-21 $\mu$m SED portions have different slopes but their average is very similar between IR-P and MM-P/IR-S2; this is confirmed by the fact that the IRAS 12 $\mu$m point (the median of the fluxes normalised to the 100 $\mu$m fluxes) is identical for the two groups of sources. Based on the IRAS fluxes alone (the lines joining the asterisks, represented in dotted lines for the 12-25 $\mu$m portion), the comparison of the SEDs shows that the main difference is in the 25 $\mu$m flux. Thus, the contribution of the IR-S2 object appears relevant only shortward of $\sim$15 $\mu$m while the SED longward of 20 $\mu$m is dominated by the MM-P object, suggesting that the [24-70] color is sufficient to disentangle an IR-P from MM-P object also if the millimeter and mid-IR counterparts in the MM-P/IR-S2 combination are not spatially separated.

In a recent evolution of Whitney's radiative transfer model, Indebetouw et al. (2006) introduced envelope clumpiness as an additional parameter. As also intuitively reasonable, the general result is that more near and mid-IR is observable for an equally massive and extended envelope, than for a model with a homogenous envelope. It may be that the sources that cannot be fitted overall with the ZAMS models, and which then we break up into an MM-P/IR-S2 combination, could instead be fitted by a clumpy envelope model. We do not think this is the case for two reasons. The model assumes density contrast ratios of more than two orders of magnitude over a spatial scale of 2-5 parsecs (or of the order of a few arcminutes for an average distance of 5 kpc, typical of many of our targets); the inhomogeneities, probably smeared by optical depth effects below 100 $\mu$m, should be visible in the submillimeter and millimeter maps even with a single-dish facility and we do not find such evidence in our millimeter maps. The second reason is that the increase in radiation escaping the envelope due to clumpiness should be visible throughout the wavelength range shortward of $\sim$60 $\mu$m. Comparing the median SEDs of the IR-P sources with the combined SEDs of MM-P/IR-S2 (see Fig. 8) it is clear instead that the MM-P/IR-S2 combination has less flux at 25 $\mu$m, and more flux for $\lambda \leq15$ $\mu$m,than an IR-P. Since the latter is fitted by a homogeneous envelope, the clumpy model could account for the raise below 15 $\mu$m toward the MM-P/IR-S2 SED, but certainly not for its simultaneous suppression in the 20-50 $\mu$m region.

5 An evolutionary sequence for massive YSOs

The difference between the SED of IR and MM sources is suggestive of different evolutionary stages, but our results offer additional tools for a more quantitative analysis. A practical diagnostic tool to follow the evolution during the embedded phase of the YSOs is the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram, expressing the relationship between the mass of the YSO envelope and its bolometric luminosity, in which objects independently known to be in different stages occupy distinct areas in the diagram. Saraceno et al. (1996) used this tool to successfully verify the evolutionary sequence of low-mass YSOs from the Class 0 to the Class II phase, and we verify its applicability to the high-mass regime. Figure 9 presents the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram for all the sources in Table 6 with heavy black symbols; light grey symbols represent for comparison the sample of Saraceno et al. (1996) containing low-mass Class 0, I and II sources. For the sources of the present work, the MM sources are indicated with empty circles; IR-P sources are indicated with filled circles, while asterisks are used for IR-S1 when they are in a field dominated by an IR-P source, and crosses for IR-S2.

Figure 9: $L_{\rm\rm bol}{-}M_{\rm env}$ diagram for sources presented in Table 6. Lines and symbols in light grey are from Saraceno et al. (1996) and depict the situation for the low-mass regime (Class 0, I and II sources are represented respectively by open circles, filled circles and crosses, with the full thin line representing the log-log linear fit to the Class I source distribution). The heavy black lines and symbols represent the sources in the present work. MM sources are represented by open circles (smaller circles for MM-S objects) while IR-P, IR-S1 and IR-S2 are indicated by filled circles, asterisks and crosses, respectively. Lines represent the models as described in the text; limited to the high-mass regime, the dot-marked full thick lines represent the accelerating accretion phase (the dots mark 104 yr intervals), the grey line is the portion of the pre-MS evolution, while the dashed diamond-marked portion (at 105 yr intervals) represents the final envelope clean-up phase. The different models (left to right) are for different initial envelope masses of 80, 140, 350, 700 and 2000 $M_{\odot }$. The full and dashed black lines are the best log-log linear fit to the IR-P and MM-P sources, respectively.

Sources that our SED analysis identifies as different based on the 10 $\mu$m-1 mm SED appearance tend to occupy different regions in the diagram. The IR-P sources define a strip which is the extension to higher luminosities and masses of the trend defined by the typical Class I source at the low-mass end of the diagram (source mass was computed from millimeter fluxes of Saraceno et al. (1996) using the same opacity we adopt in the present paper). A linear log-log fit to the distribution of IR-P objects (the heavy black line) lies at nearly the exact prolungation of the Class I analogous fit (the thin full line), with a slope of $\sim$1.27. The IR-S objects are situated to the left of the ``Class I'' strip, in analogy with the low-mass Class II objects, with the IR-S2 further to the left, on average, than the IR-S1. Below the IR-P strip we find all the MM sources, which can be considered analogues of the typical Class 0 of the low-mass regime. A linear log-log fit yields a slope $\sim$1.13, almost identical to the IR-P fit, but lying $\sim$0.6 dex (or a factor $\sim$4) below it in luminosity.

The clustered nature of star formation around massive YSOs has an impact on our analysis, in the sense that the objects we have modeled have been measured with inadequate spatial resolution both in the mid-IR and in the millimeter, and may have been erroneously assumed to be single objects. Since in principle this would be equally true for MM-P and IR-P objects, however, their relative locations in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram would not be altered.

5.1 A simple model for YSO evolution

The diagram in Fig. 9 confirms the suggestion of an evolutionary sequence coming from the analysis of the SED differences. In order to verify whether a physical model of an accreting massive YSO can explain the different populations of objects seen in the diagram, we make the assumption that the process is a scaled-up version of the traditional inside-out collapse scenario which is consolidated for the formation of low-mass stars. Convincing evidence has been collected, e.g., from the observations of outflows and disks, that accretion plays a dominant role also in the formation of massive stars (Churchwell 1997; Zhang et al. 2001; Beuther et al. 2002b). We will then adopt the model of collapse in turbulence-supported cores (McKee & Tan 2003) to predict the behaviour of luminosity and circumstellar mass as a function of time in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram.

5.1.1 The accelerating accretion phase

An important feature in these models is that accretion is a time-dependent process; its rate increases with time because the inside-out collapse wave expands, reaching out to material characterised by supersonic speeds and where higher accretion rates are possible. McKee & Tan (2003) offer a convenient analytical prescription for the accretion rate as a function of time. In particular, for a standard free-fall envelope density profile the accretion rate is expressed as

$$\dot{m}_{\star}=4.6\times 10^{-4} \left( {{ m_{\star \rm f} }... ...ar} }\over{ m_{\rm\star f} }}\right) ^{0.5}~M_{\odot}/{\rm yr}$$ (4)

where $m_{\star}$ and $m_{\rm\star f}$ are the instantaneous and the final stellar mass; $\Sigma _{\rm cl}$ is the core surface density in g cm^-2. Once $m_{\rm\star f}$ and  $\Sigma _{\rm cl}$ are fixed it is easy to compute the evolution of $m_{\star}$ as a function of time; in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram we will follow the envelope mass which, from a given initial value $M_{{\rm env}_i}$ will decrease over time both because of the material accreting onto the star at a rate  $\dot{m}_{\star}$, and because of the material expelled by the powerful molecular outflows which characterize the earliest phases of massive YSO formation. The fraction of accreting material that is expelled in the outflow is of the order of $f_{\rm out}\sim0.5$ (Matzner & McKee 2000; McKee & Tan 2003), although even higher fractions have been proposed (Churchwell 1997) to justify the large ratio between the masses of the flow and of the driving protostar observed in high-mass YSOs. An additional effect contributing to the clump envelope dispersal is the fact that massive YSOs form in clusters, and so there will be other less massive YSOs in the same cluster that will drain material from the molecular clump. Faustini et al. (subm.) have detected young stellar clusters toward many of the objects included in the present work, and a relationship is found between the mass of the most massive star in the cluster and the total stellar mass

$${\rm Log}~ [m_{\rm\star Tot}] = 0.55 + 1.41 ~ {\rm Log}~ [m_{\rm\star Max}].$$ (5)

Assuming $m_{\rm\star Max}=m_{\rm\star f}$, for each choice of $m_{\rm\star f}$ we derive an estimate of the total mass bound in stars; adopting a mean formation time of $t_{\rm cl}\sim 3\times 10^6$ years for the clusters (Faustini et al. subm.), we can define a time-averaged rate for clump dispersal into cluster stars (except the most massive one), as

$$\dot m_{\rm cl} \sim {{(m_{\rm\star Tot}/f_{\rm out}) - m_{\rm\star f}}\over {t_{\rm cl}}}$$ (6)

where the assumption that a $f_{\rm out}$ fraction of accreting mass is lost via outflows is adopted also here. The evolution of the clump mass will then be

$$M_{\rm env}(t) \sim M_{{\rm env_i}} - \left[ t~ \left( {\dot m_{\star}} + \dot m_{\rm cl}\right) \right].$$ (7)

Concerning the total luminosity emitted by the star, McKee & Tan (2003) plot in their Fig. 6 the luminosity as a function of stellar mass, taking into account the contributions from accretion, deuterium burning and core contraction, for different values of $m_{\rm\star f}$ and $\Sigma _{\rm cl}$. All tracks in their figure were digitized and interpolated to extract the luminosity for any given stellar mass in our simulations. The relevant free parameters in modelling the evolution of the massive YSO in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram are the initial mass of the core envelope and the final mass of the star $m_{\rm\star f}$ where we assume that the main accretion stops; this part corresponds to the black-dotted full line section of the evolutive tracks reported in Fig. 9 for different initial envelope masses reported in Col. 2 of Table 7.

 

 

Table 7: Parameters for the evolutive tracks.

Track	$M_{{\rm env}_i}$	$m_{\rm\star f}$	$t_{\rm\star f}$	$\dot m_{\rm\star f}$
#	($M_{\odot }$)	($M_{\odot }$)	(105 yr)	(10-4 $M_{\odot }$/yr)
1	80	6.5	4.5	0.6
2	140	8	3.7	0.9
3	350	13.5	2.7	1.9
4	700	18	2.1	3.4
5	2000	35	1.5	9

The tracks all start at very low luminosity; as the stellar mass grows the tracks cross the area occupied by the MM sources and proceed almost vertically toward the IR-P sources. Important at this point is the mass $m_{\rm\star f}$ for which we choose to terminate the main process of accelerating accretion according to Eq. (4): for the tracks in Fig. 9 we adopted the $m_{\rm\star f}$ values in Col. 3 of Table 7 to simulate a situation where the stellar formation is stopped at the position corresponding to the IR-P sources. This is clearly arbitrary, and alternative possibilities will be discussed later. Looking at the dot frequency along the tracks (each dot represents a span of 10⁴ years) we see that evolution speeds up as we go toward the rightmost track not only because of the higher  $m_{\rm\star f}$ which justifies a higher accretion rate according to Eq. (4), but also because of higher clump surface densities $\Sigma _{\rm cl}$. Indeed we assign $\Sigma _{\rm cl}$ using a relationship which we find in our MM sources between the core mass and the column density, determined by dividing the masses in Table 6 by the core extensions derived from the fits to the millimeter cores as explained in Sect. 3.1.1:

$${\rm Log}~ (\Sigma _{\rm cl}) = -1.8 + 0.64 ~ {\rm Log}~ (M_{{\rm env_i}}).$$ (8)

Based on our choice of the initial envelope masses we derived, following Eq. (8), $\Sigma _{\rm cl}=[0.3$, 0.4, 0.7, 1.1, 2.2] g cm^-3, in good agreement with the range of values measured for Galactic cores as reported by Plume et al. (1997). Columns 4 and 5 of Table 7 report the times by which the final masses $m_{\rm\star f}$ are reached in the simulations of Fig. 9, and the values of the accretion rate at that time. These rates have significant impact on the mass at the arrival onto the ZAMS. Palla & Stahler (1992) showed that the accretion flow has the important effect of delaying the time of arrival onto the Main Sequence; a YSO that would normally join the MS at 8 $M_{\odot }$ for $\dot{M}$ $_{\rm acc} =10^{-5}$ $M_{\odot }$ yr^-1, would instead join the ZAMS for M=15 $M_{\odot }$ if $\dot{M}$ were 10 times higher. Rates even one order of magnitude higher are usually postulated based on the inferred mass loss rates measured for molecular outflows in high-mass YSOs (Zhang et al. 2001, 2005); under these conditions the arrival onto the ZAMS can be deferred to higher masses provided the required ever higher accretion rates can be maintained. Using the results of Palla & Stahler (1992), although extrapolated to higher $\dot{M}$ values with respect to the range they explored, the above values of $\dot m_{\rm\star f}$ imply ZAMS masses of $M_{\rm zams}=[12.0$, 14.0, 18.7, 22.8, 33.0] $M_{\odot }$ for the five tracks, the last three within 30% of the desired final stellar masses.

5.1.2 The envelope clean-up phase

When $m_{\star}=m_{\rm\star f}$ the accretion that has been driving the protostar assembly according to Eq. (4) will stop; our choice of $m_{\rm\star f}$ is clearly arbitrary at that point, and in turbulent core models this corresponds to setting the core radius and density structure. At this time the YSO is already on the ZAMS, or very close to it, if sufficiently massive (the rightmost three tracks in Fig. 9), or it will approach it following standard PMS tracks (Palla & Stahler 1993; Behrend & Maeder 2001 [BM01]) if not. This early ZAMS phase corresponds to the stages observationally identified with Hot Cores and UCHII regions (the IR-P sources are mostly High sources with UCHII-like FIR colors, see Sect. 6.4). The very high accretion rates that we find at the end of this ascending phase in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram are compatible with the rates deduced from observations of outflows emanating from HC and UCHII regions (Beuther et al. 2002b; Zhang et al. 2005). They are also compatible with the rates required to choke-off or confine the expanding HII regions, as required to explain the well known finding that the lifetime of Hot Cores and UCHII regions as estimated from number statistics ($\sim$ $3\times 10^5$ years; Wood & Churchwell 1989) is much higher than the dynamical timescales for the expansion of HII regions after the star reaches the ZAMS. A confining/quenching mechanism is needed, and accretion is very efficient in that respect; the required value is of the order of 10^-5 $M_{\odot }$ yr^-1 for an O9 ZAMS star (adopting the stellar parameters from Thompson 1984), and decreases to a few 10^-6 $M_{\odot }$ yr^-1 for a B0 (Walmsley 1995).

After the end of the accelerating accretion phase, the evolution of the clump mass will be qualitatively similar to the first phase, and it will be dominated again by the material expelled in the molecular flows (at rates = 0.5  $\dot{m}_{\star}$) and drained by accretion onto other objects forming in the same region (see above). We neglect any additional contribution to the envelope dispersal from the expanding HII region.

Molecular outflows are known to be present down to the Ae/Be stage and since this requires disk accretion to feed the flow engine, the formed YSOs will continue to gain mass; the mean values observed for the mass loss rate $\dot{M}$ $_{\rm out}$ in flows originating from Herbig Ae/Be stars (e.g. Lorenzetti et al. 1999) is of the order of several 10^-7 $M_{\odot }$ yr^-1, and hence we will assume a nominal accretion rate of $\dot M_{\rm acc} = 10^{-6}$ $M_{\odot }$ yr^-1. A mass outflow rate decreasing with time would be more realistic, but we do not know the temporal evolution of this quantity and we prefer to keep it constant and accept that the evolutionary tracks immediately following the end of the accelerating accretion phase may be incorrect.

The luminosity radiated can be expressed as $L=L_{\star} + L_{\rm acc}$, where $L_{\star}$ is either the luminosity radiated during the approach to the ZAMS along standard PMS tracks (if the object is not already on the ZAMS), or the ZAMS luminosity. In our simulations we find that objects with $m_{\rm\star f}\la 18$ $M_{\odot }$ do not reach the ZAMS at the end of accelerating accretion phase; these objects will have a PMS phase and a well defined luminosity evolution (Palla & Stahler 1993, BM01). We adopt the tabulated tracks of BM01 for the mass value which is closer to $m_{\rm\star f_i}$; this is a first approximation, the second being that in this phase we neglect the fact that residual accretion increases the stellar mass beyond the values $m_{\rm\star f_i}$. When the objects are on the ZAMS, at the end of the PMS track or immediately at the end of the accelerating accretion, we will assign a luminosity according to the models of Bressan et al. (1993) using the analytical representation from Binney & Merrifield (1998). The accretion luminosity is computed using Eq. (74) in McKee & Tan (2003). We note that both $L_{\star}$ and $L_{\rm acc}$ will vary with time because both depend on the stellar mass, which will increase with time due to the persisting residual accretion.

This second phase of evolution described above corresponds to the dashed section of the evolutive tracks in Fig. 9, and spans the regions between IR-P and IR-S sources; the diamond symbols mark periods of 10⁵ years. This simplified model does not have the ambition to provide strict quantititative predictions on the early massive YSO evolution, but it is sufficient to qualitatively demonstrates that the different types of sources identified in this study can plausibly be part of a consistent evolutive context. The simulations are stopped when the clump has disappeared and the objects are optically visible; the timescales for the entire evolutive path for the different tracks in Fig. 9 are $t_{\rm final}\sim [2.1, 2.7, 3.5, 4.8, 4.5]\times 10^6$ years.

5.1.3 Caveats

We want to draw attention to a few specific aspects. The drop in luminosity at the end of the accelerating accretion phase is likely not real. PMS tracks start from the birthline, and its location in the H-R diagram shifts toward higher luminosities as a function of the accretion rate; BM01 tracks are computed for an accretion rate that stays at a very low level until $m_{\star}\sim 5$ $M_{\odot }$, and then rises much more quickly than what we assumed in Eq. (4). In our simulations the accretion rate is always higher than what is assumed in BM01 for the lower $M_{{\rm env}_i}$ runs, implying that the PMS section of the tracks should be considered as a lower limit. Besides, the stellar mass will increase in the PMS phase due to the residual accretion; this would imply changing the PMS track as a function of time, and we prefer to neglect this effect.

The timescales for the evolutive tracks at the end of the envelope clean-up phase are likely to be incorrect. It is not likely that the envelope dispersal agents have a constant efficiency with time. For example we have assumed a steady outflow mass loss rate in this phase, at values compatible with mean rates for Herbig Ae/Be stars (10^-6 $M_{\odot }$ yr^-1). It is more realistic to assume that, for example, the mass loss rate are still quite high at the beginning of this phase and then decrease with time; this would modify the timings, speeding-up the evolution in the portion closer to the IR-P strip. In the absence of an observational constraint, a prescription in this sense would be a free parameter. Likewise, we have neglected any contribution to the envelope dispersal that may come from the expanding HII region. This would result from a complex interplay between the expanding HII regions (the IR-P and IR-S are assumed to be on the ZAMS in these simulations) and the on-going accretion, and there are no estimates for the timescales and effectiveness of this process.

6 Discussion

   
6.1 Source statistics and formation timeline

The evolution of YSOs in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram plausibly confirms the working assumption that MM objects, along a path of rapidly rising luminosity, evolve toward the IR-P stage where the accelerating accretion stops and the YSO is already on the ZAMS if sufficiently massive, or approaches the ZAMS along PMS tracks.

From an observational point of view, our analysis presented in the last section proposes that the strip identified in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram by the IR-P sources may be the locus that marks the end of the main accretion phase. From a statistical viewpoint, the distribution of IR-P objects in a confined region would require that they spend a relatively long time there. The timescales for the different portions of our evolutionary tracks support this; considering that dots and diamonds in Fig. 9 corresponds to steps of 10⁴ and 10⁵ years respectively, an evolving object spends the larger fraction of its time in the IR-P phase. The subsequent evolution in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram is carried out in a relatively shorter time (the distance between diamond symbols increases) and this is also in agreement with the much sparser distribution of IR-S1/S2 objects. No such consideration can be made for MM objects, because we are much more likely to suffer from sensitivity effects which prevent detection of objects with low luminosity (for any given envelope mass). Indeed, an object in the earlier portions of the evolutionary tracks would have a similar envelope mass to the MM objects we detect, but a much lower luminosity that implies lower temperatures or extraordinarily high optical depth; in both cases the obvious consequence from an observational viewpoint is that the Far-IR portion of the SED in such an object would be considerably depressed, and because of our 60 $\mu$m flux threshold in the selection criteria (Palla et al. 1991) it would not be included in our sample. Thus, it is likely that the number of MM-P sources is underestimated. Unbiased, non flux-limited far-IR and submm surveys are needed to extend the source coverage to equally massive but lower luminosity objects; the InfraRed Dark Clouds revealed by MSX (Egan et al. 1998; Carey et al. 2000; Simon et al. 2006) are very promising objects to investigate for the earlier stages of massive star formation.

Out of the 42 fields examined in this work, 15 have an MM-P source and 27 have an IR-P source; taking into account the different number of High (16) and Low (26) fields included in the present sample, and the statistics of MM-P/IR-P occurrence in both groups (see Sect. 6.4), the occurrence of an IR-P object is a factor $\sim$2.3 more likely than that of an MM-P object. If, as we have been assuming throughout our analysis and interpretation, the IR-P sources are representative of the UCHII phase with an expected timescale of $\sim$10⁵ years (Churchwell 2002), then the duration of the phase traced by the MM-P sources we detect should be of the order of $\sim$ $5\times 10^4$ years. This is on average consistent (it varies with the mass) with the timescales that can be deduced from Fig. 9 for the transition from MM-P to IR-P stage. If the formation timeline represented by model tracks in Fig. 9 is correct, we should expect at least a factor of 2 more MM-P type objects with respect to the population we detect. This may be consistent with the fact, noted above, that the sensitivity bias toward luminous sources probably leads to an understimate of the number of MM-P sources.

Figure 10: Left panel: relationship between the initial envelope mass and $m_{\rm\star f}$, the stellar mass at the end of the accelerating accretion phase. Right panel: relationship between $m_{\rm\star f}$ and the accretion rate at the end of accretion phase (at $m=m_{\rm\star f}$). The asterisks represent the values for each of the tracks in Fig. 9; the thin lines are log-log linear fits whose analytical expressions are given in Eqs. (10) and (9) respectively.

The existence of massive and relatively low luminosity objects can only be established with multiband photometric observations longward of 100 $\mu$m; the availability of Herschel, Laboca@APEX, SCUBA2@JCMT and ALMA in the near future will clarify this issue.

   
6.2 The IR-P strip as the ZAMS location in the L_bol-M_env diagram

In our simulations, the IR-P strip is the location where the accelerating accretion phase ends. In the right panel of Fig. 10 we plot, for each of the tracks shown in Fig. 9, the accretion rate at the end of the accelerating accretion phase ( $t=t_{\rm\star f}$) as a function of the stellar mass at that point of the evolution ( $m=m_{\rm\star f}$). The linear log-log appearance of the IR-P strip in Fig. 9 translates into a similar feature in Fig. 10, and the first order fit (the gray line in the figure) provides

$${\rm Log}~ (\dot{m}_{\rm f}) = -5.5 + 1.6~ {\rm Log}~(m_{\rm\star f}).$$ (9)

The slope of this power law is close to the value of 1.5 suggested by Norberg & Maeder (2000) as the one that best matches the distribution of intermediate and high-mass YSOs in the H-R diagram. In that study the authors concentrate on relatively evolved objects (like Herbig Ae/Be stars) for which an estimate of the effective temperature is possible, and it is interesting that we obtain a similar result for objects still heavily embedded in their parental cocoon. Norberg & Maeder also emphasize the remarkable similarity of their adopted power law slope with the one (1.54) derived from the relationship (Churchwell 1997) between YSO bolometric luminosity and the mass loss rates measured for the related molecular outflows, under the assumptions that the mass loss is proportional to the mass accretion and that the objects are on the ZAMS.

There is however a fundamental difference. Norberg & Maeder (2000) use their parametrization of the accretion rate throughout the accretion process, while our Eq. (9) is valid only at the end of the main accretion phase when the object is very close to or on the ZAMS. As already mentioned in the previous section, however, the dependence of $\dot{m}_{\star}$ as expressed in Eq. (4) for the turbulent core model is very different from that in Norberg & Maeder. On one side the power law in Eq. (4) is much less steep as a function of the instantaneous stellar mass, and on the other hand the absolute value of the power law scales, quite realistically, with global clump parameters. The difference in power law exponent for the expression of the instaneous accretion rate has a dramatic impact on the evolution of a YSO in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram. If we adopt the Norberg & Maeder (2000) prescription in the simulations of Fig. 9 we find that the YSO would spend almost all its accretion time at a luminosity below $\sim$100 $L_{\odot}$ and then cover the final portion of the ascending tracks, between the MM and the IR-P sources, in just few 10³ years or about ten times faster than using the turbulent core formulation (independently of the initial envelope mass, i.e. for all tracks). We can rule out this possibility based on the statistics of IR-P and MM-P objects (Sect. 6.1), and considering the much longer duration of the envelope clean-up phase with respect to the accelerating accretion phase.

The consistency of Eq. (9) with the $\dot M\propto M^{1.54}$ prescription derived (Norberg & Maeder 2000; Behrend & Maeder 2001) using the $\dot M_{\rm out}{-}L_{\rm bol}$ relationship for YSOs on the ZAMS (Churchwell 1997) strenghtens our proposal that the IR-P strip in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram represents the ZAMS. At the same time the timescales implied by the distribution in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram of YSOs in different evolutionary stages suggests that that prescription is a satisfactory representation for the ZAMS phase only.

A further consideration concerns the $m_{\rm\star f}$ mass values where we stop the accelerating accretion phase (i.e. the end of the full line portion of the simulation tracks). Instead of stopping the accelerating accretion in the position corresponding to the IR-P sources, we could have done it in correspondence of the MM sources. In this case we find that the stellar mass is always at least a factor $\sim$2 below the ZAMS mass appropriate for the instantaneous accretion rate. Similar to the standard case, (Fig. 9) the PMS portion of the tracks would be well confined in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram which means that the objects should reach the ZAMS (and then start to trigger the formation of the HII region) while still well in the MM-P region. After the end of the main accretion the evolution from MM-P to IR-P objects would not only be due to an increase of stellar mass and luminosity caused by the residual accretion, but also to the envelope clearing process. However, at this stage the accretion runs at a much slower rate, and the duration of this process would be of the order of few 10⁵ years. Since the IR-P sources are mostly of the High type and, as such, typically associated with UCHII regions, these very long timescales seem unrealistic compared to the typical lifetime of the latter objects. Besides, the low values for residual accretion after $t=t_{\rm\star f}$ would not justify the high $\dot{M}$ observed for the molecular outflows emanating from these objects.

The opposite possibility that the main accretion can be extended to later times with respect to the standard scenario seems to be excluded by the fact that the MM-P and IR-P sources are distributed in a rather well-confined strip of the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram, and the tracks would end in a scarcely populated area of the plot. Since the accelerating accretion runs much faster than the subsequent residual accretion phase, it is quite unlikely that the accelerating accretion can be prolonged beyond the IR-P region.

6.3 Turbulent cores or competitive accretion?

The formulation we adopted for the various quantities used in our simple model is taken mostly from the turbulent core model of McKee & Tan (2003). Its application to a variety of sources allows verification of its consistency throughout a range of masses and luminosities. The left panel of Fig. 10 presents the relationship between the initial envelope mass and $m_{\rm\star f}$ for the various tracks of Fig. 9. The linear log-log fit provides the relationship:

$${\rm Log}~ (M_{{\rm env}}) = 0.4 + 1.9~ {\rm Log}~(m_{\rm\star f}).$$ (10)

In turbulent core models it is assumed that a significant fraction ($\sim$50%) of the mass contained in the core will end up forming the star. The mass traced by the millimeter emission in our MM and IR-P sources is more than one order of magnitude higher than the final stellar mass, clearly implying that the mass participating in the formation of the most massive star is a minor fraction of the mass we trace in our observations. The ratio between the stellar mass on the IR-P line (where the accelerating accretion is assumed to stop) and the initial adopted envelope mass decreases from $\sim$7% to $\sim$2% with increasing initial envelope mass. These numbers are very similar to the star formation efficiencies estimated on average across the Galaxy (Williams & McKee 1997), and much smaller than assumed in the turbulent core models. An obvious explanation is that the structures we trace in the millimeter continuum are not the ``cores'' but rather the ``clumps'' in McKee & Tan terminology. The decrease in efficiency as a function of initial envelope mass can also be understood considering that the clumps we detect do not form a single star, but will give origin to many more stars likely of lesser mass than the most massive one whose evolution we model in our simulations. This follows directly from the clear evidence that massive YSOs are indeed forming in clusters with other stars, and this has been included in our simple model (see Eq. (6)) allowing a contribution of other accreting YSOs to the clean-up of the clump during the massive YSO evolution. The mass of the hosting clump is related not only to the mass of the most massive star forming in the cluster, but also to the total number of stars populating the cluster. In other words more massive clumps will generate more populous clusters, implying that the relative fraction of mass going into the most massive star will be lower. Using Eq. (5) we can guess the total mass in stars for a given values of $m_{\rm\star f}$; dividing the total stellar mass by the initial envelope mass we find efficiencies of the order of 50% varying by a factor of two at most across the mass range (compared to an almost factor of 4 considering only the most massive star).

In competitive accretion models (Bonnell et al. 2003) the star initially forms by normal collapse of gravitationally bound material; soon the increasing number of stars in the young forming cluster will deepen the potential well of the system, draining material from the outer regions of the clump toward the center. Stars in the center of the cluster will compete for this material, and the winner will end up being the most massive star in the system. Since most of the stellar mass will be accreted from initially unbound clump material flowing inward under the action of gas and stellar gravity, the efficiency of accretion, and hence the accretion rate, will be higher the lower the velocity difference between the accreting star and the inflowing material (Bonnell & Bate 2006). The gas with lower velocity is located near the center of the clump, so that the stars initially forming at the clump center will be in the best conditions to accrete mass at high initial rates; the star at that point is massive enough to win the competition for the mass flowing from the outer clump regions, although at lower rates given the higher relative velocity of this gas. The result is then that a significant fraction of the mass of a star is accreted early, when the star is surrounded by low relative velocity gas. The simulations reported in Bonnell et al. (2004) indeed show that 60% of the final stellar mass is accreted in the first half of the formation time; on the contrary, in our simulations a 28 $M_{\odot }$ star accretes only 25% of its mass in the same fraction of time. Adopting this different prescription for the accretion rate would modify the time tick spacing in our evolutive tracks, increasing their density toward the end of the ascending tracks; the consequence is then that we should expect fewer MM-P objects in the low luminosity portions of the tracks compared to the turbulent core prescriptions used to generate the tracks in Fig. 9). As already anticipated (see Sect. 6.1), this issue cannot be resolved with the present data and will have to await next generation Far-IR/submm facilities (Herschel and ALMA).

Another consequence of the mechanism of competitive accretion is that the final mass of a star is unrelated to the mass of the clump in which it was initially found (Bonnell et al. 2004). The existence of the IR-P strip establishes instead a clear relationship between these two quantities (Fig. 10 left). However, the term ``clump'' in competitive accretion models is used to denote a region surrounding a forming star where the density profile is monotonically decreasing, and is generally much smaller than the size of the clumps in McKee & Tan's terminology. Yet Bonnell et al.'s simulations show that the majority of the mass that finally ends up in a 30 $M_{\odot }$ star is located more than 0.2 pc away from the star's position at the time of initial formation, while the deconvolved radii of most of the millimeter cores of MM-P objects are below 0.25 pc. We then conclude, although with caution, that the result of Figs. 9 and 10 left are not in agreement with the predictions from competitive accretion.

   
6.4 The root of the High/Low dichotomy

We can now ask the question whether the presence of an MM-P or IR-P source in a field may be associated with the High/Low nature of the fields themselves, mainly based on the [25-12 $]\lessgtr 0.57$ threshold. From the discussion of Fig. 8 it seems clear that the main difference between MM-P and IR-P objects is in their 25 $\mu$m flux relative to the 60 and 100 $\mu$m flux. Is the difference in average 25 $\mu$m flux between IR-P and MM-P sufficient to justify the [25-12] IRAS colour distribution between High and Low sources? Looking at the dotted lines in Fig. 8, this may be the case. To quantify this we plot in Fig. 11 the histogram of the fields associated with an MM-P source (full line) or an IR-P source (dashed line), as a function of the [25-12] IRAS color.

Figure 11: Distribution of fields with an MM-P (full line) and with an IR-P (dashed line) source as a function of the [25-12] IRAS color. The dash-dotted line represents the threshold value of the [25-12] color to distinguish between High and Low sources.

There is clearly no threshold-like difference between the two histograms, as could be expected because the High/Low groups, as in all classifications, should not be considered closed boxes. A difference in the distributions is however apparent, and shows that as we move toward higher [25-12] color (going from Low to High sources) the fraction of fields with IR-P sources increases significantly: the MM-P source number ratio is 54% for Low sources and only 6% for High sources. Not only does this result confirm again that the Low group contains a higher fraction of younger sources (the MM-P ones), but it suggests that the fact that the most massive YSO in a field is an IR-P or an MM-P source is indeed at the root of the difference itself between High and Low sources. The lower [25-12] value in Low sources then results from the combination of the steeper 24-70 $\mu$m SED of the MM-P compared to the IR-P, and the close proximity of the MM-P to an IR-S2 source. Clearly, MM-P sources without a close-by IR-S2 source will exist but they would not emit at 12 $\mu$m (e.g. Fig. 1a) and hence they would have not been included in our original sample that was selected using all 4 IRAS bands.

At this point a final consideration is in order concerning the selection criteria that have been applied to extract source samples from the IRAS catalogue in a variety of studies similar to the one we have been carrying out in the past decade. The selection criteria devised by Wood & Churchwell (1989) to identify sources associated with UCHII regions are at the foundation of several searches for high-mass protostellar candidates undertaken in the past decade, like Beuther et al. (2002a), Faundez et al. (2004), Hill et al. (2005); in this sense these samples are analogous to our High sources. Very young objects not associated with radio emission and with continuum emission detected only longward of 20 $\mu$m are found in our High group as well as in the samples of the above-mentioned authors. However, our study is different from all of the above because we are the only ones to include in the investigation the Low sources, and our analysis shows that the frequency of occurrence of the MM-P sources is several times higher in this group. As we have shown above, the possibility of revealing MM-P sources in IRAS selected samples seems to be at least partially due to their fortuitous proximity to a lower mass IR-S object (the IR-S2 source in our analysis). Nonetheless the conclusion remains that the Low group contains many more young sources compared to the High group and, likewise, to the various samples studied by other groups who only consider High-like sources.

6.5 A unified scenario

An important aspect to be considered in building a consistent scenario to interpret our results is the association of the dominant components with radio sources. Radio continuum in massive YSOs marks the presence of a developing HII region and therefore identifies sources where the ZAMS has already been reached. Figure 12 reports the distribution of the distances between the dominant source (highest circumstellar mass) in each field and the nearest detected radio continuum source. The radio data are taken from Molinari et al. (1998a, 2002) for sources in the northern sky, and from the Parkes 6cm survey (Wright et al. 1994) for the southern sky.

Figure 12: Distribution of fields with an MM-P (full line) and with an IR-P (dashed line) source as a function of the distance (in pc) between the dominant source in each field and the nearest radio continuum source.

Figure 13: Sketch of the proposed model for the evolution of star formation in massive star forming regions. The prototypical MM-P and IR-P fields are on the left and on the right, respectively.

The distribution of source-radio distances for IR-P sources shows a skewed distribution peaking at 0 distance, while MM-P sources are more uniformly spread over a range of separations. There are fields (mostly in the southern sky) where no radio sources were detected (the sensitivity of the Parkes 6cm survey is only 35 mJy, much higher than the $\sim$sub-mJy level reached in the dedicated VLA observations done for the northern sources); these fields however are more or less equally distributed between fields with IR-P and MM-P sources. Limiting ourselves to the northern sources for which the data have arcsec spatial resolution, the images and the annotations in the appendix show that out of 16 fields no radio sources are detected in # 38, 59, 116, 118 and 139 within 2 arcmin of the core peak; within the remaining fields 5 out of 6 sources classified as IR-P (# 8, 12, 136, 143 and 148) have an associated radio continuum source, while 4 out of 5 sources classified as MM-P (# 75, 98, 151 and 160) do not. In the latter cases the radio sources can be associated to other sources in the vicinity and visible either at 8 $\mu$m or in the near-IR.

Most (35/42) of the fields examined in this work are associated with young stellar clusters or bright single stars visibile in the $K_{\rm s}$ band; there is no difference between fields with an IR-P or MM-P source. For the sources taken from Molinari et al. (1996) the cluster association is established via dedicated near-IR observations and analysis (Faustini et al., subm.); for the sources taken from Fontani et al. (2005) we visually inspected the 2MASS images, and the clusters are clearly visible without the need for detection algorithms.

We then conclude that the differences outlined at the beginning of this paper between IRAS 23385+6053 and IRAS 20126+4104 can be assumed to be typical of two distinct stages in the evolution of a massive YSO, which our results and analysis codify in the MM-P and IR-P types of sources. Considering the evidence in its entirety we can draw a representative schematic in Fig. 13 which condenses the basic differences between the two types of sources.

The massive forming objects in all the examined fields are very young; the models that we find to satisfactorily describe our results speak of formation times that will depend on the mass of the star and the properties of the hosting clump but that are of the order of a few 10⁵ years at most. MM-P objects are likely completing their luminosity rise that will bring them to be IR-P sources, on the ZAMS, with properties typical of Hot Cores and UCHII regions; by analogy (see also Fig. 9) they represent the high-mass counterparts of the transition from Class 0 to Class I objects. Given the shape of the IMF it is reasonable to expect that the stellar cluster will also contain intermediate mass stars able to drive UCHII regions, and these may be at the origin of the radio emission that we find in the vicinity of MM-P sources, but not associated with them. This is particularly clear in cases like IRAS 05345+3157 (Sect. A.2) or IRAS 23385+6053 (A.42) where the clusters are more prominent and the radio continuum data are of better quality.

The facts that the gas clumps where massive stars form appear to be virialised, the spherical symmetry (on average) of the associated stellar clusters and their mean ages, suggest that these regions are dynamically relaxed systems existing for several dynamical timescales. The association of these objects with conspicuous young stellar clusters where the low mass members have median ages of a few 10⁶ years confirms that the most massive YSO in a star forming region is the last one to form. This cannot be explained presently by competitive accretion models; however, given that the massive star appears to characterize the last few 10⁵ years of the life of the hosting clump, the question of what mechanism prevents the massive star from forming earlier needs to be addressed by any theoretical model.

7 Summary

The 8-1200 $\mu$m SED of YSOs in 42 regions (26 Low and 16 High sources) of massive star formation has been constructed using MSX, IRAS, and submm data, newly obtained or available from previous works; the aim was to verify and quantify a possible evolutionary scheme for massive star formation based on the shape of the SEDs and the $L_{\rm\rm bol}{-}M_{\rm env}$ relationship, which are tools that proved quite effective in low mass star formation studies. The YSO SEDs were fitted to an extensive grid of envelope models with embedded ZAMS star, available from the literature. Satisfactory fits were found for a fraction of sources, which we called ``IR''; for all the others we assumed that the mid-IR and the Far-IR/submm portions of the SED come from different physical sources which happen to be positionally associated. A two-component fitting was then done, again using an embedded ZAMS model for the mid-IR emitting source (again called ``IR''), and a single-temperature optically thin greybody for the longer wavelength emitting component, called ``MM''. Sources where then subclassified as ``-P'' if they have the most massive envelope in a particular star forming region (so there can only be one ``-P'' source in each region), or ``-S'' if they are not.

The main conclusions may be summarised as follows:

-

Out of the 42 examined fields, 27 have an IR-P object and 15 have an MM-P object. MM-S and IR-S are also found in many fields. The SEDs of IR-P and MM-P objects mostly differ shortward of 100 $\mu$m, with the IR-P emitting all along the mid-IR and the MM-P only longward of 60 $\mu$m. In particular, the MIPS [24-70] colour is the most sensitive being $\sim$1 for IR-P sources and $\sim$4 for MM-P sources.

-

Similarly to the low-mass regime, we find that the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram is a valuable tool to follow the pre-MS evolution of massive YSOs in their embedded phase. The distribution of massive YSOs in the diagram clearly suggests that MM-P, IR-P and IR-S objects are the high-mass analogues of Class 0, I and II sources in the classical classification established for low-mass YSOs.

-

Using a simple model based on the prescriptions from the turbulent core model (McKee & Tan 2003), we show that MM-P sources preceed IR-P sources in the evolution toward the ZAMS. The distribution of the different classes of objects with respect to these evolutionary tracks suggests that the IR-P stage may correspond to the arrival of the YSO on the ZAMS. The MM-P objects would then represent an earlier stage.

-

The incidence of MM-P occurrence is 54% among Low sources, and only 6% in High sources. The Low sources are once again confirmed as the group with the highest incidence of the youngest massive protostellar candidates. Our sample is the only one to include such objects, while the other best studied HMPO candidates only include High-like objects. We verified that the SED shape difference between MM-P and IR-P objects shortward of 60 $\mu$m explains very well the color difference between the two classes of sources and it is then at the root of the High/Low division.

-

We find that our results can be best explained within the framework of the turbulent core model. In particular the existence of the IR-P strip in the $L_{\rm\rm bol}{-}M_{\rm env}$ diagram translates into a power-law relationship between the final mass of the massive star and the mass of the envelope it initially accreted from, and this is in contrast to the predictions from the competitive accretion models.

-

We find timescales for massive star formation of the order of 1$\div$ $4\times 10^5$ years, shortening with increasing mass. The association of these relatively recent objects with young clusters, mostly revealed in the near-IR, with a median stellar age of a few 10⁶ years, confirms that star formation lasts for several dynamical timescales of the hosting star-forming region and that the most massive object is among the last ones to form. This argues in favour of massive star formation as dominated by accretion from bound material, or a scaled-up analogue of the low-mass star formation mechanism.

Acknowledgements

We thank the anonymous referee for very useful comments which considerably improved the paper's clarity.

References

 
• André, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122  
• Behrend, R., & Maeder, A. 2001, A&A, 373, 190  
• Beltrán, M. T., Brand, J., Cesaroni, R., et al. 2006, A&A, 447, 221  
• Benjamin, R. A., Churchwell, E., Babler, B. L., et al. 2003, PASP, 115, 953  
• Beuther, H., Schilke, P., Menten, K. M., et al. 2002a, ApJ 566, 945  
• Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002b, A&A, 383, 892  
• Brand, J., & Blitz, L. 1993, A&A, 275, 67  
• Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488  
• Bonnell, I. A., Bate, M. R., & Vine, S. G. 2003, MNRAS, 343, 413  
• Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 745  
• Brand, J., Cesaroni, R., Palla, F., & Molinari, S. 2001, A&A, 370, 230  
• Bressan, A., Fagotto, F., Bertelli, G., & chiosi, C. 1993, A&AS, 100, 647  
• Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton University Press), 280  
• Carey, S. J., Feldman, P. A., Redman, R. O., et al. 2000, ApJ, 543, L157  
• Cesaroni, R., Felli, M., Testi, L., Walmsley, C. M., & Olmi, L. 1997, A&A, 325, 725  
• Cesaroni, R., Felli, M., Jenness, T., et al. 1999, A&A, 345, 949  
• Churchwell, E. 1997, ApJ, 479, L59  
• Churchwell, E. 2002, ARA&A, 40, 27  
• Egan, M. P., Shipman, R. F., Price, S. D., et al. 1998, ApJ, 494, L199  
• Evans, N. J., Shirley, Y. L., Mueller, K. E., & Knez, C. 2002, in Hot Stars Workshop III: The Earliest Stage of Massive Star Birth, ed. P. A. Crowther, ASP Conf. Proc., 267, 17  
• Faundez, S., Bronfman, L., Garay, G., et al. 2004, A&A, 426, 97  
• Faustini, F., Molinari, S., Testi, L., Brand, J., & Cesaroni, R. 2007, in preparation  
• Fontani, F., Cesaroni, R., Testi, L., et al. 2004, A&A, 414, 299  
• Fontani, F., Beltrán, M. T., Brand, J., et al. 2005, A&A, 432, 921  
• Hill, T., Burton, M. G., Minier, V., et al. 2005, MNRAS, 363, 405  
• Indebetouw, R., Whitney, B. A., Johnson, K. E., & Wood, K. 2006, ApJ, 636, 362  
• Ivezic, Z., & Elitzur, M. 1997, MNRAS, 287, 799  
• Kurtz, S., Cesaroni, R., Churchwell, E., Hofner, P., & Walmsley, C. M. 2000, in Protostars & Planets V, ed. V. Mannings, A. P. Boss, & S. S. Russell, 299  
• Lada, C. J., & Wilking, B. A. 1984, ApJ, 287, 610  
• Lorenzetti, D., Tommasi, E., Giannini, T., et al. 1999, A&A, 346, 604  
• Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364  
• McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850  
• Molinari, S. L., Brand, J., Cesaroni, R., & Palla, F. 1996, A&A, 308, 573  
• Molinari, S., Brand, J., Cesaroni, R., Palla, F., & Palumbo, G. G. C. 1998a, A&A, 336, 339  
• Molinari, S., Testi, L., Brand, J., Cesaroni, R., & Palla, F. 1998b, ApJ, 505, L39  
• Molinari, S., Brand, J., Cesaroni, R., & Palla, F. 2000, A&A, 355, 617  
• Molinari, S., Testi, L., Rodriguez, L. F., & Zhang, Q. 2002, ApJ, 570, 758  
• Norberg, P., & Maeder, A. 2000, A&A, 359, 1025  
• Palla, F., & Stahler, S. W. 1992, ApJ, 392, 667  
• Palla, F., & Stahler, S. W. 1993, ApJ, 418, 414  
• Palla, F., Brand, J., Cesaroni, R., Comoretto, G., & Felli, M. 1991, A&A, 246, 249  
• Plume, R., Jaffe, D. T., Evans, N. J., II, Martin-Pintado, J., & Gomez-Gonzalez, J. 1997, ApJ, 476, 730  
• Preibisch, Th., Ossenkopf, V., Yorke, H. W., & Henning, Th. 1993, A&A, 279, 577  
• Price, S. D., Egan, M. P., Carey, S. J., Mizuno, D. R., & Kuchar, T. A. 2001, AJ, 121, 2819  
• Ramesh, B., & Sridharan, T. K. 1997, MNRAS, 284, 1001  
• Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256  
• Saraceno, P., André, Ph., Ceccarelli, C., Griffin, M., & Molinari, S. 1996, A&A, 309, 827  
• Schnee, S. L., Ridge, N. A., Goodman, A. A., & Li, J. G. 2005, ApJ, 634, 442  
• Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23  
• Simon, R., Jackson, J. M., Rathborne, J. M., & Chambers, E. T. 2006, ApJ, 639, 227  
• Tej, A., Ojha, D. K., Ghosh, S. K., et al. 2006, A&A, 452, 203  
• Testi, L., Palla, F., & Natta, A. 1999, A&A, 342, 515  
• Thompson, R. I. 1984, ApJ, 283, 165  
• Thompson, M. A., & Macdonald, G. H. 2003, A&A, 407, 237  
• Walmsley, C. M. 1995, Rev. Mex. Astron. Astrof. Ser. Conf., 1, 137  
• Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003, ApJ, 591, 1049  
• Whitney, B. A., Indebetouw, R., Bjorkman, J. E., & Wood, K. 2004, ApJ, 617, 1190  
• Williams, J. P., & McKee, C. F. 1997, ApJ, 476, 166  
• Wood, D. O. S., & Churchwell, E. 1989, ApJ, 340, 265  
• Wright, A. E., Griffith, M. R., Burke, B. F., & Ekers, R. D. 1994, ApJS, 91, 111  
• Zhang, Q., Hunter, T. R., Brand, J., et al. 2001, ApJ, 552, L167  
• Zhang, Q., Hunter, T. R., Brand, J., et al. 2005, ApJ, 625, 864

 

  
8 Online Material

   
Appendix A: Notes and data analysis log for individual objects

This appendix summarizes the analysis steps that we carried out for each of the target fields. The detailed procedure is described in the main body of the paper; the sections below illustrates the detailed assumptions that may have been necessary for the specific fields. For each field we report in the first figure the MSX images for all 4 bands, with band A on the top-left and band E on the bottom-right, with the letter identification of the various components considered. The second figure shows in color scale the MSX band A (unless otherwise noted) image as in previous figure, with the millimeter emission overlayed (either 1.2 mm or 850 $\mu$m) in contours and the numbers indicating the different millimeter components considered. Also overlaied on this figure are, with red crosses, any radio continuum source detected at 6 or 3.6 cm in Molinari et al. (1998a, 2002), or in the Parkes 6 cm survey (Wright et al. 1994).

The second figure also shows the coordinate grid; since the color image is the same as the band A image in the first figure (unless otherwise noted), it is then easy to determine source positions in the other MSX bands.

The rest of the figures illustrate the results of the best-fitting procedure. In all plots the full thick and dashed-triple-dotted lines are used for accepted and rejected Whitney's models, respectively; the dashed lines are used to indicate the greybody models, and the dotted lines are used for the sum of Whitney and greybody models in the case of composite fits.

A.1 Mol008/IRAS 05137+3919

Figure A.1: The field of Mol008/IRAS 05137+3919 in the MSX bands A ( top-left), C ( top-right), D ( bottom-left) and E ( bottom-right); identified peaks are indicated with letters.

Figure A.2: Mol008/IRAS 05137+3919 in the MSX A band (color image) with superimposed in green a contour map of the 850 $\mu$m SCUBA emission; naming of the submm peaks as used in the paper is indicated in the figure. Red crosses indicate the presence of radio continuum sources.

There are two distinct peaks in the mid-IR (A.1) and in the submillimeter; the peaks are separated by 35$\div$40 $^{\prime \prime}$, so they both contribute to the IRAS source, and the counterparts at the different wavelengths are positionally coincident (A.2). Both sources are associated with a $K_{\rm s}$ cluster. The submillimeter images are taken with SCUBA in jiggle mode and are too small with respect to the size of the region to be sure we reach the background; it is then very difficult to estimate any contribution from extended emission. The overall emission is very well fitted with two compact peaks on a flat underlying plateau, but it is not clear if the plateau value is an extended component in the cloud or is the background level. It is then not meaningful to try and extrapolate this plateau value to the far-IR at 60 and 100 $\mu$m; however, the SES1₁₀₀=0 and CIRR2 = 2 IRAS flags are not suggestive of a relevant contamination from extended emission in this field, so we simply split F₆₀ and F₁₀₀ between the two peaks proportional to the submm fluxes. Very good ZAMS fits are found for both objects (Figs. A.3 and A.4).

Figure A.3: Observed SED for Mol008/IRAS 05137+3919 A/1. Asterisks are for MSX and submillimeter data; diamonds are for IRAS data (arrows indicate upper limits). The thick full line is the fitted Whitney's model.

Figure A.4: Same as A.3 for Mol008/IRAS 05137+3919 B/2.

   
A.2 Mol011/IRAS 05137+3919

Figure A.5: Same as Fig. A.1 for Mol011/IRAS 05345+3157.

Figure A.6: Same as Fig. A.2 for Mol011/IRAS 05345+3157.

There is an extended blob of emission in the mid-IR, called A, which is evident at all MSX bands (Fig. A.5); no attempt is made to identify and fit multiple components. $K_{\rm s}$ images show a conspicuous cluster coinciding with A. A second component, B, appears only in band E, beyond 20 $\mu$m. The submillimeter emission (Fig. A.6) shows a diffuse halo of emission encompassing the entire region that contains both A and B. The main 850 $\mu$m peak, MM1, coincides with B; there is a second millimeter core immediately S and less intense. A third MM core lies about 10 $^{\prime \prime}$ E of the peak of A and is tentatively associated with A in the SED fit. Likewise Mol008, the submm map is not sufficiently extended to estimate a potential contribution of the extended emission; the usual IRAS extended emission confusion flags SES1 and CIRR2 suggest that this component, if any, may not be significant; we will then partition the 60 and 100 $\mu$m fluxes according to the submm fluxes, because their reciprocal distance is such that all three potentially contribute to the IRAS Far-IR fluxes.

A good global ZAMS fit is found for component A/3 (Fig. A.7). No satisfactory ZAMS model can be found for B/1 (Fig. A.8) and MM2 (Fig. A.9), which are therefore successfully fit as simple greybodies.

Figure A.7: Same as Fig. A.3 for Mol011/IRAS 05345+3157 A/3.

Figure A.8: Symbols as for Fig. A.3 for Mol011/IRAS 05345+3157 B/1. The dashed-triple-dotted line is the best fit Whitney's model which is, however, rejected; the dashed line represents the greybody fit.

Figure A.9: Symbols as for Fig. A.3 for Mol011/IRAS 05345+3157 MM2. The dashed line is the greybody model fit.

A.3 Mol012/IRAS 05373+2349

Figure A.10: Same as Fig. A.1 for Mol012/IRAS 05373+2349.

Figure A.11: Same as Fig. A.2 for Mol012/IRAS 05373+2349.

There is a single mid-IR peak visible in all MSX bands (Fig. A.10) coincident with a $K_{\rm s}$ cluster and a radio continuum source. The submillimeter emission (Fig. A.11) shows a strong peak sitting over a plateau and positionally coincident with the mid-IR source. The source was also detected in the 3.4 mm continnum with interferometry (Molinari et al. 2002). The integral of the extended plateau has been extrapolated to 60 and 100 $\mu$m as specified in Sect. 3.1.3. The fit with the embedded ZAMS model is satisfactory (Fig. A.12). The 12 $\mu$m IRAS flux is assumed as an upper limit because it likely contains contributions from other stars in the cluster; the fit has a deficiency at 850 $\mu$m, although the agreement with the 3.4 mm point is very good.

Figure A.12: Same as Fig. A.3 for Mol012/IRAS 05373+2349.

   
A.4 IRAS 08211-4158

Figure A.13: Same as Fig. A.1 for IRAS 08211-4158.

Figure A.14: Same as Fig. A.2 for IRAS 08211-4158.

The MSX band A image (Fig. A.13) shows a very complex region with considerable extended structure at all scales. However, only two sources are clearly revealed in all bands, A and B. Source A coincides with a $K_{\rm s}$ cluster. The distance between the two sources is $\sim$2.5$^{\prime}$ so they likely both contribute to the IRAS 60 and 100 $\mu$m fluxes; indeed source B coincides with another IRAS source (08209-4200) which the PSC2 gives as confused with the main one beyond 25 $\mu$m and with fluxes 2.1 and 10.1 Jy at 12 and 25 $\mu$m respectively.

The analysis of the 1.2 mm image (Fig. A.14) shows a core+halo structure which we fit simultaneously. As usual, the estimated integral plateau is used to estimate an extended emission contribution at 60 and 100 $\mu$m. The mm peak is $\sim$35 $^{\prime \prime}$ from source A; given the size of the reconstruction beam in the MSX images and the SEST beam at 1.2 mm we cannot exclude that the two counterparts are associated. Fitting the entire SED of source A +MM provides a fit (Fig. A.15) with $\xi=0.14$ and, as such, it should be accepted (the dashed-triple-dotted line). However, the model has a remarkable (more than 60%) deficiency at 100 $\mu$m (which has a relatively low weight in the fit and thus $\xi$ is low); since no other important compact components are clearly visible in the 1.2 mm map, the only way to reconcile model and data is to posit that the 100 $\mu$m flux is much more contaminated by extended emission than our estimate would seem to suggest. This would be possible if the dust temperature of this diffuse component was higher than assumed; T=25 K instead of 20 K for instance, would be sufficient to increase the extended component contribution so that the 100 $\mu$m flux assigned to the source would be in agreement with the model. However, this slight temperature increase would have a dramatic impact on the 60 $\mu$m flux to be assigned to the extended component; as a result the 60 $\mu$m flux which can be assigned to the source would be 25% of what was originally assumed, missing the model by about a factor of 4. For this reason we reject the model and go ahead with a two-component fit using a ZAMS model for the mid-IR portion and a greybody for the Far-IR/submm portion of the SED; these are reported as a full line and a dashed line in Fig. A.15 (plus the sum of the two as a dotted line). The corresponding model parameters are those reported in Table 6.

Source B does not have a millimeter counterpart; considering also the smaller MSX fluxes with respect to source A, it is reasonable to assume that its contribution to the far-IR is negligible. We then fit source B with a ZAMS model assuming the 60 and 100 $\mu$m fluxes as severe upper limits (Fig. A.16).

Figure A.15: Symbols as for Fig. A.3 for IRAS 08211-4158 A - 1 sources. The dashed-triple-dotted line is the best fitting Whitney's model (rejected) to the whole SED; the thick line is the accepted Whitney's model for the mid-IR portion of the SED while the dashed line is the greybody fit to the rest of the SED. Their sum is the dotted line.

Figure A.16: Same as Fig. A.3 for IRAS 08211-4158 B.

A.5 IRAS 08470-4243

Figure A.17: Same as Fig. A.1 for IRAS 08470-4243.

Figure A.18: Same as Fig. A.2 for IRAS 08470-4243.

The MSX band A image (Fig. A.17) shows 4 sources and a diffuse halo of emission encompassing them all. A is clearly the brightest, and only B survives through band E. An additional source, E, is only visible in bands D and E and coincides with a radio continuum source (compare Figs. A.17 and A.18). The 1.2 mm emission (Fig. A.18) shows a core+halo structure which we treat as usual; it is slightly offset toward the SW with respect to source A and hence we tentatively associate them in the fit. A ZAMS model is found for the complete SED and is plotted in Fig. A.19. The slight deficiency above 60 $\mu$m is well within 50% and it is considered acceptable.

We do not attempt any fit for source B because our grid of ZAMS models do not cover situations where the mid-IR SED is flat or decreasing with $\lambda$, as is the case for this source; this source is likely a Herbig Ae/Be star.

Source E has an increasing SED with wavelength, but the absence of detectable mm emission suggests that although in principle it may contribute to the far-IR IRAS fluxes (it is $\sim$35 $^{\prime \prime}$ W-NW of source A), this contribution is probably neglegible. With only two photometric points (MSX bands D and E) we find a large set of models with a B3 and B5 central star that plausibly can be considered, and as there is insufficient information to select one or the other there will be no entry for this source in Table 6.

Figure A.19: Same as Fig. A.3 for IRAS 08470-4243 A/1 source.

A.6 IRAS 08477-4359

Figure A.20: Same as Fig. A.1 for IRAS 08477-4359.

Figure A.21: Same as Fig. A.2 for IRAS 08477-4359.

This field presents one dominant peak (Fig. A.20) (A, coincident with a $K_{\rm s}$ cluster) plus two fainter sources which are only visible in the MSX A band. The 1.2 mm image (Fig. A.21) shows a complex morphology with one main main peak about 1.5$^{\prime}$ SE of source A, and two fainter peaks about 20 $^{\prime \prime}$ away. The latter peaks are immersed in a general halo and they do not stand out sufficiently brightly to plausibly constrain the JMFIT fitting. MM1 is instead fitted with a Gaussian+plateau; the total contribution of the extended emission is obtained computing the integrated emission of the entire structure and subtracting the integrated flux of MM1. The extrapolation of this contribution to 60 and 100 $\mu$m is done as usual.

Given the reciprocal distance we can exclude the association of A with MM1, while the question of the association might be raised for MM2 or MM3. We fit the SED of source A (Fig. A.22) using the IRAS 60 and 100 $\mu$m fluxes as upper limits and a millimeter flux of 0.3 Jy, which is the value corresponding to MM2 and MM3. We cannot find a ZAMS model that can reconcile the almost flat mid-IR portion of the spectrum with the assumed 1.2 mm flux. The model we present in Fig. A.22 is well below the assumed millimeter data point and argue against association of either MM2 or 3 with source A; the model deficiency in MSX band A, also in this case, could be due to an excess flux from PAH contamination.

Source MM1 is then fitted with a greybody (Fig. A.23), where IRAS 60 and 100 $\mu$m fluxes corrected for the extended emission have been used; the SED fitted to source A has been used as an additional contraint in the fit since the reciprocal distance is sufficient to expect potential contamination in the far-IR IRAS fluxes. The fit is shown in Fig. A.23, where in mid-IR we have used as additional contraints the upper limits estimated from the MSX maps corresponding to MM1. The fit misses the 60 $\mu$m point, and we cannot explain this; if we force the parameters so that the 60 $\mu$m is also fit, we obtain mid-IR fluxes largely in excess of the upper limits.

Figure A.22: Same as Fig. A.3 for IRAS 08477-4359 A source.

Figure A.23: Same as Fig. A.9 for IRAS 08477-4359 MM1 source.

A.7 IRAS 08563-4225

Figure A.24: Same as Fig. A.1 for IRAS 08563-4225.

Figure A.25: Same as Fig. A.2 for IRAS 08563-4225.

There is one bright peak in the mid-IR (Fig. A.24) and it coincides with a $K_{\rm s}$ cluster. The 1.2 mm image (Fig. A.25) shows a complex structure extending NE-SW, with a dominant peak only slightly SW of the mid-IR peak. Cuts through the mm peak axis shows clearly a main compact peak overlaid on a shallower component, and as such they are fitted with JMFIT. It is however impossible to plausibly estimate the contribution of the extended emission in IRAS far-IR because of the extension of the mm emission. We tried to integrate the emission in a 5$^{\prime}$ beam centered on the main mm peak, but the usual extrapolation to 100 $\mu$m yields a value equal to the entire IRAS flux and this is of course not acceptable. On the other hand extended emission at 100 $\mu$m should be present as suggested by the SES1 = 1 and CIRR2 = 4 flags, so in the fit we will assume the 60 and 100 $\mu$m fluxes as upper limits. A ZAMS model is indeed found (Fig. A.26).

Figure A.26: Same as Fig. A.3 for IRAS 08563-4225 A/1 sources.

A.8 IRAS 08589-4714

Figure A.27: Same as Fig. A.1 for IRAS 08589-4714.

Figure A.28: Same as Fig. A.2 for IRAS 08589-4714.

The MSX band A image (Fig. A.27) shows one peak (A), seen also in the other bands, sitting in the convexity of an arc-shaped structure; there are other peaks but only one (named B) shows up marginally at 21.4 $\mu$m. In the 1.2 mm image (Fig. A.28) there is one peak (MM1) which is slightly shifted SW with respect to A, and which sits in a small oblate halo. Cuts through the axis of the structure allow us to estimate the size of the peak, and we use these values as contraints in the JMFIT of the peak+plateau fit. The difference between the integrated flux of the structure and the integrated flux of the fitted peak provides the contribution of the extended halo and, with the usual procedure, we extrapolate this to 60 and 100 $\mu$m and correct the IRAS fluxes to be used in the source SED fit. Another fainter peak (MM2) is located south of source B; clearly the dominant contribution to the far-IR fluxes comes from the A+MM1 combination, so that any attempt to fit B and/or MM2 would be poorly constrained.

The ZAMS model fit to the A+MM1 SED (we keep the IRAS 12 $\mu$m flux as an upper limit due to the considerable extended emission visible in the MSX C band map) provides a satisfactory fit (Fig. A.29) apart from a deficiency at 8 $\mu$m and 100 $\mu$m which are both acceptable given the uncertainties in the 100 $\mu$m flux assignment and the potential level of PAH contamination of the MSX A band flux.

Figure A.29: Same as Fig. A.3 for IRAS 08589-4714 A - 1 sources.

A.9 IRAS 09131-4723

Figure A.30: Same as Fig. A.1 for IRAS 09131-4723.

Figure A.31: Same as Fig. A.2 for IRAS 09131-4723.

Four sources are visible in the MSX band A images (Fig. A.30), but only two, A and B, are visible also at 21.4 $\mu$m. Each of these objects is associated with a cluster of $K_{\rm s}$ sources. Source C is also associated with a bright near-IR star and with another fainter IRAS source which is separated form our main IRAS source only at 12 $\mu$m. The 1.2 mm image shows (Fig. A.31) a three-peaked structure immersed in a larger halo. MM1 and 2 are locate respectively E ($\sim$15 $^{\prime \prime}$) and W ($\sim$35 $^{\prime \prime}$) of source A, while M3 is halfway between A and B. We perform brightness cuts in various directions across the peaks to estimate the size of the cores, and use this information to improve the JMFIT we performed to derive the Gaussian fits to the peaks. There is plateau of emission and we estimate its contribution subtracting the integrated flux of the three cores from the integrated flux of the entire emission patch. This value is estrapolated to 60 and 100 $\mu$m as usual, and the results are subtracted from the IRAS fluxes. the remaining fluxes are split among the three millimeter peaks proportionally to the 1.2 mm flux.

We assume that A and MM1 are associated and attempt a ZAMS fit to the assembled SED and a satisfactory fit is indeed found (Fig. A.32) with departures from the data that are within acceptable limits (see discussions for other sources above and the assumptions discussed in the text (Sects. 3.1.3 and 3.1.2). MM2 and MM3 are fitted with greybodies and good fits are found (Figs. A.33 and A.34); their emission shortward of 20 $\mu$m is well below the dominant contribution from source A/1.

We fit source B (Fig. A.35) taking the IRAS fluxes as upper limits; the source is at the boundary of the extended 1.2 mm emission and we will take the local value of the 1.2 mm emission as an additional upper limit constraint. A good fit is found.

Figure A.32: Same as Fig. A.3 for IRAS 09131-4723 A - 1 sources.

Figure A.33: Same as Fig. A.9 for IRAS 09131-4723 MM 2 source.

Figure A.34: Same as Fig. A.9 for IRAS 09131-4723 MM3 source.

Figure A.35: Same as Fig. A.3 for IRAS 09131-4723 B.

A.10 IRAS 10019-5712

Figure A.36: Same as Fig. A.1 for IRAS 10019-5712.

Figure A.37: Same as Fig. A.2 for IRAS 10019-5712.

There are two peaks in the mid-IR (Fig. A.36), and the 1.2 mm image (Fig. A.37) shows an elongated structure embedding them, which plausibly results from the blending of two cores; source B is coincident with a very bright $K_{\rm s}$ source. We estimate the size of the cores using cross cuts in the brightness profiles; the fits to these cuts suggests that the cores are lying on a diffuse halo. The estimated sizes are used as constraints in JMFIT. An additional contraint is assumed for the position of the two cores; the first core is taken at the position of the mm peak, while the second is constrained at an offset equal to the offset between A and B. In this way we can obtain a comprehensive fit of the two core+plateau; the extended emission contribution is estimated by subtracting the integrated flux of the two cores from the integrated flux of the entire emission patch, and extrapolated as usual to the IRAS far-IR wavelengths. IRAS fluxes are then corrected and split between the two cores in proportion to the 1.2 mm fluxes.

We associated A with MM1 and B with MM2, and obtain very good fits with ZAMS models for both objects (Figs. A.38 and A.39). In both cases the IRAS 12 and 25 $\mu$m fluxes have been taken as upper limits because they are blended in the IRAS beam.

Figure A.38: Same as Fig. A.3 for IRAS 10019-5712 A/1 source.

Figure A.39: Same as Fig. A.3 for IRAS 10019-5712 B/2 source.

   
A.11 IRAS 10095-5843

Figure A.40: Same as Fig. A.1 for IRAS 10095-5843.

Figure A.41: Same as Fig. A.2 for IRAS 10095-5843.

There are two objects, A and B (Fig. A.40); A is coincident with a $K_{\rm s}$ cluster and brightens with wavelength, while B fades and almost disappears at 21.4 $\mu$m. The 1.2 mm image (Fig. A.41) shows a peak+halo structure, with the peak situated about 15 $^{\prime \prime}$ away from A. We carry out the usual procedure with JMFIT to separate the core form the halo; the integrated flux of the core is subtracted from the entire emission patch and extrapolated as usual to far-IR IRAS wavelengths; these extrapolations are then subtracted from the IRAS fluxes and the residuals are used in the model fit of the source. A ZAMS fit is found with $\xi=0.14$ with model departure at 100 $\mu$m at the limit of 50%. The examination of the 10 best-$\xi$ models provides an alternative model with $\xi=0.17$ which better fits the 100 $\mu$m while still fitting the mid-IR portion within acceptable levels. This model is presented in Fig. A.42 and its parameters are reported in Table 6.

Figure A.42: Same as Fig. A.3 for IRAS 10095-5843 A/1 source.

A.12 IRAS 10184-5748

Figure A.43: Same as Fig. A.1 for IRAS 10184-5748.

Figure A.44: Same as Fig. A.2 for IRAS 10184-5748.

There are several peaks visible in the MSX A band (Fig. A.43), with considerable extended emission embedding all of them. Only 4 sources survive through the 21.4 $\mu$m band and these are the ones that we will try to fit. The 1.2 mm image (Fig. A.44) is equally complex and we see a complex emission patch extending N-S and embedding sources A, B, C and G. Sources F and D coincide with very bright $K_{\rm s}$ objects, while no evidence of a cluster is seen toward source A.

The strongest millimeter peak (MM1) is coincident with A, and there is a prominence toward the N corresponding to B; source C is in a local minimum of millimeter emission, while G corresponds to a fainter and relatively isolated local maximum of emission (MM2). Given the complexity of the emission we do not attemp a comprehensive fit of the peaks; cross-cuts of the emission profiles are done to verify the core-halo structure of the emission and to measure the extent of the cores. MM1 is fitted simultaneously with an underlying variable plateau. We cannot find a plausible JMFIT solution for the mm emission on top of B, so we will take the local value of the emission as an indication. The same is done for the mm emission on top of source C; in that case, however, the emission will be an upper limit because there is no apparent core corresponding to C. The flux of MM2 is not on a plateau, and we integrate the entire emission patch. The procedure to estimate the contribution of extended emission at far-IR fluxes is quite dubious in this case; we integrate the entire emission in a 5$^{\prime}$ beam centered on A, and subtract the integrated flux of MM1. The result is extrapolated to the IRAS bands and subtracted from the IRAS fluxes to perform the fit to source A+MM1 (Fig. A.45). For sources B, C and G (Figs. A.46-A.48), we will assume the IRAS fluxes as strict upper limits.

We find good fits with ZAMS models for all fitted sources.

Figure A.45: Same as Fig. A.3 for IRAS 10184-5748 A source.

Figure A.46: Same as Fig. A.3 for IRAS 10184-5748 B source.

Figure A.47: Same as Fig. A.3 for IRAS 10184-5748 C source.

Figure A.48: Same as Fig. A.3 for IRAS 10184-5748 G source.

   
A.13 IRAS 10501-5556

Figure A.49: Same as Fig. A.1 for IRAS 10501-5556.

Figure A.50: Same as Fig. A.2 for IRAS 10501-5556.

Quite simple field; there is only one peak visible in all 4 MSX bands (Fig. A.49) and it coincides with a bright $K_{\rm s}$ source with a cluster. The mm image (Fig. A.50) shows a peak on a halo; the analysis is performed as usual and a good ZAMS model is found (Fig. A.51); the 8 $\mu$m data excess is within acceptable limits for PAHs contamination.

Figure A.51: Same as Fig. A.3 for IRAS 10501-5556 source.

A.14 IRAS 10555-6242

Figure A.52: Same as Fig. A.1 for IRAS 10555-6242.

Figure A.53: Same as Fig. A.2 for IRAS 10555-6242.

Two sources (A and B) (Fig. A.52), but B is only visible at 8 $\mu$m; there is perhaps a very small $K_{\rm s}$ cluster toward source A, while source B seems to be associated with a bright $K_{\rm s}$ source. Peak+halo structure in the 1.2 mm image (Fig. A.53); the peak is slightly offset with respect to A. The usual analysis is performed and a good fit of A+MM1 with a ZAMS model is found; the model deficiency at 1.2 mm is well within 50%, while the 8 $\mu$m departure (Fig. A.54) is within acceptable limits for PAHs contamination.

Figure A.54: Same as Fig. A.3 for IRAS 10555-6242 source.

A.15 IRAS 11380-6311

Figure A.55: Same as Fig. A.1 for IRAS 11380-6311.

Figure A.56: Same as Fig. A.2 for IRAS 11380-6311.

There is one mid-IR peak in a complex pattern of extended emission and fainter peaks (Fig. A.55). Only the main central peak is visible throughout the mid-IR, and is coincident with a $K_{\rm s}$ peak. The 1.2 mm image (Fig. A.56) has a strong peak about 30 $^{\prime \prime}$ E of the mid-IR peak, with an elongation toward the NE. The analysis is carried out in the usual way. We find a very good fit (Fig. A.57); the 8 $\mu$m departure is at the limit of acceptability but we accept the model since the rest of the SED is fitted very well.

Figure A.57: Same as Fig. A.3 for IRAS 11380-6311 source.

A.16 IRAS 12063-6259

Figure A.58: Same as Fig. A.1 for IRAS 12063-6259.

Figure A.59: Same as Fig. A.2 for IRAS 12063-6259.

Strong mid-IR peak (coincident with a radio continuum source) with some other components that disappear beyond 8 $\mu$m (Fig. A.58). The source is associated with a $K_{\rm s}$ star+nebulosity; there might also be a small cluster. The mm image (Fig. A.59) shows the usual peak-halo structure; the usual analysis is performed for this quite simple field. Excellent fit to a ZAMS model (Fig. A.60).

Figure A.60: Same as Fig. A.3 for IRAS 12063-6259 source.

A.17 IRAS 12127-6244

Figure A.61: Same as Fig. A.1 for IRAS 12127-6244.

Figure A.62: Same as Fig. A.2 for IRAS 12127-6244.

There is one strong peak in all mid-IR bands (Fig. A.61), coinciding with a radio continuum source, in a complex emission pattern; it is associated with a $K_{\rm s}$ source with a small nebulosity. There is an elongated emission pattern east of the main peak. The 1.2 mm image (Fig. A.62) has a strong peak (MM1) very close to source A, in an extended structure that is resolved in fainter peaks (MM2 is the brightest) toward the east, corresponding to the dark region in the 21.4 $\mu$m image between source A and the elongated mid-IR structure. We perform cross-cuts across the two peaks and we use the estimate of the core size in the JMFIT analysis. The rest of the emission is attributed to extended emission and we carry out the usual analysis to extrapolate to far-IR IRAS bands and derived corrected fluxes for the SED fit to source A+MM1. The fit of MM1+A to a ZAMS model gives a very good match (Fig. A.63); MM2 is fitted with a simple greybody (Fig. A.64), and its contribution shortward of 25 $\mu$m is negligible with respect to the main source.

Figure A.63: Same as Fig. A.3 for IRAS 12127-6244 A/1 source.

Figure A.64: Same as Fig. A.9 for IRAS 12127-6244 MM2 source.

   
A.18 IRAS 14000-6104

Figure A.65: Same as Fig. A.1 for IRAS 14000-6104.

Figure A.66: Same as Fig. A.2 for IRAS 14000-6104.

The mid-IR shows a very complex region (Fig. A.65) with many peaks in extended emission. There a number of peaks that are visible throught the four MSX bands. Several peaks are also found in the 1.2 mm image (Fig. A.66) and some of them can be positionally associated with mid-IR sources. Given the complexity of the region we prefer not to attempt any estimate of the extended emission; considering that MM cores 1 to 4 should all contribute the most to the IRAS 60 and 100 $\mu$m fluxes, we split those fluxes among the 4 cores in proportion to the 1.2 mm fluxes.

Source A+MM1 cannot be fit by a ZAMS model (dashed-triple-dotted line in Fig. A.67). We perform a two-component fit as usual in these cases: the ZAMS model is a factor of 2 below the data at 8 $\mu$m, although the background subtraction for the mid-IR flux (especially at 8 $\mu$m, see Fig. A.65) is critical in this complex region.

We try to fit source B together with MM5. There is another IRAS source (14004-6104) midway between sources B and C; the fluxes at 60 and 100um are confused with the main IRAS source, so that we will use the IRAS fluxes as strict upper limits. In spite of the low value of $\xi$ the fit (Fig. A.68) significantly departs from the data at 8 $\mu$m. We choose to accept the model anyway because it would not make sense to try a two-component fit without a reliable estimate of the FIR fluxes; besides, the error has no impact on the conclusions since most of them are dependent on the dominant sources in each field (which is not the case for the present source).

To fit source C we take the local value of the millimeter emission as an upper limit. The IRAS fluxes are those of the IRAS source between B and C. We find a good ZAMS model (Fig. A.69).

Source D coincides with yet another IRAS source (13595-6100) whose 60 and 100 $\mu$m fluxes are given at FQUAL quality 1, and so we treat them as upper limits in the fit. D is positionally associated with MM7, and we indeed find few ZAMS models that fit the data (Fig. A.70).

Source H may be associated with MM3 and as such we find a ZAMS model fit which is not entirely satisfactory beyond at 1.2 mm and 8 $\mu$m, but we accept it for the same arguments as for source B above (Fig. A.71).

Source I may be associated with MM4. The fit is very difficult in the sense that the estimated fluxes in the far-IR are affected by very large uncertainty. We find a ZAMS model (Fig. A.72) which is not satisfactory, similarly to source H above.

Since F and G do not have a mm counterpart we neglect them.

Figure A.67: Same as Fig. A.15 for IRAS 14000-6104 A - 1 source.

Figure A.68: Same as Fig. A.3 for IRAS 14000-6104 B - 5 source.

Figure A.69: Same as Fig. A.3 for IRAS 14000-6104 C source.

Figure A.70: Same as Fig. A.3 for IRAS 14000-6104 D/7 source.

Figure A.71: Same as Fig. A.3 for IRAS 14000-6104 H/3 source.

Figure A.72: Same as Fig. A.3 for IRAS 14000-6104 I/4 source.

A.19 IRAS 14057-6032

Figure A.73: Same as Fig. A.1 for IRAS 14057-6032.

Figure A.74: Same as Fig. A.2 for IRAS 14057-6032.

There are several sources in the 8 $\mu$m map, but only one, A, is visible through 21.4 $\mu$m (Fig. A.73); it is also associated with a bright $K_{\rm s}$ star. The 1.2 mm image (Fig. A.74) shows a single peak with no clear presence of diffuse halo judging from the brightness profile analysis. We then make no assignment for the extended emission contribution. We find an excellent fit with a ZAMS model (Fig. A.75).

Figure A.75: Same as Fig. A.3 for IRAS 14057-6032 A/1 source.

A.20 IRAS 15068-5733

Figure A.76: Same as Fig. A.1 for IRAS 15068-5733.

Figure A.77: Same as Fig. A.2 for IRAS 15068-5733.

There are 4 peaks visible at 8 $\mu$m, plus considerable extended emission, and a fifth peak (E) that appears mostly in MSX bands D and E (Fig. A.76). Source A is associated with a $K_{\rm s}$ star. The 1.2 mm image (Fig. A.77) shows only one peak centered on source A; the brightness profile analysis shows negligible halo structure so we assign the entire flux to the source. A good fit to the ZAMS model is found for the far distance (Fig. A.78), with overall model departures from the data well within acceptable levels.

Figure A.78: Same as Fig. A.3 for IRAS 15068-5733 A/1 source.

A.21 IRAS 15454-5335

Figure A.79: Same as Fig. A.1 for IRAS 15454-5335.

Figure A.80: Same as Fig. A.2 for IRAS 15454-5335.

There are two peaks in the mid IR images at all wavelengths (Fig. A.79), but only source A seems to be associated with a 1.2 mm couterpart (Fig. A.80). Source A also seems to be associated with a bright $K_{\rm s}$ source. The mm peak sits on a plateau that as usual we use to estimate the contribution of extended emission, extrapolate it to far-IR wavelengths and correct IRAS fluxes for the SED fit. Source A+MM1 is fitted by a single ZAMS model within acceptable levels (Fig. A.81).

To fit source B we keep the IRAS fluxes as upper limits since they are dominated by source A, and a good fit is easily found (Fig. A.82).

Figure A.81: Same as Fig. A.3 for IRAS 15454-5335 A - 1 source.

Figure A.82: Same as Fig. A.3 for IRAS 15454-5335 B source.

A.22 IRAS 15519-5430

Figure A.83: Same as Fig. A.1 for IRAS 15519-5430.

Figure A.84: Same as Fig. A.2 for IRAS 15519-5430.

There are several peaks in the mid-IR (Fig. A.83), but only two, A and B, are visible in all four MSX bands; additionally, one peak, G, is only visible in bands D and E. The IRAS source is associated with A, which is also associated with a $K_{\rm s}$ cluster, while corresponding to B and G there two other IRAS sources. Source B is associated with IRAS 15514-5427, with fluxes 10.6, 9.7, 45.2 (with FQUAL=1 and hence an upper limit) and 100 $\mu$m flux confused with the main IRAS source. Source G is associated with IRAS 15517-5432 with fluxes of 2.5, 11.1 in the first two bands, and confused with the main IRAS source at 60 and 100 $\mu$m. The 1.2 mm image (Fig. A.84) shows emission associated with A only (about 30 $^{\prime \prime}$ SW of A). The brightness profile analysis reveals a plateau underlying the main peak, and we proceed with the usual analysis.

To fit source A we take the 12 and 25 $\mu$m fluxes as upper limits given the high level of extended emission seen in all the MSX images towards A. There is no ZAMS model that can fit the entire SED, so we use the two-component fit (Fig. A.85).

Source B cannot be fit by any of our ZAMS models; it either requires a later spectral class or a fainter envelope, therefore being similar to a classical Ae/Be star.

Source G can instead be fitted by a ZAMS model (Fig. A.86).

Figure A.85: Same as Fig. A.15 for IRAS 15519-5430 A - 1 source.

Figure A.86: Same as Fig. A.3 for IRAS 15519-5430 G source.

A.23 IRAS 15557-5337

Figure A.87: Same as Fig. A.1 for IRAS 15557-5337.

Figure A.88: Same as Fig. A.2 for IRAS 15557-5337.

There are 4 sources plus some extended emission in the mid-IR (Fig. A.87); two sources, A and B, are much brighter and have a strong $K_{\rm s}$ counterpart. Given the strength and the vicinity of the 4 sources we fit their fluxes with JMFIT importing them into AIPS. The 1.2 mm image (Fig. A.88) shows two bright peaks coinciding with A and B and a small underlying plateau. We then carry out the usual analysis to obtain the source SED fit, assuming the totality of the far-IR IRAS flux coming from A and B (split in proportion to the 1.2 mm fluxes). We assume the near distance, and find a good fit with ZAMS models for all 4 objects (Figs. A.89 to A.92); for the fit of C and D the IRAS fluxes are taken as strict upper limits.

Figure A.89: Same as Fig. A.3 for IRAS 15557-5337 A/1 source.

Figure A.90: Same as Fig. A.3 for IRAS 15557-5337 B/2 source.

Figure A.91: Same as Fig. A.3 for IRAS 15557-5337 C source.

Figure A.92: Same as Fig. A.3 for IRAS 15557-5337 D source.

A.24 IRAS 16148-5011

Figure A.93: Same as Fig. A.1 for IRAS 16148-5011.

Figure A.94: Same as Fig. A.2 for IRAS 16148-5011.

There are 5 sources in the 0.2$^{\circ}$ $\times$ 0.2$^{\circ}$ MSX images that are clearly visible throughout 21.4 $\mu$m (Fig. A.93). The 1.2 mm image (Fig. A.94) shows a two-peaked structure corresponding to source A, which also coincides with a $K_{\rm s}$ cluster. The two peaks (the faint emission feature on the E side is neglected) are well fitted with Gaussians using JMFIT, and the brightness profile analysis does not show a significant extended component; the IRAS far-IR fluxes will be split between the two mm peaks in proportion to the 1.2 mm fluxes. There is another very compact peak (MM3) coinciding with source C. Unfortunately the millimeter map does not cover source B.

Source A is midway between MM1 and MM2 at an almost equal distance of 30 $^{\prime \prime}$ from both. There are no problems of image registration or astrometric accuracy because source C and MM3 are exactly coincident, so we will assume that A is not associated with either MM1 or MM2. The two component fit is performed for MM1 and 2 (Figs. A.95, A.96), associating each of them in turn to source A merely to contrain the mid-IR portion of the greybody fit. We point out that since in this case we have a strong positional argument to rule out the association of mid-IR and mm counterparts, the figures will not show a ``closest match'' global fitting attempt (the dashed-triple-dotted line of other figures).

There is another IRAS source midway between sources B and C (16153-5016 with fluxes 102 89 727 1839); from the comparison of the mid-IR MSX fluxes of sources B and C we see that while source C has a SED rising to 21.4 $\mu$m, source B peaks at 12 $\mu$m and then consistently declines at longer wavelengths; we then chose to assign the totality of the 60 and 100 $\mu$m fluxes to source C. Source C is fitted with a ZAMS model (Fig. A.97).

Figure A.95: Same as Fig. A.15 for IRAS 16148-5011 A - 1 source. No fit to the overall SED has been attempted with Whitney's models since the mid-IR and submm components are spatially separated (see above text).

Figure A.96: Same as Fig. A.95 but for source A and source 2.

Figure A.97: Same as Fig. A.3 for IRAS 16148-5011 C/3 source.

A.25 IRAS 16232-4917

Figure A.98: Same as Fig. A.1 for IRAS 16232-4917.

Figure A.99: Same as Fig. A.2 for IRAS 16232-4917.

There is one dominant peak in the mid-IR (Fig. A.98), which coincides with a $K_{\rm s}$ cluster and with the peak of the 1.2 mm emission (Fig. A.99). The latter shows no evidence of an extended component. We successfully fit a ZAMS model to the whole SED (Fig. A.100); the model deficiency at 100 $\mu$m is within 50% and hence tolerable.

Figure A.100: Same as Fig. A.3 for IRAS 16232-4917 A source.

   
A.26 IRAS 16501-4314

Figure A.101: Same as Fig. A.1 for IRAS 16501-4314.

Figure A.102: Same as Fig. A.2 for IRAS 16501-4314.

There are three main peaks at 8 $\mu$m (plus two filamentary structures going from A to SE and SW), but only peak A survives through 21.4 $\mu$m (Fig. A.101). The 1.2 mm image (Fig. A.102) shows two peaks on each side of source A (MM1 SE and MM2 W of A); cross cuts through the peaks are done to estimate the size of the peaks and show that there is an underlying extended component, which is treated as usual to correct the IRAS 60 and 100 $\mu$m fluxes for the extended emission. The resulting values are split among MM1 and MM2 in proportion to the 1.2 mm fluxes. The two peaks plus the underlying extended component are immersed in a larger emission patch that is superimposed on the filamentary structures mostly seen at 8 $\mu$m.

Fitting the global SED of A+MM1 (or equivalently MM2, since the two cores have basically identical parameters) does not provide any good fit either for the near or for the far distance. Adopting the near distance the closest match is shown as usual as the dashed-triple-dotted line in Fig. A.103; the 100 $\mu$m flux is missed by slightly more than 50%, and the MSX bands A and C fluxes are underestimated by about an order of magnitude. The 100 $\mu$m flux assigned to the source could be decreased if we increase the contribution assigned to the extended emission component, but the same arguments discussed for source IRAS 08211-4158 A (Sect. A.4) apply here. As for source IRAS 10095-5843 A (Sect. A.11), we examined the 10 best-$\xi$ models and we could find alternatives with $\xi<0.3$ which provide a better fit in the mid-IR, but the amount by which they miss the mid-IR part is even more than an order of magnitude.

We then use a two-component model, shown with the usual line coding. An identical situation is obtained for peak MM2 since it is identical to MM1.

Figure A.103: Same as Fig. A.15 for IRAS 16501-4314 A - 1 source.

   
A.27 IRAS 16535-4300

Figure A.104: Same as Fig. A.1 for IRAS 16535-4300.

Figure A.105: Same as Fig. A.2 for IRAS 16535-4300.

We see three sources in the MSX 8 $\mu$m image; however A and B fade considerably going toward 21.4 $\mu$m (Fig. A.104), while source C remains strong throughout the mid-IR wavelength range. Curiously, the 1.2 mm image (Fig. A.105) shows a peak with no or little evidence for an extended component. Source B is associated with a $K_{\rm s}$ cluster, while for source A this evidence is dubious (in the sense that it is not immediately visible to the eye, as in all cases above where an association with a $K_{\rm s}$ cluster has been assigned).

It is plausible to anticipate, given the mid-IR SED of source A, that the 1.2 mm peak is not associated with A, and indeed no satisfactory ZAMS model is found (Fig. A.106). A two-component model is then used, although the fading of the mid-IR going toward MSX band E is such that no satisfactory model is found for the mid-IR; this is because probably later spectral classes than the ones spanned by our model grid should be attempted.

Source B cannot be fit for the same reasons.

Source C is close to IRAS 16536-4305: fluxes are 4.2, 14.9, 58.5 and confused with the main IRAS source at 100 $\mu$m. However the 60 $\mu$m flux has FQUAL = 1 so it will be assumed as an upper limit. A good ZAMS model fit is found (Fig. A.107).

Figure A.106: Same as Fig. A.15 for IRAS 16535-4300 A - 1 source.

Figure A.107: Same as Fig. A.3 for IRAS 16535-4300 C source.

A.28 IRAS 17082-4114

Figure A.108: Same as Fig. A.1 for IRAS 17082-4114.

Figure A.109: Same as Fig. A.2 for IRAS 17082-4114.

We identify 5 main peaks in the mid-IR (Fig. A.108), but only A and B are prominent throught the entire MSX range. The $K_{\rm s}$ images show faint nebulosities coincident with these two objects, while very bright $K_{\rm s}$ sources are found coinciding with both sources C and D. The 1.2 mm image (Fig. A.109) shows a complex three-peaked structure, where none of the peaks is coincident with either source A or B. We carefully estimate the size of the three peaks using brightness profile analysis in several directions, and this information can then be used to constrain the JMFIT fitting. The contribution of the extended emission is not prominent.

We try to associate A with any of the MM peaks (Figs. A.110 to A.112) (B is similar), but no ZAMS models are found to fit the entire SED (dashed-triple-dotted line is the closest match). We then try the two-component fit constraining each time MM1, 2 and 3 to the ZAMS model found for A in the mid-IR. We could have also contrained the fit to the mid-IR ZAMS model of source B and the results would have been the same.

Figure A.110: Same as Fig. A.15 for IRAS 17082-4114 A - 1 source.

Figure A.111: Same as Fig. A.110 but trying to associate source A with source 2.

Figure A.112: Same as Fig. A.110 for IRAS 17082-4114 A - 3 source.

A.29 IRAS 17425-3017

Figure A.113: Same as Fig. A.1 for IRAS 17425-3017.

Figure A.114: Same as Fig. A.2 for IRAS 17425-3017.

There are several sources in the MSX 8 $\mu$m image, but only source A appears throughout 21.4 $\mu$m (Fig. A.113). Source B is associated with a bright $K_{\rm s}$ cluster. The 1.2 mm image (Fig. A.114) is analysed with cross-cuts through the peak and a peak+core halo structure is visible. The usual analysis is done. The peak is offset about 20$\div$30 $^{\prime \prime}$ north of source A.

The fit to source A does not provide a ZAMS model for the entire mid-IR/mm SED, and the two-component approach is used (Fig. A.115).

Figure A.115: Same as Fig. A.15 for IRAS 17425-3017 A - 1 source.

A.30 Mol038/IRAS 18024-2119

Figure A.116: Same as Fig. A.1 for Mol038/IRAS 18024-2119.

Figure A.117: Same as Fig. A.2 for Mol038/IRAS 18024-2119.

All sources appearing at 8 $\mu$m fade away at longer wavelengths. There is only one faint 21.4 $\mu$m source (Fig. A.116), E, which seems to coincide with a bright peak at 450 and 850 $\mu$m (Fig. A.117). There seems to be an extended submm component but the SCUBA observations were done only in a small field around the IRAS source position and the estimate of the real background is problematic. We assume the IRAS 100 $\mu$m flux as an upper limit in the fit.

Source G is the only one we try to fit because is the only one visible at 21.4 $\mu$m. We do not find a ZAMS model fitting the entire SED (Fig. A.118); since there is only one data point in the mid-IR we cannot go ahead with the usual two-component approach, so we fit the entire SED with a greybody model and find a very good fit from 21 $\mu$m to 850 $\mu$m; in addition to 100 $\mu$m, upper limits have been estimated for source G in the first three MSX images. The IRAS 12 $\mu$m flux has also been assumed as an upper limit due to the contamination of the nearby objects.

Figure A.118: Same as Fig. A.8 for Mol038/IRAS 18024-2119 G/1.

A.31 Mol045/IRAS 18144-1723

Figure A.119: Same as Fig. A.1 for Mol045/IRAS 18144-1723.

Figure A.120: Same as Fig. A.2 for Mol045/IRAS 18144-1723.

There are several peaks in the mid-IR (Fig. A.119), but only source A survives through 21.4 $\mu$m. There is a radio continuum source coinciding with source E.

The 450 and 850 $\mu$m images (Fig. A.120) show a bright compact peak on a diffuse halo. Cross-cuts through the peak are performed to estimate the size of the peak and better contrain the JMFIT analysis. The extent of the halo is less than the extent of the map, so that a plausible estimate of the extended emission at 60 and 100 $\mu$m is possible.

We find a ZAMS model fitting the entire SED (Fig. A.121) of A and MM1; the fit is not very good at 850 $\mu$m, but the departure is within the 50% limit so we accept the fit.

Figure A.121: Same as Fig. A.3 for Mol045/IRAS 18144-1723 A - 1.

A.32 Mol059/IRAS 18278-1009

Figure A.122: Same as Fig. A.1 for Mol059/IRAS 18278-1009.

Figure A.123: Same as Fig. A.2 for Mol059/IRAS 18278-1009.

There are several mid-IR peaks that are visibile up to 21.4 $\mu$m, but A is the brightest (Fig. A.122). It coincides with a submm peak visible both at 450 and 850 $\mu$m (Fig. A.123). These images show the presence of a peak+halo structure, but the maps show large regions with negative signals, suggesting that there may have been problems in the background subtraction in the SCUBA jiggle map. The integrity of the peak flux should be preserved because it is estimated on top of the background, but the absolute value of the background is highly uncertain. We therefore will not attempt any estimate of the extended emission and assume the 100 $\mu$m flux as an upper limit.

The fit to source A+MM1 provide a satisfactory fit with a ZAMS model; the departures in the submm are within tolerances (Fig. A.124).

We do not attempt a fit with sources B, C and D because their mid-IR peaks around 12 $\mu$m and then fades with wavelength, and we have already seen that our model grids fails in finding models for this type of SED (see Sect. A.27).

Figure A.124: Same as Fig. A.3 for Mol059/IRAS 18278-1009 A - 1.

   
A.33 Mol075/IRAS 18511+0146

Figure A.125: Same as Fig. A.1 for Mol075/IRAS 18511+0146.

Figure A.126: Mol075/IRAS 18511+0146. ISOCAM 15 $\mu$m emission image with superimposed contour plot of the 450 $\mu$m emission.

Figure A.127: Mol075/IRAS 18511+0146. ISOCAM 15 $\mu$m emission image with superimposed contour plot of the 850 $\mu$m emission.

The main source A is bright throughout the MSX wavelength range (Fig. A.125). ISOCAM images at 6.75 and 15 $\mu$m show that source A is split into a bright component (A1) and a much fainter component (A2) immediately S-SE. Both the 450 (Fig. A.126) and 850 $\mu$m (Fig. A.127) images show a very bright peak almost coinciding with A1, and an elogation S-SE close to A2. The spatial resolution of the submm images is not sufficient to uniquely deblend two cores, so we will force JMFIT to find two cores constraining the relative distance between the two cores equal to the relative distance between A1 and A2. An additional constraint is established assigning the FWHM of the peaks as estimated from cross-cuts through the brightness profile analysis. Clearly we are then associating A1 with MM1 and A2 with MM2.

The usual procedure of estimating the contribution of extended emission to the IRAS 60 and 100 $\mu$m bands, and splitting the residuals between the two sources in proportion to the submm fluxes, cannot be reliably done in this case; the flux ratio between the two sources is indeed very different at 450 and 850 $\mu$m. We then choose to keep the nominal IRAS fluxes in the fit for both objects as upper limits.

When we try to fit A1+MM1 with a ZAMS model (Fig. A.128) we find a model with $\xi=0.17$ (dashed-triple-dotted line in Fig. A.128 which fits the SED very well except for a departure at 450 $\mu$m of more than a factor of 3, which we deem not acceptable. A two-component model is then used.

The SED of A2+MM2 is instead very well fitted by a ZAMS model (Fig. A.129).

There is a radio continuum source associated with the faint source located about 2$^{\prime}$ S of A in the 8 $\mu$m MSX image.

Figure A.128: Same as Fig. A.15 for Mol075/IRAS 18511+0146 A1/1.

Figure A.129: Same as Fig. A.3 for Mol075/IRAS 18511+0146 A2.

A.34 Mol098/IRAS 19092+0841

Figure A.130: Same as Fig. A.1 for Mol098/IRAS 19092+0841.

Figure A.131: Mol098/IRAS 19092+0841. ISOCAM 15 $\mu$m emission image with superimposed contour plot of the 450 $\mu$m emission.

Figure A.132: Mol098/IRAS 19092+0841. ISOCAM 15 $\mu$m emission image with superimposed contour plot of the 850 $\mu$m emission.

The brightest peak in the MSX images is A (Fig. A.130), which is further resolved in the ISOCAM images (A.131 and A.132). Both the 450 and 850 images show a bright peak SE of A with a very faint halo structure. Since also in this case the size of the halo is less than the size of the map, and the background level of the entire image seems to be close to 0 and relatively constant, we estimate the far-IR contribution of the extended emission in the usual way: the peak integrated flux is subtracted from the integral of the whole structure at 850 $\mu$m and then extrapolated at 60 and 100 $\mu$m.

There is no ZAMS model that can fit A+MM1 (see Fig. A.133 for the closest match represented by the dashed-triple-dotted line), so a two-components model is used.

Figure A.133: Same as Fig. A.15 for Mol098/IRAS 19092+0841 A - 1.

A.35 Mol116/IRAS 20062+3550

Figure A.134: Same as Fig. A.1 for Mol116/IRAS 20062+3550.

Figure A.135: Same as Fig. A.2 for Mol116/IRAS 20062+3550.

Two peaks are visible at 8 $\mu$m, but only A survives up to 21.4 $\mu$m (Fig. A.134). The 850 $\mu$m image (Fig. A.135) shows a bright peak coinciding with A and a faint underlying halo. Since this halo is smaller than the map extent, we can plausibly estimate its contribution to 60 and 100 $\mu$m in the usual way.

The ZAMS fit is very good (Fig. A.136); there is a slight deficiency at 850 $\mu$m, which however is within acceptable limits.

Figure A.136: Same as Fig. A.3 for Mol116/IRAS 20062+3550 A - 1.

   
A.36 Mol118/IRAS 20106+3545

Figure A.137: Same as Fig. A.1 for Mol118/IRAS 20106+3545.

Figure A.138: Same as Fig. A.2 for Mol118/IRAS 20106+3545.

There are two peaks in the 8 $\mu$m image (Fig. A.137). The northern one (B) is large and disappears at 21.4 $\mu$m; the 2MASS images show some stars and a certain level of diffuse emission in $K_{\rm s}$ so B is probably associated with a small cluster of relatively evolved stars. Source A coincides with a bright and slightly elongated 850 $\mu$m peak, which is resolved into two peaks (MM1 and MM2) at 450 $\mu$m. It is relatively easy to separate the two peaks at 450 $\mu$m; at 850 $\mu$m (Fig. A.138) we perform the JMFIT constraining the relative distance between the two cores as seen at 450 $\mu$m. The procedure is clearly uncertain and the results should be taken with care. The extended emission extends to the edge of the map and hence we cannot be sure of the correct background subtraction in the SCUBA jiggle map, so no estimate of the potential contribution of extended emission will be made. The IRAS fluxes are then split between the two sources in proportion to the 850 $\mu$m flux (it is basically the same if we use the 450 $\mu$m flux), keeping in mind that since there is no correction for extended emission contamination (which appears to be present in the submm images) the assigned IRAS fluxes should not be exceeded.

As in the case of IRAS 10095-5843 A (see Sect. A.11), the best ($\xi=0.09$) model fit (Fig. A.139) does a good job overall, but misses the 100 $\mu$m point by more than a factor of 2. In principle this may be acceptable given the unknown contribution of extended emission to the 100 $\mu$m flux; however, since this is unknown, we prefer to verify whether an alternative model can be found where the 100 $\mu$m departure from the data is more limited. Examining the 10 best $\xi$ models we indeed find one, with $\xi=0.11$, nearer to the 100 $\mu$m flux, while still being well within tolerances in the mid-IR.

MM2 is successfully fitted with a greybody (Fig. A.140); the fluxes shortward of 60 $\mu$m are well below the fluxes of source A.

Figure A.139: Same as Fig. A.3 for Mol118/IRAS 20106+3545 A/1.

Figure A.140: Same as Fig. A.9 for Mol118/IRAS 20106+3545 MM2.

A.37 Mol136/IRAS 21307+5049

Figure A.141: Same as Fig. A.1 for Mol136/IRAS 21307+5049.

Figure A.142: Same as Fig. A.132 for Mol136/IRAS 21307+5049.

The MSX images show one peak with a bar-like structure to the NE and running SE-NW (Fig. A.141). There are radio sources within a distance of 20$\div$30 $^{\prime \prime}$. The ISOCAM higher resolution map (Fig. A.142) clarifies the point-like nature of the mid-IR peak. The 850 $\mu$m image (contours in A.142) shows a bright peak nearly coincident with the mid-IR with an underlying halo which we use as usual to estimate the contribution of the extended emission to far-IR wavelengths.

The fit to the entire SED assuming the mid-IR and the submm peaks as associated gives a very good ZAMS model (Fig. A.143). The IRAS 12 $\mu$m flux is considered an upper limit due to the contribution of the extended features visible in the MSX and ISOCAM images.

Figure A.143: Same as Fig. A.3 for Mol136/IRAS 21307+5049 A/1.

A.38 Mol139/IRAS 21519+5613

Figure A.144: Same as Fig. A.1 for Mol139/IRAS 21519+5613.

Figure A.145: Same as Fig. A.2 for Mol139/IRAS 21519+5613.

There are three sources visible at 8 $\mu$m, but only A is also visible throughout 21.4 $\mu$m (Fig. A.144). The 850 $\mu$m image (Fig. A.145) shows a bright peak coincident with A. The integrated flux is an estimate assuming a peak+background combination in JMFIT. The structure of the background in the image is not regular and the extent of the map is not sufficient to obtain a reliable estimate of the extended emission, if any (we then keep the 100 $\mu$m as an upper limit in the SED fit); we then allow for a variable (i.e. not flat) background in the JMFIT analysis. As a result the peak integrated flux estimate is highly uncertain.

The fit of a ZAMS models to A+MM1 provides a good fit (Fig. A.146), with a deficiency $\sim$50% the submm; we deem this acceptable also given the uncertainties in the peak flux estimate.

Figure A.146: Same as Fig. A.3 for Mol139/IRAS 21519+5613 A/1.

A.39 Mol143/IRAS 22172+5549

Figure A.147: Same as Fig. A.1 for Mol143/IRAS 22172+5549.

Figure A.148: Same as Fig. A.2 for Mol143/IRAS 22172+5549.

There is one peak in the mid-IR (Fig. A.147), superimposed on a filamentary structure oriented NE-SW, which makes the source flux estimates difficult in MSX bands A and C, and makes us assume the IRAS 12 $\mu$m flux as an upper limit in the SED fit. The mid-IR source is very close to a radio continuum source, and coincides with an 850 $\mu$m peak (Fig. A.148) superimposed on a patch of extended emission that overlaps the mid-IR filament. The overall background in the submm image has negative values areas; probably there were subtraction problems in the SCUBA jiggle map and the whole absolute level of the image is quite uncertain (the integrated peak flux should be reliable, though). Since we cannot estimate the extended emission contribution from the submm map, we will assume the IRAS 100 $\mu$m flux as an upper limit in the SED fit.

We can find a ZAMS model that reasonably well fits the mid-IR and submm counterparts (Fig. A.149).

Figure A.149: Same as Fig. A.3 for Mol143/IRAS 22172+5549 A/1.

A.40 Mol148/IRAS 22305+5803

Figure A.150: Same as Fig. A.1 for Mol148/IRAS 22305+5803.

Figure A.151: Same as Fig. A.2 for Mol148/IRAS 22305+5803. the millimeter image is at 1.2 mm taken with the IRAM 30 m telescope.

There is one mid-IR peak visible at all wavelengths (Fig. A.150); it coincides with a $K_{\rm s}$ cluster and a radio continuum source. For this field the millimeter image was taken at 1.2 mm (Fig. A.151) using the IRAM 30 m telescope. The mm peak, coinciding with the mid-IR source, lies on top of a halo and the extent of the map is such that a reliable background subtraction is possible; we can then reliably estimate the contribution to far-IR wavelengths of the extended emission and subtract it from the IRAS fluxes.

The fit to the entire SED is very good with a ZAMS model (Fig. A.152).

Figure A.152: Same as Fig. A.3 for Mol148/IRAS 22305+5803 A/1.

A.41 Mol151/IRAS 22506+5944

Figure A.153: Same as Fig. A.1 for Mol151/IRAS 22506+5944.

Figure A.154: Same as Fig. A.2 for Mol151/IRAS 22506+5944.

There is one mid-IR peak (Fig. A.153) nearly coinciding with a strong 850 (Fig. A.154) and 450 $\mu$m peak (MM1). The brightness profile analysis via cross-cuts through the peak reveals the presence of a halo and the core probably is made up of two components. The resolution is however not sufficient to deblend them, so we fit it with a single Gaussian with a plateau. We find here the usual difficulty, as with most of the SCUBA images, that the background level is not stable and reflects off-source subtraction difficulties in the jiggle maps. We prefer not to make any extended emission correction to the far-IR fluxes. There is another peak (MM2) to the SE but it is very faint and irrelevant for the present analysis.

In spite of all the uncertainties we cannot find a ZAMS model that fits the mid-IR and MM1 counterparts together; the closest match we find is reported in Fig. A.155 with the line coding as usual. The two component model is successful. In principle the fitted greybody shown as a dashed line in Fig. A.155 is an upper limit to the real one, because of no-estimated contamination of the extended emission to the 100 $\mu$m IRAS flux, implying that the mass and the luminosity of the MM1 core will be upper limits. This, however, has no substantial impact on the conclusions based on Fig. 9, because M and L will scale with the absolute greybody level in the same way; the MM1 core placement on that diagram may be different from the one implied by the fit in Fig. A.155, but still it will be basically along the strip identified by the other MM sources.

Figure A.155: Same as Fig. A.15 for Mol151/IRAS 22506+5944 A/1.

   
A.42 Mol160/IRAS 23385+6053

Figure A.156: Same as Fig. A.1 for Mol160/IRAS 23385+6053.

Figure A.157: Same as Fig. A.132 for Mol160/IRAS 23385+6053.

Figure A.158: Same as Fig. A.8 for Mol160/IRAS 23385+6053 MM.

This source has been object of detailed studies. The MSX (Fig. A.156) and the ISOCAM images (Fig. A.157) show a complex pattern of emission with several components that tend to faint at 21.4 $\mu$m. The mid-IR patches are coincident (over the entire extent) with a large stellar cluster seen in the near-IR, and with two radio continuum sources. The 850 $\mu$m image (the contours in A.157) shows a strong peak on a halo, which is positionally anti-correlated with the mid-IR emission; the mm peak is also detected in interferometry at 3.4 mm and does not show evidence of radio continuum.

The only source we try to fit is the cold submm peak; the far-IR fluxes used are the IRAS ones corrected for the contribution of the extended emission estimated in the usual way (for this source the background is regular and around 0 so it is credible. We cannot find a ZAMS model that fits the SED (constrained in the mid-IR by the upper limits estimated in ISOCAM images at the mm peak location); the closest match is indicated as usual as a dashed-triple-dotted line in Fig. A.158. We successfully fit a single greybody to the $\lambda\geq60$ $\mu$m fluxes.