YOUNG STELLAR OBJECTS IN THE GOULD BELT

The Spitzer c2d survey (PI: N. J. Evans) imaged seven large, nearby molecular clouds in the GB and approximately 100 isolated dense molecular cores. Evans et al. (2003) review both the primary science objectives of c2d and its basic observation plan, and Evans et al. (2009) summarize the results in the large clouds. The Spitzer GB survey (PI: L. E. Allen) is a follow-up program that imaged an additional 11 molecular clouds in the GB, completing most of the remaining clouds in the GB except for the Taurus and Orion molecular clouds (see Section 1). Both surveys imaged their respective targets at 3.6–8.0 μmwith the Spitzer Infrared Array Camera (IRAC; Fazio et al. 2004), and at 24–160 μm images with the Multiband Imaging Photometer (MIPS; Rieke et al. 2004). A data pipeline developed by c2d was applied to the observations obtained by both programs. This pipeline includes data reduction, source extraction, and band-merging the results from each wavelength to produce final images and source catalogs, and has been described in detail elsewhere (Harvey et al. 2006; Evans et al. 2007).

For each cloud targeted by c2d or GB, Table 1 lists the parent survey, the assumed distance and reference for the distance, the total area mapped, the total cloud gas mass, and references to in-depth studies of each cloud based on the c2d or GB data, when available. The Spitzer Astronomical Observation Requests (AORs) for each cloud can be found in these references to individual cloud studies; Appendix A lists the AORs for the clouds without published data references. The isolated cores surveyed by c2d are not considered here except for a few cases where the cores are actually parts of larger cloud complexes and included in the results for those clouds (namely, Cepheus and Ophiuchus North; Kirk et al. 2009; Hatchell et al. 2012). Note that the distances to these clouds are not all well determined, and some cloud distances are still under debate. One such example is the uncertainty over the distance to Aquila, with two distances (260 and 429 pc) commonly adopted in the literature (see, e.g., Harvey et al. 2006; Gutermuth et al. 2008; Dzib et al. 2011; Maury et al. 2011, for details). For this reason we do not list formal distance uncertainties and instead refer the reader to the references given in Table 1 for full discussions of the distances to each cloud.

The total cloud areas mapped and total cloud gas masses contained within these areas are taken from Heiderman et al. (2010), with the masses determined from extinction maps (see Heiderman et al. for details). For most clouds, the mapped areas were chosen to map the entire cloud above a minimum visual extinction in the range 2 < A_V < 3, with the exact value for each cloud depending on the total area that could be mapped within the allocated time. For Serpens and Aquila, only the cloud areas above A_V = 6 were mapped due to time constraints and difficulty in separating the cloud from the general extinction close to the Galactic midplane.

The data reduction pipeline used by both the c2d and GB surveys creates band-merged source catalogs where the sources extracted at each Spitzer wavelength are matched with each other and with the 2MASS catalog, creating a final catalog of 1.25–70 μm photometry for each cloud. Candidate YSOs are identified using a method developed by Harvey et al. (2007a) and summarized in all of the publications focused on individual clouds that are listed in Table 1. In this method, the Spitzer SWIRE Legacy survey of the ELAIS N1 extragalactic field (Lonsdale et al. 2003) is processed to simulate the same sensitivity and extinction distribution of the c2d and GB clouds, and then used to locate background galaxies observed through the target clouds in various infrared color–magnitude diagrams. Using these results to define the expected colors and magnitudes of background galaxies, each source extracted in the cloud maps with infrared colors that can not be explained as an extincted background star is assigned an unnormalized "probability" of being a galaxy or YSO based on its position in each color–magnitude diagram, its morphology in the two shortest Spitzer IRAC bands (3.6 and 4.5 μm), and its 24 and 70 μm flux densities. Since it is unnormalized, this quantity is not a true probability but a relative parameter expressing the likelihood of being a background galaxy (higher values indicate higher likelihood of being galaxies; see Harvey et al. 2007a, for details). The final boundary between candidate YSO and candidate galaxy in this quantity is set to provide a nearly complete elimination of SWIRE sources, likely at the cost of excluding YSOs with overlapping colors and magnitudes (e.g., Hsieh & Lai 2013).

One potential drawback to the c2d+GB method of identifying YSOs is that it requires detections in all four IRAC bands (3.6–8.0 μm) as well as the first MIPS band (24 μm). As the 24 μm MIPS band has lower resolution and sensitivity, some of the faintest and/or most clustered YSOs may be missed (e.g., Gutermuth et al. 2009, R.A. Gutermuth et al. 2015, in preparation). We note that several alternative classification methods for separating background galaxies and YSOs have been presented in the literature (e.g., Gutermuth et al. 2009; Rebull et al. 2010; Kryukova et al. 2012; Hsieh & Lai 2013). While they are all based on a similar philosophy of separating sources in various colors and magnitudes, they differ in their detailed implementations of this philosophy. A full comparison of these methods and ultimate consensus on best practices for extracting YSOs from infrared surveys of star-forming clouds is left for future work.

