LUMINOSITY FUNCTIONS OF SPITZER-IDENTIFIED PROTOSTARS IN NINE NEARBY MOLECULAR CLOUDS

Aperture photometry of the IRAC data (Program IDs (PID) 43, 50, 30641, 50070) in Orion was taken from S. T. Megeath et al. (2012, in preparation). The IRAC aperture photometry from Cep OB3 and Mon R2 (PID 20403) was taken from R. A. Gutermuth et al. (2012, in preparation) with no modification. MIPS data for Orion (PID 43, 47, 58, 30641, 50070), and Cep OB3 and Mon R2 (PID 20403) were reduced, calibrated, and mosaicked using the MIPS instrument team's Data Analysis Tool (Gordon et al. 2005). An additional MIPS field in the Cep OB3 cloud containing an extension of the Cep OB3b cluster (PID 40147) was reduced using Cluster Grinder (Gutermuth et al. 2009). PSF-fitting photometry was used for the MIPS 24 μm sources in each of these regions due to the lower angular resolution in the wavelength band, the resulting confusion with other sources, and also nebulosity present at this wavelength. The PSF-fitting photometry was performed directly on the 24 μm data mosaic. The steps of the PSF extraction are described below.

First, MIPS 24 μm sources were identified using the source finding routine in PhotVis (Gutermuth et al. 2008). The photometry was then extracted by PhotVis using an aperture size of 5 pixels and sky annulus between 12 and 15 pixels, with each pixel having a size of 125. Point-like sources were automatically selected by PhotVis, but careful visual inspection of the images added a few sources per image, mostly in nebulous regions. A zero-point magnitude of 16.05 was adopted; this included the aperture correction from 5 pixels to infinity of 1.70. The images were in units of DN s⁻¹ per detector pixel. The sources were detected in PhotVis above a 6σ threshold, then each image was carefully inspected for missing sources.

Next, we created two PSFs using a sample of bright but not saturated stars; these were selected to be relatively isolated and in regions with little nebulosity. Typically 15 sources were used for a given mosaic, but in a subset of mosaics there were fewer than 15 suitable sources available. We used two PSFs: one that would be used to extract photometry from the unsaturated sources and one used for the saturated sources, the latter of which was larger in order to encompass the extended PSF wings of the saturated sources. The larger PSF used for saturated sources has a radius of 35 pixels and a fitting radius of 10 pixels, while the PSF used on unsaturated sources has a radius of 25 pixels and fitting radius of 2 pixels.

We then fit the appropriate PSF to the sources using the IDL implementation of NSTAR from the IDL astronomy library (Landsman 1993). To minimize the effect of crowding on the photometry, the sources were first grouped together in clusters of potentially overlapping source PSFs using the routine GROUP. This routine groups together sources that were within 8 pixels of one of the other group members; the sources in the group were then fit simultaneously using NSTAR. Sources that contained pixels with signals greater than the saturation limit were flagged. To fit the saturated sources, we modified NSTAR to ignore saturated pixels, which were flagged by an "NaN" value in the image. Unsaturated sources were fit with the smaller PSF while saturated stars were fit with the PSF with the larger size and fitting radius to ensure that there were enough unsaturated pixels to provide a robust fit. Any unsaturated sources that were grouped together with saturated sources were then treated like saturated sources for PSF fitting and fit using the larger PSF. We then used NSTAR to extract fluxes by first extracting fluxes for grouped unsaturated sources and then for grouped saturated sources. Afterward, the PSFs were scaled by the photometry and subtracted out with the SUBSTAR routine; this process was done twice, once first for the unsaturated and once for the saturated sources. The residual image was then inspected for sources which had not been extracted. These sources, typically one or fewer for each of the mosaics, were added into the photometry. The process of grouping, fitting, and extraction was iterated until we found no more new sources.

Very saturated sources (those with many saturated pixels and for which PhotVis had difficulty estimating fluxes) may not be easily fit by NSTAR. These sources were then given estimated magnitudes and run through NSTAR using these estimations. After flux extraction, a residual image was created by SUBSTAR, and visually inspected for poorly subtracted PSFs. Saturated sources were typically oversubtracted and were apparent in the residual image. The input magnitudes initially from PhotVis aperture photometry or from our initial estimation were manually adjusted until the residual in the image was minimized. Photometric uncertainties for these sources are due to the varied environments of these sources, including areas of crowding or nebulosity, and were then found by determining the change in magnitude needed to create a noticeable over- or undersubtracted residual. This constrains the source magnitude to within typically ±0.25 mag. We fit six protostar candidates in Orion in this way. For some sources, positions needed to be adjusted after NSTAR fitting as well, which was done so that the center of the PSF fell directly over the center of the source. Since a formal least-squares fit and uncertainty could not be obtained, we use the results of this "fit by eye" for these sources. Three protostar candidates in Orion were given positions in this way. The process of grouping, fitting, and extraction was iterated while we found more new sources.

Since each MIPS frame is the subtraction of readouts from the beginning and end of the integration, some very saturated pixels had values below our adopted threshold for saturation. These saturated pixels can be identified by comparing the original image with a "fake" image. We built the "fake" image by first creating an image of the background nebulosity with the stellar sources subtracted. This was done by first running SUBSTAR and smoothing the resulting image to minimize the artifacts from the subtraction, including the reduction of "holes" in the image due to the oversubtraction of saturated sources. The subtracted stars were then added back to the smoothed image using one of our two PSFs scaled to the best-fit photometry. Where the resulting simulated data exceeded the saturation limit we masked out pixels. The entire photometric extraction process was then repeated with the newly identified saturated pixels masked. This process was repeated until no new saturated pixels were found.

The uncertainties listed in Table 1 for Orion, Cep OB3, and Mon R2 are the relative uncertainties of the data taken into account the photon noise and variations of the signal in the sky annulus. The 24 μm uncertainties are those returned by the NSTAR program. For the highly saturated stars that are "fit by eye," the uncertainties are given by the change in the magnitude to create a noticeable change in the residual and are noted in Table 1 for Orion. The reported uncertainties do not include absolute uncertainties due to calibration errors; these uncertainties are 5% for the IRAC and MIPS instruments.

Photometry for Serpens, Perseus, Ophiuchus, Chamaeleon, and Lupus were obtained from the Cores to Disks (c2d) Spitzer Legacy project 4th high-reliability data release.⁷ We used only the fluxes that were not band-filled,⁸ thus requiring that the source was detected in a specific band before adopting the tabulated photometry. These sources were identified in the c2d photometry tables as having an entry of "−2" in the "image type" column in the high-reliability catalog. In addition to the high-reliability catalog, we inspected the 24 μm images of these clouds to search for sources that were not included in the c2d catalog due to saturation. For this purpose, we used the 24 μm image of Serpens from Winston et al. (2007) and c2d images⁷ for the Perseus, Ophiuchus, Lupus, and Chamaeleon II clouds. By comparing source list positions with the images, we identified sources that do not appear in the c2d catalog and appeared saturated in the images. We then attempted to fit these sources using the PSF technique described in Section 2.1. We were successful in obtaining 24 μm fluxes for nine saturated sources across all of the c2d clouds using the modified version of NSTAR. We used a sample of unsaturated sources from the Ophiuchus cloud to compare our photometry extraction method with that of the c2d survey. We found in a comparison of 72 sources, the mean difference between the m₂₄ that we extracted and the m₂₄ from the c2d data release is 0.0061 mag with a standard deviation of 0.12 mag. Similarly, we compared our extracted Serpens m₂₄ with that from the c2d data release, and found the mean difference among 168 sources to be 0.051 mag with a standard deviation of 0.13 mag. These differences are well below the uncertainties in the flux calibration and the uncertainties in the measurements.

The positions and 24 μm photometry for the saturated sources are listed along with the rest of the photometry in Table 1 for Ophiuchus, Taurus, Chamaeleon II, Lupus I, Perseus, and Serpens. In Table 1, the uncertainties listed for these locations for the Two Micron All Sky Survey (2MASS), IRAC, and unsaturated MIPS sources are adopted from the c2d catalog.

We used Taurus photometry from the Spitzer Legacy Taurus Molecular Cloud Survey data release S14.⁹ We add to our sample 13 sources in the L1551 region from Gutermuth et al. (2009); L1551 is a small satellite of the Taurus region with several well-known protostars that was not covered by the Spitzer Legacy map of Taurus (Ungerechts & Thaddeus 1987). L1551 contains some of the most luminous known protostars in Taurus and thus we desired inclusion of these sources in our Taurus photometry sample. For two protostars from Gutermuth et al. (2009) in L1551 and six from the Legacy survey that do not have 24 μm detections due to saturation, we used IRAS fluxes at 25 μm from the Faint Source Catalog except for L1551 IRS5, which is from the IRAS Point Source Catalog. We then convert from 25 μm flux to 24 μm flux using the approximation F₂₄ = K × F₂₅ where K = 0.986 for a cool blackbody (T = 70 K; MIPS Handbook, version 2.0). The sources receiving their 24 μm photometry in this way are L1551 IRS5, IC 2087 IR, 04365 + 2535, HL Tau, GV Tau, 04239 + 2436, 04278 + 2253AB, and 04361 + 2547AB.

