WIDE-FIELD INFRARED SURVEY EXPLORER OBSERVATIONS OF THE EVOLUTION OF MASSIVE STAR-FORMING REGIONS

We have developed a scheme using WISE photometry to identify young stars from their mid-IR colors and magnitudes, in the same way as has been performed in studies using Spitzer photometry (for example, Gutermuth et al. 2008).

At ages 1–2 Myr, 60%–70% of young stars possess optically thick disks which boost their near- and mid-IR emission above the photospheric level (Strom et al. 1989; Haisch et al. 2001). By an age of 5 Myr, the fraction of low-mass stars (0.1–1.2 M_☉) with optically thick disks drops to around 20%, while those with mass >1.2 M_☉ exhibit almost no excess emission at wavelengths λ < 16 μm (Carpenter et al. 2006; Dahm & Hillenbrand 2007). No more than 10% of stars older than 10 Myr exhibit any near-IR and 3–8 μm excess emission from circumstellar disks (Strom et al. 1989).

Clearing of the inner parts of the disk by photoevaporation by the star or by planet formation produces stars that exhibit excess emission only at wavelengths >20 μm, the so-called transition disks. The absence or low level of excess emission from 1 to 10 μm and large excess beyond 20 μm has been modeled as arising from a truncated optically thick outer disk and a dissipated inner disk or a radial gap in the inner disk (for example, Calvet et al. 2005; D'Alessio et al. 2005). These "transition disks" may be created by binary systems; however, recent work by Pott et al. (2010) has shown that these may be rare occurrences and planet formation or brown dwarf companions may be responsible for the inner disk clearing (Andrews et al. 2011). The ages of transition disks are uncertain. Their appearance has been noted in young regions of the Taurus star-forming cloud by Luhman et al. (2010). Giant planet formation may be able to produce disk clearing on timescales <5 Myr. Photoevaporation of the disk may take a longer time to operate (for example, Alexander & Armitage 2009).

At a still later stage of evolution, debris disks appear. These are likely the result of disks replenished in small dust by collisions between planet-size bodies, and are thus formed by second generation dust, while transition disks are evolved primordial disks that have been cleaned in their inner parts (Currie et al. 2008; Kenyon & Bromley 2004). However, these disks have only weak excesses and are unlikely to be detected at the large distances (>1 kpc) surveyed in this paper as WISE cannot see optically thin debris disks for stars where the photosphere could not be detected by band 4 at 22 μm.

These properties of young stars make them possible to distinguish from field stars—which will have bluer spectral energy distributions (SEDs) in the infrared wavelength range—and crudely classify them by their colors. Even though young stars that have completely lost their circumstellar disks will be missed with this technique, clusters and stellar groups younger than 10 Myr are still identifiable in the mid-infrared by their disk frequency which is larger than in the old field population.

Classification of sources is performed using a multi-phase mid-IR colors method. The basis for this scheme is the work of Allen et al. (2004) and Hartmann et al. (2005) who first investigated the distributions of young stars in Spitzer color space. The scheme for identification of young stars is similar to those in Spitzer studies that use mid-IR colors to distinguish YSOs and mitigate contamination (e.g., Gutermuth et al. 2008, 2009). In the Appendix, we detail our adapted WISE method in full. In summary, we use the numerous available flux ratios, or colors, to identify and classify candidate YSOs and, to the extent possible, mitigate the effects of contamination and reddening. The scheme makes use of combined WISE and Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) data. We match the extracted WISE source lists for each field to the 2MASS JHK_S point-source catalog. We match using a 3'' search radius and remove 2MASS objects with photometric quality flags D, E, F, or X in any band. These flags refer to the following: D—the signal-to-noise ratio <5; E—profile-fit photometry was very poor or did not converge; F—a reliable estimate of the photometric error could not be determined; and X—no valid brightness estimate could be extracted using any algorithm.

Contamination arises from non-YSO sources, including star-forming galaxies, broad-line active galactic nuclei (AGNs), unresolved knots of shock emission from outflows colliding with cold cloud material, planetary nebulae, and asymptotic giant branch (AGB) stars (both carbon-rich and oxygen-rich; Robitaille et al. 2008).

Once separated from the diskless photospheres (Class III pre-main-sequence stars and unrelated field stars) and the aforementioned contaminants, the YSOs are classified into the canonical categories of Class I YSOs (protostars with infalling envelopes, including flat spectrum objects) and Class II YSOs (pre-main-sequence stars with optically thick disks), with supplemental categories of "deeply embedded sources" (added to the Class I tallies) and "transition disks." Class I YSOs have rising or flat SEDs from 3.3 to 22 μm, while Class II YSOs have a falling SED slope (Lada 1987; Greene et al. 1994). Transition disks are objects with optically thick excess emission at long wavelengths and little to no excess at short wavelengths (Strom et al. 1989). As described in the Appendix, we require detections in WISE bands 1, 2, and 3 or 2MASS K_S plus WISE bands 1 and 2 to eliminate background galaxies and classify sources as Class I or II YSOs. The transition disk category requires a detection in WISE bands 1, 2, and 4. We note that Class II YSOs viewed at high inclination will have similar SEDs to Class I YSOs. Edge-on Class I YSOs can have characteristics of Class 0 YSOs (an earlier evolutionary stage with a more massive infalling envelope and high internal extinction; André et al. 1993; Robitaille et al. 2006). We emphasize that all Class I, II, and transition disk objects found in this paper are YSO candidates, but hereafter refer to them simply as YSOs.