Applying these selection criteria results in the identification of 3239 candidate YSOs in the 18 c2d and GB clouds. We note here that we only consider sources that are located in the overlap regions imaged by both the IRAC and MIPS instruments, as these are the sources that can be most reliably classified using the method summarized above. For technical reasons, the MIPS observations often covered much larger areas (see the references listed in Table 1 for more details on the observations of each individual cloud). Several of the individual cloud studies listed in Table 1 also present candidate YSO lists in the MIPS-only regions using alternative identification and classification techniques, but such sources are not considered here. We visually inspected each of the 3239 candidate YSOs to remove resolved galaxies and image artifacts that are sometimes misidentified (see Evans et al. 2009, for details). Additionally, a few known YSOs that were saturated or undetected at one or more wavelengths were missing from the initial list and were added by hand. In the end, we were left with a final list of 2966 YSOs, and since all have passed visual inspection, we follow Evans et al. (2009) and drop the prefix "candidate" at this point.

Finally, we note that initial results from Herschel surveys of the c2d and GB clouds suggest that a modest number of the most deeply embedded protostars are undetected by Spitzer, increasing the number of Class 0 protostars¹⁹ by 15%–30% (e.g., Harvey et al. 2013; Stutz et al. 2013; Sadavoy et al. 2014). Revisions to the YSO lists presented here will thus be necessary once complete Herschel results are available for each cloud.

Our final sample of YSOs is likely contaminated by background giant stars with infrared excesses. Oliveira et al. (2009) performed optical spectroscopy toward 78 YSOs in Serpens and identified 20 (26%) as background AGB stars. While the overall contamination rate is relatively low, Figure 1 shows that the contamination fraction strongly depends on the infrared spectral index α.²⁰ Starting from positive values, as α decreases, the contamination fraction remains close to zero until the Class II/III boundary is reached at α = −1.6 (see Section 3.2 for a full definition and discussion of YSO classes), and then steeply rises to >50% for more negative values of α. On average, 62% of the Class III sources examined by Oliveira et al. are actually background AGB stars, whereas only 5% of the Class II sources they examined were found to be contaminants.

Very high rates of AGB star contamination among Class III sources are also found by Romero et al. (2012), also with optical spectra. They found Class III contamination rates of 100% based on observations of 30 objects in Lupus V and VI. While the fact that they only observed objects that met initial selection criteria for transition disks may bias their results, they observed 45% and 39% of the total Lupus V and VI Class III populations, respectively; thus, the true contamination fractions are at least ∼40% in each cloud.

The true contamination fraction is likely a strong function of galactic latitude, and the specific clouds considered above (Serpens, Lupus V, and Lupus VI) all lie within 2°–10° of the galactic plane. Indeed, Romero et al. (2012) found a lower contamination rate of 25% for Ophiuchus (b = 18°), consistent with the Ophiuchus contamination rate of 20% measured by Cieza et al. (2010) using similar methods. To further investigate the contamination rate in our full sample of YSOs, Figure 2 plots a Spitzer color–color diagram consisting of [3.6]–[24] versus [3.6]–[4.5], using the observed photometry. Inspection of Figure 2 shows that the Class III YSO population can be divided into a population of objects clustered at lower (bluer) values of both colors, and a second, more distributed population extending to higher (redder) colors. Using optical spectra of candidate transition disks, Cieza et al. (2010, 2012) and Romero et al. (2012) found that background AGB stars are predominately found at low values of [3.6]–[24], suggesting that the clustered population in Figure 2 may be dominated by AGB stars whereas the more distributed population may be dominated by YSOs. Indeed, the above studies found that all of their objects with [3.6]–[24] ≤ 1.5 were in fact AGB stars, with decreasing contamination rates at higher colors. We find that 298 of the objects in our final YSO sample have [3.6]–[24] ≤ 1.5, all but three of which are classified as Class III objects. This represents 25% of the Class III YSO sample and 10% of the total YSO sample. We consider these objects as likely AGB stars and mark them as such throughout the remainder of this paper. However, we do not remove them from the sample since the studies cited above only targeted a subset of objects that were identified as transition disk candidates and thus there is no guarantee that their finding of 100% contamination at the lowest values of [3.6]–[24] applies to our full sample. Confirmation of the status of these objects requires either spectroscopic or proper motion studies. Such work should be considered as a high priority for future investigations seeking to refine the YSO sample presented here.

Zoom In Zoom Out Reset image size

The 25% Class III contamination by AGB stars found above is almost certainly a lower limit to the full number, since some contaminating AGB stars are found to have [3.6]–[24] > 1.5 (Cieza et al. 2010, 2012; Romero et al. 2012). However, these studies target samples that are too small and too biased to determine contamination rates as a function of [3.6]–[24] color. To obtain an upper limit to the possible contamination fraction, Figure 3 compares the distributions of α for our YSOs with those identified in young stellar clusters by Gutermuth et al. (2009), using a selection method that is also based on infrared colors and magnitudes but is different in its details (see Gutermuth et al. for details). The c2d+GB sample contains a much higher fraction of Class III objects. While a full comparison of different methods for identifying YSOs is beyond the scope of this work, under the assumptions that the Gutermuth et al. method perfectly separates out background AGB stars and our regions contain the same fractions of sources in each Class, 87% of the Class III objects in the c2d+GB sample would need to be removed to match the Gutermuth et al. (2009) Class III fraction. We note here that the Gutermuth et al. survey targeted smaller regions centered on dense clusters whereas the c2d and GB surveys mapped much larger areas, including those more distant from cluster centers. As the Class II and III YSO populations are generally more distributed than earlier Classes, the intrinsic fraction of Class III sources is likely higher in the c2d and GB surveys than in Gutermuth et al. (2009). For this reason, our above assumption that the surveys contain the same fractions of sources in each Class results in an upper limit to the Class III contamination fraction.