In Table 1, we list the photometry we have adopted for the Taurus protostar candidates. For most sources the photometry is from the Taurus Legacy Survey; we note in this table the sources for which we use photometry and adopted uncertainties from Gutermuth et al. (2009), the IRAS Faint Source Catalog, and the IRAS Point Source Catalog.

Taurus is a relatively well-studied cloud and it is important to ensure that known Taurus protostars are included or accounted for in our sample. We compare our list with a list of known protostars from Furlan et al. (2008). We find that all but two of the 28 sources listed in their Table 1 are also in our photometry sample; the remaining sources fall off the Taurus Survey map and L1551 field. We compare the results of our protostar selection criteria with the protostars tabulated in Furlan et al. (2008) in Section 3.1.

The Spitzer Space Telescope, with its ability to rapidly map entire molecular clouds with ⩽5'' resolution and high sensitivity between 3.6 and 24 μm, has dramatically increased the number of known protostars within 1 kpc (e.g., Enoch et al. 2009; Evans et al. 2009; Gutermuth et al. 2009). Protostars can be readily identified by their flat or rising SEDs over the Spitzer wavelength range (Greene et al. 1994). Several schemes for identifying protostars in Spitzer surveys have been developed; these rely either on fitting models to the SED (ranging from simple power laws to detailed radiative models) or on measuring mid-IR colors that depend on the SED slope (Allen et al. 2004; Megeath et al. 2004; Muzerolle et al. 2004; Whitney et al. 2004; Robitaille et al. 2006; Harvey et al. 2007; Winston et al. 2007; Gutermuth et al. 2009). In this paper, we identify protostars using criteria first developed by Megeath et al. (2009) and S. T. Megeath et al. (2012, in preparation) for protostars with 24 μm detections. These criteria are satisfied by SEDs that have a power-law slope (defined as α = dlog (λF_λ)/dlog (λ)) between the 4.5 μm and 24 μm bands greater than α  =  −0.3. This corresponds to the lowest α for flat-spectrum sources in Greene et al. (1994). We sought sources where the emission at 24 μm was dominated by thermal emission from infalling envelopes (Whitney et al. 2003). These include flat-spectrum sources, Class I sources, and Class 0 sources (Lada & Wilking 1984; Calvet et al. 1994; Winston et al. 2007; Enoch et al. 2009).

Before applying the criteria designed to identify protostar candidates, several detection and signal-to-noise criteria were applied to the point sources recovered by our photometry. We began by restricting ourselves to the set of protostars with detection in the MIPS 24 μm band, which detects thermal emission from the envelope. In contrast, the shorter wavelength IRAC bands are typically dominated by light from the inner disk, which is scattered by dust in the envelope (Whitney et al. 2003). In the next section, we will show that the 24 μm band is needed to estimate luminosities of the protostars; thus, likely protostars without 24 μm detections are ignored. Protostars that are located in saturated areas of our 24 μm maps, i.e., the Orion Nebula or NGC 2024, were not included in this survey. To minimize contamination from galaxies, we required that protostar candidates be projected on regions of the clouds where A_V > 3, as determined from extinction maps of the clouds (these maps are discussed at the end of this section). The contamination from galaxies was further reduced by adopting a maximum 24 μm magnitude (Megeath et al. 2009); as discussed later, a different limit was set for each cloud. Finally, in each set of criteria, we required that the utilized photometry had uncertainties in J, H, K_s, 3.6, and 4.5 μm ⩽ 0.1 mag and uncertainties in 5.8 and 8.0 μm ⩽ 0.15 mag, and a detection at 24 μm.

The adopted protostar selection criteria were first described in Megeath et al. (2009) with a few additional modifications taken from S. T. Megeath et al. (2012, in preparation). These selection criteria are initially based on 4.5 and 24 μm photometry alone since these are typically the most sensitive bands for detecting protostars; at 5.8 and 8.0 μm the IRAC detectors are less sensitive and near 8 μm the protostellar SED typically dips (Whitney et al. 2003; Megeath et al. 2009; S. T. Megeath et al. 2012, in preparation). The most deeply embedded sources are sometimes not detected at 3.6 μm; in these cases, we assigned a photometric upper limit magnitude at 3.6 μm of 15.5 mag (a conservative limit equal to the ∼80% completeness level determined from Harvey et al. 2006; a similar limit is used by S. T. Megeath et al. 2012, in preparation) to determine lower limits for the colors for these sources. For sources with uncertainty in m_4.5 ⩽ 0.1 mag and with 24 μm detections, we required protostars to fulfill the following criteria:

$$\begin{eqnarray} && \qquad [4.5] \hbox{--} [24] \ge 4.761 \nonumber \\ &&\qquad \qquad \mbox{ and } \nonumber \\ && [3.6] \hbox{--} [4.5] \ge 0.752 + \sigma _{3.6,4.5} \end{eqnarray} \tag{ 1 }$$

or

$$\begin{eqnarray} && \qquad [4.5] \hbox{--} [24] \ge 5.303 \nonumber \\ && \qquad\qquad \mbox{ and } \nonumber \\ && [3.6] \hbox{--} [4.5] \ge 0.652 + \sigma _{3.6,4.5}, \end{eqnarray} \tag{ 2 }$$

where for any two photometric bands a and b, σ_{a, b} = $\sqrt{\sigma ^{2}_{a} + \sigma ^{2}_{b}}$. The colors of [3.6]–[4.5] = 0.6520 and [4.5]–[24] = 4.761 correspond to a power-law SED with of α = −0.3. Note that we require a minimum color of [4.5]–[24] = 5.303, corresponding to α = 0 for sources where [3.6]–[4.5] = 0.652. This was done to reduce contamination from reddened Class II objects (i.e., pre-main-sequence stars with disks) based on the examination of the color–color diagrams in Figure 1 as well as those of regions with A_V < 3 where few protostars are expected.

Zoom In Zoom Out Reset image size

Examination of the data led to the identification of protostellar objects with highly reddened [4.5]–[24] colors but relatively blue [3.6]–[4.5] colors. These sources showed resolved yet compact scattered light nebulae in I-band images (J. Bally 2008, private communication); indicating that they were protostars where the blue colors were due to a strong contribution from scattered light. Based on the colors of these sources, we established the following criteria to ensure that such sources were included in our protostar catalog. For sources with detections at 4.5, 5.8, and 8 μm and uncertainties of σ_4.5 ⩽ 0.1 mag, σ_5.8 ⩽ 0.15 mag, and σ_8.0 ⩽ 0.15 mag, we applied the following criteria:

$$\begin{eqnarray} & [4.5] \hbox{--} [24] \ge 7\nonumber \\ &\qquad \mbox{ and } \\ & [5.8] \hbox{--} [8.0] \le 1.75,\nonumber \end{eqnarray} \tag{ 3 }$$

where the [5.8]–[8.0] criterion is used to eliminate SF galaxies with strong polycyclic aromatic hydrocarbon emission in the 8 μm band.

Sources without a 4.5 μm detection may still be protostars (see, for example, Fischer et al. 2010). To identify such sources, we used detections in either the 5.8 or 8.0 μm band. Sources with detections in the 5.8 μm band with σ_5.8 ⩽ 0.15 were identified as protostar candidates if they satisfied the color criteria

$$\begin{eqnarray} & [5.8] \hbox{--} [24] \ge 4.117 \nonumber\\ &\qquad \mbox{ and } \\ & [3.6] \hbox{--} [5.8] \ge 1.296 + \sigma _{3.6,4.5},\nonumber \end{eqnarray} \tag{ 4 }$$

while for sources with 8 μm detections and σ_8.0 ⩽ 0.15, we identified protostar candidates with the criteria

$$\begin{eqnarray} & [8.0] \hbox{--} [24] \ge 3.20 \nonumber\\ &\qquad \mbox{ and } \\ & [3.6] \hbox{--} [8.0] \ge 2.12 + \sigma _{3.6,4.5}.\nonumber \end{eqnarray} \tag{ 5 }$$

In both cases, the colors corresponded to those of a power-law SED where α = −0.3.