To calculate the rate of contamination in our YSO samples, we used three non-Galactic plane fields (the Boötes SDWFS field of Ashby et al. 2009 and two polar fields, all of ≈10 deg²). The combined range of these fields, including Poisson errors, is 1.4–2.9 Class I sources, 2.2–5.4 Class II sources, and 0.9–2.7 transition disk sources deg⁻². This derived contamination rate is generally applicable, but likely overestimates the amount of contamination in regions of high near-IR extinction or heavy cirrus confusion. Without these foreground issues, extragalactic source confusion should be isotropic. We note that our scheme will not generally exclude all galaxies, but at the same time it will also overcompensate for contamination because faint YSOs will also be removed from consideration. A further step of looking for "extended source" status in the 2MASS catalog (which we do not perform here) will exclude most resolved galaxies, but those point-like at 2MASS resolution will remain contaminants. These typically consist of centrally concentrated star-forming galaxies at a redshift of 0.1–0.4 (D. L. Padgett et al. 2011, in preparation; Rebull et al. 2010).

Field Galactic AGB stars and Classical Be stars (CBe; Porter & Rivinius 2003) can also mimic the properties of young stars in surveys such as this one (Robitaille et al. 2008). With their weak infrared excesses above the photosphere, they often resemble "anemic" disks (Lada et al. 2006) with excess at all infrared wavelengths but at a level below that of optically thick T Tauri stars or Herbig AeBe star accretion disks. CBe stars are most commonly in the spectral type range O9 to B3 and thus will be bright but relatively rare. Similarly, Robitaille et al. (2008) estimated a rate of a few AGB stars per square degree at most from the GLIMPSE survey region in the inner Galaxy. Since our fields are in the outer Galaxy, the density of AGB stars will be lower. We expect ∼few contaminating AGB and CBe stars in each cluster studied.

In Table 2, we list the number of YSOs classified as Class I, II, and transition disks within a radius R_max, defined as the by-eye radius that encloses the visible extent of the cluster (for NGC 659, 663, and 957) or the outer edge of the expanding shells as traced out by the bright 12 and 22 μm cirrus emission in the younger regions. The location and size of these radii R_max are shown in Figure 1 as white dashed-line circles. The centers of these circles are adjusted by eye to match the observed mid-IR morphologies at 12 and 22 μm and are discussed in detail in Section 4. The W 5 region is divided into the clusters IC 1848 and W 5 East, with the coordinates given in the SIMBAD online database for HD 18326 for W 5 East. The W 4 region is divided into a region roughly centered on the cluster IC 1805, and another to its south to trace out the apparent second bubble of emission in this location. The Bk 59 region is divided into an eastern and a western bubble.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass completeness limit for our survey is dictated by our photometric completeness and the wide span of YSO mid-IR colors. The primary constraint in completeness for the samples of YSOs is the 12 μm detection limit, since this band has low sensitivity relative to WISE bands W1 and W2 and very bright background contamination from polycyclic aromatic hydrocarbon (PAH) feature emission. The additional YSOs found with the K_S − [3.4] − [4.6] color–color diagram (see the Appendix) helps to mitigate the loss of the faintest YSOs. We find from matching the sample of detected YSOs in the W 5 region to that of Koenig et al. (2008), the WISE sample of YSOs is 90% complete to a magnitude of 12 in the 2MASS K_S band. We are able to recover roughly 30% of the YSOs (Class I, Class II, and transition disk) found with Spitzer and listed by Koenig et al. (2008) with our WISE data. Using the evolutionary models of Siess et al. (2000), at age 2 Myr, at 2 kpc with extinction A_V = 2 (characteristic of the H ii region cavity in W5; Hillwig et al. 2006), K_S ∼12 equates to a limiting mass 3.75 M_☉. However, since we are dealing in the majority with T Tauri and Herbig AeBe-like stars that possess infrared excess at near- and mid-infrared wavelengths, we must take into account the excess emission produced by the accretion disks around these stars. Meyer et al. (1997) demonstrate that T Tauri stars typically have a K-band excess between 0.1 and 1.1 mag. Taking a median value excess of ΔK = 0.65, the completeness limit converts to a mass of 2.0 M_☉. We carry out the same estimation for each region in our sample and summarize the results in Table 2.

As a check on these results we also compared our data for IC 1396 to the work of Barentsen et al. (2011) and Sicilia-Aguilar et al. (2006), and our data for IC 1795 to the work of Roccatagliata et al. (2011). In IC 1396, we recover all of Barentsen et al.'s YSOs down to a computed mass of 1 M_☉, and we are 90% complete in matching to Sicilia-Aguilar et al.'s sample to a K_S magnitude of ≈13. With extinction A_V = 2, this similarly converts to a mass ≈1 M_☉ at 4 Myr in the models of Siess et al. (2000). In IC 1795 our 90% completeness limit (after applying the same WISE and 2MASS photometry requirements that we use to find YSOs) is at a Spitzer [3.6] magnitude of 10.5 when compared to the sample of Roccatagliata et al. Assuming a typical YSO K_S − [3.6] excess of 1 mag and A_V = 2, this limit converts to a mass between 2 and 3 M_☉ at an age of 1 Myr, dependent on the level of K_S-band excess. At an age of 3 Myr—the lower limit age derived by Roccatagliata et al.—the limit becomes 6 M_☉, much higher than our quoted value. However, their sample included many diskless YSOs found through their Chandra X-ray emission, to which we are not sensitive in this work. We thus note that the completeness limits in Table 2 should only be used for a relative comparison between regions in the survey, given the strong positional dependence within regions of bright nebular background, extinction, and age, all of which influence the completeness calculation.