Zoom In Zoom Out Reset image size

Synthesizing all of the above information, our sample of Class III YSOs suffers from between 25% and 90% contamination by AGB stars. These values are lower and upper limits, respectively, for the reasons listed above, and the most likely value is somewhere in between. Tighter constraints on the true contamination rate will require large optical spectroscopic campaigns and/or proper motion studies coupled with revisions to the selection criteria used to identify YSOs in infrared surveys. These contamination rates must be kept in mind when considering the results presented here and in other studies based on the c2d+GB YSO list (such as those discussed in Section 1).

Following the method first summarized by Evans et al. (2009) and updated by Dunham et al. (2013), we compile as complete an SED as possible for each of the 2966 YSOs. We refer to Dunham et al. for a full description of this process, but give a summary here. To construct the SEDs, we included the following.

• 1.  
  Selected ground-based optical and infrared photometry as compiled by the authors of the individual cloud papers listed in Table 1.
• 2.  
  Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) 12 and 22 μm photometry from the all-sky catalog.²¹
• 3.  
  Spitzer 160 μm photometry for sources that are detected but not located in saturated or confused regions. Since the c2d pipeline does not extract 160 μm sources, we measured flux densities using aperture photometry (using the IDL procedure aper) and standard aperture corrections as given by the MIPS Instrument Handbook.²²
• 4.  
  Selected 350 μm photometry from a targeted survey of protostars with SHARC-II²³ (Wu et al. 2007, Suresh et al. 2015).
• 5.  
  Selected 450 and 850 μm photometry from the SCUBA²⁴ Legacy Catalog (Di Francesco et al. 2008).
• 6.  
  Other submillimeter and millimeter photometry from the literature, where available from large surveys of nearby molecular clouds (Young et al. 2005, 2006; Enoch et al. 2006, 2007; Belloche et al. 2011a, 2011b; Maury et al. 2011).

Regarding the last three items, we have access to complete submillimeter or millimeter surveys for only 6 out of the 18 clouds (Chamaeleon I, Chamaeleon II, Chamaeleon III, Ophiuchus, Perseus, and Serpens), plus a partial survey of Aquila and piecemeal coverage of other clouds from the SCUBA Legacy Catalog (Di Francesco et al. 2008) and our own 350 μm observations (Wu et al. 2007, A. Suresh et al. 2015, submitted). We consider a Spitzer source to match a submillimeter or millimeter detection if it is located within one beam of the detection, using the appropriate beam size for each survey. If multiple Spitzer sources match a single submillimeter or millimeter detection, we do not attempt to split the flux density and simply match it to the closest Spitzer source unless the SED clearly indicates a better match is obtained with a different Spitzer source.

Finally, we correct all photometry for foreground extinction following the procedure outlined by Evans et al. (2009). We only wish to correct for the foregound extinction from both the clouds and the material between us and the clouds, and not the local extinction from dense cores surrounding protostars (this extincted emission is reprocessed to longer wavelengths and thus included in our compiled SEDs). Thus, we adopt extinction values as follows.

• 1.  
  Whenever possible, we adopt extinction values from the literature for Class II and III YSOs (classified via infrared spectral index; see Section 3.2), as determined from optical or near-infrared spectroscopic studies of each cloud that provide reliable spectral type (and thus extinction) measurements. We refer to the data references listed in Table 1 for more details on these studies.
• 2.  
  For all Class II and III YSOs with no published extinction values, we de-redden to the intrinsic near-infrared colors of an assumed spectral type of K7. This particular spectral type is chosen because it is found to be representative of the YSOs in the c2d clouds (Oliveira et al. 2009, 2010, see also Evans et al. 2009 for details).
• 3.  
  We calculate the mean extinction toward all Class II YSOs in each cloud, and then adopt that mean value for each of the Class 0 + I and Flat-spectrum YSOs in that cloud.

With extinctions assigned as above, we de-redden the SED of each YSO using the Weingartner & Draine (2001) extinction law for R_V = 5.5. Whereas Evans et al. (2009) adopted a version of the Weingartner & Draine (2001) R_V = 5.5 extinction law that only included dust absorption, here we update the procedure to correct for both dust absorption and scattering. The choice of this particular extinction law is motivated by results showing that it is appropriate for the dense molecular clouds in which stars form (e.g., Chapman et al. 2009). While we do caution that our approach only provides approximate extinction corrections, especially for the Class 0 + I and Flat-spectrum sources where the corrections only account for the mean cloud extinction, it is nevertheless the best that can currently be done and is significantly more reliable than ignoring the effects of extinction altogether.