Additionally, we identified sources with less infrared excess and/or fainter m₂₄ that are likely more evolved YSOs without envelopes using the criteria for Class II sources from Gutermuth et al. (2009). Color–color and color–magnitude diagrams for protostar candidates and identified Class II sources are shown in Figures 1 and 3, respectively. We refer to these objects as protostar candidates because we will find that a fraction of these sources may be galaxies, reddened Class II objects, and edge-on disks. For each protostar candidate, we calculated the spectral slope, α, by using a best-fit line to the plot of log (λF_λ) versus log (λ) over the IRAC and MIPS 24 μm detections. In Figures 1 and 3, we distinguish the flat-spectrum −0.3 ⩽ α ⩽ 0.3 and the rising-spectrum α ⩾ 0.3 protostars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A 24 μm magnitude cutoff was imposed to minimize background galaxy contamination (Megeath et al. 2009). At magnitudes fainter than this cutoff, the contamination rises quickly and surpasses the number of protostar candidates. To determine the optimal value of the cutoff, a histogram of the 24 μm magnitudes for sources satisfying the above color criteria was created for each cloud (Figure 2). Galaxy contamination was estimated by applying our protostar candidate identification criteria to the SWIRE¹⁰ sample assuming that the density of galaxies in the SWIRE field is similar to that toward our sample of molecular clouds. Our sample of molecular clouds is distributed over a broad range in Galactic coordinates from ∼30° < l < ∼350° and ∼−22° < b < ∼20°. Since the m₂₄ cutoff is primarily concerned with extragalactic contamination, which is distributed approximately evenly over the sky, the SWIRE field is appropriate for estimating contamination due to galaxies. Extragalactic wide field surveys such as SWIRE have been used by a number of investigators to estimate galactic contamination (Harvey et al. 2006; Gutermuth et al. 2008). This number of contaminants was scaled to the size of the A_V > 3 areas of the coverage maps for each region, namely, 4.25, 8.88, 0.83, 0.45, 2.46, 0.86, 5.49, 1.37, and 0.32 deg² for Ophiuchus, Taurus, Lupus, Chamaeleon II, Perseus, Serpens, Orion, Cep OB3, and Mon R2, respectively. These are plotted on the same histogram as the candidate protostars toward the clouds (Figure 2). The adopted m₂₄ cutoff is shown for each region. The cutoff is also shown in the color–magnitude diagram shown in Figure 3. In the following analysis, we only include protostars that are brighter than the adopted cutoff magnitude.

For the c2d clouds, we compared our sample of protostars with the sample of 132 c2d sources with envelopes detected at 1.3 mm (Enoch et al. 2009). We classified 88 of these sources as protostars and 14 as Class II sources. Of the remaining 30 c2d envelope sources, we rejected 17 as protostars because of missing or band-filled 24 μm detections or because m₂₄ exceeded the cutoff magnitude, four are rejected for uncertainties above our requirements for protostar selection, three satisfy our criteria but are rejected because they fall off the available A_V maps or do not have an A_V > 3, 3 sources are rejected because they are not consistent with our criteria, and three sources are found only in the full c2d source catalog and not found in the c2d high-reliability catalog. Given the relatively low angular resolution of the Bolocam data (30''), we suggest that some of the envelope sources with Class II colors may be due to chance coincidences between pre-main-sequence stars and dense cloud material. In other cases, the 1.3 mm measurements may be detecting massive disks. Additionally, we identified 52 protostars not identified as likely envelope sources since they do not have 1.3 mm detections of envelopes; we expect these sources to have envelope masses less than 0.5 M_☉ (Evans et al. 2009).

In the Taurus region, we compared the candidate protostars identified using our criteria with the protostars studied in Furlan et al. (2008). We find that of the 26 protostars in Furlan et al. (2008), 18 sources also satisfy our criteria, although one source has A_V < 3. Of the eight remaining Furlan et al. (2008) protostars six are identified as Class II sources, one has a m₂₄ below our cutoff, and one is not detected in m₂₄. We note that Furlan et al. (2008) also find that several of their objects may be sources transitioning between the protostellar and Class II phase.

Figure 4 shows the extinction maps and spatial distributions of all protostar candidates found in our study using our selection criteria as well as the sources that fit the selection criteria but reside in A_V < 3 regions, again separating the flat and rising SED sources. We took the A_V maps for the c2d clouds from the c2d Legacy project 4th high-reliability data release,¹¹ the map of Orion is from S. T. Megeath et al. (2012, in preparation), the maps of Cep OB3 and Mon R2 clouds are from R. A. Gutermuth et al. (2012, in preparation), and the Taurus map is from Lombardi et al. (2010). Also shown in this figure is an A_V > 3 contour on each of the maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (11 panels)

Tarball of figure set images

The majority of protostellar luminosity is emitted at wavelengths longer than 24 μm, wavelengths at which data are not available for most sources in our study. Since the total bolometric luminosity, L_bol, is dominated by the flux at these longer wavelengths, we need a method for estimating L_bol using the luminosities integrated over the available wavelength bands. We present here a method for using the SED slope and the total mid-infrared luminosity in the Spitzer 3.6–24 μm bands, L_MIR, to estimate L_bol. We first discuss the calculation of the SED slope and of L_MIR and then the technique used to estimate bolometric luminosity.

We created SEDs for each of our protostars and used the Spitzer 3.6–24 μm bands to determine the SED slope. We calculated α by using a best-fit line to the plot of log(λF_λ) versus log(λ) over IRAC and MIPS 24 μm detections. Mid-IR luminosities were calculated by integrating the SED over the available fluxes from J, H, K_s, and IRAC bands. Detections were converted to fluxes F_λ, using zero-point fluxes of 1594, 1024, 666.8 Jy for J, H, and K_s (Cohen et al. 2003) and 280.9, 179.7, 115.0, 64.13 (IRAC Instrument Handbook, version 1.0), and 7.17 Jy (MIPS Instrument Handbook, version 2.0) for F_3.6, F_4.5, F_5.8, F_8.0, and F₂₄, respectively. We divide F₂₄ by a color correction of 0.967; for −1 < α < 2 the color correction ranges from 0.960 to 0.967 (Stansberry et al. 2007). The adopted IRAC fluxes are for a flat-spectrum source, the color correction ranges from 1 to 1.03 for sources with α = −1 to α = 2 (IRAC Instrument Handbook, version 1.0). We use the central wavelength for the J, H, and K_s bands and approximate bandwidths from Cohen et al. (2003). We estimated the bandwidths to be 1.073–1.397 μm for F_J, 1.411–1.913 μm for F_H, and 1.897–2.420 μm for $F_{K_s}$. We used estimated bandwidths from the IRAC and MIPS handbooks of 3.175–3.925 μm for F_3.6, 3.986–5.000 μm for F_4.5, 5.019–6.443 μm for F_5.8, 6.420–9.324 μm for F_8.0, and 20.800–26.100 μm for F₂₄. We used rectangular integration over each band by summing the product of the flux and bandwidth from each band and converting to luminosity using the following equation:

$$\begin{eqnarray} L_{{\rm MIR}} &=& [19.79 F_{\lambda }(J) + 16.96 F_{\lambda }(H) + 10.49 F_{\lambda }(K_s)\nonumber\\ && + 5.50 F_{\lambda }(3.6) + 4.68 F_{\lambda }(4.5) + 4.01 F_{\lambda }(5.8) \nonumber \\ & & + 4.31 F_{\lambda }(8.0) + 0.81 F_{\lambda }(24)] \times 10^{-6} \times d^{2} \,[{\rm pc}] \,L_{\odot },\nonumber \\ \end{eqnarray} \tag{ 6 }$$

where d is the distance to the cloud in pc and fluxes F_λ are in Jy. We list the distances in Table 2.

Dunham et al. (2008) compared the protostellar flux, νF_ν, to the internal luminosity, L_int, for a grid of radiative transfer models and protostars with complete SEDs from the c2d survey. For the 70 μm band, they found a linear correspondence between νF_ν and L_bol. An approximately linear relationship was also evident at shorter wavelengths, particularly at 24 μm, but the scatter was much higher, with an order of magnitude variation in luminosity for a given 24 μm flux. Because 70 μm photometry is unavailable for the majority of our sources, we attempt here to reduce the scatter between the mid-IR fluxes and L_bol. Since the conversion factor between the mid-IR luminosity and L_bol must depend on how rapidly the SED rises at wavelengths longer than 24 μm, we examined the dependence of the ratio of L_MIR/L_bol (which will be constant if there is a linear dependence between L_MIR and L_bol) and the slope of the SED between 3.6 and 24 μm.

To establish a relationship between the slope and L_MIR/L_bol, we used protostars selected from the c2d program with well-established L_bol. This sample of c2d identified protostars spans the range of colors and magnitudes characteristic of protostars as shown in Figures 1 and 3. Evans et al. (2009) determined bolometric luminosities for YSOs in the c2d catalog using available photometry between 0.36 μm and 1300 μm and the method described in Dunham et al. (2008). We scaled these luminosities to the distances adopted in our study, listed in Table 2. We used 87 YSOs from Serpens, Perseus, Chamaeleon II, Lupus, and Ophiuchus, which are identified as protostar candidates using our criteria, had L_bol flagged as "good" in Evans et al. (2009), and had 70 μm detections. We then used the photometry from this sample to determine the slope of the SED and L_MIR/L_bol.

We found a relationship between the SED slope and the ratio of L_MIR/L_c2d using 66 rising-spectrum c2d protostars. We chose the slope calculated from 3.6 to 24 μm as it shows a stronger trend than the slope calculated using only 4.5 and 24 μm, which are preferentially used in determining protostar candidate status, or the entire range of available wavelengths for all clouds (J band through 24 μm). The 24 μm flux is useful because it is less affected by inclination effects or the geometry of the outflow cavities, but we found that using L_MIR gives a better fit than using the 24 μm flux alone. We found the smallest residuals using a linear relationship between $\ssty\sqrt{\frac{L_{{\rm MIR}}}{L_{{\it c2d}}}}$ and log (α). We also attempted a polynomial fit to this relationship, but again found that the linear fit has smallest residuals. Thus, we found that $\ssty\sqrt{\frac{L_{{\rm MIR}}}{L_{{\it c2d}}}}$ goes linearly with log (α). We expect some intrinsic scatter in the relationship, which is not due to the uncertainties of the measurements. Figure 5 shows the best fit for this relationship, and a plot relating log (α) and L_c2d/L_MIR. We found that this relationship is