The resolving power of the WISE telescope means that as we sample more and more distant regions, the likelihood that there will be multiple sources in the beam increases. At 1 kpc (the distance to the nearest regions in the sample), the 6'' resolution of bands 1–3 corresponds to a distance of about 0.03 pc, while the 12'' resolution at band 4 corresponds to 0.06 pc. In the Orion Nebula Cluster (ONC), with a volume density of 2 × 10⁴ pc⁻³, the peak spacing between cluster members is about 35, or about 0.007 pc, with a median nearest neighbor distance of 0.09 pc (using the full compilation of Rebull et al. 2006, supplied by the authors). This result implies that we are likely confused in the densest clustered regions in each field. However, the ONC sample extends down to stellar masses ∼0.04 M_☉, while at best our survey completeness is 0.7 M_☉ as given in Table 2, which will mitigate this effect somewhat. Of the 2500 sources with a listed mass in the data of Rebull et al. in the inner part of the ONC (R.A. = 835–842 and decl. = −64 to −45), we find that considering only sources with mass >0.7 M_☉, the peak nearest neighbor distance is 0.16 pc. Considering only sources with mass >1.0 M_☉, the peak nearest neighbor distance is 0.26 pc.

We urge caution prior to follow-up of the YSO lists presented in this paper. An additional source of contamination is from confusion caused by the change of resolution from short to long wavelengths (i.e., what was star A at W1 is star B at W4 22 μm). Star/galaxy superposition is also significant, as is star/nebula superposition. A star on top of a nebula can have a "fake" excess if the nebula gets steadily brighter in proximity to the star (since the background subtraction will then be wrong). For example, an object resembling a transition disk (a stellar SED at short wavelengths but an excess at 22 μm) could be mimicked by the combination of a foreground or cluster member stellar source and a background galaxy that only appears at the longest wavelength. A direct comparison of the WISE detected source list with the high-resolution (<1'') r'z' imaging data of Koenig & Allen (2011) in the W 5 region shows that at 2 kpc, 10% of sources are actually multiple. We also tested the results of our classification scheme by comparing the YSO catalog of Koenig et al. (2008) for the W 5 region from Spitzer with that produced by our WISE algorithm. Of 852 candidate YSOs identified by WISE that matched with Spitzer sources, 141 (so roughly 16%) appeared to be spuriously produced by the effects described above. In particular, the list of WISE transition disks was 90% contaminated by false detections. We thus consider all transition disks in regions beyond 1 kpc to be unreliable; however, due to their low numbers in each region this fact does not alter our analysis.

In Table 2, we use the YSO counts to derive a star formation rate for the star-forming regions within 2 kpc. In each case we first extrapolate the number of YSOs to a full 0.1–100 M_☉ initial mass function (IMF) using a standard Salpeter IMF (Salpeter 1955). We assume that the YSOs we have detected range in mass from our 90% completeness limit up to 8 M_☉ (an approximate upper limit to the masses of Herbig AeBe stars; Hernández et al. 2006). Since we are only detecting stars with disks and are missing those young stars that have already lost their disks, this calculation gives the amount of recent star formation. Following the example of Heiderman et al. (2010) and Lada et al. (2010), we estimate the star formation rate by adopting an average mass of 0.5 M_☉ (see, for example, Chabrier 2003) and a timescale of 2 Myr (the time between the formation and end of the Class II phase; see Evans et al. 2009). Thus, the star formation rate is

$$\begin{equation} \hbox{SFR} = \frac{N_{{\rm YSO}} \times 0.5\, M_\odot }{2\hbox{ Myr }} \hbox{ (}M_\odot \hbox{ yr$^{-1}$).} \end{equation} \tag{ 1 }$$

Although we could plausibly derive this parameter for all the regions surveyed, in the more distant clouds the lack of completeness combined with the likely different timescale for disk erosion in the more massive young stars makes the calculation more uncertain.

We describe the regions we have surveyed in order of increasing age as listed in Table 1. In Figure 1, we show the results of the survey. For each region we plot two figures: a color composite of the four WISE bands, overlaid with the CO data from Leisawitz et al. (1989) and a grayscale WISE 12 μm image overlaid with the YSOs that we have classified for each region and known O and B stars from the literature.

The key observations in this sample are (1) a rough trend of decreasing Class I YSO content and decreasing CO emission in the central portions of each region, (2) ubiquitous presence of bright arcs and pillars of dust emission coincident with condensations of ¹²CO emission at the periphery, (3) clusters of Class II YSOs around the central O stars in the largest regions, and (4) smaller embedded clusters of YSOs in the peripheral dust arcs and pillars.

3.5.1. NGC 2175

3.5.2. W 3 (IC 1795)

3.5.3. NGC 1893

3.5.4. Berkeley 59

3.5.5. NGC 7380

3.5.6. W 4 (IC 1805)

3.5.7. NGC 281

3.5.8. IC 1396

3.5.9. Roslund 4

3.5.10. NGC 1624

3.5.11. W 5 (IC 1848 and W 5 East)

3.5.12. NGC 957

3.5.13. NGC 659