For each of the 2966 YSOs, Table 2 lists a running index, the cloud in which the YSO is found, the Spitzer source name (which also gives the source coordinates), the visual extinction used for extinction corrections, and the infrared spectral index (α), the bolometric temperature (T_bol), and the bolometric luminosity (L_bol). These last three quantities are calculated using both the observed and extinction corrected photometry, and are denoted with primes in the latter case (α', ${T}_{\mathrm{bol}}^{\prime }$, and ${L}_{\mathrm{bol}}^{\prime }$). The second-to-last column of Table 2 indicates whether or not each object has [3.6]–[24] ≤ 1.5 and is thus likely a background AGB star, as described above in Section 2.3.

Table 2.  Properties of the c2d+GB YSOs

		Spitzer		Observed			Extinction Corrected
		Source Name	AV		Tbol	Lbol		${T}_{\mathrm{bol}}^{\prime }$	${L}_{\mathrm{bol}}^{\prime }$	Likely
Index	Cloud	(SSTc2d or SSTgb +)	(mag)	α	(K)	(L⊙)	α'	(K)	(L⊙)	AGB?	Core?
1	Aquila	J180144.8–044941	8.2	−2.16	1700	3.0	−2.52	2000	8.2	N	...
2	Aquila	J180145.8–045902	8.0	−2.14	1700	1.4	−2.56	2000	3.7	N	...
3	Aquila	J180154.2–043753	12.4	−2.57	2100	0.52	−3.22	2600	4.7	Y	...
4	Aquila	J180229.2–051552	8.3	−2.03	1700	3.6	−2.38	2000	9.8	N	...
5	Aquila	J180302.0–044233	12.4	−2.45	2100	0.11	−3.15	2500	0.92	N	...
6	Aquila	J180308.0–043523	8.0	−2.18	1700	0.24	−2.52	2000	0.64	Y	...
7	Aquila	J180313.8–042845	6.9	−2.22	1800	0.72	−2.53	2000	1.7	Y	...
8	Aquila	J180423.9–042248	9.0	−2.15	1700	2.2	−2.44	2000	6.3	Y	...
9	Aquila	J180435.7–044700	8.1	−2.26	1700	1.9	−2.65	2000	5.3	Y	...
10	Aquila	J180444.5–043706	12.4	0.53	540	0.055	−0.10	1000	0.11	N	...

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The final column of Table 2 indicates whether or not each object is associated with a dense core as traced by submillimeter or millimeter continuum emission at λ ≥ 350 μm. As a lack of such an association does not necessarily imply that no core is present (see below), to avoid any confusion in the interpretation of this column we list either y or no entry rather than y or n. Dunham et al. (2013) used such associations to define the sample of protostars in the c2d+GB clouds, with the rationale behind this definition being that detections in existing submillimeter and millimeter surveys of star-forming regions trace dense cores but not circumstellar disks due to the relatively low spatial resolutions and mass sensitivities of these surveys (see Dunham et al. 2013, for details). By requiring a dense core but making no assumptions about the infrared colors of protostars, this definition ensures a reliable sample of protostars that includes those viewed edge-on down outflow cavities, but comes at the cost of incompleteness to those protostars that are embedded in low-mass cores or located in clouds with incomplete or no available submillimeter or millimeter continuum surveys (see Dunham et al. 2014, for details). Associations with dense cores and the resulting sample of protostars must be revisited once the James Clerk Maxwell telescope (JCMT) SCUBA-2 survey of the GB is complete, the first results of which are just starting to become available (e.g., Pattle et al. 2015).

The infrared spectral index (α, α') is defined as the slope of the infrared SED in log(λ S_λ) versus log(λ), where λ is the wavelength and S_λ is the flux density at that wavelength:

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{d\;\mathrm{log}(\lambda {S}_{\lambda })}{d\;\mathrm{log}(\lambda )}.\end{eqnarray} \tag{ 1 }$$

In practice, we calculate α with a linear least squares fit to all available 2MASS and Spitzer photometry between 2 and 24 μm. The bolometric temperature of a source is defined to be the temperature of a blackbody with the same flux-weighted mean frequency (Myers & Ladd 1993). It is in essence a protostellar equivalent of stellar effective temperature; it starts at very low values (∼10 K) for cold, starless cores and eventually increases to the effective temperature of the central (proto)star once the core and disk have fully dissipated. L_bol, ${L}_{\mathrm{bol}}^{\prime }$, T_bol, and ${T}_{\mathrm{bol}}^{\prime }$ are calculated by integrating over the full SEDs:

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}=4\pi {d}^{2}{\displaystyle \int }_{0}^{\infty }{S}_{\nu }d\nu ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{bol}}=1.25\times {10}^{-11}\;\displaystyle \frac{{\displaystyle \int }_{0}^{\infty }\nu {S}_{\nu }d\nu }{{\displaystyle \int }_{0}^{\infty }{S}_{\nu }d\nu }\quad {\rm{K}}.\end{eqnarray} \tag{ 3 }$$

In practice, we calculate these integrals using the trapezoid rule to integrate over the finitely sampled SEDs.

Figure 4 plots ${L}_{\mathrm{bol}}^{\prime }$ versus ${T}_{\mathrm{bol}}^{\prime }$ (a "BLT" diagram; Myers & Ladd 1993) for all 2966 YSOs identified in the c2d and GB clouds. The luminosities of Class 0 and Class I sources are distributed over greater than 4 orders of magnitude and extend to much lower luminosites than predicted by simple evolutionary models. A full discussion of this "protostellar luminosity problem" can be found in Offner & McKee (2011), Dunham et al. (2013, 2014), but to summarize, mass-dependent and/or time variable mass accretion is required to explain the observed luminosities of protostars. Class II and Class III objects exhibit both a very narrow range of bolometric temperatures and a bifurcation in ${T}_{\mathrm{bol}}^{\prime }$ into two groups. Both effects are artifacts introduced by our extinction corrections and are discussed in Appendix C. To summarize, the first effect is because the bolometric temperature of a Class II or III source depends more on the spectral type of the object than the amount and characteristics of circumstellar material, and for most objects a constant spectral type of K7 is assumed (see above). The second effect (bifurcation) is a consequence of optical data only being available for a subset of our YSOs.