$$\begin{eqnarray} & \frac{L_{{\rm MIR}}}{L_{{\rm bol}}} = (-0.466\pm 0.014 \times \log (\alpha) + 0.337\pm 0.053)^2. \end{eqnarray} \tag{ 7 }$$

The flat-spectrum sources did not follow the same trend. We chose not to fit the flat-spectrum sources (−0.3 < α < 0.3) because there were so few in our sample; instead we assumed that flat-spectrum sources exhibit a constant ratio in bolometric to mid-IR luminosity. For these sources, we adopted the ratio of L_MIR/L_bol given in Equation (7) for α = 0.3:

$$\begin{eqnarray} & \frac{L_{{\rm MIR}}}{L_{{\rm bol}}} = 0.338. \end{eqnarray} \tag{ 8 }$$

We find that the mean value of L_MIR/L_c2d for the 21 flat-spectrum sources not used to create the fit is 0.297 with a standard deviation of 0.11, which is consistent with the adopted ratio of Equation (8).

Zoom In Zoom Out Reset image size

We note that in addition to the internal heating of protostars, the contribution from external heating, L_ext may also affect the measured luminosity. Evans et al. (2001) estimate that the effect of L_ext is typically of order of 0.1 L_☉, which is a small contribution for all but the faintest sources. Therefore, we neglected the effect of L_ext on L_bol.

In Figure 6, we compare the luminosities derived using Equations (7) and (8) with the bolometric luminosities from Evans et al. (2009). The protostars from Evans et al. (2009) used to determine Equation (7) range in bolometric luminosity from 0.03 L_☉ ⩽ L_c2d ⩽ 19.6 L_☉. The left panel of Figure 6 compares the luminosity function constructed from our estimated luminosities to that constructed from the bolometric luminosities in Evans et al. (2009). The luminosity functions are similar; a Kolmogorov–Smirnov (K-S) test comparing the two distribution results in a probability of 0.83 that they are from the same parent distribution. In the right panel, we plot the estimated luminosities versus the luminosities for these sources from Evans et al. (2009). In Figure 7, we show a histogram of the difference between our estimated log (L_bol/L_MIR) and log (L_bol/L_MIR) derived from the Evans et al. luminosities. The standard deviation of the difference is 0.35. From this analysis, we conclude that we can recover the luminosity function of the protostar candidates with reasonable fidelity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can test whether the protostar candidates we identify in Taurus follow a similar relationship. We compared the luminosities determined with model fits by Furlan et al. (2008) with the luminosities determined with our relationship. These sources are shown in Figure 5. We find that the Taurus protostar candidates with well-determined bolometric luminosities agree well with our fit; the K-S probability that our Taurus luminosity function is from the same parent distribution as the luminosity function using the model fit luminosities of Table 1 in Furlan et al. (2008) is 0.93. The assumption we will make for the remainder of the paper is that the protostar candidates in more distant clouds forming massive stars will show the same relationship.

Values of α and L_bol for each protostar candidate are listed in Table 1.

We have used the method above to calculate luminosities for each of the protostars in the nine clouds. In Figure 8, we show the resulting luminosity functions. The protostars of Taurus, Lupus, and Chameleon are again combined because of the similar distances, the dispersed star formation of these clouds, and the low number of protostar candidates in these regions. The individual luminosity functions of these regions and the combined luminosity function are also displayed in Figure 9. We also show the combined luminosity function for the clouds forming high-mass stars (Orion, Cep OB3, and Mon R2) and for the clouds forming low- to intermediate-mass stars (Perseus, Ophiuchus, Taurus, Lupus, and Chameleon).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Galaxies can have colors very similar to those of YSOs (Harvey et al. 2006; Gutermuth et al. 2009). Although we have minimized contamination from galaxies by applying the m₂₄ cutoff to each of our clouds (Section 3.1), a small number of extragalactic sources brighter than this limit is expected. To quantify the contamination from the remaining galaxies on our luminosity functions, we took a sample of SWIRE galaxies from the Elias N2 region scaled to the angular coverage of each cloud, and subjected them to our protostar selection criteria including the m₂₄ cutoff. Any of these known galaxies which were identified as protostars are considered to be contaminants. Using Equations (7) and (8), we determined faux bolometric luminosities for the galaxy contaminants (i.e., the luminosities they would have if they were protostars in the observed cloud) and created luminosity functions for the contaminants. These are shown in Figure 8. Although present, galaxies are the smallest source of contamination, comprising only 2.5% of our protostar candidate sample.

Pre-main-sequence stars seen through their flared disk may have a rising SED and can be misidentified as protostar candidates using our criteria. The fraction of Class IIs seen through their disk is difficult to estimate theoretically, since it depends on poorly constrained properties of the disks, including the amount of flaring and the outer radius of the disk. Instead, we employed an empirical estimate using the technique of Gutermuth et al. (2009). We first identified a large cluster of young stars in a region where the gas has been dispersed and the extinction is low. In our survey, the best example is the Cep OB3b cluster (Allen et al. 2012). Although the gas has been cleared by the OB stars in the cluster, objects with protostellar-like colors were detected in this low extinction cavity. Following Gutermuth et al. (2009), we assumed that YSOs with protostellar-like colors in the low extinction regions were edge-on or nearly edge-on disks (hereafter we use "edge-on disks" to refer to disks that are close enough to an edge-on inclination that they are observed through their disks). We then calculated the ratio of protostars to Class II objects and multiplied this ratio by the number of Class II objects in each cloud to calculate the number of expected edge-on sources. Since some of the sources in Cep OB3b may be actual protostars that have survived gas dispersal, this assumption gave us an upper limit to the number of edge-on disk sources that have colors similar to protostar candidates. In the following analysis, we will set the number of edge-on disks equal to this number.

We identified a region within the Cep OB3b cluster where the total extinction is A_V < 3; this low extinction region contained 34 protostar candidates and 568 Class II sources and is shown in the map of Cep OB3 in Figure 4. The number of Class IIs, however, was for the entire range of magnitudes and did not include a cutoff in m₂₄. Furthermore, these sources were not corrected for possible extragalactic contamination. To calculate the appropriate ratio, R, which imposed a cutoff at m₂₄, we calculated the ratio of edge-on disks to Class II objects using the following equation:

$$\begin{eqnarray} & R = \frac{N_{{\rm Cep \;OB3b}}^{p}-N_{{\rm gal}}^{p}}{N_{{\rm Cep \;OB3b}}^{d}-N_{{\rm gal}}^{d}} , \end{eqnarray} \tag{ 9 }$$

where N^p_Cep OB3b is the number of sources with protostellar-like colors in Cep OB3b with m₂₄ < cutoff, N^d_Cep OB3b is the number of Class II sources with m₂₄ brighter than the cutoff, and N^p_gal and N^d_gal are the expected contamination from galaxies with protostellar-like and disk-like colors, respectively. Although R is an upper limit, we treat R as the actual fraction of edge-on disks in the remainder of this paper and thus we may overestimate the contamination by edge-on disks. The number of edge-on disks in the A_V > 3 region was calculated as

$$\begin{eqnarray} & N_{{\rm edge}} = R \times \big[N_{p}+N_{{\rm CII}}-N_{{\rm gal}}^{d}-N_{{\rm gal}}^{p}\big], \end{eqnarray} \tag{ 10 }$$

where R is calculated for Cep OB3b using the m₂₄ cutoff of that cloud; but N_CII and N_p are the number of Class II sources and protostars brighter than the m₂₄ cutoff determined for the A_V > 3 region of the cloud. We included the number of protostars (N_p) since a fraction of our protostars may have been edge-on disks. This overestimates slightly the number of edge-on disks; however, not including N_p would have resulted in slight underestimation. The edge-on disk sample was corrected for galaxy contamination by removing edge-on disk candidates that have luminosity within 0.2 log (L) of a source from the galaxy sample from Section 4.1 scaled to the size of the Cep OB3b A_V < 3 region.

To determine the number of Class II sources (N_CII), we used the criteria from Gutermuth et al. (2009). We counted all Class II sources with de-reddened m₂₄ brighter than the cutoff. The relatively large size of the Taurus map (44 deg²) along with the proximity of the region resulted in more contamination to the Class II sample. To minimize the contamination, we modified the criteria for this cloud as well as the Lupus and Chamaeleon clouds, which were combined with the Taurus region in our analysis. The active galactic nucleus (AGN) identification criteria for Taurus, Lupus, and Chamaeleon were taken from the criteria in the Appendix of Gutermuth et al. (2008) and adjusted by lowering the m_4.5 threshold magnitude by 2.5 mag.

For each of the candidate protostars in the A_V < 3 region of Cep OB3b (i.e., our sample of likely edge-on disks), we estimated a faux bolometric luminosity using Equation (7). Then, for each expected edge-on disk, we randomly selected one of the luminosities from this sample of likely edge-on disks. We repeated this N_edge times until we have constructed the luminosity function of edge-on disks.