The regions we have assembled in this study suggest an evolutionary sequence by their morphology. At the young end are regions where the YSO clusters are dominated by Class I sources (like those in W 3, NGC 1624, NGC 281, and NGC 2175). Intermediate cases show a low-density cavity/H ii region that has been swept out by the first generation of star formation (as in W 4, W 5, NGC 7380, Roslund 4, Bk 59, NGC 1893, and IC 1396), with YSOs in pillars at the edges of shells of interstellar gas. The final end state of massive star-forming regions is practically devoid of YSOs (as in the regions around open clusters like NGC 957, NGC 663, or NGC 659). The fraction of YSOs in the Class I category appears to follow this morphological scheme with NGC 2175, NGC 281, IC 1795, and NGC 1624 having 30.2%, 23.1%, 20.4%, and 14.3% of their YSOs as Class I. The ages listed for the regions in Table 2 do not exactly agree with this morphological scheme, with NGC 281 and NGC 1624 at 3.5 and 4 Myr, respectively. However, these more distant regions will suffer more strongly from incompleteness in the Class II sample than in the Class I population (since the latter are brighter in the mid-IR) and so may have an artificially high Class I:Class II ratio. In Figure 2, we plot the fraction of stars in Class I relative to Class I and Class II combined (dropping the transition disk candidates which we consider unreliable) as a function of age. The error bars represent simple $\sqrt{N}$ counting statistics and do not account for completeness or possible contamination. A trend of declining ratio with age is visible by eye, however, each region has multiple episodes of star formation within it (and so an age range may be more appropriate rather than a single value), and the ratio is affected by the level of bright background which would likely lower the ratio in the younger regions like NGC 2175 (the leftmost point in the figure). The large intrinsic age spreads combined with the large uncertainties in determining stellar ages also mean that the trend as presented here is not statistically robust. Continuously occurring star formation and the timing of the appearance of secondary generations of star formation will also alter the relationship between the ratio and (central cluster) age.

Zoom In Zoom Out Reset image size

The later episodes of star formation around the edges of the young massive regions are frequently seen in pillar or mountain-like structures that point back toward the ionizing sources that created the cleared out H ii regions. However, they are not all the same sizes. Within the same region, a star-forming clump at the edge of the H ii region bubble may contain a cluster with many young stars forming within it, or instead a solitary star or very small cluster. The outcome of these secondary generations of star formation is likely driven by the initial properties of the molecular gas at these locations, prior to their formation in the process described by Thompson et al. (2002). We show a sample of image cutouts around the peripheral clusters in the regions of this survey in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Froebrich et al. (2005) discovered that the mass of gaseous globules detected within their extinction map of IC 1396 increased with increasing distance from the central O star HD 206267. They also tentatively claim that the observed star density (number of stars versus the projected area of the globules) falls off with increasing distance from the central star or stars. We find in IC 1396 and other regions like Bk 59 that stars are forming in both large and small clumps at the same projected distance from the central stars. For example, the Bk 59 panels in Figure 3 show that large structures can form clusters while at the same time, at likely the same distance from the central O star, small globules are forming small clusters or single stars. This observation suggests that star formation at the boundaries of H ii regions may be strongly influenced by the mass and perhaps the density structure in the molecular core at the time it encounters the O star wind and pressure from the adjacent freshly ionized gas.

Triggered star formation is a process that has been proposed and discussed at some length in the literature but still remains only a plausible hypothetical process (e.g., Elmegreen 1998). From the observed distributions of young stars in this sample of star-forming regions, we attempt to assess how well the ensemble average distribution follows a pattern consistent with passive, stochastic star formation, or actively triggered star formation.

To investigate these two cases, we measure the distribution of YSOs found in the survey in each region except the three older open clusters (NGC 659, 663, and 957). In each case, we derive the surface density histogram as a function of radial distance from the cluster centers as listed in Table 2. We calculate the distribution out to the radius 1.5× R_max, where R_max represents the visually identified boundary of the expanding H ii region shell. The 12 and 22 μm emission that traces these arcs and rings is likely produced by dust and/or PAH grains illuminated on the surface of cold, dense gas. We interpret these circles as expanding bubbles and investigate the distribution of young stars picked out by our selection criteria within them. In the case of W 3/W 4/W 5 there are five circle centers (IC 1795, IC 1805, IC 1805 South, IC 1848, and W 5 East). In W 4 and Bk 59 there are two. To calculate the surface density histograms, we measure the radial distance to each star from the assumed center, normalized to R_max, and then derive the surface density of stars in annuli around the cluster center of width 0.1 dimensionless distance units. We use only the lists of Class I and Class II sources, since the transition disk list is more unreliable and only contributes a small number of objects.

In Figure 4, we show the results for the individual clusters. In this sample, the young regions present similar configurations of star formation, in that more embedded clusters of stars often coincide with the edges of the bubbles we have traced out. This distribution can be seen in the peak-up in space density at normalized radii ∼1 in NGC 7380, NGC 2175, IC 1805 South, Roslund 4, and Bk 59 East. However, several regions, NGC 281, NGC 1893, IC 1396, Bk 59 West, IC 1805, IC 1848, and W5 West, also contain large, roughly central clusters of YSOs that dominate the density distributions at small and intermediate radii.

Zoom In Zoom Out Reset image size

Next we construct an averaged surface density histogram to sum over the peculiarities of this sample of regions and arrive at a more general measurement of the configuration of young stars in massive star-forming regions. The average plot shown in Figure 5 is produced by adding all young stars from all regions at their normalized radii from the cluster centers and calculating a new surface density distribution from that. The act of combining the distributions of these several regions that are likely at different ages and different stages of star formation produces a time-average of the bubble expansion and triggering process.

Zoom In Zoom Out Reset image size

