THE SPITZER c2d LEGACY RESULTS: STAR-FORMATION RATES AND EFFICIENCIES; EVOLUTION AND LIFETIMES

Five large clouds were selected for the c2d project: Serpens (Eiroa et al. 2008), Perseus (Bally et al. 2008), Ophiuchus (Wilking et al. 2008), Lupus (Comerón 2008), and Chamaeleon (Luhman 2008). More specifically, we targeted Lupus I, III, and IV and Chamaeleon II (hereafter Cha II). These were chosen to lie within about 300 pc of the Sun, to span a range of previously-known star-formation activity, and to complement observations of smaller regions in these clouds obtained by Guaranteed Time Observers (GTOs). We have included data from the GTO observations with our data for a complete picture of the clouds.

The clouds in this study are listed in Table 1, with the adopted distances, solid angles and areas with both IRAC and MIPS data, a measure of turbulence, the total mass in the mapped area, and the crossing time. Some of these clouds were targets of ISO surveys of smaller regions: 0.7 deg² in Ophiuchus (Bontemps et al. 2001); 0.13 deg² in Serpens (Kaas et al. 2004); and 0.2 deg² in Cha II (Persi et al. 2003). The results presented here are based on detailed studies of the individual clouds: Serpens (Harvey et al. 2006, 2007a, 2007b); Perseus (Jørgensen et al. 2006; Rebull et al. 2007; S.-P. Lai et al. 2008, in preparation); Ophiuchus (L. Allen et al. 2008, in preparation; Padgett et al. 2008); Lupus (Chapman et al. 2007; Merín et al. 2008); and Cha II (Young et al. 2005; Porras et al. 2007; Alcalá et al. 2008; Spezzi et al. 2007, 2008). Numbers in those papers have been updated based on a consistent analysis of the combined data set.

None of these clouds are forming very massive stars; the most massive star in the close vicinity of a cloud mapped by c2d is HD147889, a B2 (III/IV) star lying within our map of Ophiuchus. The σ Sco group is farther away (projected distance from the cloud of 4 pc), consisting of an O9V star and a B2III star and at least two other B stars (Pigulski 1992). The star-formation patterns in Ophiuchus may have been affected by interaction with these stars (Loren 1989; Nutter et al. 2006). For Perseus, there are hints of an interaction with 40 Per, a B0.5 star in the Per OB Association, which is 26 pc from the L1451 core, as discussed by Walawender et al. (2004) and Kirk et al. (2006). Serpens contains a Herbig Ae star, VV Ser, which produces a large nebulosity (green object toward the south in Figure 1) discussed by Pontoppidan et al. (2007a, 2007b). In addition, Serpens contain three other A-type stars and a B8 star (Oliveira et al. 2009). The most massive star in the Cha II cloud is the F-type star DK Cha (Hughes & Hartigan 1992); all the rest have spectral type later than K2 (Spezzi et al. 2008). The Lupus clouds are in the vicinity of the Scorpius–Centaurus association, which flanks the Lupus clouds at a distance of about 5 degrees, or 17 pc at a distance of 200 pc. The sub-groups of the association are called Upper Scorpius (5–6 Myr) and Upper Centaurus–Lupus (14 Myr; de Geus et al. 1989). The strong high-energy radiation and supernova remnants from those OB stars might have had an important role in the formation and evolution of the Lupus clouds (Tachihara et al. 2001). Specifically, the highly fragmented cloud structure is thought to be related to the effects of the nearby OB association.

Zoom In Zoom Out Reset image size

There are distance uncertainties for all five clouds, which are discussed in the detailed papers on those clouds. Here we supply our adopted uncertainties and propagate those into the areas and masses. Recent parallax observations of radio emission from young stars in Ophiuchus yield a distance of 120 ± 4.5 pc (Loinard et al. 2008). By combining extinction maps with parallaxes from Hipparcos and Tycho, Mamajek (2008) found a distance of 135 ± 8 pc, while Lombardi et al. (2008) have derived a distance of 119 ± 6 pc. These are all consistent with our assumed distance and uncertainty of 125 ± 25 pc. Because it is not clear that the whole cloud is at the same distance (Loinard et al. 2008; Lombardi et al. 2008), we retain our value and uncertainty in the analysis. A recent astrometric measurement of water masers in NGC1333, in Perseus, provides a distance of 235 ± 18 pc (Hirota et al. 2008). Distance estimates for IC348 are larger, so we retain our standard distance and uncertainty of 250 ± 50 for the Perseus cloud as a whole. Lombardi et al. (2008) also derived a distance to the Lupus complex of 155 ± 8 pc, with some evidence that individual Lupus clouds were at different distances, both consistent with our assumed values of 150 ± 20 pc for Lupus I and IV and 200 ± 20 pc for Lupus III. We use a distance to Serpens of 260 ± 10 pc (Straizys et al. 1996), as was used in all our previous papers on this cloud, but we note that there is some recent evidence for a distance of 230 ± 20 pc (Eiroa et al. 2008).

Seeking a common measure of the turbulence of these clouds, we chose to give the FWHM line width of the ¹³CO J = 1 → 0 line, ideally averaged over a full map of the cloud. This was not available in all cases, so we have resorted to accepting statements like "a typical line width is..." (Hara et al. 1999) for Lupus. For Perseus, Serpens, and Ophiuchus, J. E. Pineda (2008, private communication) has provided the mean and median line widths for the ensemble of spectra in the COMPLETE maps (Ridge et al. 2006), and the average and difference of the median and mean are given in Table 1 as the value and uncertainty. N. F. H. Tothill et al. (2008, in preparation) have provided the same kind of information for the J = 2 → 1 line of ¹³CO in Lupus, based on maps with AST/RO. The line widths are based on dividing the integrated intensity by the peak intensity and correcting for the channel width of 0.91 km s⁻¹, assuming Gaussian lines. The results are substantially larger than the estimates from the ¹³CO J = 1 → 0 lines. For Lupus I, 〈Δv〉 = 2.17 ± 0.05; for Lupus III, 〈Δv〉 = 2.11 ± 0.05; and for Lupus IV, 〈Δv〉 = 1.53 ± 0.03. It would be valuable to have more consistent and sophisticated measures of the turbulence in these clouds.

The cloud mass (gas and dust) has been derived from extinction maps made from our c2d data, together with 2MASS data, toward background stars (Evans et al. 2007). We used the extinction law with R_V = 5.5 (Weingartner & Draine 2001), which reproduces reasonably well the data for molecular clouds (Flaherty et al. 2007; Chapman et al. 2009), to determine A_V. The mass was then calculated using the relation N_H/A_V = (1.086C_ext(V))⁻¹ (Draine 2003) and the value of C_ext(V) from the online tables¹⁵ for R_V = 5.5. This grain model results in a conversion from extinction to hydrogen column density of 1.37 × 10²¹ cm⁻² mag⁻¹, instead of the usual 1.87 × 10²¹, established for diffuse interstellar medium (ISM) gas (Bohlin et al. 1978). In our earlier papers (Harvey et al. 2007a; Alcalá et al. 2008; Chapman et al. 2007; Merín et al. 2008), the value for the diffuse gas was used, resulting in cloud masses overestimated by a factor of 1.4 relative to those in Table 1.

The mass includes all cloud mass above an extinction contour with A_V = 2 mag, with no assumptions needed about geometry. The mass derived from our extinction maps refers to the same area covered in our surveys, making it ideal for calculations of efficiency, etc. The two exceptions to this statement are Serpens, for which our survey was designed to cover completely down to A_V = 6, and Ophiuchus, for which we tried to cover completely down to A_V = 3. The coverages were based on extinction maps by Cambrésy (1999), but some areas down to A_V = 2 were covered incompletely. The mass inside the A_V = 6 contour for Serpens is 1532 M_☉, 76% of the mass within A_V = 2. The mass inside the A_V = 3 contour for Ophiuchus is 1914 M_☉, 88% of the mass within A_V = 2. We use the value for A_V = 2 in what follows, but the difference will not be large for these clouds. Table 1 also gives the cloud crossing time in Myr, calculated from the area and the mean speed of turbulent motions, approximated by 〈v〉 = 0.68Δv, with Δv the ¹³CO J = 1 → 0 line width.

With the sensitivity of Spitzer, we can detect sources with luminosity as low as 10⁻³ L_☉ at 350 pc, but contamination by other kinds of sources ultimately limits our sensitivity, as discussed below. A similar study of the Taurus cloud has been done by D. L. Padgett et al. (2008, in preparation), and other clouds in the Gould Belt are being mapped with the same techniques (L. Allen et al. 2008, in preparation). When those studies are completed, the analysis in this paper can be extended to cover most star formation in large clouds within about 300 pc of the Sun.

The main challenge is to separate YSOs from contaminants. By far, most contaminants are stars, mostly background, without infrared excess. These have colors very close to zero in a color–color diagram using the IRAC bands (Allen et al. 2004) and are easily removed. The primary remaining contaminants are then background galaxies with active star formation; these may have colors very similar to those of embedded young objects, so we must use magnitude information as well as colors.

The automated criteria described by Harvey et al. (2007a) were used for all five clouds to produce catalogs of YSOc, which are part of our final data delivery. Compared to earlier catalogs, advances in distinguishing YSOs from background galaxies (Harvey et al. 2007a) have provided much cleaner samples. These criteria assign to each source with sufficient information an unnormalized probability (P(Galc)) that the source is a galaxy. The calculation of the probability is complex and fully described in Harvey et al. (2007a) and Evans et al. (2007), so we summarize it here. The probability depends on the location of the source in three color–magnitude diagrams, whether the source is extended, and whether the source is detected above a threshold at 70 μm. The assignment of probabilities is entirely empirical, based on excluding objects in a suitable sample of extragalactic observations.

By comparing to a set of sources from the SWIRE survey of the ELAIS N1 extragalactic field (Surace et al. 2004), suitably processed to simulate the c2d sensitivity and the extinction distribution of each cloud, we define a threshold in P(Galc) below which a source becomes a YSOc. For the catalogs of the five large clouds, the aggregate plot is shown in Figure 3. As was true for the individual clouds, most YSOc separate nicely from galaxy candidates (Galc), but there are some sources in a lower probability tail from the big peak in the Galc sector that are ambiguous. Thus, confusion with star-forming galaxies can be a problem for both contamination and completeness. Sources with logP(Galc) < −1.47 were assigned the YSOc moniker in the c2d catalogs, based on the analysis of Serpens by Harvey et al. (2007a). From the surface density of objects from the degraded SWIRE data that would be misclassified as YSOc for each cloud, multiplied by the surface area of the cloud, we estimate that there could be as many as 51 contaminating galaxies in our total 15.5 deg². That would still be a small fraction of the total of 1086 YSOc (see below). We next describe attempts to minimize the number of contaminants.

Zoom In Zoom Out Reset image size