Highly reddened Class II YSOs can have colors and 24 μm magnitudes that satisfy the classification criteria for protostars. To find the amount of contamination due to these sources and estimate the luminosities that would be derived for these contaminants, we used a Monte Carlo simulation. This simulation randomly applied a realistic distribution of extinctions to the observed colors and magnitudes of a fiducial sample of Class II objects and then selects the objects that have artificially reddened photometry that fits the protostar selection criteria. Since the reddening applied to each source is determined from the probability distribution of extinctions, we repeat this simulation 1000 times to sample the range of likely contamination.

The fiducial sample was constructed from Class II objects identified in low A_V regions. We first identified Class II objects in the low extinction (A_V < 3) region of our extinction map. The selected Class II sources were then de-reddened using the reddening law from Flaherty et al. (2007). Sources with de-reddened m₂₄ brighter than the cutoff magnitude were the selected for the fiducial sample. We used these objects as a fiducial sample of Class II objects, with the accompanying assumption that the colors and magnitudes of this sample are representative for all Class II objects in the cloud.

We then identified the full sample of Class II objects and protostar candidates in the A_V > 3 region. We included protostar candidates since a fraction of the protostar candidates may be the reddened Class II sources. Furthermore, protostar candidates are in more highly reddened locations and failure to include them would bias our distribution of A_V to lower values. We sorted the Class II objects and protostar candidates (hereafter YSOs) into regions of higher and lower stellar density by using the distance to the 4th nearest YSO neighbor: sources with 4th nearest-neighbor distances less than the median value of the protostar candidate sample in a particular cloud were considered to exist in regions of higher stellar density ("high stellar density") and those with distances greater than the median value were considered to exist in regions of lower stellar density ("low stellar density"). We then executed the following analysis for the high and low stellar density sources independently. This was done because high stellar density sources may be in regions of systematically higher gas column density (and reddening) compared to sources in regions of lower stellar density.

We extracted the A_V values coincident with each YSO using the A_V maps; this gave us two A_V distributions, one for the high-stellar density sources, and one for the low stellar density sources. These were the maximum A_V values, A_V(max) at the position of a particular YSO. Since the YSOs are embedded in the cloud, the A_V value to that YSO is between 0 and A_V(max). For each YSO with a de-reddened m₂₄ brighter than the cutoff for that cloud, we randomly selected an A_V(max) value and set the A_V to a value drawn from a uniform distribution of extinctions between 0 and A_V(max). We then randomly selected a Class II object from our fiducial sample and applied the A_V using the reddening law from Flaherty et al. (2007). To estimate the fraction of (high or low stellar density) reddened disks masquerading as protostar candidates, we then tested to see how many of the reddened Class IIs were identified as protostar candidates using our criteria. This process was repeated for 1000 iterations for the high and low stellar density sources independently. In each of the 1000 iterations, we estimated the faux bolometric luminosities of the contaminants using the same technique used for the protostar candidates using Equations (7) and (8). For each cloud, we were left with two distributions of contaminant luminosities; one for the sources in regions of high stellar density and one for the sources in regions of low stellar density.

Once we estimated the luminosities for the contaminating sources, we sought to remove protostar candidates with similar luminosities from our sample. Instead of binning the data, and subtracting the luminosity function of the contaminants from that of the candidates (see Figure 8), we eliminated individual protostar candidates from the sample using the following method. This method allows us to create contamination-subtracted cumulative distributions of the luminosities without binning.

First, we divided the sources in each cloud into high and low stellar density protostar candidates using the method described in Section 4.3. We note that this division between high and low stellar density YSOs made a difference for the reddened disk contamination, since the high stellar density regions have systematically higher extinctions, and for the galaxy contamination because likely background galaxies are assumed to be found among the lower stellar density sources, which cover a larger region of the sky.

For each high and low stellar density sample, we generated 1000 trials of estimated contaminant luminosities for each type of contaminant: galaxies, edge-on disks, and reddened Class II sources. In each of the 1000 trials, the protostar candidate with the closest luminosity to each contaminant luminosity, such that the difference was δlog (L) < 0.2, was flagged as contamination and removed from the (high stellar density or low stellar density) protostar candidate sample. If there was no protostar candidate that satisfied these criteria, no source was removed. The final product is 1000 realizations of the contamination-subtracted protostar sample.

These realizations will be used to generate luminosity functions and do statistical analyses in Section 5.1. In each realization, different candidate protostars are removed as contamination. Thus, each of the 1000 realizations will be used in the analyses that follow to take into account the uncertainty in the contamination removal.

For each of the clouds we created contamination-subtracted luminosity functions, which were generated from the 1000 realizations of the protostar candidate sample. After each of the 1000 iterations, we determined the average number of remaining protostars per bin to generate the luminosity functions for the high and low stellar density protostars; these in turn were combined to generate one protostellar luminosity function per cloud. We found a luminosity cutoff, L_cut, which is the m₂₄ cutoff magnitude translated to luminosity. We calculated L_cut by assuming an SED slope steeper (redder) than 90% of all sources in the region and a mid-IR luminosity equal to the 24 μm luminosity at the cutoff magnitude and the relationship in Equation (7). We note that the luminosity cutoff is found from the m₂₄ cutoff we used to reduce the number of galaxies (Figures 2 and 3) and is above the sensitivity cutoff. Thus, we expect our samples to be complete down to the cutoff, except in regions with very bright nebulosity. A more detailed description of incompleteness is given in the Appendix and in Section 5.5.

The protostellar luminosity functions and L_cut are shown in Figure 9. We also show the combined luminosity functions for the clouds with high-mass star formation (Orion, Cep OB3, and Mon R2, hereafter the high-mass SF clouds) and for the clouds forming low- to intermediate-mass stars (Serpens, Perseus, Ophiuchus, Taurus, Lupus I and III, and Chameleon II, hereafter the low-mass SF clouds). Properties of the luminosity functions, including the number of sources comprising each luminosity function, the peak, mean, median, 1σ values, and L_cut, are listed for each cloud as well as the combined low-mass SF clouds and high-mass SF clouds in Table 3.

The protostellar luminosity functions of the high-mass SF clouds have peaks near 1 L_☉ and tails extending toward higher luminosities upward of 100 L_☉. A combined luminosity function for these regions shows a similar peak and tail. The median luminosity in each of the high-mass SF clouds is ∼1 L_☉. Cep OB3 has the highest mean luminosity of the high-mass SF clouds, at 12.06 L_☉. In Orion we expect that we are missing some of the most luminous sources in the saturated regions of the Orion nebula and NGC 2024. In all three massive clouds the most massive objects are missed due to spatially extended regions of saturation.

The low-mass SF cloud luminosity functions do not show a consistent trend. Perseus does not show a peak in the luminosity function above L_cut, but instead rises toward lower luminosities down to L_cut. Ophiuchus shows a marginal peak in its luminosity function, but with small number statistics this peak is not significant. The luminosity function of Tau/Lup/Cha peaks near 1 L_☉, similar to the high-mass SF clouds. The Serpens luminosity function peaks at the highest luminosity, near 2.60 L_☉. The median protostar luminosity for most of the low-mass SF clouds is below 1 L_☉ except in Serpens, which has a median luminosity of 3.07 L_☉. The mean protostar luminosity ranges from 0.45 L_☉ in Ophiuchus to 5.20 L_☉ in Serpens. The low-mass SF clouds do not contain protostars at luminosities at or above 1000 L_☉, and do not exhibit a distinct tail near 100 L_☉ as in the high-mass SF luminosity functions. The combined luminosity function does not show a peak, but rises toward lower luminosities down to L_cut, a feature akin to the Perseus luminosity function.

To establish whether the observed differences in the luminosity functions are statistically significant, we perform a K-S test on each of the 1000 realizations of the contamination-subtracted luminosity functions. We compare realization 1 of cloud 1 with realization 1 of cloud 2, realization 2 of cloud 1 with realization 2 of cloud 2, and so on (in each realization we have combined the high and low stellar density protostars into a single luminosity function). All 1000 realizations are used to take into account the uncertainty in the contamination removal. Table 4 gives the median K-S probabilities for each combination. We choose a threshold probability of 0.0027, equivalent to significance at the 3σ level in Gaussian statistics, to determine whether we can rule out the possibility of two realizations coming from the same parent distribution. We find that among the high-mass SF clouds, we cannot rule out that they are drawn from the same parent distribution: the median probability in the comparison of Orion and Cep OB3 is 0.0095, of Orion and Mon R2 is 0.0653, and of Cep OB3 and Mon R2 is 0.1168. Comparison with the lower mass SF regions shows that the Serpens luminosity function is not likely from the same distribution as the Ophiuchus luminosity function. The Ophiuchus luminosity function is not likely from the same distribution as the Serpens luminosity function, or the luminosity functions of any of the high-mass SF clouds. The Perseus luminosity function is not likely from the same parent distribution as any of the luminosity functions of the high-mass SF clouds. We also note that the comparison of the Ophiuchus and Perseus luminosity functions shows a probability close to our threshold. The Tau/Lup/Cha luminosity function cannot be distinguished from any of the distributions; this may result from the small number of protostars in these regions. We conclude that there are significant differences between the luminosity functions of the observed clouds, with some clouds showing relatively similar luminosity functions and others showing distinctly different luminosity functions. Of particular interest is the differences between the high-mass SF clouds and the low-mass SF clouds. The luminosity functions of high-mass SF clouds peak at and extend to higher luminosities than the low-mass SF clouds, which do not show a common, distinct peak in their luminosity functions. In Figure 10, we show the resulting distribution of K-S probabilities for the combined high-mass SF cloud luminosity function compared with the luminosity function from the combined low-mass SF clouds. The median probability is log (prob) = −4.61. The observed differences in the protostellar luminosity functions suggests a real difference between the properties of the protostars in these two cloud environments.