Zoom In Zoom Out Reset image size

Tables 3 and 4 list, for each YSO, the same running index as in Table 2, followed by flux density and flux density uncertainty pairs for the 2MASS (1.25, 1.65, and 2.17 μm) and Spitzer (3.6, 4.5, 5.8, 8.0, 24, and 70 μm) wavelengths. All flux densities and flux density uncertainties are listed in mJy and rounded to two significant digits. Table 3 lists the observed values while Table 4 lists the extinction corrected values. The observed and extinction corrected flux densities that we have compiled at all other wavelengths in order to construct SEDS, as described above, are tabulated in Appendix D.

Over the past few decades a standard model for low-mass star formation has emerged. This model, developed by Adams et al. (1987) and summarized by Shu et al. (1987), resulted from the merger of an empirical classification system based on the infrared spectral indices of YSOs (Lada & Wilking 1984) with theoretical studies of the stages of the collapse of a dense, rotating core (Shu 1977; Terebey et al. 1984). Recent reviews on classification of young stars and the associations between observed SED classes and physical evolutionary stages are given by Evans et al. (2009) and Dunham et al. (2014).

Including the revisions by Greene et al. (1994), these empirical Classes are as follows:

• 1.  
  Class 0 + I: α' ≥ 0.3;
• 2.  
  Flat-spectrum: −0.3 ≤ α' < 0.3;
• 3.  
  Class II: −1.6 ≤ α' < −0.3;
• 4.  
  Class III: α' < −1.6.

We denote the infrared spectral index as α' to emphasize that we only consider extinction corrected values in this study, and we refer to the first SED class as Class 0 + I instead of the more commonly used Class I to emphasize that Class 0 and I objects are both included in this definition (see below).

Table 5 lists, for each cloud, the numbers of YSOs in each Class using the extinction corrected infrared spectral indices, both with and without the likely background AGB stars removed. Summed over all clouds, out of the 2966 total YSOs, 326 (11%) are classified as Class 0 + I, 210 (7%) are classified as Flat-spectrum, 1248 (42%) are classified as Class II, and 1182 (40%) are classified as Class III. If we first remove the 298 sources identified as likely background AGB stars, out of the 2668 remaining YSOs, 326 (12%) are classified as Class 0 + I, 210 (8%) are classified as Flat-spectrum, 1245 (47%) are classified as Class II, and 887 (33%) are classified as Class III. We emphasize that there are large variations between clouds in the relative numbers in each Class. To reinforce this point visually, Figures 5 and 6 plot, for each cloud, four-color bars where each color corresponds to a Class and the fractional length of each color represents the fraction of YSOs with the corresponding Class. Figure 5 displays the results including all 2966 YSOs whereas Figure 6 displays the results after removing the 298 likely background AGB stars. Finally, we also emphasize that our method of YSO identification is only capable of identifying Class III sources with detectable infrared excesses and is thus incomplete to the full populations of pre-main sequence stars in these clouds.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As noted above, we refer to the first Class as Class 0 + I, but the original definition of this Class was simply Class I (Lada & Wilking 1984; Greene et al. 1994). Class 0 objects were later added as an earlier Class for protostars too deeply embedded to detect in the infrared but inferred through other means (i.e., outflow presence; Andre et al. 1993). They are defined as objects with a ratio of submillimeter to bolometric luminosity (L_smm/L_bol) greater than 0.5% (Andre et al. 1993), where L_smm is calculated for λ ≥ 350 μm. While they are defined observationally, the intended corresponding physical definition of Class 0 objects are protostars with at least 50% of their total mass still residing in the surrounding, infalling core (Andre et al. 1993). An alternative observational definition is given in terms of the bolometric temperature as objects with T_bol < 70 K (Myers & Ladd 1993; Chen et al. 1995). While L_smm/L_bol is generally a better discriminator between Class 0 and Class I objects than T_bol, and is less sensitive to source geometry (e.g., Young & Evans 2005; Dunham et al. 2010, 2014), it cannot be calculated without 350 μm photometry. As we do not have 350 μm photometry available for many of our YSOs, we only consider T_bol in this study.

With the sensitivity of Spitzer, Class 0 sources are easily detected in the infrared (e.g., Noriega-Crespo et al. 2004; Young et al. 2004; Dunham et al. 2006, 2008; Jørgensen et al. 2006; Tobin et al. 2010). Figure 7 plots α' versus ${T}_{\mathrm{bol}}^{\prime }$ for all 2966 YSOs. In general, the two quantities are anti-correlated, although there are some notable exceptions. There are a small number of SED Class 0 + I sources (plotted in red in Figure 7) with sufficiently high ${T}_{\mathrm{bol}}^{\prime }$ values to be classified as Class II or III by ${T}_{\mathrm{bol}}^{\prime }$. None of these objects are associated with dense cores and are likely simply more evolved YSOs viewed through large extinctions that are not fully extinction-corrected (recall that, for SED Class 0 + I sources, only the average cloud extinction is accounted for). There are also a few SED Class II YSOs (plotted in blue in Figure 7) with sufficiently low ${T}_{\mathrm{bol}}^{\prime }$ values to be classified as Class 0 by ${T}_{\mathrm{bol}}^{\prime }$. These are likely true protostars viewed at nearly pole-on inclinations through outflow cavities (e.g., Whitney et al. 2003; Robitaille et al. 2006; Crapsi et al. 2008), although detailed follow-up study is required to rule out chance alignments between more evolved YSOs and dense cores.