Certain categories of sources were checked by eye for all the clouds in an attempt to eliminate residual contaminants. Typically these were sources that had been "band-filled" at 24 μm or sources with peculiar SEDs. Sources that were well detected at some wavelengths received flux estimates at other wavelengths, based on the known position, in the process of band filling. Some of these flux estimates turned out to be emission from the wings of nearby bright sources and other artifacts. In addition, objects with the "PAH_em" suffix received extra scrutiny because this feature is common in star-forming galaxies. Objects that had galactic morphologies were removed from the final YSO list; however, many star-forming galaxies are point-like to Spitzer. A total of 91 YSOc were eliminated by this process. We do not know how many point-like galaxies were removed by this process, but the numbers are consistent with removing most of them.

While stars without infrared excess are readily removed from our sample via color criteria, post-main-sequence stars with circumstellar shells can masquerade as YSOs. Harvey et al. (2007a) were able to identify and remove four of these behind the Serpens cloud. Follow-up optical spectroscopy toward Serpens (Oliveira et al. 2009) have provided spectral types for 78 objects, 58 of which had been identified as YSOs. Based on their analysis, we rejected 11 of those 58 sources (about 20%) as background giants with infrared excesses (Oliveira et al. 2009). All 11 are classified as either Class II or Class III sources, with nine out of 11 (about 80%) classified as Class III. We have removed these 11 sources from the final sample. We so far lack the data required to identify and remove background giants in Lupus, Perseus, and Ophiuchus. However, Serpens should be the worst case because it lies at low galactic latitude and longitude. Lupus and Ophiuchus may also be affected, though less so than Serpens. Additionally, a spectroscopic study of Cha II found that 96% of the YSOs identified in our catalogs were true cloud members (Spezzi et al. 2008). These points, combined with the fact that most of the background giants (80%) are Class III sources, to which we are incomplete anyway (Section 3.2), lead us to conclude that contamination by background giants will not significantly affect our main results.

Adding these 11 background giants in Serpens to the 91 sources removed as described above gives a total of 102 sources removed from the initial sample of 1086 YSOc (Table 2).

Since the vast majority of sources in the catalogs are background stars, we require an infrared excess for a source to be considered a candidate YSO. Consequently, our sample is missing Pre-Main-Sequence Stars (PMS) that no longer have infrared excess but may have Hα emission, X-ray emission, etc. Complete surveys for such objects are needed to complete the sample of PMS objects, so we concentrate on YSOs, defined here to have an infrared excess. More complete samples of PMS objects exist for some clouds that have the needed complementary data, but not for all. For example, Alcalá et al. (2008) found 51 certain and 62 likely PMS stars in Cha II, more than double the number (26) of YSOs. The ratio of PMS to YSOs in our catalog is 2 for Class II sources and 4.8 for Class III sources. Cha II may be a worst case because many of the PMS stars are not in the region studied by c2d. Restricting attention to that region, only four PMS stars are added to the YSOs found by c2d. A similar study in Lupus would add 19 PMS stars without infrared excess to the 94 YSOs (Merín et al. 2008).

While objects of very low luminosity may still be lost among the galaxy background, we estimate that our YSO sample is 90% complete down to a luminosity integrated from 1 to 30 μm of 0.05 L_☉ and 50% complete down to 0.01 L_☉ (Harvey et al. 2007a) for our most distant cloud (Serpens at 260 pc, which also has the highest density of background stars). For comparison, the survey of Serpens with ISO reached a limit of 0.08 L_☉, but had only two bands and covered only 0.13 deg² (Kaas et al. 2004). They found 61 YSOs compared to our 227, mostly because of their smaller areal coverage. The 1–30 μm wavelength range covers the bulk of emission from YSOs with infrared excesses arising from circumstellar disks; thus these completeness limits are a good proxy for our completeness to such objects. A 2 Myr old object at the mass boundary between stars and brown dwarfs has a luminosity of about 0.01 L_☉. Indeed some objects believed to be brown dwarfs with disks based on more complete analysis do not make our YSOc list, either because their fluxes are not of high enough quality or because they have logP(Galc) ⩾ −1.47. For example, two of the four added PMS objects in Cha II are very low-mass brown dwarfs with excesses (Allers et al. 2007).

For younger YSOs still embedded within their dense cores, this 1–30 μm wavelength range is less appropriate since the bulk of the emission is reprocessed by the surrounding envelope to the far-infrared. A separate search for embedded objects with luminosities less than 1 L_☉ found that we reach a similar completeness limit for the internal luminosity of these objects (L_int of 4 × 10⁻³ L_☉ at 140 pc, or ∼0.014L_☉ at 260 pc; Dunham et al. 2008), where this limit is set by the sensitivity of the 70 μm MIPS-2 observations. The internal luminosity measures the contribution from the embedded source, after correction for the effects of heating by the interstellar radiation field.

We require good quality detections in all four IRAC bands and MIPS-1 for a source to be classified as a YSOc. This requirement can result in our missing two types of sources: deeply embedded sources that were not detected in all bands and strong sources that saturated the detectors. Note that the latter can include both embedded objects and more evolved YSOs no longer embedded within their dense cores. A total of 40 sources that are clearly YSOs but which did not make the automatically-generated YSOc list for one of the above two reasons were added by hand to our final sample, bringing the final sample size to 1024.

Of these 40 sources, 36 were added by comparison to the searches for embedded objects presented by Dunham et al. (2008), Enoch et al. (2009), and Jørgensen et al. (2007, 2008). Two were added by being well-known YSOs that were saturated in one or more of the Spitzer bands but not included in the samples of embedded objects compiled by the above authors. For all saturated sources, data from other telescopes were substituted for the saturated Spitzer data. Three of the added sources actually have source types of Galc, one each in Lupus, Perseus, and Serpens. The source in Lupus was added by Merín et al. (2008) because it lies essentially on the border between YSOc and Galc, and it is candidate to be a brown dwarf with a disk (Allers et al. 2006). The source in Perseus was added because it is contained in the sample of embedded objects compiled by Jørgensen et al. (2007). The Serpens source is associated with an outflow.

Turning the luminosity completeness limit into a limit on stellar mass requires further analysis. For the early, embedded stages, the luminosity depends on the product of stellar mass and mass accretion rate. As discussed in Section 9.2, accretion at the mean rate expected in a Shu-type model would predict a luminosity of 1.6 L_☉ for a central object of 0.08 M_☉. Our luminosity limit of 0.014 L_☉ would then translate into a mass limit of 7 × 10⁻⁴ M_☉. However, there is strong evidence for highly variable accretion rates (Section 9.2), so this limit is highly suspect. For the later stages when accretion luminosity is negligible, spectral types and accurate extinction corrections are necessary to determine masses. For the most distant cloud (Serpens), a main-sequence luminosity of 0.05 L_☉ implies a mass of 0.08 M_☉, while 0.01 L_☉ corresponds to 0.04 M_☉ (Chabrier et al. 2000). At young ages, the mass limits would be lower. In Serpens, the masses determined by Oliveira et al. (2009) from spectral types range from 0.2 to 3.0 M_☉. In Lupus, masses of YSOs are complete down to 0.1 M_☉ (Merín et al. 2008), and in Cha II, masses extend well below the stellar limit to 0.015 M_☉ (Spezzi et al. 2008). We do clearly miss some substellar objects with disks, as discussed above. Until spectroscopy is available for a larger fraction of the young objects, we can only estimate our mass completeness to be near the stellar/substellar boundary.

For each source with sufficient data, the c2d catalogs provide a least-squares fit to all photometry between 2 μm and 24 μm in the c2d catalog to determine the spectral index α:

$$\begin{equation} \alpha = {d\log (\lambda S(\lambda)) \over d\log (\lambda)}, \end{equation} \tag{ 1 }$$

where λ is the wavelength and S(λ) is the flux density at that wavelength. These were used to classify objects into the four classes defined by Greene et al. (1994), as described in detail in Section 5.1: Class I, Flat, Class II, and Class III.

An alternative classification scheme, used especially for more embedded objects, employs the bolometric temperature (T_bol). We will also need the bolometric luminosity (L_bol) for later analysis. These were computed, using both Spitzer data and auxiliary data at both shorter and longer wavelengths, by calculating the first two moments of the SED (Chen et al. 1995) using integration methods described in more detail by Dunham et al. (2008). Uncertainties in L_bol and T_bol are dominated for early phases by incomplete sampling of the SED at wavelengths from 70 to 1000 μm; they were estimated to be 20–60% (Enoch et al. 2009; Dunham et al. 2008). The poor sampling at far-infrared wavelengths arises in part because of the large beam of Spitzer, causing confusion in crowded regions, and from the lack of data at 350 μm. A recent analysis of the Serpens B cluster using the fine-scale 70 μm mode and adding 350 μm data (Harvey & Dunham 2009) was able to separate sources that were confused in the c2d data. They find that the values of L_bol and T_bol change by about the amount estimated above. Nonetheless, only one source out of 18 would change classification (from Class I to Class 0).

For some of the analysis, it is desirable to correct the flux densities for extinction. Based only on the c2d data, such corrections would be highly uncertain. Optical follow-up studies have provided spectral types for the Class II and Class III populations in Cha II (Spezzi et al. 2007, 2008), Lupus (Comerón et al. 2009; Merín et al. 2008), and part of Serpens (Oliveira et al. 2009). For Ophiuchus and Perseus, along with the remainder of Serpens, we assume a spectral type of K7 for the Class II and Class III sources. We then used the spectral type (known or assumed), the near-infrared photometry (J, H, and K in order of preference), and the Weingartner & Draine (2001) extinction law for R_V = 5.5 to correct the photometry for extinction. The extinction law was chosen to match that used for determining our extinction maps and cloud masses (Section 1.1). For Class I and Flat objects, we assume the mean extinction to all the Class II objects in the same cloud. The idea is that this correction removes foreground extinction, but not local extinction from the surrounding envelope, which will be reradiated in the far-infrared. The mean extinctions to Class II sources for the five clouds are A_V = 3.95 mag in Cha II, A_V = 2.91 mag in Lupus, A_V = 5.92 mag in Perseus, A_V = 9.57 mag in Serpens, and A_V = 9.76 mag in Ophiuchus. Extinction values calculated for Class II and III sources without spectral type information should be regarded as highly uncertain. Extinction values for Class I and Flat sources should be regarded as averages for each cloud only and not as actual extinctions toward each object.