Zoom In Zoom Out Reset image size

Protostars are often found in extended regions of high extinction that can further redden the protostars already reddened by their infalling envelopes. Gutermuth (2005) devised a scheme for de-reddening pre-main-sequence stars that extended the approach of Meyer et al. (1997) to the Spitzer wavelength bands. This scheme takes into account sources with infrared excesses due to disks by de-reddening the stars onto the CTTS locus established in the combined 2MASS J, H, K_s, and IRAC 3.6 and 4.5 μm color space. Although any object detected in two of the five bands can be de-reddened, this procedure may not be appropriate for protostars for two reasons. First, we cannot distinguish between foreground extinction and the extinction by the envelope. Second, much of the light from the protostars at 1–4 μm is scattered by dust in the envelope, thus making the objects appear more blue. Both of these effects will result in an overcorrection for reddening and an artificially high luminosity.

Since we cannot reliably de-reddened protostars using their observed colors, we instead calculate the reddening from the Class II objects in the vicinity of the protostars. This approach, which was also adopted by the c2d team, assumes that the Class II objects are reddened by the same foreground reddening as the protostars, although Class II sources may not be as deeply embedded as the protostars. This foreground reddening includes both that from the extended molecular gas exterior to the protostellar envelope, and that from the interstellar medium between the observer and the molecular cloud. Although the reddening may vary on smaller scales than accounted for in this technique, we can still assess the magnitude of the effect reddening has on the luminosity function. The Class II objects can be de-reddened on the basis of the J, H, K_s, 3.6, and 4.5 μm photometry using the technique described by Gutermuth et al. (2005, 2009). To take into account variations in the foreground extinction due to the structure in the parental molecular cloud, we adopt the median extinction of the five Class II objects nearest to each protostar. We then use the reddening law of Flaherty et al. (2007) to de-redden the photometry of the protostar candidates. After de-reddening, we re-calculate the SED slope and L_MIR and find bolometric luminosities. The distribution of the ratio of de-reddened bolometric luminosity to uncorrected bolometric luminosity for all clouds is shown in Figure 11. We create luminosity functions using the de-reddened luminosities; the peak, median, mean, and 1σ are listed in parenthesis in Table 3. K-S tests give the likelihood that the uncorrected and de-reddened luminosity functions are from the same parent distribution with probabilities that range from 0.08 in Orion to 0.72 in Mon R2. This shows that although in certain regions it can be important, in most cases the reddening does not seem to have a large impact. Since the reddening corrections for the protostars are uncertain, we use both the uncorrected and de-reddened luminosity functions in the following analysis. We determine a de-reddened L_cut by adding the median increased luminosity (from uncorrected to de-reddened) to the uncorrected L_cut for each cloud. These de-reddened L_cut values are (in log (L/L_☉)) −1.82, −1.59, −1.21, −1.24, −1.33, −1.37, and −1.26 for Ophiuchus, Tau/Lup/Cha, Perseus, Serpens, Orion, Cep OB3, and Mon R2, respectively.

Zoom In Zoom Out Reset image size

We identified protostar candidates as sources with a flat or rising spectrum in the mid-IR and determine the bolometric luminosity of these sources using Equation (8). There have been some question of the relationship between flat-spectrum sources and rising-spectrum sources. Some flat-spectrum sources may be rising-spectrum sources observed from a face-on orientation through which emission from the warm inner layers can escape through the outflow cavity (Calvet et al. 1994; Whitney et al. 2003). Alternatively, protostars resulting from the collapse of a flattened sheet-like cloud can also give a flat SED (Hartmann et al. 2006). Finally, sources with tenuous envelopes can give a flat spectrum, in part because of the backwarming of the envelope (Natta 1993); such sources may be protostars at the later stages of envelope dissipation. Winston et al. (2007) also find that some flat-spectrum sources may be reddened disks, although we have accounted for those in our sample.

We define flat-spectrum sources as protostars with α between −0.3 and 0.3. The fraction of Spitzer protostars that are flat-spectrum ranges from 48% in Ophiuchus to 23% in Perseus. We compare the luminosity functions of the rising- and flat-spectrum sources in each cloud by using K-S tests to determine the probability that the luminosity functions are drawn from the same parent distribution. In each case, we only compare the distribution above the luminosity cutoff set for each region. This analysis yields the following probabilities for each region: Orion (0.95), Cep OB3 (0.01), Mon R2 (0.33), Serpens (0.88), Perseus (0.81), Ophiuchus (0.34), and Tau/Lup/Cha (0.02). Thus, we find the distributions for the flat and rising SED protostars statistically indistinguishable. Given the similarity of the luminosity functions, and given that the sample of flat-spectrum sources may contain many rising-spectrum sources observed at a face-on orientation, we only analyze the combined luminosity function for each set of protostars.

We compare the luminosity function of the 120 c2d protostars in this work with the observed protostellar luminosity function of 112 c2d protostars presented in Dunham et al. (2010). For this comparison, our Serpens, Perseus, Ophiuchus, Lupus, and Chamaeleon protostar luminosity functions were combined into a single c2d cloud protostellar luminosity function. The combined luminosity function for the c2d clouds from this work does not show a significant peak but increases toward lower luminosities and does not extend above 100 L_☉. Similar to our luminosity function, the observed protostellar luminosity function of the c2d protostars from Dunham et al. (2010) also does not extend above 100 L_☉. However, in contrast to our luminosity function, it shows a distinct peak above 1 L_☉. Thus, we find more low luminosity protostars in the c2d clouds than Dunham et al. (2010).

Offner & McKee (2011) give dimensionless quantities for the observed luminosity function from Dunham et al. (2010) in their Table 3, particularly the ratio of the median bolometric luminosity to the mean bolometric luminosity, L^med_bol/L_bol^mean = 0.3, and the standard deviation of log (L_bol), σ(logL) = 0.7. For our combined c2d luminosity function, L^med_bol/L_bol^mean = 0.22 and σ(logL) = 0.72. These values are the same as those derived for the Dunham et al. (2010) luminosity function. De-reddening the photometry using the method described in Section 5.2 brings σ(logL) to 0.77, which is consistent with the value for the Dunham et al. (2010) luminosity function, and L^med_bol/L_bol^mean to 0.15, which is less than the value for the Dunham et al. (2010) luminosity function. Thus, we find that on the basis of the quantitative statistics used by Offner & McKee (2011), our c2d luminosity function has a similar width as the luminosity function from Dunham et al. (2010), but when de-reddened, our luminosity function has a lower L^med_bol/L_bol^mean. The main difference is that we find an excess of faint objects and no peak. The reason for the difference may be that Dunham et al. (2010) require 1.3 mm envelope detections for their protostars; this may remove very low mass protostars with envelope masses less than 0.5 M_☉, from being included in their protostar sample (Evans et al. 2009).

Molecular clouds host a variety of SF environments, including regions of high and low stellar density YSOs. Although it has been shown that OB stars are typically found in clusters (Testi et al. 1999), it is unclear whether this is because of the fact that massive stars are rare and thus are only likely to be found in groups of low-mass stars (Bonnell & Bate 2006), or because these stars preferentially form in clusters (Bontemps et al. 2010). Since the luminosity of a protostar is a combination of accretion luminosity and intrinsic luminosity, we cannot determine the masses of our protostar candidates; however the most luminous protostars tend to be the most massive protostars (McKee & Tan 2003). We now compare the luminosity functions of high and low stellar density protostars to determine if the luminous protostars are found preferentially in high stellar density environments.

To provide a measure of the YSO density around each protostar candidate, we computed nearest-neighbor distances. The nearest-neighbor (nn) distance is the distance to the nth nearest Class II or protostar candidate. We chose n = 4 after considering n = 2 and n = 10 distances; n = 4 gives a better indication of clustering in both low and high stellar density SF clouds while n = 2 is dominated by random fluctuations (Casertano & Hut 1985) and n = 10 is not sensitive to clustering in smaller groups (hereafter the distance to the 4th nearest YSO will be denoted as "nn4 distance").

We show in Figure 12 the bolometric luminosity versus nn4 distance for all protostar candidates. For the three high-mass SF clouds, Orion, Cep OB3, and Mon R2, a trend can be noted between the nn4 distance and luminosity of the most luminous protostar at that nn4 distance. In these clouds, as the nn4 distance decreases the luminosity of the most luminous protostar increases. This can also be seen in the combined plot for all three high-mass SF clouds. Thus, the most luminous sources (L > 10 L_☉) are typically found with nn4 < 0.50 pc.

Zoom In Zoom Out Reset image size