Zoom In Zoom Out Reset image size

While there is a general anti-correlation between α' and ${T}_{\mathrm{bol}}^{\prime }$, this trend weakens for objects classified as Class 0 + I by α' (objects with α' ≥ 0.3). With substantial overlap in α' for objects classified as Class 0 or Class I by ${T}_{\mathrm{bol}}^{\prime }$, it is clear that these two Classes overlap in their infrared spectral indices. These results confirm those presented previously by Enoch et al. (2009) and Evans et al. (2009) and indicate that mid-infrared data alone are insufficient to separate Class 0 and Class I objects; instead, full SEDs are required.

Using our compiled SEDs, and restricting ourselves only to those YSOs associated with a dense core (both to restrict ourselves to confirmed protostars and to ensure available submillimeter or millimeter detections for reliable measurements of ${T}_{\mathrm{bol}}^{\prime }$), our ${T}_{\mathrm{bol}}^{\prime }$ measurements result in 65 Class 0 protostars and 129 Class I protostars. Thus, we find that approximately one-third (65/[65 + 129] = 0.335) of protostars are classified as Class 0, in agreement with previous results based on the smaller c2d-only sample (Enoch et al. 2009; Evans et al. 2009). While future work must revisit this ratio once complete, sensitive submillimeter coverage is available for all clouds through the Herschel and JCMT GB surveys, preliminary Herschel results suggest only modest increases to the number of Class 0 objects (generally 15%–30%, Harvey et al. 2013; Stutz et al. 2013; Sadavoy et al. 2014), leading to only small upward revisions in the ratio of Class 0 to Class I protostars.

Figure 8 presents an extinction corrected infrared color–color diagram of all YSOs identified in this work constructed from the four IRAC wavelengths: [36] – [45] plotted versus [58] – [80]. The different Classes of YSOs (classified via α' ) are generally well-separated in this diagram, as first noted by Allen et al. (2004). In cases where photometry is only available at 3.6–8.0 μm and thus the infrared spectral index cannot be reliably calculated, this diagram can be used to approximately classify YSOs. Indeed, Allen et al. (2004) used the following two conditions to define a box where Class II sources are located:

• 1.  
  0.4 ≤ [5.8] – [8.0] ≤ 1.1, and
• 2.  
  0.0 ≤ [3.6] – [4.5] ≤ 0.8.

The original definition of this "Class II box" by Allen et al. (2004) used observed photometry. Inspection of Figure 8 reveals that minor revisions are required when using extinction corrected photometry. While we acknowledge that polygons with greater than four sides could possibly give better results, for simplicity we only consider rectangular regions. To define a new Class II box, we allow all four boundaries to vary and select those which simultaneously maximize the probability that a Class II YSO is located within the box (defined as the number of Class II YSOs located within the box divided by the total number of Class II YSOs) and the probability that a YSO located within the box is classified as Class II (defined as the number of Class II YSOs located within the box divided by the total number of YSOs located within the box). In practice, we implement these definitions by selecting the boundaries which maximize the sum of these two probabilities. Our proposed new Class II box, plotted in panel (b) of Figure 8, has the following boundaries:

• 1.  
  0.45 ≤ [58] – [80] ≤ 1.55, and
• 2.  
  0.00 ≤ [36] – [45] ≤ 0.70.

We note that the shape of this revised box is qualitatively similar to that proposed by Gutermuth (2005), who followed a similar procedure. With these revised boundaries, the fraction of Class II YSOs located within the box increases from 84.5% to 89.9%, and the fraction of YSOs in the box that are Class II increases from 78.8% to 82.4%.

Zoom In Zoom Out Reset image size

Figure 9 presents a second extinction corrected infrared color–color diagram of all YSOs which uses both the IRAC and MIPS observations by plotting [36] – [58] versus [80] – [24]' . Again the different Classes of YSOs (classified via α' ) are generally well-separated, as first noted by Muzerolle et al. (2004). Indeed, Muzerolle et al. defined three regions in this color–color diagram to separate Class III/stellar, Class II, and Class 0 + I YSOs; these regions are marked in panel (a) of Figure 9 and are defined in Table 6. As above, Muzerolle et al. used observed photometry, and we find that minor revisions are required when using extinction corrected photometry. Furthermore, the original boundaries defined by Muzerolle et al. did not include a region for Flat-spectrum objects. Thus, we follow the same procedure as above, again only considering rectangular regions, and maximize the sum of the eight probabilities (two for each Class), subject to the following constraints: (1) the minimum x and y Class III boundaries be held fixed at the Muzerolle et al. values, (2) there are no maximum x and y Class 0 + I boundaries, and (3) the regions can not overlap. Our revised boundaries are plotted in panel (b) of Figure 9 and listed in Table 6. Table 6 also lists each of the two probabilities considered for both the original and revised boundaries.

Zoom In Zoom Out Reset image size

Ultimately, no region in either color–color diagram will classify YSOs with a perfect one-to-one correspondence to the infrared spectral index, but in cases where insufficient photometry at 2 and/or 24 μm is available, approximate YSO classification is possible via infrared color–color diagrams.

The standard method for determining the time spent in each YSO Class is to calculate the ratio of the number of sources in that Class to the number in a reference Class, and then multiply by the duration of the reference Class (e.g., Wilking et al. 1989; Evans et al. 2009; Dunham et al. 2014). Implicit in this method are the assumptions that star formation is continuous over at least the duration of the reference class and that time is the only variable. Regarding the first assumption, we strive to average out cloud-to-cloud variations in the recent star formation history by averaging over many clouds. Regarding the second, our results are best thought of as an average over all other relevant variables (e.g., environment, stellar mass; see Evans et al. 2009). In this section, we investigate the effects of different assumptions for the reference class and its duration on the timescales of Class 0 + I and Flat-spectrum YSOs.

We first follow the method of Evans et al. (2009) and Dunham et al. (2014), who adopt Class II as the reference Class and assume a Class II duration of 2 Myr. This duration is taken from analyses of the fraction of cluster members with disks versus cluster age for various young clusters, which generally show a disk disperal time of a few Myr (Wyatt 2008; Ribas et al. 2014, 2015). As discussed extensively by Evans et al. (2009), this time is best considered as a disk half-life rather than an absolute duration. The first row of Table 7 presents durations for the Class 0 + I and Flat-spectrum Classes calculated following this method, using the extinction corrected numbers: 0.52 and 0.34 Myr for Class 0 + I and Flat-spectrum objects, respectively.

Table 7.  Timescales for Young Stellar Objects

	τrefa	τ0 + Ib	τFc	τ0d	τIe
Reference Class	(Myr)	(Myr)	(Myr)	(Myr)	(Myr)
All Class II	2.0	0.52	0.34	0.17	0.35
All Class II	3.0	0.78	0.50	0.26	0.52
All Class II + All Class III	3.0	0.40	0.26	0.13	0.27
All Class II + 75% of Class III	3.0	0.46	0.30	0.15	0.31
All Class II + 10% of Class III	3.0	0.72	0.46	0.24	0.48

Notes.

^aAssumed duration of the reference Class. ^bCalculated duration of Class 0 + I sources (${\tau }_{\mathrm{ref}}\frac{{N}_{0+{\rm{I}}}}{{N}_{\mathrm{ref}}}$). ^cCalculated duration of Flat-spectrum sources (${\tau }_{\mathrm{ref}}\frac{{N}_{{\rm{F}}}}{{N}_{\mathrm{ref}}}$). ^dCalculated duration of Class 0 sources ($0.335\times {\tau }_{0+{\rm{I}}}$); see Section 4.2 for details. ^eCalculated duration of Class I sources ($0.665\times {\tau }_{0+{\rm{I}}}$); see Section 4.2 for details.

Download table as:  ASCIITypeset image

Table 9.  Observed 0.36–0.96 μm Flux Densities, in mJy, of the c2d+GB YSOs

	U	Δ U	B	ΔB	V	Δ V	R	ΔR	Hα12	Δ Hα12	I	Δ I	m914	${\rm{\Delta }}m914$	z	Δ z
Index	0.36 μm	0.36 μm	0.44 μm	0.44 μm	0.55 μm	0.55 μm	0.64 μm	0.64 μm	0.665 μm	0.665 μm	0.79 μm	0.79 μm	0.915 μm	0.915 μm	0.96 μm	0.96 μm
1	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
4	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
5	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
6	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
7	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
8	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
9	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
10	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Although the ages of young stars are highly uncertain, recent studies suggest a revised half-life for disk survival of 2.5–4 Myr on both observational (e.g., Kraus et al. 2012; Bell et al. 2013; Soderblom et al. 2014) and theoretical (e.g., Alexander & Armitage 2009) grounds. Very recent results by Ribas et al. (2014, 2015), using spectroscopically confirmed members of various young clusters, show that this duration depends on the wavelength where the infrared excess indicative of the presence of a disk first appears, ranging from 1.9 Myr for λ < 5 μm to 3.0 Myr for λ ≤ 24 μm (note that these times have been converted from the 1/e times presented by Ribas et al. to half-life times). As our Spitzer data are sensitive to infrared excesses out to 24 μm, and our Spitzerc2d+GB data have comparable mid-infrared sensitivities to the Ribas et al. Spitzer and WISE data used to determine disk lifetimes, their derived lifetimes suggest that 3 Myr is a more appropriate choice for the duration of the reference Class in this work. The second row of Table 7 lists the revised durations: 0.78 and 0.50 Myr for Class 0 + I and Flat-spectrum objects, respectively.

Table 10.  Extinction Corrected 0.36–0.96 μm Flux Densities, in mJy, of the c2d+GB YSOs

	U	Δ U	B	ΔB	V	Δ V	R	ΔR	Hα12	Δ Hα12	I	Δ I	m914	${\rm{\Delta }}m914$	z	Δ z
Index	0.36 μm	0.36 μm	0.44 μm	0.44 μm	0.55 μm	0.55 μm	0.64 μm	0.64 μm	0.665 μm	0.665 μm	0.79 μm	0.79 μm	0.915 μm	0.915 μm	0.96 μm	0.96 μm
1	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
4	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
5	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
6	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
7	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
8	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
9	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
10	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

While it has become somewhat standard to adopt Class II as the reference Class, this choice may result in a logical inconsistency. The duration of either 2 or 3 Myr based on disk fraction versus cluster age is derived by determining the ratio of cluster members with infrared excesses to all cluster members (in other words, the ratio of disk sources to disk + no disk sources). By adopting Class II as the reference Class, we omit all of the Class III YSOs that still have disks, even though these objects are included when other studies (such as those by Ribas et al. ) derive disk lifetimes. Since infrared colors indicative of dust are the fundamental ingredient in our identifcation of YSOs (see Section 2.2), by definition the Class III YSOs identified in this work have disks (if they are all in fact YSOs; the issue of Class III contamination will be addressed in the following paragraph). Since the recent work by Ribas et al. (2014, 2015) is based on mid-infrared data from Spitzer and WISE with comparable sensitivity to mid-infrared excesses as our c2d+GB data, and since we use disk lifetimes derived by Ribas et al. to set the duration of the reference Class, and further given that all Class II and Class III sources identified here have disks, we thus argue that it is more appropriate to adopt the sum of Class II and Class III as the reference Class than Class II alone. Thus, the third row of Table 7 lists the revised durations assuming a duration of 3 Myr for the sum of Class II and III sources: 0.40 and 0.26 Myr for Class 0 + I and Flat-spectrum objects, respectively.

The one major complication with the above analysis is the fact that our Class III sample suffers from 25% to 90% contamination, as discussed in Section 2.3. On the other hand, contamination is not a problem in the Ribas et al. (2014, 2015) samples since they only consider spectroscopically confirmed cluster members. As a result, the above calculations of Class durations likely overcount the number of objects in the reference Class and thus underestimate the Class 0 + I and Flat-spectrum durations. To account for this range of contamination rates, the last two columns of Table 7 list the calculated durations of the Class 0 + I and Flat-spectrum Classes with the Class II + Class III reference Class corrected for contamination and a reference Class duration of 3 Myr. We consider these values to be our most accurate and logically consistent lifetimes: 0.46–0.72 Myr for Class 0 + I sources, and 0.30–0.46 Myr for Flat-spectrum sources.

In Section 3.2, we reported that 33.5% of the objects associated with dense cores (and thus considered protostars according to the definition adopted in this paper) are classified as Class 0 by ${T}_{\mathrm{bol}}^{\prime }$, with the remainder classified as Class I. Under the same assumptions as above (continuous star formation and time is the only variable) and our definition of a protostar as a YSO associated with a dense core, our results imply that 33.5% of the total protostellar duration is spent in the Class 0 phase, with the remaining 66.5% of the time spent in the Class I phase. If we assume that the lifetime of the SED Class 0 + I objects is a good proxy for the lifetime of the protostellar stage of evolution (see Section 4.3 below for justification of this assumption), the range of SED Class 0 + I durations listed in Table 7 (0.40–0.78 Myr) imply resulting durations of 0.13–0.26 Myr for Class 0 protostars and 0.27–0.52 Myr for Class I protostars (see the last two columns of Table 7). If we only consider the last two rows of Table 7 (including Class III YSOs in the reference Class, correcting for estimated contamination rates in the Class III sample, and adopting a reference Class duration of 3 Myr), the resulting durations are 0.15–0.24 Myr for Class 0 and 0.31–0.48 Myr for Class I.

As noted in Section 3.2, we generally find the same ratios of SED Class 0 to SED Class 0 + I YSOs as in previous studies (Evans et al. 2009; Enoch et al. 2009; Dunham et al. 2014). While those studies derived a Class 0 duration of 0.15–0.17 Myr, our expanded range of possible durations is due to our expanded analysis of the contents and duration of the reference Class.

As noted in Section 3.2, the standard picture of low-mass star formation resulted from the merger of an empirical classification system based on the infrared spectral indices of YSOs (Lada & Wilking 1984) with theoretical studies of the stages of the collapse of a dense, rotating core (Shu 1977; Terebey et al. 1984). In this picture, SED Class 0 + I YSOs are protostars still embedded in and accreting from dense cores, SED Class II YSOs are T Tauri stars with optically thick, accreting disks, and SED Class III YSOs are pre-main sequence stars with little or no disks left (all Spitzer-identified Class III YSOs, by definition, still have infrared excesses, thus the population of Class III YSOs without disks is not identified in the work presented here). Flat-spectrum YSOs have no clear link to a distinct physical stage, and in reality the correlation between observational SED Class and physical Stage is weakened by the effects of geometry, extinction, and accretion history (e.g., Whitney et al. 2003; Young & Evans 2005; Robitaille et al. 2006; Crapsi et al. 2008; Dunham et al. 2010).

While our study provides the statistics necessary to measure durations of the observed SED Classes, it does not provide the necessary information to measure durations of physical Stages. In a separate, complementary study, Heiderman & Evans (2015) surveyed nearly all of the SED Class 0 + I and Flat-spectrum YSOs in the c2d and GB surveys for emission from the dense gas tracer HCO⁺ (3–2). Using criteria based on the detection and strength of this line, they sorted objects into Stage 0 + I (protostar) and Stage II (disk) sources. They found that the Flat-spectrum YSOs are approximately evenly divided between protostar and disk sources, with no evidence that they occupy a physically distinct evolutionary stage. They also found that the number of Flat-spectrum YSOs classified as protostars is balanced by the number of SED Class 0 + I YSOs classified as disks, coincidentally resulting in a nearly identical number of protostars (defined in their work based on the presence of dense gas traced by HCO⁺ line emission) and SED Class 0 + I YSOs (326 Class 0 + I YSOs identified here, 335 protostars identified by Heiderman & Evans). Thus, based on their results, the range of possible protostellar stage durations is the same as the range of possible SED Class 0 + I durations listed in Table 7.