After these extinction corrections were made, the values of spectral index, T_bol, and L_bol were recomputed and these values are distinguished from the observed values by primes (α', T'_bol, and L'_bol). Both sets of values are given, along with the value of A_V used for the extinction correction, in Table 12. When wavelength coverage is insufficient to obtain reliable values for T'_bol and L'_bol, a flag in the table indicates this and appropriate upper or lower limits are given (see the Appendix for more details).

Combining all the clouds, we have the following statistics based on the final c2d data release. The full catalogs contain a total of 4.26 × 10⁶ entries, of which 6.14 × 10⁵ are also in the high reliability catalog, which requires detection in at least one band with signal-to-noise ratio (S/N) ⩾7 and a second band, if available at that position, with S/N ⩾5. A total of 6.77 × 10⁵ have detections in at least three bands and can be further classified. This group provides the parent sample from which our YSOc are drawn. The great majority are classified as stars (3.32 × 10⁵) or other (2.55 × 10⁵). The "other" category contains sources that do not fit any particular template; the great majority are probably galaxies. Of those that remain, 2965 are galaxy candidates and 1086 are YSO candidates. The statistics for individual clouds are given in Table 2.

The processes described in Sections 3.1 and 3.2 resulted in the list of YSOs used in this paper. There are undoubtedly still some contaminants in this list, but we believe that they are a small fraction. Overall the ratio of YSOs to YSOc is 0.94. The ratio ranges from 0.85 in the Lupus clouds to 0.99 in Perseus. In both Perseus and Ophiuchus, the fraction is high because of the addition of deeply embedded and saturated sources not in the original YSOc list.

In the end, there are a total of 1024 YSOs in our sample. This number represents an order of magnitude increase over previous samples, and they have been selected in a uniform way from data with very similar sensitivity.

Given the mass surface density of the cloud from our extinction maps, one can predict the star-formation rate surface density from the relations employed for other galaxies using the formula from Kennicutt (1998):

$$\begin{eqnarray} \Sigma ({\rm SF})({\rm \mbox{$M_\odot $}\ yr^{-1} kpc^{-2}}) & & \nonumber \\ =(2.5\pm 0.7) \times 10^{-4} (\Sigma ({\rm gas})/1 \mbox{$M_\odot $}\ {\rm pc^{-2}})^{1.4\pm 0.15}. & & \end{eqnarray} \tag{ 2 }$$

Taking all clouds together, the surface density is 64 M_☉pc^-2, and the predicted star-formation rate would be 0.08 M_☉ Myr⁻¹ pc⁻², a factor of 20 below the observed value of 1.6. As shown in Figure 4, all clouds lie well above the prediction of Equation (2), even though these clouds are forming only low-mass stars, which would be largely invisible to the kinds of tracers, such as Hα emission, used to establish Equation (2).

Zoom In Zoom Out Reset image size

This difference is not surprising, but it reminds us that the Kennicutt relation applies to averages over much larger regions than individual clouds. The coefficient and exponent in such relations are clearly functions of the mean density of the region being averaged over. A census of all star formation within 500 pc should come from the completion of many Spitzer observations; a first crude estimate, along with information on the surface density of gas, gave rough agreement with the predictions of Kennicutt's relation (Evans et al. 2008). Note that 85% of the gas in this region is atomic hydrogen. Blitz & Rosolowsky (2006) give a more detailed discussion of these points.

Studies of star formation in dense gas have shown a linear relation between star-formation rate and the amount of gas traced by HCN emission (Gao & Solomon 2004a; Gao & Solomon 2004b; Wu et al. 2005). As discussed below, we find that star formation is sharply confined to the dense gas, so it is two steps removed from the scales on which Equation (2) applies: the formation of molecular clouds from atomic gas; and the formation of dense cores from molecular clouds (Evans et al. 2008). At which of these steps does the nonlinearity (an exponent that is greater than unity) of Equation (2) enter?

To address that question, studies that resolve individual clouds in other galaxies are needed. Studies of M51 with 0.5–2 kpc resolution continue to show a nonlinear relation with exponent around 1.4 (Kennicutt et al. 2007). Recent observations with sub-kpc resolution of both molecular and atomic gas in other nearby galaxies (Bigiel et al. 2008) reveal a linear relation between star-formation rate and molecular gas density over a range from 3 to 50 M_☉ pc^-2:

$$\begin{eqnarray} \Sigma ({\rm SF})({\rm \mbox{$M_\odot $}\ yr^{-1} kpc^{-2}}) & & \nonumber \\ =10^{-2.1\pm 0.2} (\Sigma ({\rm H}_2)/10 \mbox{$M_\odot $}\ {\rm pc^{-2}})^{1.0\pm 0.2}. & & \end{eqnarray} \tag{ 3 }$$

(Note that this relation is normalized to a different surface density.) This relation is also shown in Figure 4 over the range where it was established and extrapolated with a dotted line. Equation (3) would predict 0.051 M_☉ yr⁻¹kpc⁻², much less than our observed value of 1.6. Bigiel et al. (2008) note that they are still measuring the filling factor of clouds rather than resolving structure within molecular clouds. This comparison suggests that clouds with the properties we study are filling only about 3% of their beams. As resolution improves further, and as more studies become available within our own galaxy, especially of regions forming high-mass stars, further comparisons will be important.

We also plot the relation found for dense gas, as traced both in other galaxies and in massive dense cores in our Galaxy by HCN, by Wu et al. (2005). This relation extrapolated to lower surface densities comes closest to our observed points.

Following Elmegreen (2002), a linear relation in molecular gas would correspond to a threshold core density for star formation of 10⁵ cm⁻³, in good agreement with typical densities in the dense cores. However, he also predicts that the fraction of total gas in dense cores would be 10⁻⁴, whereas the fraction in the three clouds with the relevant data is 4.6%. Allowing for the fact that 85% of the gas in the local kpc is not molecular gets the fraction down to 6 × 10⁻³, still on the high side. The nonlinear relation that seems to apply on larger scales was interpreted by Elmegreen (2002) to mean that the critical criterion is a ratio of core density to mean ISM density of 10⁵. Caution is still warranted as the extragalactic relations need further analysis, and other factors, such as shear and bars, will affect star formation, especially in circum-nuclear regions (e.g., Jogee et al. 2005).

Comparison of the mass in YSOs to the cloud mass gives a measure of the current efficiency. Since we are only sensitive to YSOs with infrared excess, this efficiency represents an average over the last 2 Myr. We give the star-formation efficiency, defined as

$$\begin{equation} {\rm SFE} = {\mbox{$M_{\star }$}\over \mbox{$M_{\star }$}+ \mbox{$M({\rm cloud})$}}, \end{equation} \tag{ 4 }$$

in Table 4. The values range from 3% to 6% over the various clouds, with a value of 4.8% taking all clouds together. Given the star-formation rate $(\dot{M}_{\star })$ from Table 3, we can compute a depletion time for the cloud:

$$\begin{equation} \mbox{$t_{\rm dep}$}= \mbox{$M({\rm cloud})$}/\mbox{$\dot{M}_{\star }$}. \end{equation} \tag{ 5 }$$

Table 4. Efficiencies and Depletion Timescales

Cloud	${M_*\over (M({\rm cloud}) + M_*)}$	M*/M(dense)	tdep(cloud)	tdep(dense)	SFRff	Notes
			(Myr)	(Myr)
Cha II	0.030		66	⋅⋅⋅	0.028	Alcalá et al. (2008)
Lupusa	0.054		35	⋅⋅⋅	0.050	Merín et al. (2008)
Perseus	0.038	0.69	50	2.9	0.049	S.-P. Lai et al. (2008, in preparation)
Serpens	0.053	1.2	35	1.6	0.036	Harvey et al. (2007a)
Ophiuchus	0.063	3.3	30	0.6	0.064	L. Allen et al. (2008, in preparation)
All Cloudsb	0.048	1.2	40	1.8	0.040	⋅⋅⋅

Notes. ^aA sum over the three Lupus clouds of the values for each cloud. ^bFor all but SFR_ff, this number is calculated by adding all clouds with the relevant data together; for SFR_ff, it is the average over all clouds of the individual values.

Download table as:  ASCIITypeset image

As shown in Table 4, these t_dep are 30–66 Myr. Estimates of cloud lifetimes range from 10 to 40 Myr (McKee & Ostriker 2007). If clouds produce stars at the current rates for 10 Myr, the final efficiency when star formation has ended in these clouds could be as high as 15–30%; as discussed in Section 4.4, achieving such high efficiencies would require continued conversion of cloud material into dense cores.

Krumholz & Tan (2007) have emphasized that star formation, even in molecular clouds, is slow in the sense that the star-formation rate per unit free-fall time is low. Krumholz & McKee (2005) define the star-formation rate per unit free-fall time (SFR_ff) to be the fraction of an object's mass that turns into stars in a free-fall time at the object's mean density. Using our cloud masses from the extinction maps covering the same areas where we surveyed for YSOs, we can quantify this, again where we mean star formation over the last 2 Myr since we identify only the sources with infrared excess. Then we have

$$\begin{equation} \mbox{${\rm SFR}_{{\rm ff}}$}= \mbox{$\dot{M}_{\star }$}\mbox{$t_{{\rm ff}}$}/ \mbox{$M({\rm cloud})$}= \mbox{$t_{{\rm ff}}$}/\mbox{$t_{\rm dep}$}, \end{equation} \tag{ 6 }$$

where t_ff is the free-fall time for the mean density of the cloud, calculated from

$$\begin{equation} \mbox{$t_{{\rm ff}}$}= 34 {\rm Myr}/\sqrt{\mbox{$\langle n\rangle $}}. \end{equation} \tag{ 7 }$$

The density 〈n〉 in that calculation is the density of all particles, and we have assumed a mean molecular mass of 2.3 amu. Mean densities for the clouds, computed from the mass and surface area by assuming spherical clouds, are given in Table 1. The average over all clouds is 〈n〉 = 390 cm⁻³. Then, we find SFR_ff ranging from 0.028 to 0.064, with an average over all clouds of 0.040. This value is on the high side of the range of SFR_ff inferred from more global considerations by Krumholz & Tan (2007) and the theory of Krumholz & McKee (2005).

Figure 2 shows a strong correlation of YSOs, especially the younger ones, with regions of high extinction. For the clouds with millimeter continuum maps, a strong correlation of the youngest objects with dense cores identified in the continuum maps is apparent (Enoch et al. 2007). As discussed by Enoch et al. (2007), the cores identified in this way have minimum mean densities of about 2 × 10⁴ cm⁻³ and more typically more than 10⁵ cm⁻³, and almost all seem to be gravitationally bound (Enoch et al. 2008). Their mean densities are 50–200 times the mean density of the cloud; they thus represent quite distinct entities, rather than just modest peaks in column density. Thus, it is also interesting to calculate the efficiencies and rates using the total mass of gas in dense cores for the three clouds with Bolocam data, using the data from Enoch et al. (2007). The total mass in YSOs in the cloud is quite comparable to the mass in dense cores (Table 4); taking all three clouds together the ratio is 1.3. We compute the time to deplete the current stock of mass in discrete dense cores by dividing the SFR by the mass in dense cores. We find 0.6 Myr for Ophiuchus, 1.6 Myr for Serpens, and 2.9 Myr for Perseus. Taking all clouds together, t_dep(dense) = 1.8 Myr. These times are consistent with the 2 Myr timespan for detectability of YSOs that we use to calculate lifetimes and with plausible spreads of formation times in clusters. While star formation is faster and more efficient in the dense gas probed by Bolocam, it is still slow compared to a free-fall time. Using a mean density of 5 × 10⁴ cm⁻³ to calculate the free-fall time for the typical core leads to SFR_ff(dense) = 0.05 for Perseus to 0.25 for Ophiuchus.

These values for SFR_ff(dense) are somewhat higher than SFR_ff for the whole cloud, but they neglect the fact that the mass of stars is now comparable to the mass of dense gas. The depletion times also assume that all the gas winds up in stars once it reaches the density of the cores identified by Enoch et al. (2007). We can make a second calculation in which we assume that star formation began at some time in the past with a larger reservoir of dense gas and that a fraction of the dense core mass winds up in a star, while a fraction 1 − is lost to the star-forming region by outflows. By comparing the core mass function to the IMF of stars, Alves et al. (2007) estimated that = 0.3 and Enoch et al. (2008) find >0.25. At any given time, $\mbox{$t_{\rm dep}$}= \mbox{$M({\rm dense})$}/\mbox{$\dot{M}_{\star }$}$ and the mass of dense gas decreases exponentially, with a time constant of t_dep. At time t after the start of star formation, the expression for t_dep becomes

$$\begin{equation} \mbox{$t_{\rm dep}$}= {t\over \epsilon \ln (1+\eta /\epsilon)}, \end{equation} \tag{ 8 }$$

where η = M_⋆/M(dense) at time t. As usual, we take t = 2 Myr as the timescale over which we have complete statistics, assume = 0.3, and use the observed value of η. This expression then yields values for t_dep that are more like 3–6 Myr and SFR_ff(dense) that are 0.03–0.06, comparable to those for the cloud as a whole, as Krumholz & Tan (2007) would predict.

The current working model for star formation arose in the 1980s with the merger (Adams et al. 1987) of an empirical classification scheme, based on the slope of the SED between 2 and 20 μm (Lada & Wilking 1984), with a theoretical picture of star formation involving the collapse of an isolated rotating dense core (Terebey et al. 1984) forming a star and disk (Adams & Shu 1986). The essential stages of the theoretical model were three: the collapse of the envelope, forming the protostar and disk; the continued accretion of disk material onto the forming star; and the dissipation of the disk by planet formation, evaporation, etc. Each of these stages became identified with an empirical SED class as follows: the envelope collapse as Class I; the accretion disk and star as Class II; and the dissipation of the disk during Class III. With further study, other significant events could be distinguished theoretically, and further refinements of the classification system were suggested. We follow the suggestion of Robitaille et al. (2006) in referring to the physical arrangement as a "Stage" and the SED characteristic as a "Class."

We now focus on the evolution of the empirical classification system. Following the recognition by Lada & Wilking (1984) that SEDs were falling into natural groups, Lada (1987) first codified the tripartite class system using the spectral index. The original boundaries were as follows:

• I 0 < α ≲ 3;
• II −2 ≲ α ⩽ 0;
• III −3 < α ≲ −2.

Note that Lada and coworkers consistently use the symbol a, but we have standardized on α. These were supplemented by descriptions; for example, Class III sources could have some "mid-infrared excess," but "no or little excess near-infrared emission." Class II and III objects were all visible, but Class I objects were all invisible for λ < 1 μm. Class II sources were identified as T Tauri stars. It is worth noting that this system and many of the later elaborations were based on the study of the cluster of sources in the L1688 (Ophiuchus) cloud.

By 1989, Wilking et al. (1989) had noted that some SEDs had "double humps," one around 1 μm and one in the far-infrared. These were put into subclasses, labeled II-D and III-D. Class II-D sources were T Tauri stars with far-infrared excesses and were thought to be intermediate between Classes I and II. Class III-D sources looked like reddened photospheres plus a far-infrared excess (Lada 1991). By comparing the numbers of objects in different classes, Wilking et al. (1989) also assigned relative lifetimes to the different classes (discussed below).

The next major development was the introduction of a class "before" Class I by André et al. (1993). Since zero was not part of the original Roman number system, the Arabic symbol for zero was used. Since the Class 0 objects could not then be observed at the wavelengths originally used for classification, André et al. (2000) listed three criteria for a Class 0 source:

• Indirect evidence for a central YSO, as indicated by, e.g., the detection of a compact centimeter radio continuum source, a collimated CO outflow, or an internal heating source;
• Centrally peaked but extended submillimeter continuum emission tracing the presence of a spheroidal circumstellar dust envelope (as opposed to just a disk);
• High ratio of submillimeter to bolometric luminosity, suggesting that the envelope mass exceeds the central stellar mass: L_smm/L_bol > 0.5%, where L_smm is measured longward of 350 μm. In practice, this often means an SED resembling a single-temperature blackbody at T ∼ 15–30 K.

The boundary between Class 0 and Class I was associated with a new distinction between physical stages (André et al. 1993): the point at which the masses of the protostar and the remaining envelope were about equal.

The last major development came in 1994, when Greene et al. (1994) formalized the 4-class system,¹⁷ again based on the L1688 cluster, as follows:

• I 0.3 ⩽ α;
• Flat −0.3 ⩽ α < 0.3;
• II −1.6 ⩽ α < −0.3;
• III α < −1.6.

Greene et al. (1994) described the evolutionary status of sources in the new Flat Class as "uncertain," but Calvet et al. (1994) showed that these could be interpreted as infalling envelopes. Greene et al. (1994) noted that they found no sources with −1.9 < α < −1.6 or −0.3 < α < 0. These classifications were based on the slope between 2 and 10 μm, though they claimed that comparison to classes based on slopes between 2 and 20 μm showed that this did not matter. Greene et al. (1994) did not list Class 0 in their list, but they did acknowledge André & Montmerle (1994), who pointed out that millimeter wavelength emission became much weaker for α < −1.5, as a rationale for the revised boundary between Classes II and III.

Improving sensitivity at millimeter wavelengths led to the detection of some starless cores (Ward-Thompson et al. 1994), and these were called "pre-protostellar cores" (PPCs) or later "prestellar cores." To fit them into the class system, one might use Class −1, but Boss & Yorke (1995) argued that they should instead be Class −2 to make room for a theoretically important event, the formation of the first (molecular) hydrostatic core, which would then be called Class −1. These terms have not caught on, and the recent definitive review by di Francesco et al. (2007) simply distinguishes "prestellar" cores among the larger set of starless cores as being gravitationally bound.

With the increasingly baroque nature of the class system, it was natural to seek a continuous variable to capture the transitions. By analogy with the effective temperature of stars and spectral classes, Myers & Ladd (1993) suggested the use of a bolometric temperature, defined as the temperature of a blackbody with the same flux-weighted mean frequency as the actual SED, considering all wavelengths with available data. Chen et al. (1995) showed that the traditional classes, including 0, could be associated with certain ranges of T_bol, as follows:

• O T_bol < 70;
• I 70 ⩽ T_bol ⩽ 650;
• II 650 < T_bol ⩽ 2800.

For very early stages, data at short wavelengths were lacking, and T_bol could be quite sensitive to whether or not such data existed. Detections at shorter wavelengths could increase T_bol enough to move a source from Class I to II. Orientation of aspherical envelopes and disks can also affect T_bol substantially. At the other wavelength extreme, improved availability of submillimeter data considerably increased the number of Class 0 sources (Visser et al. 2002; Young et al. 2003).

For distinguishing Class 0 from Class I, the ratio of luminosities (L_bol/L_smm) may be most useful. In this scheme, a ratio exceeding 200, the inverse of criterion 3 of André et al. (2000), marks the transition from 0 to I. This form of the ratio has the virtue of increasing with the evolutionary progression, as does T_bol, and Young & Evans (2005) found that this ratio was a more robust indicator than T_bol of the ratio of mass in the star to mass in the envelope, based on models of collapsing cores. Both Visser et al. (2002) and Young et al. (2003) found that quite a few objects would be classified as 0 by the ratio, but as I by T_bol. A drawback to using L_bol/L_smm is the lack of complete data at 350 μm, the shortest wavelength included in L_smm; we lack such data for many of our sources in this study. Clearly, conclusions about evolution will depend on establishing a clear connection of such parameters to physically relevant stages of evolution.

It remains difficult to capture the increasing amount of information in any single parameter. With the greatly increased and more uniform data set available from the c2d project, we will re-evaluate the different tracers to learn what, if any, well-defined criteria can trace evolution through all stages. The distinction between physical stages of evolution and SED classes, however defined, has been usefully emphasized recently by Robitaille et al. (2006, 2007), who presented large grids of model SEDs from two-dimensional radiative transfer calculations, following earlier work by Whitney et al. (2003), which demonstrated the importance of inclination effects. Similarly, Crapsi et al. (2008) have discussed the possibility of confusion between stages and classes in the context of a grid of two-dimensional models.

One of the goals of research into star formation is to constrain the lifetimes associated with different stages of evolution. As noted above, Wilking et al. (1989) used the numbers of objects in various classes to establish relative lifetimes. This method assumes that the census is complete and that star formation in the sample has been a continuous process at a steady rate for longer than the dwell time in the slowest phase considered. When extrapolated from a single region, it also assumes that other variables are irrelevant, such as the total mass available, the conditions in the region, such as turbulence, etc.

With these caveats in mind, Wilking et al. (1989) used their data on the Ophiuchus cluster to suggest lifetimes. They estimated the ages of the Class II sources, finding an average of 0.39 ± 0.17 Myr, with the oldest at about 1.5 Myr. They found roughly equal number of Class I and Class II sources, suggesting equal lifetimes. With some other considerations, they suggested a duration for the Class I phase of 0.2–0.4 Myr. Later, Greene et al. (1994) found that Class I plus Flat SED lifetimes in Ophiuchus were 75% of the Class II lifetime, which they took to be 0.4 Myr. A somewhat older population in Ophiuchus, with ages around 2 Myr, was later established by Wilking et al. (2005).

Kenyon & Hartmann (1995) analyzed data from the Taurus cloud and found about 10 times as many Class II plus Class III sources as Class I sources, very different from the situation in Ophiuchus. They found a similar number of Flat SED sources and suggested lifetimes of 0.1–0.2 Myr for each of the Class I and Flat classes, based on a duration for Class II plus Class III of 1–2 Myr, which they found from comparison to evolutionary tracks.

More recent analyses of the Perseus cloud find a lifetime for the Class 0/I SED of 0.25–0.67 Myr within 95% confidence level (Hatchell et al. 2007). These estimates include Spitzer data on embedded sources and analysis of the age spread of the IC348 cluster from Muench et al. (2007). Early estimates of 0.01 Myr for the lifetime of the Class 0 phase were based on the small number of Class 0 sources in Ophiuchus (André et al. 1993). Because the first examples of Class 0 sources found were quite luminous with powerful, collimated outflows, the short lifetime suggested a phase of rapid accretion in which at least half of the final stellar mass was accreted (André et al. 2000).

Most of these estimates were based on samples of 50–100 sources and in particular clouds. They are subject to small number statistics and to possible differences from cloud to cloud, especially for the earlier classes. An exception was a compilation of 95 Class 0/I objects (Froebrich 2005). The positions of these sources in the T_bol–L_bol plane was compared to grids of models, and Froebrich et al. (2006) favored a Class 0 lifetime of 0.02–0.06 Myr.

Table 5 presents the numbers of YSOs in each Greene et al. (1994) class as determined by 2MASS–Spitzer 24 μm photometry. They are also shown graphically in Figure 5. We remind the reader that Greene et al. (1994) used photometry only out to 10 μm, but found no significant differences in α for sources with data out to 20 μm. We have not yet separated Class 0 sources from Class I sources in Table 5; this separation will be discussed in Section 7.1. In addition, Class III is certainly missing sources, as discussed in Section 3. So at this point, we restrict comparison to Class I, Flat, and Class II sources.

Zoom In Zoom Out Reset image size

With 1024 YSOs and more than 100 in each class, the uncertainties from small number statistics are decreased to less than 10% effect in any class, almost certainly smaller than other sources of uncertainty. We find that 60% of our YSOs are in Class II, with 16% in Class I, 12% in Flat, and 12% in Class III. The results with classes defined by α' (after extinction corrections) are also shown in Table 5 for each cloud and summarized in Table 6. The main effect of extinction corrections is to decrease the number of Class I (to 14%) and Flat (to 9%), while increasing the number of Class II (to 64%) and Class III sources (but still at 13%).

With or without extinction corrections, our result is quite unlike that of Wilking et al. (1989) who found roughly equal number of Class I and Class II objects in Ophiuchus. We find many more Class II objects, as also found by Wilking et al. (2005). The variations from cloud to cloud are substantial, as can be judged from Figure 5, though the effects of small number statistics are substantial for Cha II and Lupus (see Table 5). Among the three clouds with substantial numbers of sources in early classes, Perseus stands out as particularly rich in Class I sources.

If star formation has been continuous over a period longer than the age of Class II sources and if we can average all our clouds, we can obtain relative lifetime estimates for each phase by taking the ratios of number counts in each class, and multiplying by the lifetime for Class II. Recall that this method assumes that all objects of all masses behave identically and that no other variables enter importantly. Since these assumptions are not very realistic, caution is needed in interpreting these, or any previous, estimates for lifetimes.

The lifetime of the Class II phase is still uncertain. In the classic study of 83 cTTS in Taurus, Strom et al. (1989) found that 60% of those with ages less than 3 Myr had K-band excesses. The sample of cTTS would probably be biased toward those with long-lived disks. More recently, Haisch et al. (2001) studied young clusters of different ages and showed that half the stars had lost their disks, as indicated by an L-band excess, in ⩽3 Myr. A study of disk frequency in wTTs (Cieza et al. 2007) shows that half have lost their disks within about 1 Myr; since the wTTS sample is undoubtedly biased toward stars that lose disks early, this is probably a lower limit for the disk half-life. The median age of stars in IC348, assuming a constant formation rate and the distance we have adopted, is 3.0 Myr (Muench et al. 2007). The ratio of wTTS to cTTS in IC348 is 1.5, and somewhere between 30% (Lada et al. 2006) and 22% (Cieza & Baliber 2006) of the wTTs have disks. However, many of these would not be Class II objects. These numbers suggest an upper limit of 3 Myr for the Class II phase. For Cha II, Spezzi et al. (2008) find a mean age for the YSOs of about 2 Myr. Both Bontemps et al. (2001) and Kaas et al. (2004) also assumed 2 Myr in their analysis of ISO data. Based on all these pieces of information, it seems that the best estimate for the lifetime of the Class II phase is 2 ± 1 Myr. The uncertainty is dominated by uncertainties in stellar ages, which are uncertain by factors of 2 at these early times (e.g., Haisch et al. 2001; Hillenbrand 2008). Since the marker used to define the transition was that half the stars lacked infrared excess, this lifetime may be best thought of as a half-life rather than a lifetime which is the same for all objects (see also Hillenbrand 2008). In particular, there is some evidence for longer lifetimes for infrared excesses around lower mass stars or brown dwarfs (Allers et al. 2007).

If we take 2 Myr for the duration of the Class II phase, then the lifetime of the Class I phase, t(I) = 0.54 Myr, with an additional t(F) = 0.40 Myr for the Flat SED phase. Using the numbers after extinction corrections, we obtain t'(I) = 0.44 Myr and t'(F) = 0.35 Myr. With or without extinction corrections, these estimates are substantially longer than some recent estimates for the Class I lifetime of 0.1–0.2 Myr (Greene et al. 1994), but this estimate depends on the assumption of an age for the Class II phase. Interestingly, our estimate agrees better with the original estimate of 0.39 Myr for the Class I stage by Wilking et al. (1989) because they assumed an age of 0.4 Myr for the Class II stage in the L1688 cluster and found equal numbers. So we obtain the same answer but with totally different data and assumptions.

The uncertainties in these numbers are clearly not dominated by small number statistics any longer. While the original definition of classes did not include extinction corrections, it appears that an approximate correction leads to 10% and 20% shorter lifetimes for the Flat and Class I phases. The limitations of the underlying assumptions dominate the uncertainties. In particular, the continuous formation assumption is dubious. If star formation in a particular region decreased or stopped within the last 2 Myr, we would underestimate the lifetime for the Class I and Flat classes. Conversely, if star formation began or increased within the last 2 Myr, we would overestimate those lifetimes. Evidence for these effects can be seen in the lifetimes calculated cloud by cloud in Table 5. For Class I, these vary from t(I) = 0.2 Myr (Cha II and Lupus) to t(I) = 0.77 Myr in Perseus, which supplies over half the Class I objects and nearly 40% of all YSOs in our sample. Since Perseus clearly had significant star formation in the IC348 region 3 Myr ago, it does not obviously violate the continuous formation assumption, but there may have been a recent burst. If we take an extreme position and leave Perseus out of the average, 〈t(I)〉 = 0.40 Myr.

Our only defense against these issues is to hope that averaging over all the clouds provides some cancellation of biases. Since we selected our clouds to have a range of known star-formation activity, the sample should not be too heavily biased, but it does contain three of the most active nearby regions. Once the data from the Taurus cloud and the Gould Belt Survey are available in the same form as the c2d data, these numbers should be recomputed.

One possible systematic bias toward overestimating Class I lifetimes would arise if low-mass objects are found in our sample when they are young and accreting but fall into the galaxy confusion region during the Class II stage. We start to lose Class II sources around the stellar/substellar boundary, based on a few examples. If we assume that we miss all substellar objects, which is too extreme, we would miss 10–30% of Class II sources (e.g., Andersen et al. 2006; Luhman et al. 2007), which would introduce a comparable percentage error in the Class I lifetime.

The original definition (Lada & Wilking 1984) used wavelengths between 2 and 20 μm to define α. We use a least-squared fit to any data in Table 8 between 2 and 24 μm. More recently, Lada et al. (2006) have used a fit to only the four IRAC bands (3.6–8.0 μm) in a study of disks in IC348. How much difference does the choice of wavelength range make? We classified sources both ways (Table 6). Using only IRAC bands moves sources from earlier to later classes in such a way that we get more Class III (by 47%), fewer Class I (by 8%), fewer Flat SEDs (by 15%), and fewer Class II (by 8%). Reddening will have a larger effect at shorter wavelengths, which biases α upward, producing "Flat" SEDs that are clearly reddened Class II sources when the full SED is examined. Thus, different definitions of α can result in 10–15% difference in inferred lifetimes for the earlier SED phases. We did another experiment, using sources in Cha II, for which we have spectral types (Spezzi et al. 2008), allowing photometry to be corrected for extinction. Indeed, the values of α calculated after extinction corrections were smaller, with a mean difference of 0.2. This was not enough to affect the classification of sources in Cha II. With the approximate extinction corrections used for the full data set, we find somewhat larger effects, but still 10–20%. These effects are probably smaller than other sources of uncertainty, inherent in assumptions such as the continuous formation assumption (Section 5.2).

Values for L_bol and T_bol were calculated as described in Section 3.3 and are given in Table 12. The observed values of L_bol and T_bol are easily compared to most other data on embedded sources in the literature, but in the paper connecting T_bol to the classes defined by α, Chen et al. (1995) corrected the observed flux densities for extinction before computing T_bol. For less embedded objects with A_V determinations from optical data, they used those values. If a source had no measured A_V but had a known spectral type, they used the spectral type and V − R to calculate A_V. If no spectral type was known, they assumed an M0 spectral type. For sources with no optical data, they assumed an average A_V = 1 for Taurus and off the core region in Ophiuchus, and A_V = 10 for the core sample in Ophiuchus. Our extinction corrections are similar to what they assumed in Ophiuchus and substantially larger in other clouds (see Section 3.3). They analyzed the effects of extinction; these could be substantial in both T_bol and L_bol, primarily for the sources with higher intrinsic T_bol, where most of the energy emerged at shorter wavelengths, and these sources generally had reasonably well determined values for A_V. They considered the uncertainties introduced by extinction uncertainties and found a factor of two for high T_bol sources and 10% for low T_bol sources (Chen et al. 1995). Most subsequent application of T_bol has been to deeply embedded sources, so no corrections for extinction are typically made; the energy absorbed at short wavelengths is assumed to be re-radiated at longer wavelengths.

Our method of extinction corrections attempts to capture the essence of the Chen et al. (1995) approach. If we compare the numbers in classes defined by T'_bol, the effects are more dramatic than was the case for classes defined by α. The observed fluxes produce only one Class III source with the T_bol dividing line from Chen et al. (1995), while using T'_bol give 485 Class III sources! Clearly, the boundary in T_bol between Classes II and III is a sensitive function of extinction.

Class 0 objects cannot be separated from Class I objects by using α. Previously this statement was true because Class 0 objects could not be detected in the mid-infrared (André et al. 2000). With IRAC on Spitzer, we now detect many, though not all (Jørgensen et al. 2005), Class 0 sources. However, the resulting values of α from the c2d catalog scatter widely and do not correlate with other discriminants such as T_bol for the early classes (Enoch et al. 2009). We show the full sample of sources with enough data to constrain both α and T_bol in Figure 6. While there is some tendency for Class 0 sources, as identified by T_bol, to have larger values for α than Class I sources, the scatter precludes any Class 0 criterion based only on α. These conclusions do not change using the photometry after extinction corrections to calculate T_bol (right panel of Figure 6). The range of α may be due in part to the geometry of the envelope–disk interface, the inner radius of the envelope (Jørgensen et al. 2005), different amounts of scattered light, emission from jets in the 4.5 μm band, and deep ice features (see, e.g., Huard et al. 2006; Dunham et al. 2006; Bourke et al. 2006; Boogert et al. 2008). The fluxes do not increase monotonically from 3.6 to 24 μm in some Class 0 sources.

Zoom In Zoom Out Reset image size

Conversely, there was no definition of the Flat SED class in terms of T_bol. Figure 6 shows a strong concentration around 650 K, the Class I/II boundary defined by Chen et al. (1995), with some outliers. Histograms of the number of sources in linear bins of T_bol and T'_bol are shown in Figure 7. The different classes according to α and α' are color-coded. The mean value of T_bol for Flat SED sources is 649 K. We suggest a range of T_bol from 350 to 950 K for Flat SED sources; this range includes 79% of the Flat SED sources and contains only 23% of the Class I sources (by α) and 22% of the Class II/III sources. Enoch et al. (2009) have also suggested a boundary between Stage I and Stage II sources in the 400–500 K range.

Zoom In Zoom Out Reset image size

The distribution of Flat SED sources in extinction-corrected T'_bol is somewhat broader and shifted to higher T'_bol, with a mean value of 844 K (right panel of Figure 6). A range of 500–1450 K contains 77% of the Flat SED sources, while only 14% of Class I and 18% of Class III sources lie in this range. Figure 7 suggests that Flat SED sources do form a distinct class of objects. However, this seeming homogeneity can in part be caused by observing a flat SED over a restricted wavelength range; T_bol has little freedom to vary without data at shorter and longer wavelengths.

Class II and III sources have barely distinct peaks in Figure 7, but substantial overlap. This figure illustrates the sensitivity of the Class II/III boundary to extinction corrections that was noted above. It also suggests that T_bol is not particularly good at distinguishing Class II from Class III objects and α is preferred.

We can also compare the source types in the c2d catalogs, as described in Section 3, to classifications based on α or T_bol. Figure 8 shows the source types versus T_bol, both with and without extinction corrections, with the color code showing the α SED class. The categories such as YSOc_red, red, and rising are almost exclusively Class 0/I or Flat. The source type YSOc_star+dust(IR1) contains a broad range of classes; an excess starting at 3.6 μm can continue to increase (Class I), stay constant (Flat), or fall (Class II or Class III). When the excess begins at longer wavelengths, such as YSOc_star+dust(MP1), almost all are Class III, with a few exceptions. The exceptions may include the "cold disks," such as those studied by Brown et al. (2007). The generic YSOc class contains a heterogeneous sample of classes. Similar patterns are seen in Figure 9, which shows the source type versus α.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Color–color diagrams constructed from IRAC photometry (Allen et al. 2004) or from combined photometry from IRAC and MIPS (Muzerolle et al. 2004) have also proved useful in separating sources in different phases of evolution. These are shown for our YSO data in Figures 10 and 11, both with data after extinction corrections in the right panel. The color codes correspond to the classification by α except that Class 0 and I sources are separated by T_bol. The overall correspondence is good, with Class II sources heavily concentrated in the expected locations in the two diagrams, as indicated by the boxes, taken from Allen et al. (2004) and Muzerolle et al. (2004). Flat SED sources were not separated by Allen et al. (2004) or Muzerolle et al. (2004), but they tend to lie between the Class II and Class 0/I regions (see also Allen et al. 2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Class 0/I sources scatter more widely but generally have redder colors than the Flat SED sources. The Class 0 and I sources are not clearly separated in these diagrams, similar to what we found in Figure 6, but on average, Class 0 sources have redder colors. There are exceptions, such as Class 0 sources with [5.8]–[8.0] colors near zero. Also, there are some Class II sources rather far from their home box, especially in [3.6]–[4.5], many of which find their way home in the figure corrected for extinction.

Our large sample suggests that some adjustments to the box locations and perhaps shapes would be advisable for better segregation of classes in the [3.6]–[5.8] versus [8.0]–[24] diagram. In particular, the divisions between Class II, Flat, and Class 0/I may be captured better by diagonal lines in Figure 11. The bottom two panels of Figure 11 show the divisions between the simulated SEDs for sources in different stages (Robitaille et al. 2006), which do capture the divisions seen in the data well, especially for the de-reddened data. This general congruence indicates that class definitions generally agree with the physical stages, despite exceptions. We consider some of the remaining issues and extend the analysis to classes not separated in Figure 11 (prestellar, 0, and Flat) in the following section.

As discussed previously, distinguishing Class 0 from Class I sources can be ambiguous without very well-sampled SEDs through the submillimeter region. Lacking data at 350 μm for many of our sources, we will use T_bol = 70 K to define the Class boundary (Chen et al. 1995).

To include still earlier phases, prestellar cores, we rely on surveys in the millimeter continuum. We have complete surveys at 1.1 mm with reasonably uniform sensitivity for Perseus, Serpens, and Ophiuchus and a less sensitive survey of Cha II. Focusing on the first three clouds, Enoch et al. (2008) report on a total of 201 cores (defined as sources of millimeter continuum emission above some level). By comparing these cores to the Spitzer data, Enoch et al. (2009) separate these into cores with and without central luminosity sources. The latter are assigned to prestellar cores. Cores with central luminosity sources are divided into Class 0, I, and II using the boundaries in T_bol listed in Section 5.1, based on Chen et al. (1995). There are only a handful of Class II sources associated with dense cores, and these are likely due to projection effects. The results are given in Table 6.

To get lifetimes for these classes, we must bootstrap, using the lifetime derived above for Class I (now re-interpreted as the lifetime for Class 0 and Class I combined). Using the same continuous formation assumptions and a lifetime for Class 0/I of 0.54 Myr, we would infer a mean lifetime for Class 0 of 0.16 Myr for all three clouds combined. If we use the same sample, but use T'_bol, calculated after extinction corrections, the resulting t'(0) = 0.10 Myr. The 60% decrease in lifetime results partly from the shorter lifetime for Class I, but primarily from the smaller number of Class 0 sources after extinction corrections (Table 6). Doing the calculation cloud by cloud, with the values of t(I) for that cloud, yields a wide range, from Ophiuchus with t(0) = 0.043 Myr to Perseus with 0.32 Myr. The Perseus result is consistent with the results of Hatchell et al. (2007). While the Ophiuchus lifetime is 3 times higher than the original estimate of André et al. (1993), it is much shorter than for the other clouds. In addition, it could be overestimated, as it rests on three Class 0 sources, two of which are borderline Class 0/I sources. Either Class 0 sources in Ophiuchus do evolve very quickly into Class I sources, or the continuous formation approximation breaks down there. A decrease in the star-formation rate beginning 0.16 Myr ago would result in very few Class 0 sources now. The variation from cloud to cloud is less if the same value of t(I) of 0.54 Myr is used for all clouds, as in Enoch et al. (2009). Note that these results use the ratio of Class 0 to Class I sources with both required to have millimeter emission.

To get a lifetime for the prestellar (t(pre)) cores that are seen in the Bolocam surveys, Enoch et al. (2008) compared the number of cores that remain starless to the number of embedded protostars. The latter number is slightly less than the number of Class 0/I sources because a small number of Class I objects do not have distinct Bolocam cores. Since this number is small, we still use the lifetime for the combined Class 0 and I phase to represent the lifetime of such protostellar cores. The result, using all cores in all clouds, is t(pre) = 0.46 Myr. In this case, extinction corrections remove enough Class 0/I sources to balance the decrease in Class I lifetime, so the extinction-corrected result is t'(pre) = 0.47 Myr. These numbers differ very slightly from that (0.45 Myr) given by Enoch et al. (2008) because of small differences in the source list. One must immediately caution that this is the lifetime for cores remaining starless after they become dense enough to be detected with Bolocam. Based on careful analysis of the instrumental limitations, Enoch et al. (2007) concluded that the minimum mean density is about 2 × 10⁴ cm⁻³. The variation in t(pre) from cloud to cloud is 0.23 (Serpens) to 0.78 (Perseus) Myr if we use the t(I) for each cloud; Ophiuchus is not unusual in this category with a value (0.37 Myr) near the mean. However, the total mass of the cores in Ophiuchus is small; it may be a cloud that has just finished a burst of star formation that has exhausted much of its dense gas. Again, the range is less if the same t(I) is used for all clouds, as in Enoch et al. (2008).

Comparison to a separate study using SCUBA data can provide some idea of the uncertainties in these estimates. Jørgensen et al. (2007) compared SCUBA data in Perseus to c2d infrared data, and Jørgensen et al. (2008) have extended this study to Ophiuchus. Compared to the Bolocam data, the SCUBA data are taken at a slightly shorter wavelength (0.85 mm vs. 1.1 mm), have a smaller beam (15'' vs. 30''), and are more sensitive to low mass cores (∼0.5M_☉ vs. 0.8 M_☉ for Perseus), but they do not cover the clouds as completely. The sensitivity to density is similar, but they may be probing objects of somewhat lower mean density. Other differences exist in detail in the methods of selecting sources, so comparison to the Bolocam results will give a good idea of the current systematic uncertainties. Jørgensen et al. (2008) find more starless cores in Ophiuchus and fewer in Perseus than does Enoch et al. (2008). Calculating the lifetimes in the same way, we find a prestellar core lifetime of 0.55 Myr for Perseus and 0.75 Myr for Ophiuchus. Taking the two clouds together yields a value of 0.61 Myr, somewhat larger than the value of 0.46 Myr, obtained by Enoch et al. (2009).

Considering the two studies, a reasonable estimate for the lifetime of prestellar cores above a mean density of about 2 × 10⁴ cm⁻³ would be about 0.5 ± 0.1(±0.3) Myr, where the uncertainty in parentheses includes the range over different clouds. These results are in general agreement with a third analysis of Perseus by Hatchell et al. (2007).

In comparing our results to previous results, summarized in Figure 2 of Ward-Thompson et al. (2007), our lifetime would appear to be lower than what they found if plotted at our minimum density of n = 2 × 10⁴ cm⁻³. However, most cores are denser than the minimum value, and Enoch et al. (2008) have done a more careful analysis, dividing the prestellar cores into two bins of mean density greater than or less than n = 2 × 10⁵ cm⁻³. The inferred lifetimes of 0.32 Myr and 0.11 Myr are quite consistent with values found by Kirk et al. (2005) and with the overall relation found by Ward-Thompson et al. (2007) indicating lifetimes decreasing toward, but not yet reaching, the free-fall time (t_ff) as the density increases, as originally suggested by Jessop & Ward-Thompson (2000). At this point, we should recall the assumptions required for these estimates. In particular, we assume no dependence on core mass in deriving a lifetime. Clark et al. (2007) suggested that lower mass cores collapse more quickly, while Hatchell & Fuller (2008) note that, to the contrary, a higher fraction of the more massive cores in Perseus contain embedded objects, suggesting faster evolution for more massive cores. The core mass distribution for Bolocam cores with embedded sources is broader and flatter than the prestellar core distribution as well (Enoch et al. 2008), also suggesting faster evolution for the more massive cores, though other explanations are possible. The issues connecting the core mass function and the IMF of stars are explored in a number of papers (Enoch et al. 2008; Hatchell & Fuller 2008; Clark et al. 2007; Goodwin et al. 2008; André et al. 2007; Swift & Williams 2008).

Note that we have referred to these starless cores as prestellar cores because almost all with suitable molecular line data were found to be gravitationally bound (Enoch et al. 2008) and hence worthy of the "prestellar" moniker, as defined by di Francesco et al. (2007) and Ward-Thompson et al. (2007).

With the large sample of embedded objects now available, Enoch et al. (2009) have subdivided the Class 0 and Class I SEDs into early and late versions. While the SED evolution appears to be continuous, suitable choices of T_bol can produce distinctive SEDs in the mean. By counting the number in each bin, they found a variation in T_bol that increases with time: T_bol = 25K + C(t/10⁵yr)^1.8, with C = 51 K, though the evolution must slow down to match the Class II lifetime.

The boundary between Classes I and II is also problematic. As discussed in Section 5.1, we will use the terminology of "stage" to describe the physical configuration and "class" to describe the SED. Models with two-dimensional radiative transfer by Robitaille et al. (2006) and Crapsi et al. (2008) show that Stage II objects with large inclination angles can have Class I SEDs. Stage I and II objects occupy largely distinct regions of the various color–color diagrams but there is also significant overlap in some parts of these diagrams (Robitaille et al. 2007). Some Stage I objects viewed relatively face-on will be classified as Class II objects using popular color–color diagrams such as Figure 10. Longer wavelengths are helpful in distinguishing these situations; in particular, extended millimeter continuum emission or tracers of dense gas, peaked on the source in either case, unambiguously indicate the presence of an envelope.

Based on his grid of models, Crapsi et al. (2008) argues that 34% of Stage II (star+disk) sources are misclassified as earlier types because of inclination effects. Potentially, all Flat spectrum sources could be nearly edge-on Stage II objects (Crapsi et al. 2008). If correct, the effect on lifetimes of physical stages could be substantial; further comparison of the properties of Flat spectrum sources to known edge-on disks is needed to evaluate this suggestion. In addition, Crapsi et al. (2008) argue that half of the objects with Class I SEDs could be part of the edge-on Stage II population. If correct, the lifetime for envelope dissipation would be less than half that of the combined lifetime for Class 0 and Class I derived above. Crapsi et al. (2008) considered a number of methods for determining the SED class, including α, T_bol, and the absorption or emission in the silicate feature; all were subject to confusion by inclination effects.

The most reliable indicator of how much envelope remains is the ratio of single dish to interferometer flux density at millimeter wavelengths; these probe the envelope and disk, respectively (Crapsi et al. 2008). When such observations become available for the large number of Class I sources in our sample, Stage I lifetimes can be refined. For the three clouds with complete Bolocam data, we can already determine what fraction of objects we classify as Class I based on α can be associated with millimeter continuum emission. Overall, 65% can be associated, larger than predicted by Crapsi et al. (2008). If we assign all Flat SED class and 33% of the Class I objects to Stage II, the lifetime of Stage I would be 0.28 Myr. Broken down by cloud, the percentage that lie within one FWHM of a dense core is 69% for Perseus, 75% for Serpens, and 49% for Ophiuchus. Again, these numbers differ slightly from those in Enoch et al. (2008), owing to slight differences in Class I source definition. A random distribution of infrared sources would produce about five more-evolved objects projected by chance on a millimeter core for all three clouds taken together (Enoch et al. 2008). A similar comparison for the SCUBA data finds that 66% and 53% of Class I sources in Perseus and Ophiuchus lie within 30'' of the center of a SCUBA core, and one false identification is predicted among the 66 SCUBA cores in Ophiuchus (Jørgensen et al. 2008).

The large fraction of the Class I sources without 1 mm emission in Ophiuchus may result from the generally high extinction in front of the YSOs in that cloud. As discussed earlier, high extinction will increase the value of α. This problem has been explored by van Kempen et al. (2009), who mapped objects in Ophiuchus that had been previously classified as Class I using HCO⁺ J = 4 → 3 and C¹⁸O J = 3 → 2 emission. They found that half (11 of 22) of the objects assigned to Class I by their c2d value of α had no HCO⁺ emission peaked on the source. They found a similar result (6 of 13) for the Flat SED sources. They reclassified those objects as disks. The very high extinction toward many sources in Ophiuchus had increased the value of α derived from observations, by amounts as high as 1.1. As can be seen from Table 5, extinction corrections decrease the number of Class I sources in Ophiuchus from 35 to 27, using α. There is a smaller decrease for Flat sources, from 47 to 44. The core of Ophiuchus has even higher extinction than our assumed value of A_V = 9.76 mag. Calculating separately the extinction toward Class II sources inside the contour of A_V = 15 in the extinction maps, the mean A_V is 15.7 mag inside the contour and 5.7 mag outside. Using extinction corrections based on those numbers for Class I and Flat SED sources inside or outside that contour, the number of Class I sources is unchanged, but the number of Flat SED sources drops to 28. Thus, very high extinction may substantially inflate the number of Flat (by 68%) SED and Class I (by 30%) sources as classified by observed α. Ophiuchus seems to be exceptional in this regard, as can be seen in Figure 12. There are many Flat SED sources in regions of very high extinction, as measured by the cloud extinction maps. These would de-redden to Class II sources.

Zoom In Zoom Out Reset image size

The objects classified as Class I by α but lacking dense gas tracers tended to have characteristic SEDs which peaked at 24 μm and fell toward longer wavelengths. Similarly, Class I sources without millimeter continuum emission tended to have such SEDs (Enoch et al. 2009; Dunham et al. 2008). This characteristic SED may be a way to identify suspicious Class I sources. Another approach to sorting out these borderline sources may be studies of veiling of absorption lines (Casali & Matthews 1992); a few sources in Ophiuchus without evidence of dense gas do have veilings that are generally associated with Class I objects (Doppmann et al. 2005).

Finally, the simple division between Class II and Class III sources may be obsolete. First, note that there is no minimum value for α in the Greene et al. (1994) scheme, as stars with no infrared excess, but other evidence that they are pre-main-sequence, can be included in Class III. Values as low as α = −3 are possible for the Rayleigh–Jeans tail of a hot star. For most of our clouds, we lack a complete sample of such objects, because we require a significant infrared excess to identify a source as a YSO, thus effectively setting a lower limit on α. Thus, our counts of Class III sources, on this definition, are incomplete.

A second point is that the classical idea that disk evolution leads to a more negative value of α is based on all the excess fluxes decreasing roughly together. In fact, we see a wide variety of SED shapes around the Class II/III boundary, leading us to suggest a new, two-dimensional classification system.

As discussed in Cieza et al. (2007) and in more detail by B. Merín et al. (2008, in preparation), the diversity of SEDs for star and disk systems can be captured better in a two-dimensional system. Many stars show no excess below a certain wavelength. We define λ_turn-off as the last wavelength with no detectable excess over a photosphere, as fitted to the photometry at shorter wavelengths. Then we define α_excess as the slope of the SED beyond λ_turn-off. Thus, α_excess characterizes the properties of the excess emission, rather than mixing the stellar and disk properties.

As shown by Cieza et al. (2007) and B. Merín et al. (2008, in preparation), classical T Tauri stars have excesses that start short (λ_turn-off ⩽ 2 μm) and a tight distribution of α_excess, around −1. In contrast, weak-line T Tauri stars with excesses usually have λ_turn-off ⩾ 2 μm, with a wider range of α_excess. The spread in α_excessvalues grows with increasing λ_turn-off until it covers the range from −3 to 1, similar to the range seen for debris disks (Chen et al. 2005).

The increase in λ_turn-off values may be caused by several things, including unresolved binaries, either stellar or substellar in mass, and inside-out photo-evaporation. The c2d sample includes a substantial number of "cold" disks, with large λ_turn-off and positive α_excess, but they remain a small fraction of the sample, suggesting that the duration of this phase is short. These issues will be explored by B. Merín et al. (2008, in preparation).

The five clouds in this study span a small range in level of turbulence, as crudely probed by the mean line width of the ¹³CO J = 1 → 0 line. We have searched for correlations between line width and cloud mass, star-formation rate, depletion time, and the lifetime for Class I sources. None are apparent, probably because the range of turbulence over our cloud sample is small.

Schmeja et al. (2008) have used an earlier c2d delivery to count and classify YSOs, and they compared the nature of the clusters in Perseus, Serpens, and Ophiuchus. They concluded that the most centrally condensed cluster (Serpens) was associated with the lowest turbulence, as measured by the NH₃ emission lines from the region of the cluster. In contrast, Table 1 indicates that Serpens has the largest line width averaged over the cloud, as measured in ¹³CO J = 1 → 0. Clearly, the characterization of the turbulence can change conclusions. In addition, the earlier delivery was considerably more contaminated with background galaxies. Inclusion of background galaxies would tend to make earlier classes look less tightly clustered.

Some pictures of star formation invoke a very dynamic picture with protostars or dense cores competing for gas (e.g., Bonnell et al. 2004). We can assess the relevance of such pictures to the regions we study by computing the timescale for collisions between dense cores from t_coll = (n_coreπr²v)⁻¹. For the Serpens cluster A (Harvey et al. 2007a), the data from Enoch et al. (2007) on n_core and r, combined with a velocity dispersion among cores of 1 km s⁻¹ (Williams & Myers 2000) yields t_coll = 3.1 Myr, 6 times the Class I lifetime and 3 times the combined Class I and prestellar core lifetimes spent above a density of about 2 × 10⁴ cm⁻³. For this cluster, we conclude that collisions of dense cores at the stage that we observe them are not dominant. Rather the cores we observe are more likely to evolve individually to form stars. André et al. (2007) reached a similar conclusion for the main cluster in Ophiuchus, and even for very compact subclusters within the main cluster. Also, Walsh et al. (2004) found very small (⩽0.1 km s⁻¹) motions of dense cores relative to surrounding less dense gas in a larger survey of dense cores. We caution, however, that the answer may be very different in regions forming still richer clusters, such as Orion.

We can make some comparisons of the lifetimes derived from the relative numbers of objects in different classes or stages in Section 5.3 to models of star formation. As noted by Enoch et al. (2008), the timescale for prestellar cores to begin forming stars, once their mean densities are sufficient for detection by Bolocam surveys, is a few times the free-fall time at those densities, consistent with the earlier work summarized by Ward-Thompson et al. (2007). Such lifetimes are broadly consistent with the evolution of magnetically controlled models (Tassis & Mouschovias 2005; Adams & Shu 2007) once the density exceeds the threshold for detection by Bolocam. The figure from Ward-Thompson et al. (2007) that summarizes lifetimes as a function of mean density suggests that they get shorter, relative to a free-fall time, as the density increases, also broadly consistent with magnetic models.

The lifetime for the Class 0 stage averaged over all clouds is about 0.16 Myr, but only 0.043 Myr in Ophiuchus. The very short lifetime often quoted for the Class 0 phase (0.01 Myr) was originally derived from observations in Ophiuchus (André & Montmerle 1994), but this now seems anomalous. If star formation in Ophiuchus has declined from rates that were higher about 0.1 Myr ago, the number of Class 0 sources would be much lower compared to Class I sources. Such a scenario is broadly consistent with the fact that Ophiuchus also has the highest ratio of mass in YSOs to mass in dense gas of the three clouds with Bolocam data (Table 4). Another factor that may have contributed somewhat to the short Class 0 lifetime in Ophiuchus is an overestimate of the number of Class I objects (Section 7.2) because of very high extinction in the Ophiuchus core.

The longer lifetime for the Class 0 phase in most clouds is more consistent with the luminosity distribution shown in Figure 13 and Enoch et al. (2009). In particular, there is no clear evidence in our data for an early phase of very rapid mass accretion, which would produce luminosities higher than the mean values we observe. In fact, the luminosities seen in Class 0 sources are mostly similar to or less than predictions (Young & Evans 2005) of the evolution in L_bol–T_bol space based on simple inside-out collapse models (Shu 1977), with constant accretion rates and unit efficiency. Models with mean accretion rates onto the star faster than the Shu (1977) rate are not supported by these observations.

For Class I sources, most of the observed L_bol fall well below the predictions by Young & Evans (2005), suggesting that accretion occurs at a rate even lower than that of the Shu picture, which is the smallest among star-formation models. The mass in the envelope, on the other hand, decreases about as predicted by this model (Enoch et al. 2009). This result suggests that infall occurs at about the rate expected in the Shu picture, but that accretion onto the star occurs at a much lower rate most of the time. The lower rate may be partially explained if only a fraction of the available mass in a dense core winds up on the star. As discussed earlier, comparison of the core mass function and the IMF suggest ≈ 0.3. This helps explain the low luminosities, but it is insufficient to explain the very low values that we see. In addition, the spread in L_bol of 2–3 orders of magnitude at a given T_bol range practically demand a picture in which accretion onto the star is not steady. While no Class 0 sources have L_bol < 0.1 L_☉ in Figure 13, this can be misleading; heating of a substantial envelope just by the interstellar radiation field can produce an apparent luminosity at this level (Evans et al. 2001). The effect decreases as the envelope dissipates, but is still present during the Class I phase.

The distribution of luminosities for all Class 0 or I sources with millimeter emission is shown in Figure 14. The panels on the bottom cover 1–20 L_☉; there are three sources more luminous, with the most extreme at 70 L_☉. However, the vast majority of sources have L_bol ⩽ 3 L_☉. As shown in the top panels, the trend of increasing numbers at lower luminosity continues down to 0.1 L_☉, where selection effects come into play. A more thorough study (Dunham et al. 2008) of the low-luminosity end that dealt with these effects found 15 objects with internal luminosity, L_int ⩽ 0.1 L_☉. Dunham et al. (2008) was able to rule out a continued rise in the luminosity distribution to lower luminosities, which could have called into question whether any cores were truly lacking in internal sources. In agreement with an earlier analysis (Kirk et al. 2007), Dunham et al. (2008) found that 75–85% of cores thought to be starless before Spitzer remain so down to a limit of L_int > 4 × 10^-3(d/140 pc)² L_☉. As discussed by Dunham et al. (2008), the luminosity distribution in Figure 14 is quite inconsistent with constant mass accretion, which would dominate the luminosity in the early stages. For an object at the stellar/brown dwarf boundary of 0.08 M_☉ and a radius of 3 R_☉, the standard accretion rate from a Shu model $(\mbox{$\dot{M}_{\star }$}= 2\mbox{${} \times 10^{-6}$} M_{\odot}\,{\rm yr}^{-1})$ would produce L_int = 1.6  L_☉. For the sources in Figure 14, 59% have L_bol < 1.6L_☉; L_bol generally provides an upper limit to L_int for deeply embedded sources. Our results reinforce and aggravate the "luminosity problem" first noted by Kenyon et al. (1990). As suggested by Kenyon et al. (1990) and further developed by Kenyon et al. (1994) and Kenyon & Hartmann (1995), non-steady accretion is the obvious solution, and it is nearly required by our results. Furthermore, it appears to begin in the very early, Class 0, phase. Such a scenario also invalidates any attempt to infer stellar mass from luminosity in phases where accretion luminosity dominates stellar luminosity.

Zoom In Zoom Out Reset image size

Non-steady accretion could be caused by a magnetic wall that forms in some models, resulting in "spasmodic" accretion as the wall periodically fails (Tassis & Mouschovias 2005). Presumably, turbulence in the infalling envelope could also result in variations in infall rates as well. Simulations of star formation in turbulent clumps show variations in mass accretion rates of about 1 order of magnitude (e.g., Offner et al. 2008). Because almost all the luminosity of accretion is released when material falls onto the star, relatively constant infall can store mass in a disk, producing very little luminosity until a disk instability dumps a lot of mass onto the star, as in FU Orionis events (Kenyon & Hartmann 1995). Statistical evidence for episodicity extends into the Class 0 stage (Dunham et al. 2008) and there is direct evidence in one case for a current accretion rate much lower than a time average rate (Dunham et al. 2006).

Episodic accretion at peak rates higher than predicted by a Shu model would explain the few sources that lie above the tracks from Young & Evans (2005) in Figure 13. Short periods of rapid accretion would presumably also drive strong outflows. If such episodes are more frequent in the Class 0 period, that could explain the higher ratio of CO outflow luminosity to L_bol seen in Class 0 objects (André et al. 2000). Some evidence for that picture might be seen in Figure 13, if one accounts for the lower luminosity predicted in the Class 0 phase for constant accretion models.

We can also apply our derived lifetimes to simple models in a consistency check. Considering a Shu model (Shu 1977), the wave of infall propagates outward at the sound speed until it reaches the outer bound of the region that will collapse. This should occupy roughly half the time spent in the Class I phase, as the other half is occupied by the last shell falling inward. In this oversimplified picture, the feeding zone for a protostar would be given by

$$\begin{equation} r_f = 0.5 a t(I), \end{equation} \tag{ 9 }$$

where a is the sound speed and t(I) is the lifetime for the Class I objects. For a purely thermal sound speed at T_K = 10 K, r_f = 0.096t(I) pc with t(I) in Myr, or 0.052 pc if t(I) = 0.54 Myr, or 0.042 pc using t'(I) = 0.44 Myr. These radii are comparable to the typical size of cores (Enoch et al. 2008), and the modal spacing between YSOs in a study of a large number of clusters is also about 0.05 pc (R. Gutermuth et al. 2008, in preparation). The age we derive for the Class I phase, together with a Shu model, is consistent with the data at this level. However, core sizes are defined more by observational techniques than by evidence of a sharp boundary. In clustered regions, the core size is set by the overlap with nearby cores, as seen in projection. Sizes of cores taken from CLUMPFIND operating on the SCUBA data give somewhat smaller sizes (0.036 pc in Perseus and 0.018 pc in Ophiuchus), suggesting that resolution effects are important (Jørgensen et al. 2008). In much more crowded environments, such as the Orion cluster, the feeding zone could be even smaller.

How much mass would accumulate in this time? For the inside-out collapse model (Shu 1977), the mass infall rate is given by

$$\begin{equation} \dot{M} = m_0 a^3/G, \end{equation} \tag{ 10 }$$

with m₀ = 0.975 and a is the sound speed. As discussed by Shu et al. (1987), the form of this expression is general, but $\dot{M}$ is constant only for certain initial density distributions. For a temperature of 10 K, we have $\dot{M} = 1.6\mbox{${} \times 10^{-6}$} M_{\odot}\,{\rm yr}^{-1}$. With an efficiency of for matter to wind up in the star, $\mbox{$\dot{M}_{\star }$}= \epsilon 1.6\mbox{${} \times 10^{-6}$} M_{\odot}\,{\rm yr}^{-1}$, and the final stellar mass would be M_⋆ = 0.85 M_☉ if t(I) is 0.54 Myr. For = 0.3, suggested by Alves et al. (2007), this produces a star of 0.25 M_☉ (0.20 M_☉ using t'(I) = 0.44 Myr), near the modal mass of the IMF. If the Flat SED objects correspond to the infall of the remnant envelope, one could add another 0.20 M_☉ in the time attributed to that class, but we suspect that many Flat sources do not have any substantial envelope, as discussed above. Previous, shorter lifetimes for the Class I phase would have predicted lower final masses, so it is not trivial to understand how one builds any stars of mass exceeding 0.5 M_☉. One could make more massive stars by increasing or the sound speed, but the latter solution would predict even larger L_bol, and we already have a luminosity problem, illustrated in Figure 13. Models with faster mean infall rates will have to explain the low mean values and broad spread of luminosities in this figure, probably with recourse to spasmodic or episodic accretion.

To test whether or not a picture of episodic mass accretion could result in a typical mass star, we consider a very simple model. The observed L_bol for each source classified as Class 0 or Class I according to T_bol and associated with an envelope (112 sources) is turned into a mass accretion rate, $\dot{M}_{\star }$, by

$$\begin{equation} \mbox{$\dot{M}_{\star }$}= \frac{\mbox{$L_{\rm bol}$}R}{G \mbox{$M_{\star }$}}. \end{equation} \tag{ 11 }$$

Here, R is the stellar radius, assumed to be constant at 3 R_☉, and M_⋆ is the protostellar mass, assumed to be 0.25 M_☉ for this calculation. This value of M_⋆ is chosen because it is half of the mean mass of 0.5 M_☉ assumed in Section 4. In reality, M_⋆ will grow as the protostar accumulates mass. However, if we assume that the average mass accretion rate stays constant throughout the embedded phase, the under- and over-estimates in $\dot{M}_{\star }$ introduced by including sources with masses above and below 0.25 M_☉ should roughly cancel out. The final stellar mass accumulated, M_⋆(t_f), is then

$$\begin{equation} \mbox{$M_{\star }$}(t_f) = \displaystyle \sum _{\mbox{$\dot{M}_{\star }$}} \mbox{$\dot{M}_{\star }$}t(\mbox{$\dot{M}_{\star }$}), \end{equation} \tag{ 12 }$$

where $t(\mbox{$\dot{M}_{\star }$})$ is the amount of time spent in each accretion state and is given by the number of sources at each $\dot{M}_{\star }$ divided by the total number of sources, multiplied by the Class I lifetime of 0.54 yr. With a bin size in $\dot{M}_{\star }$ of 1.0 × 10⁻⁸M_☉ yr⁻¹, we calculate M_⋆(t_f) = 0.7 M_☉. Smaller bin sizes result in identical results, while larger bin sizes overestimate M_⋆(t_f) since very low mass accretion rates are lumped together in one bin. Even though 92 of the 112 sources (82%) have L_bol implying $\mbox{$\dot{M}_{\star }$}\le 1.6 \times 10^{-6} M_{\odot}\,{\rm yr}^{-1}$, the end result is still a typical mass star. In this picture, much of the mass of stars accumulates during brief periods of high accretion. For this specific model and data sample, half the final mass is contributed during the 0.039 Myr (7% of the Class I lifetime) represented by the eight highest luminosity sources. Since the duration of the events with the highest accretion rates is the shortest, even larger samples may be needed to find the rarest accretion events.

Myers (2008) has pointed out that the dense cores do not have sharp boundaries, but instead inhabit "clumps" of lower, but still substantial, density. In this case, the final stellar mass may be determined by a competition between free fall onto the protostar and dispersal of the clump gas. More massive stars could be formed in this way. If the infalling material from the clump has sufficient opacity to make the SED fall into Class I, the time limit on this process would still be the lifetimes derived above, and those would limit the final stellar mass to values similar to those derived above.

To compare more specifically to theoretical models, it is necessary to translate the models into simulated observations. Young & Evans (2005) did that for the standard Shu model (Shu 1977), assuming unit efficiency and spherical radiative transfer. They predicted that the observed Class 0 phase would be short, about 0.05 Myr, terminated by the collapse of the first hydrostatic core as the central source shrank to 3 R_☉ and accretion luminosity became important, and significant emission in the mid-infrared began to be produced. That lifetime was very uncertain because it depended on assumptions about the first hydrostatic core, so the time may be longer; alternatively, episodic accretion could move sources back and forth across the Class 0/I boundary. The Class 0 sources with high L_bol would need to have thicker envelopes or be viewed edge-on through a non-spherical envelope. The transition from Class I to Class II in the models of Young & Evans (2005) depended strongly on the mass of the initial envelope, ranging from 0.1 to 0.4 Myr for masses of 0.3 to 3 M_☉. This transition occurred near, but not at, the end of infall. The time spent in the range of T_bol we might attribute to the Flat stage was very short, reaching 0.1 Myr only for the model with a mass of 3 M_☉. Clearly, the effects of non-spherical models and non-steady accretion will have to be considered to bring observations and models into better congruence.