This trend is more clear for the Orion clouds, in which we have the largest sample of protostar candidates. However, the most luminous protostar candidate in the Orion sample has an nn4 distance of 0.52 pc, greater than the median nn4 distance. In addition, the fifth most luminous protostar candidate in Orion has an nn4 distance of 0.38 pc, greater than the typical distance for sources of comparable luminosity in the Orion cloud. These sources, Reipurth 50 and V883 Ori, are thought to be undergoing FU Ori eruptive events (Strom & Strom 1993). Thus, these sources appear to be low-mass protostars undergoing luminous outbursts in which the accretion luminosity dominates the intrinsic luminosity (Hartmann & Kenyon 1996). Additionally, an outlier is seen in the third most luminous protostar candidate in Cep OB3 at an nn4 distance of 0.47 pc. We speculate that this may be an outburst source as well.

The low-mass SF clouds do not show the same trend. Instead, we find no convincing correlation between the most luminous protostar candidates and nn4 distance. We perform a two-dimensional K-S test comparing the distributions of nn4 distance and luminosity for protostar candidates with luminosities above the most conservative L_cut in the high-mass SF clouds with those in the low-mass SF clouds. We find a probability of 0.0023 that these distributions might be from the same parent distribution, which is below our threshold. The low-mass SF clouds are different from high-mass SF clouds because they contain few sources above 10 L_☉. In addition, the broad range of nn4 distances found in the high-mass SF clouds is not apparent in most of the low-mass SF clouds.

Is the trend that we see in Figure 12 the result of a decreasing sample size at larger nn4 (and hence fewer of the rare luminous protostars) or is it due to a real change in the luminosity function? To test this, we separate the Orion protostars into two equal-sized samples based on nn4 distance. The two samples are separated using the median protostar candidate nn4 distance, D_c, such that one sample has nn4 distances less than D_c, and the other sample has nn4 distances greater than D_c. Then, we compare the luminosity functions of protostars in regions of higher stellar density (nn4 distances <D_c) and luminosity functions of protostars in regions of lower stellar density (nn4 distances >D_c). In Figure 13, we show the luminosity functions of the protostars in regions of high and low stellar density in Orion. The left panel shows the distributions using the m₂₄ cutoff implemented in this work; the right panel emphasizes the small effect a 1 mag difference in the m₂₄ cutoff would have on completeness and is discussed in further detail below. The luminosity function of the protostars in regions of high stellar density peaks above 1 L_☉, while the luminosity function of the protostars in regions of lower stellar density peaks at 1 L_☉. We see for the contamination-subtracted luminosity functions, as in Figure 12 for all protostar candidates, that the majority of sources above 10 L_☉ are in regions of higher stellar density, though again we see that the most luminous protostar (the potential outbursting source) is indeed in a region of low stellar density.

Zoom In Zoom Out Reset image size

To test whether the differences between the high and low stellar density luminosity populations are significant, we perform 1000 K-S tests comparing the contamination-subtracted luminosity functions of the high and low stellar density regions in each cloud using the same method as the cloud-to-cloud comparison. The results are shown in Table 5 for both the uncorrected and de-reddened luminosity functions. Listed are the median K-S probabilities for the 1000 trials that high and low stellar density protostar luminosity functions are from the same parent distribution. For Orion, all 1000 probabilities are shown in Figure 14. We again use the threshold of 0.0027 to determine whether we can rule out the possibility of the high and low stellar density luminosity functions coming from the same parent distribution. We find that for both the uncorrected and de-reddened Orion luminosity functions we rule out the possibility that the luminosity functions in the high and low stellar density regions come from the same parent distribution.

Zoom In Zoom Out Reset image size

Could the difference observed between the high and low stellar density populations in the Orion clouds result from biases due to saturation and/or incompleteness in the 24 μm band? For many of the saturated sources, we are able to recover photometry using the method described in Section 2.1. We are missing sources in the regions of extended saturation in the Orion nebula and NGC 2024 nebula. It is in these high stellar density regions where the most luminous sources and densest clusters are found. We expect that the inclusion of these regions would enhance the differences we see between high and low stellar density environments in the Orion clouds.

However, incompleteness may still affect the result. There is a bias toward recovering fainter sources in less nebulous regions than in more nebulous regions because of the mid-IR nebulosity that exists most prominently in the MIPS 24 μm images. Source crowding will affect our sample less since we have used PSF photometry extraction at 24 μm. The Appendix discusses the completeness on our sample. This shows that in both the high and low stellar density regions we do not expect to detect all protostars down to our luminosity cutoff. The bright nebulosity and sources common in high stellar density regions may alter the luminosity function in these regions by preferentially hiding the faintest protostars in those regions.

One way of addressing the issue of incompleteness is to change the limiting m₂₄ used in the Orion luminosity function. Since the incompleteness comes primarily at the faintest magnitudes (the Appendix), changing the magnitude limit should reduce the effect of incompleteness. In Figure 13, we show the average of 1000 iterations of low and high stellar density luminosity functions after adopting a limiting magnitude of m₂₄ < 6 mag and m₂₄ < 7 mag. We find the median K-S probability that the high and low stellar density protostar luminosity functions are from the same parent distribution to be 0.0013 for the original m₂₄ < 7 mag cutoff. If we use a cutoff of m₂₄ < 6 mag, the median probability is 0.0015. Thus, the result that the high and low stellar density luminosity functions are significantly different does not change after raising the limiting m₂₄ by 1 mag.

Another way to assess the incompleteness is to compare the distribution of m₂₄ for Class II sources in high and low stellar density regions in Orion. If our sample is incomplete in regions of high or low stellar density, we expect to find differences between distributions of m₂₄ of Class II sources in these two regions. We use m₂₄ since we cannot use our conversion to luminosity for Class II objects. We identify Class II sources in Orion which are within D_c of a protostar candidate, then classify the Class II source as being in a region of high stellar density if the nearest protostar candidate is in a high stellar density region or classify the Class II source as being in a region of low stellar density if the nearest protostar candidate is in a low stellar density region. Figure 15 shows the m₂₄ distribution for high and low stellar density regions. A K-S test gives a 0.61 probability that the distributions of Class II sources in regions of high and low stellar density are from the same parent distribution. We note that the m₂₄ histogram for the protostars in high and low stellar density regions of Orion shows the same difference as the luminosity functions: the histogram for the high stellar density regions extends to brighter magnitudes with only a probability of 0.0089 that the two distributions are drawn from the same parent population. If we remove the two outbursting sources, this probability drops to 0.0051. Although these are above our 0.0027 threshold, they still correspond to the probabilities equivalent to a 2.6σ detection in Gaussian statistics. We conclude that the differences we see between m₂₄ distributions of protostar candidates in regions of high and low stellar density are not due to incompleteness in m₂₄ because we do not see differences in the m₂₄ distributions of the sample of Class II sources in regions of high and low stellar density.

Zoom In Zoom Out Reset image size

The same trend of increasing luminosity with decreasing nn4 distance is also apparent for Mon R2. The de-reddened Mon R2 luminosity functions show a probability that the low and high stellar density regions come from the same parent distribution equal to our threshold; this is further evidence that the luminosity functions of the high and low stellar density regions differ. The trend is less distinct for Cep OB3. The Cep OB3 low- and high-density luminosity functions are statistically indistinguishable. This may be due to the lower number of protostars at high stellar density (as evidenced by the D_c). Nevertheless, this cloud contains several more luminous objects that appear to be widely isolated; these objects should be followed up in future studies to determine whether they are undergoing outbursts.

The sample of protostellar luminosity functions presented in this paper provides the means to examine protostellar evolution in a diverse set of nearby clouds. The luminosity of a protostar is a combination of the intrinsic luminosity of the central protostar and the luminosity generated by accretion. Since much of the total protostellar luminosity is derived from accretion, the luminosity functions can put constraints on the rate of accretion onto the central star. This in turn can constrain models that determine the rate of infall of material on the central disk, and the subsequent accretion of matter from the disk onto the star. For an initial look as to how our data may constrain these models, we compare our results to recent models of protostellar luminosity functions in the literature.

Dunham et al. (2010) compared the c2d luminosity function with models of protostars based on a singular isothermal sphere collapse. This comparison includes five models that extend the analysis of Young & Evans (2005) by incorporating the following: scattering, an axisymmetric disk, an axisymmetric envelope, outflows and mass loss, and episodic accretion (EA). In each case, except the EA model, the resulting model luminosity function peaks above the observed luminosities presented in Dunham et al. (2010) and in our sample. One explanation for the lower observed luminosities is that the infalling material builds on protostellar accretion disks, and that the material on the disk is episodically dumped onto the protostar creating a large, but brief, jump in the source luminosity (Kenyon & Hartmann 1995; Dunham et al. 2010). Dunham et al. (2010) find that models which include EA can reproduce the range of luminosities observed in their c2d luminosity function.

Figure 16 shows our combined low-mass SF cloud luminosity function and the luminosity function from the EA model. This model has been designed for the low-mass SF c2d regions with a mass limit of 3 M_☉, and thus we do not compare this model with the high-mass SF clouds. The range of luminosity produced by the EA model matches well with our low-mass SF cloud luminosity function. However, the Dunham et al. (2010) model has a peak near 5.5 L_☉, while the low-mass SF cloud luminosity function from this work has no significant peak above L_cut. De-reddening (see Section 5.2) does not produce a significant peak in the low-mass SF cloud luminosity function nor does it improve the match between the low-mass SF cloud luminosity function and the EA model. We thus find that our luminosity function does not exhibit the peak evident in the EA model from Dunham et al. (2010).