In the case of passively occurring star formation, although new young stars will not uniformly fill space as they form inside a giant molecular cloud, averaged over many regions, multiple bursts of star formation arranged randomly might approximate a uniformly filled sphere. Although an individual cloud's internal structure might not be uniform and smooth (consisting of filaments or cores, etc.), it could, to a zeroth-order approximation, be uniformly filled. The average region initially might have contained a giant molecular cloud whose size is comparable to that of the regions studied (∼50 pc or so). We calculate the expected surface density distribution of a uniform filled sphere of size R_max at impact parameter x:

$$\begin{equation} \hbox{Surf. Dens. } \propto \sqrt{1 - \left(\frac{x}{R_{{\rm max}}}\right)^{2}}. \end{equation} \tag{ 2 }$$

We also derive the distribution produced by a uniform empty shell of size R_max with a thickness 2% of its radius in bins of size Δx:

$$\begin{eqnarray} \hbox{Surf. Dens.} &\propto& \frac{R_{{\rm max}}}{x \Delta x}\left(\sqrt{R_{{\rm max}}^2 - x^2 - \left(\frac{\Delta x}{2}\right)^2 + x\Delta x}\right.\nonumber\\ &&\left. - \sqrt{R_{{\rm max}}^2 - x^2 - \left(\frac{\Delta x}{2}\right)^2 - x\Delta x} \right). \end{eqnarray} \tag{ 3 }$$

We overplot these two simple cases in the left panel of Figure 5. Neither the filled nor the empty shell model comes close to describing the measured average YSO distribution.

At small radii (≪R_max), the central clusters in the assembled regions plausibly dominate the measured average surface density distribution. Thus, we also compare the measured distribution with a King profile (King 1962). The full expression for a King profile is

$$\begin{equation} \hbox{Surf. Dens. } \propto \left[ \frac{1}{ \sqrt{1 + \left(\frac{x}{r_c}\right)^2}} - \frac{1}{ \sqrt{1 + \left(\frac{r_t}{r_c} \right)^2}} \right] ^2, \end{equation} \tag{ 4 }$$

where in Figure 5 (right-hand panel) we have set r_c, the scale factor or "core radius," to 0.35 and r_t the limiting or tidal radius to 1.3. Only the central portion of the average YSO density distribution follows the shape of a King profile; a King-like distribution drops off far too rapidly in density with distance from the center to match the data. The King profile was originally developed to match the density distribution of gravitationally relaxed stellar groups such as globular clusters. The lack of a match in the case of the present survey is not surprising and simply reflects that these central clusters are not dynamically relaxed associations of stars. Lada et al. (1984) found that clusters that remain bound after gas dispersal undergo significant expansion as they emerge from a cloud and in the process may lose 10%–80% of their stars. Such evaporating stars would contribute to a smooth, flatter-than-King surface density profile. However, evaporating stars would not tend to cluster in the molecular globules and elephant trunks that we see at the periphery of the studied regions.

Since the actual density appears to fall off smoothly, we also investigated power-law models for the volume density distribution of YSOs, with density ρ(r) = ρ₀(r/R_max)ⁿ, where r is the three-dimensional distance from the cluster center. We found that an r⁻¹ law corresponds to a surface density distribution similar to the measured distribution, as shown in the center panel of Figure 5. The analytical form of the profile is

$$\begin{equation} \hbox{Surf. Dens. } = 2 \rho _0 R_{{\rm max}} \ln \left[ \frac{R_{{\rm max}}}{x} + \sqrt{ \left(\frac{R_{{\rm max}}}{x}\right)^{2} - 1}\right]. \end{equation} \tag{ 5 }$$

Given that ρ∝r⁻¹ approximates the YSO volume density distribution, we can integrate to derive the total count of YSOs, C = 2πρ₀R³_max, and the count within radius r₁, C(r₁) = 2πρ₀R_maxr²₁. The fraction of stars that lie within the inner half-radius r₁ = 0.5 R_max is thus 1/4. In other words, three-fourths of the YSOs in the spherical distribution lie at a radial distance r > R_max/2.

We note that the cleared out H ii region cavities in these massive star-forming regions will have much lower line-of-sight extinction than the peripheral embedded clusters and so the rate of contamination by background galaxies and stars will be much higher there (radii 0.4–0.8 in Figure 5). Since the peripheral new embedded clusters form in regions of dense gas and thus high extinction, the young stars that we find there are probably less affected by contaminants. Accounting for this effect would likely accentuate the peak-up in stellar density at radii ∼1. At all radii our survey is incomplete, since our YSO identification technique misses cluster members that no longer possess an infrared excess and is compromised on regions of bright nebulosity. Thus, the true averaged distribution may not exactly resemble the ρ∝r⁻¹ trend.

We now consider what this result suggests in regard to models of triggered star formation. Two major paradigms are relevant here. The first is the "collect and collapse" process of triggered star formation, where the initial generation of stars forms and then material is swept away from it by the massive stars' winds and ionizing radiation (see Elmegreen 1998 for a summary of the relevant literature). Once a critical density of material has accumulated in the swept-up gas, new bursts of star formation occur. The varying initial distributions of molecular gas in the molecular cloud will likely produce broken rings or arcs of swept up material rather than complete shells. More importantly for our analysis, this process would be likely to produce an average distribution with two components to it: the interior, initial cluster and a separate exterior shell of the second, triggered generation of stars. The second mechanism is the so-called RDI of pre-existing overdensities in the molecular cloud. As sketched out by Reipurth (1983) and developed analytically by Bertoldi (1989) and Bertoldi & McKee (1990), regions of higher density are eroded more slowly as the initial burst of massive star formation ionizes and pushes away the molecular cloud, leaving them exposed. The pressure of the ionized H ii region compresses the revealed clump, causing new stars to collapse therein, as it is pushed farther away. If we assume that the distribution produced by the collect and collapse process will consist of a central peak in density, a lower density surrounding region and a second peak up at the cluster boundary, while the RDI process will not result in such a two-component distribution, we argue that the fact that our data is approximated well by a simple density law more strongly supports RDI as the dominant process of triggering in massive star-forming regions.

Recent simulations of massive star formation have been presented by Dale & Bonnell (2011). The Dale et al. simulations propose that the collapse of a giant molecular cloud will create a filamentary structure of dense gas within it, which collapses still further to form stars. The evacuated bubbles and peripheral embedded clusters in pillars and elephant trunks seen in massive regions—such as those in this survey—represent the influence of the central massive star clusters on pre-existing structures in the gas (for example, RDI) rather than an accumulation and collapse process.

We argue that we are seeing the relic of the initial filamentary structure of the gas in the distributions of the YSOs in these regions. The r⁻¹ distribution that is seen in the averaged distribution plot is simply the result of combining individual irregularities in the many different clouds we have surveyed. The central, most massive group of stars that forms—containing the most massive stars—erodes the cloud around it, and while stars may already be forming in the dense filaments, the actions of RDI will enhance this star formation by compressing the collapsing structures as the H ii region bubble expands.

Since most of the gas is in the filaments, there is not enough material in the less dense regions to be gathered up as in the collect and collapse process model, although this mechanism cannot be ruled out completely. For example, the expanding bubble could impinge on a filament tangential to its surface and pile up additional material there, producing new star formation.

The distributions of YSOs that we see in this sample of regions suggest that star formation may in general be triggered around massive star-forming clusters rather than passively formed, although it is certainly possible that episodes of spontaneous star formation occur around these regions and even alongside triggered star formation. The most crucial tests of this theory will be a precise and accurate determination of the relative ages of the central and surrounding clusters, a measurement of their three-dimensional proper motions, and establishing a YSO population free of external contamination.

If the r⁻¹ is indeed a good fit to the distribution of stars, we infer that triggered star formation is self-propagating because each "original" young star yields something like three "new" young stars in its neighborhood; YSOs at r > 0.5 R_max outnumber those at smaller radial distances by 3 to 1. If an initial generation of star formation triggers a large enough secondary generation to populate the upper end of the IMF, we postulate a "fireworks hypothesis" whereby some of the new stars formed will be massive and lead to subsequent generations of star formation. In an analogy to a nuclear chain reaction, the process presumably will continue until the fuel source (nearby molecular cloud material) is exhausted or dispersed or until the new generation of stars does not contain enough massive stars to continue the propagation of this process. We note that since our study is only sensitive to objects with infrared excess, it is possible that since the triggered population is younger, and it likely has a higher disk fraction than the central, triggering stars. As we have noted, the exact function of the distribution of stars may not be r⁻¹, and thus the number of new YSOs may be less than three times the initial generation.

This process may be rapid. Assuming that the central clusters we have detected are the first generation to form, the peripheral clusters are evidently forming very soon after that the first stars have not had time to disperse most of their optically thick, mid-IR emitting accretion disks. Whether or not there is a causal relationship between the interior and peripheral populations of star formation, this observation implies that the time difference between the two events is short, given that the half-life for the thick-disk population (i.e., disks that are bright from 3 to 24 μm) is roughly 2 Myr (Hernández et al. 2008). If the bubble expansion ages are also short, this would lend further weight to this picture of rapid successions of star formation. Assuming an expansion velocity of 10 km s⁻¹, the apparent bubble radii of 3–24 pc in our sample would only require at most 2.4 Myr to form. In the center panel of Figure 5, we have overplotted the averaged radial distribution (in red) of the Class I sources by themselves. As can be seen, their distribution is flatter: the ratio between the central density and the density at r ≈1 is about 1.5 in the Class Is versus 3.8 in the full YSO sample. Since Class I objects are likely younger on average than Class II YSOs, this result supports the idea that the YSOs evolve on approximately the shell expansion timescale. However, the bubble expansion times may not be quite so short. Barentsen et al. (2011) find an age difference in the YSOs in IC 1396 of 1 Myr over a distance of 5 pc, and the age differences found within other galactic and between separate extragalactic star-forming regions are more like 3 Myr over a typical distance of 10 pc (Elmegreen 2000). Since we do not know the actual disk fraction of the central clusters, if the fraction of stars with disks is 40% or 30%, the implied central cluster age would be more like 3–4 Myr, bringing the timescale closer to these slower measurements.

We emphasize that the fireworks hypothesis process described here should be placed in the context of much larger superbubble star formation. For example, Patel et al. (1998) describe a model of sequential massive star formation to explain the ∼50 pc radius Cep OB2 bubble, which contains IC 1396 as one of its secondary generations at the periphery. Similarly, Sato et al. (2008) find that the space motions of NGC 281 out of the Galactic Plane are consistent with it being at the edge of an even larger expanding bubble of 130 pc radius (Megeath et al. 2003). It is unclear at present how the two regimes differ in the role of triggering and the mechanisms which operate on these scales. Further expansion of the study to more regions to test our conclusions about RDI and to larger regions to compare with the process of star formation on superbubble scales will be vital to improve and refine our understanding of the progression of star formation in the Galaxy.

The sample of YSOs that we derive from the WISE data comes primarily from the matched WISE bands 1, 2, and 3 catalogs (3.4, 4.6, and 12 μm) that we produce for each field, requiring photometric uncertainty <0.2 mag in all three bands. The first step in extracting young stars is to remove objects from the merged catalog that are most likely to be unresolved external galaxies.