Zoom In Zoom Out Reset image size

Offner & McKee (2011) compared a variety of accretion models to the observed luminosity functions of the protostars from Evans et al. (2009), including isothermal sphere, turbulent core (TC), two-component turbulent core (2CTC), competitive accretion (CA), and two-component CA models with both accelerating and non-accelerating star formation rate, and tapered or untapered accretion. We compare our observed luminosity functions to the Offner & McKee models in three ways. First, we again use the ratio L^med_bol/L_bol^mean and the standard deviation σ(logL) (see Section 5.4) to compare our observed luminosity functions with the corresponding values for the models from Offner & McKee (2011). The values for the low-mass SF clouds come from Table 3 of Offner & McKee (2011), which use an upper mass limit of 3 M_☉; values for the high-mass SF clouds use an upper limit of 10 M_☉ and are listed in Table 6 (S. Offner 2012, private communication). Second, we present a visual comparison of the data and models in Figure 17. Finally, in Table 7, we show the K-S test result of the comparison of the models to the data. We do not provide a comparison to the isothermal sphere model to our luminosity functions since it clearly does not match any of our observed luminosity functions (Offner & McKee 2011).

Zoom In Zoom Out Reset image size

The CA model shows a good match to L^med_bol/L_bol^mean and σ(logL) for the low- and high-mass SF clouds. The CA model also provides a reasonable match to the high- and low-mass SF cloud luminosity functions in Figure 17. Compared to the other models, the peak of the CA model is better matched to the broad plateau of the low-mass SF regions and it is well matched to the luminosity of the peak in the high-mass SF regions, particularly for the de-reddened luminosity functions, although the model underpredicts the amplitude of the peak observed toward the high-mass regions. The K-S test probabilities show that the observed high- and low-mass luminosity functions have the highest probability of being drawn from the CA model, although there is still a low probability that the high-mass luminosity function, particularly the uncorrected luminosity function, is drawn from the CA model.

The TC model provides a reasonable match for the L^med_bol/L_bol^mean and σ(logL) of the high-mass SF clouds, but overestimates L^med_bol/L_bol^mean by a factor of two for the low-mass regions. This is apparent in Figure 17, where this model shows a strong peak which is not apparent in the low-mass SF regions. On the other hand, it provides a good match to the high-mass regions, although the peak is at a too high luminosity. The observed high- and low-mass luminosity functions show a significantly lower probability of being drawn from the TC model than from the CA model.

The 2CTC model provides a poor match to L^med_bol/L_bol^mean and σ(logL) of the high- and low-mass SF cloud sample. The luminosity distributions of the 2CTC models peak at higher luminosities than the high- and low-mass regions, and the peak of this model is much sharper than apparent in the low-mass SF clouds. This is reflected in the K-S probabilities, which show a very low probability that our data are drawn from this model luminosity function.

Another two-component model is analyzed by Myers (2011), particularly a model similar to the 2CTC model of McKee & Tan (2003) and the two-component thermal and non-thermal model of Myers & Fuller (1992). Myers (2011) find that this model reproduces a realistic IMF between 0.1 and 10 M_☉. The model luminosity function shows a peak near 1 L_☉, akin to the luminosity functions from our sample of clouds which form high-mass stars. P. C. Myers (2012, in preparation) compares the resulting luminosity function with our Orion protostellar luminosity function and finds good agreement.

In conclusion, the shape of the luminosity functions of our high-mass and low-mass SF clouds are best matched by those produced by the CA luminosity function of Offner & McKee (2011). We note that there are still significant differences between this luminosity function and the observed luminosity functions. We will further discuss the implications of this in Section 6.3.

Of particular interest is whether the different luminosity functions observed in high and low stellar density regions in the Orion clouds imply a difference in the distribution of masses for emerging stars. This would be a form of primordial mass segregation, with potentially more massive stars forming preferentially in denser regions. We stress again that the bolometric luminosities we estimate are the sum of the intrinsic luminosity of the protostar and the accretion luminosity resulting from gas falling onto the protostar. In low-mass protostars, the accretion luminosity probably dominates, but for higher mass objects the intrinsic luminosity becomes increasingly important (McKee & Tan 2003; Offner & McKee 2011). We make no attempt to estimate the contribution of each of these luminosity sources, and thus, it is not clear whether the most luminous protostars are those with the highest intrinsic luminosity or whether they have the highest accretion rates. If they are sources with a higher intrinsic luminosity, then these sources probably are higher in mass, particularly if they follow a stellar birth line (Hartmann et al. 1997; Palla & Stahler 1991). In the case that the luminosity is dominated by accretion, then a higher luminosity might imply a higher protostellar mass or a higher accretion rate (Young & Evans 2005; Myers 2011). This might also imply a higher outcome mass for the protostar; however, if the accretion rate is higher and the accretion time is shorter, then sources in clusters may exhibit higher luminosities but have the same outcome IMF. Thus, although one interpretation of the different luminosities is that the high stellar density regions contain protostars that are more massive or that will form stars of higher mass, this is not a unique interpretation.

If we are indeed seeing primordial mass segregation, then this would suggest that the presence of massive stars in the dense center of clusters is due to the environment found in these crowded regions, and not the statistical sampling of a constant IMF (Elmegreen 1999; Bonnell & Clarke 1999). The CA model predicts a correlation between the density of stars in a cluster and the mass of the most massive member (Bonnell et al. 2004). This correlation is due to the large mass of gas drawn in by the gravity of the entire cluster of stars and accreted onto the most massive members of the cluster. The result is a higher mass accretion rate for the protostars in dense, clustered environments, particularly the most massive stars.

Although we observe a dependence of the luminosity function on stellar density, this may in effect be a dependence between the luminosity function and the gas column density. Gutermuth et al. (2011) presented observational evidence for a Schmidt-like star formation law in molecular clouds where the star formation rate per unit surface area is proportional to the gas density squared. In this case, the stellar density may be just tracing the overall column density of the natal gas. Consequentially the potential dependence of the stellar masses on stellar density may in fact be a dependence of the stellar masses on gas column density.

A dependence between stellar density and gas density is predicted by the CA model, but could also result from Jeans fragmentation in sheet-like clouds (Gutermuth et al. 2011). However, although Jeans fragmentation might explain the increasing stellar density with increasing gas density, it cannot explain the increasing stellar mass with increasing gas and stellar density since the Jeans mass decreases with increasing gas density. An increasing mass accretion rate with increasing gas column density is predicted by TC models, but it is not clear whether such models could also explain the higher stellar densities in the vicinity of higher mass stars (Cunningham et al. 2011). Thus, the competitive model is attractive in that it can explain the observed protostellar luminosity functions and their dependence on stellar density.

Based on the observational study presented in this paper, we can draw the following conclusions: that the peak of the luminosity function of clouds forming high-mass stars is near 1 L_☉ while for the clouds not forming high-mass stars the luminosity functions show a broad plateau that peaks below 1 L_☉, that the luminosity functions of the high- and low-mass SF clouds are different, and finally, that the luminosity functions drawn from low- and high-density regions in Orion are different. Our comparison of the luminosity functions in this study with the current models in the literature shows that none of the available models provide an excellent match to the observed luminosity functions. This is not surprising. The Offner & McKee models do not take into account details of protostellar evolution such as the clearing of the envelope. On the other hand, Dunham et al. (2010) only consider models starting with the collapse of a singular isothermal sphere. We expect that these results, as well as future luminosity functions derived from Herschel data, will drive the development of realistic models of protostellar evolution.

Within the available models, the CA model provides the best match. It best reproduces the shape of the low-mass SF cloud and high-mass SF cloud luminosity functions. We note that this model works well since it invokes an accretion rate that is dependent on both the instantaneous and final mass of the star; this is different than the models of isothermal sphere collapse where the infall rate is constant and depends only on the sound speed in the gas (Offner & McKee 2011). The different luminosity functions in the high- and low-mass SF clouds can be approximately reproduced with luminosity functions of an ensemble of protostars forming a Chabrier (2005) IMF truncated at upper mass limits of 3 M_☉ and 10 M_☉ for the high- and low-mass SF clouds, respectively (note that we have excluded the regions of high-mass star formation from our study of high-mass SF clouds due to saturation). Furthermore, the variations of the protostellar luminosity function with the spatial density of YSOs is also a prediction of the CA model, where the most massive stars are formed preferentially in the centers of dense clusters of stars (Bonnell et al. 2004). For these two reasons, we currently favor the CA model. The protostellar luminosity function, however, does not provide a definitive test for models of protostellar evolution; these models must be tested on other grounds as well, such as the measured velocities of the protostars and the density of the surrounding gas (Krumholz et al. 2005; Ayliffe et al. 2007; Winston et al. 2007) and the distribution of sources in a bolometric luminosity and temperature diagram (Dunham et al. 2010). Only though comprehensive observations that examine in detail the many facets of protostars and their environments, can we both test models of protostars and drive the refinement of those models to incorporate all the essential elements of protostellar evolution.