The Spitzer Deep, Wide-Field Survey catalog (Ashby et al. 2009) has been matched to WISE data for the same 8.5 deg² area of Boötes. We developed our criteria by plotting these data, as they will likely be dominated by external galaxies in WISE color space. For comparison, we also extracted WISE photometry from the north and south equatorial poles (decl.  >  +8822 and also decl. < −8822) which we assume will also consist mostly of galaxies.

In Figure 6 we plot a color–color diagram of the objects detected in the north polar field, requiring photometric error <0.2 in WISE bands 1, 2, and 3. For comparison, we overplot tracks showing the location of the galactic SED templates of Assef et al. (2010). Although labeled as E (elliptical), Sbc (spiral), Im (Magellanic irregular), and AGN, these templates are not built to look like specific galaxy types, but are used to model galaxies as composites of different populations. The templates are not orthogonal by design, but are intended to reproduce the typical components of galaxies: old (E), intermediate (Sbc), and young stellar populations (Im), as well as nuclear activity (AGN).

Zoom In Zoom Out Reset image size

Galaxies with elevated star formation activity exhibit increased PAH-feature emission which gives them red colors in the Spitzer 5.8 and 8 μm bands and this property is similarly observed in the WISE [4.6] − [12] color. We thus label objects following all the constraints below as PAH/star-forming galaxies. These equations hold to delimit the region drawn in Figure 7:

$$\begin{displaymath} [3.4] - [4.6] < 0.46\times \left([4.6]-[12] - 1.7\right) \end{displaymath}$$

$$\begin{displaymath} [3.4] - [4.6] > -0.06\times \left([4.6]-[12] - 4.67\right) \end{displaymath}$$

$$\begin{displaymath} [3.4] - [4.6] < -1.0\times \left([4.6]-[12] - 5.1\right) \end{displaymath}$$

$$\begin{displaymath} [3.4] - [4.6] > 0.48\times \left([4.6]-[12] - 4.1\right) \end{displaymath}$$

$$\begin{displaymath} [4.6] > 12 \end{displaymath}$$

$$\begin{displaymath} [4.6] - [12] > 2.3. \end{displaymath}$$

Zoom In Zoom Out Reset image size

Unresolved broad-line AGNs possess mid-IR colors very similar to young stars. However, on average they will be fainter than a typical young star in Galactic regions closer than ∼5 kpc, thus we use a brightness cut in WISE band 2 as the primary discriminant. We consider as candidate AGN those sources following either set of conditions below (see Figure 8):

$$\begin{displaymath} [4.6] > 1.9\times \left([4.6]-[12] + 3.16\right) \end{displaymath}$$

$$\begin{displaymath} [4.6] > -1.4\times \left([4.6]-[12] - 11.93\right) \end{displaymath}$$

$$\begin{displaymath} [4.6] > 13.5 \end{displaymath}$$

or

$$\begin{displaymath} [3.4] > 1.9\times \left([3.4]-[12] + 2.55\right) \end{displaymath}$$

$$\begin{displaymath} [3.4] > 14.0. \end{displaymath}$$

Zoom In Zoom Out Reset image size

Both star-forming galaxies and AGNs are removed from further consideration as YSOs in the following treatment. As discussed in Gutermuth et al. (2009), Koenig et al. (2008), and Rebull et al. (2011), these criteria are imperfect as they misclassify some fraction of the extragalactic objects and also remove some fraction of true YSOs from the catalog, since faint YSOs will overlap the AGN and star-forming galaxy color/magnitude boundaries. We estimate from the Boötes field sample that 4.7 ± 0.7 galaxies deg⁻² will appear in the YSO list as Class II sources, 2.4 ± 0.5 deg⁻² will appear as Class I sources, and 1.8 ± 0.5 will be false transition disk objects. The north and south pole fields result in a slightly lower number of fake YSOs: the combined range of these fields, including Poisson errors, is 1.4–2.9 Class I sources, 2.2–5.4 Class II sources, and 0.9–2.7 transition disk sources per square degree.

These values are sensitive to the relative depths of WISE coverage (which is deeper at the poles than in the ecliptic, due to the orbit of the WISE telescope), any bright extended emission in the image and natural variability in the space density of stars and galaxies. For example, in the study of the wider environment around the Taurus star-forming region by Rebull et al. (2011), 724 of their 1014 YSO candidates were found to be likely contaminant galaxies, foreground or background stars, or other objects. Their sample required signal-to-noise ratio >7 in WISE bands 1, 2, 3, and 4, and resulted in 2.75 total fake YSOs per square degree after comparison with available spectroscopy and existing catalogs. Extrapolating this result to a WISE band 1, 2, and 3 sample in line with this paper, we still find only 4.8 fake YSOs per square degree. Thus, the quoted level of contamination that we present is likely an overestimate in most Galactic star-forming regions.

We next remove two classes of image contaminants from the remaining catalog. As discussed by Gutermuth et al. (2009), the first are resolved shock emission knots—in Spitzer images these are prominent in the 4.5 μm band. The same is likely true of the WISE 4.6 μm band. The second class is resolved structured PAH emission; in WISE data, these are mainly fake detections at 12 μm due to the comparable point source-response-function of WISE at 12 μm and size scale of structure in the copious PAH nebulosity that is prevalent in star-forming regions.

Shock objects are those with the following colors:

$$\begin{displaymath} [3.4] - [4.6] > 1.0 \end{displaymath}$$

and

$$\begin{displaymath} [4.6] - [12] < 2.0. \end{displaymath}$$

Resolved PAH emission objects require either

$$\begin{displaymath} [3.4] - [4.6] < 1.0 \end{displaymath}$$

and

$$\begin{displaymath} [4.6] - [12] > 4.9 \end{displaymath}$$

or

$$\begin{displaymath} [3.4] - [4.6] < 0.25 \end{displaymath}$$

and

$$\begin{displaymath} [4.6] - [12] > 4.75. \end{displaymath}$$

Having removed the previously defined contaminants, we next identify YSOs using the list of objects with good three-band WISE detections, i.e., those possessing photometric uncertainty <0.2 mag in WISE bands 1, 2, and 3. We have tested the following scheme by comparing the various color criteria with the photometry of Taurus region YSOs listed in Rebull et al. (2010) that are also detected by the WISE telescope. See Figure 9 for the following steps. Class I YSOs (candidate protostars) are the reddest objects, and are selected if their colors match:

$$\begin{displaymath} [3.4] - [4.6] > 1.0 \end{displaymath}$$

and

$$\begin{displaymath} [4.6] - [12] > 2.0. \end{displaymath}$$

Zoom In Zoom Out Reset image size

Class II YSOs (candidate T Tauri stars) are slightly less red objects, and are selected with colors:

$$\begin{displaymath} [3.4] - [4.6] - \sigma _1 > 0.25 \end{displaymath}$$

and

$$\begin{displaymath} [4.6] - [12] - \sigma _2 > 1.0, \end{displaymath}$$

where

$$\begin{displaymath} \sigma _1 = \sigma \left([3.4] - [4.6]\right) \end{displaymath}$$

and

$$\begin{displaymath} \sigma _2 = \sigma \left([4.6] - [12]\right), \end{displaymath}$$

where σ(...) indicates a combined error, added in quadrature.

Many objects visible in WISE bands 1 and 2 will lack a reliable band 3 or 4 detection, due to the bright background emission present at these longer wavelengths and the drop in sensitivity. To make up for this, we match the WISE catalog to the 2MASS JHK_S point-source catalog (Skrutskie et al. 2006). We match using a 3'' search radius, removing objects from the 2MASS catalog with photometric quality flags D, E, F, or X in any band. This stage of the scheme requires a measurement of the intrinsic colors of objects, removing line-of-sight extinction to each source. We directly adopt the technique described by Gutermuth et al. (2005) and use the 2MASS JHK_S point-source catalog (cleaned for photometric quality as above) to generate an extinction map for each region with 36'' pixels. For each source in the merged WISE+2MASS source list, we take the median extinction value in a 3 × 3 pixel grid around the nearest pixel in the extinction map. Each source's WISE+2MASS photometry is then dereddened using the extinction law presented in Flaherty et al. (2007).

Using these dereddened colors, we find additional Class II objects requiring the following dereddened colors from the remaining unclassified objects (see Figure 10):

$$\begin{displaymath} \left[[3.4] - [4.6]\right]_0 - \sigma _1 > 0.101 \end{displaymath}$$

and

$$\begin{displaymath} \left[K_S - [3.4]\right]_0 - \sigma _3 > 0.0 \end{displaymath}$$

and

$$\begin{eqnarray*} &&\left[K_S - [3.4]\right]_0 - \sigma _3 > -2.85714\times \left(\left[[3.4] - [4.6]\right]_0 - 0.101\right)\\ &&\quad +\, 1.05 \end{eqnarray*}$$

and

$$\begin{displaymath} [3.4]_0 < 13.8, \end{displaymath}$$

where

$$\begin{displaymath} \sigma _1 = \sigma \left([3.4] - [4.6]\right) \end{displaymath}$$

$$\begin{displaymath} \sigma _3 = \sigma \left(K_S - [3.4]\right), \end{displaymath}$$

where σ(...) indicates a combined error, added in quadrature.

Zoom In Zoom Out Reset image size

Among these new Class II objects, we classify the reddest objects as Class I objects with

$$\begin{eqnarray*} &&\left[K_S - [3.4]\right]_0 - \sigma _3 > -2.85714\times \left(\left[[3.4] - [4.6]\right]_0 - 0.401\right)\nonumber\\ &&\quad +\, 1.9. \end{eqnarray*}$$

In the same way that Gutermuth et al. (2008) used Spitzer MIPS 24 μm data to identify evolved so-called transition disks and confirm that Class I sources SEDs continue to rise at long wavelengths.

We first identify as "transition disks" (disks with photospheric colors between 3.4 and 12 μm but an excess at 22 μm) objects with photometric error <0.2 in WISE bands 1, 2, and 4 and

$$\begin{displaymath} [4.6] - [22] > 2.5 \end{displaymath}$$

and

$$\begin{displaymath} [3.4] < 14. \end{displaymath}$$

Figure 11 shows the location of the Class I criteria specified above and the region in which we find candidate transition disk objects.

Zoom In Zoom Out Reset image size

We check that previously classified Class I stars that possess photometric error <0.2 in WISE band 4 have rising SEDs at 22 μm. If they do not, they are considered reddened Class II sources and re-classified as such if

$$\begin{displaymath} [4.6] - [22] < 4.0. \end{displaymath}$$

We also eliminate Class II stars that possess excessively blue colors. While these may possibly be stars with disks, we conservatively reject them as T Tauri stars and place them back in the unclassified pool if

$$\begin{displaymath} [3.4] - [12] < -1.7 \times \left([12]-[22]\right) + 4.3. \end{displaymath}$$