The Protostellar Luminosity Function
Abstract
The protostellar luminosity function (PLF) is the presentday luminosity function of the protostars in a region of star formation. It is determined using the protostellar mass function (PMF) in combination with a stellar evolutionary model that provides the luminosity as a function of instantaneous and final stellar mass. As in McKee & Offner (2010), we consider three main accretion models: the Isothermal Sphere model, the Turbulent Core model, and an approximation of the Competitive Accretion model. We also consider the effect of an accretion rate that tapers off linearly in time and an accelerating star formation rate. For each model, we characterize the luminosity distribution using the mean, median, maximum, ratio of the median to the mean, standard deviation of the logarithm of the luminosity, and the fraction of very low luminosity objects. We compare the models with bolometric luminosities observed in local star forming regions and find that models with an approximately constant accretion time, such as the Turbulent Core and Competitive Accretion models, appear to agree better with observation than those with a constant accretion rate, such as the Isothermal Sphere model. We show that observations of the mean protostellar luminosity in these nearby regions of lowmass star formation suggest a mean star formation time of 0.30.1 Myr. Such a timescale, together with some accretion that occurs nonradiatively and some that occurs in highaccretion, episodic bursts, resolves the classical “luminosity problem” in lowmass star formation, in which observed protostellar luminosities are significantly less than predicted. An accelerating star formation rate is one possible way of reconciling the observed star formation time and mean luminosity. Future observations will place tighter constraints on the observed luminosities, star formation time, and episodic accretion, enabling better discrimination between star formation models and clarifying the influence of variable accretion on the PLF.
Subject headings:
stars: formation stars: luminosity function, mass function1. Introduction
Protostars are born in dense, compact molecular cloud cores (McKee & Ostriker, 2007). During their earliest evolution, protostars are both dim and heavily obscured by a dusty envelope. It is an unavoidable observational reality that the majority of the star formation process occurs while protostars are deeply embedded within their natal gas. Consequently, high extinction and significant radiation reprocessing render the details of protostellar evolution, lifetimes, and the accretion process extremely uncertain (White et al., 2007).
Recent surveys of the nearby starforming regions have successfully obtained large statistical samples of young protostars with reasonable completeness down to luminosities of (e.g., Evans et al. 2009; Enoch et al. 2009; Dunham et al. 2008). Highresolution millimeter emission maps tracing dusty envelopes, provide a measure of core masses (Enoch et al., 2008). Combined with mid to farinfrared data, the available wavelength range is sufficient to trace the spectral energy distribution (SED), from which the total luminosity may be estimated.
Using the infrared spectral slope, observed sources can be divided into four classes that can be approximately mapped to evolutionary stages (Andre et al., 2000). Class 0 protostars are heavily obscured by a dusty envelope, such that most of the radiation falls in the submillimeter band. During the Class I phase the protostar, while still embedded, becomes less obscured and may be surrounded by a thick circumstellar accretion disk. By the Class II phase, the now premain sequence star has accreted or expelled most of the initial envelope mass and produces little submillimeter emission. The remaining gas lies in thin accretion disk surrounding the star. Signatures of outflows may be apparent during both the Class I and Class II phases. During Class III, the disk dissipates and the source approaches the main sequence.
Despite this straightforward picture, cataloging sources and definitively mapping them to a physical stage remains complicated. Geometric effects shift objects over the Class 0/I boundary, distorting the correlation between physical stage and SED characteristics. Edgeon Class II protostars with higher extinction may be misclassified as Class I, while Class I sources viewed down the outflow cavity resemble Class II sources (Masunaga & Inutsuka, 2000; Robitaille et al., 2006). Variability in the accretion rate may cause the protostar to oscillate across class boundaries (Dunham et al., 2010). Measurements of the millimeter emission used as a proxy for the envelope mass can be used to distinguish between embedded, i.e. Class 0 and Class I objects, and nonembedded, Class II objects (Enoch et al., 2008), although objects with thick disks along the lineofsight may still be misclassified.
During the Class 0 and I phases, luminosity due to accretion likely dominates the total radiative output (Evans et al., 2009; Dunham et al., 2010). Consequently, upper limits for the accretion rates may be inferred from the luminosity (Enoch et al., 2009). This information gives clues about the formation timescale and the accretion process during which the protostars are deeply embedded and cannot be directly imaged. However, the large distribution of observed luminosities, and in particular the significant number of dim protostars, creates a picture that is largely inconsistent with predictions of some star formation models.
1.1. The Luminosity Problem
The accretion rates inferred and the star formation times predicted by theoretical models combine to suggest that protostars are on average too dim (e.g., Kenyon et al. 1990; Young & Evans 2005; Enoch et al. 2009). The typical accretion luminosity of a protostar
(2)  
exceeds the observed average luminosity of about and significantly exceeds the observed median luminosity of (Enoch et al., 2009). Here we have adopted a typical mass accretion rate of yr (McKee & Ostriker, 2007) and a typical protostellar radius of (Stahler, 1988; Hartmann et al., 1997). The factor is the fraction of the gravitational potential energy that is radiated away. Note that is the total accretion luminosity, including any emission from an accretion disk around the star. There are three factors that can reduce the accretion efficiency below unity: First, some of the energy of the gas that accretes onto the protostar can be extracted by a hydromagnetic wind (e.g., Ostriker & Shu, 1995). If half the energy lost by the disk is mechanical instead of radiative, then, since half the total potential energy is lost by the disk, this would give . Second, some of the energy in the infalling gas can be absorbed by the star; Hartmann et al. (1997) argue that this is not significant for accretion rates yr, and this has been confirmed by Commerçon et al. (2011) so long as the accretion flow is transparent to optical radiation. Finally, some of the accretion energy is used in dissociating and then ionizing the accreting gas (Tan & McKee, 2004), but this is significant only for very low masses. Thus, of these three effects, the mechanical loss of energy is the only one of significance. As we shall see in Section 3.4, episodic accretion can reduce the observed mean luminosity below the true mean and therefore contribute an effective reduction in .
1.2. Review of Past Work
A number of mechanisms have been suggested to resolve the luminosity problem. A successful solution must account for both the mean luminosity of the distribution and the spread of luminosities over several orders of magnitude. The first indication of the discrepancy between accretion and luminosity was observed in Taurus by Kenyon et al. (1990); Kenyon & Hartmann (1995). At that time, observations of only a handful of young protostars existed, and the authors speculated that the problem could be resolved by, among others, short, unobserved periods of high accretion or increased ages of T Tauri stars.
Fletcher & Stahler (1994a, b) performed the first derivations of protostellar mass and luminosity functions. Assuming constant accretion rates and a fixed time period during which stars formed with a constant star formation rate, they derived the timedependent luminosity function for a cluster of evolving protostars. Their model included the luminosity contribution from Class I to mainsequence stars. For steadystate star formation, our constant accretionrate model produces a protostellar protostellar luminosity function similar to that derived by (Fletcher & Stahler, 1994b) when premain sequence stars are excluded (see §3).
Myers et al. (1998) used a semianalytic model to describe the evolution of the protostar, disk, and envelope system. With simple radiative transfer, they generated bolometric temperature and luminosity tracks from Class 0 to the zeroage main sequence stage. In order to represent the decline in accretion luminosities with time, the models assumed an exponentially decreasing accretion rate. For some parameters, the tracks were able to reproduce the mean protostellar luminosities. However, the breadth of the distribution and the relative number of lowluminosity sources were not addressed.
MonteCarlo radiation transfer modeling by Whitney et al. (2003) underscored the the dependence of the SED on core geometry. In particular, they showed that the inferred bolometric luminosity may vary by 50% from the true luminosity depending upon the orientation of the disk and outflow cavity relative to the lineofsight. For a sufficient number of observed protostars sampling all geometric projections, the orientation can only alter the distribution of the luminosities, not the mean.
Young & Evans (2005) performed 1D radiative transfer observations of the insideout collapse model developed by Shu (1977). These calculations improved upon the Myers et al. (1998) work by following the hydrodynamic evolution of the protostellar core. They included six contributions to the total luminosity with the result that their tracks agreed only with the brightest protostars, further illuminating the discrepancy between observed luminosities and predictions of lowmass star formation.
Recent extensive surveys of five local starforming regions have provided more comprehensive statistics and observational data of the earliest stages of protostellar evolution (Enoch et al., 2009; Evans et al., 2009) This work demonstrated that Class 0 and Class I sources have comparable mean luminosities, suggesting that protostars may accrete significantly for longer periods of time than previously assumed. Extending the accretion epoch from 0.1 Myr, as first assumed by Kenyon et al. (1990), to 0.5 Myr presents one possible solution to the luminosity problem.
Dunham et al. (2010) revisited the Young & Evans (2005) methods and explored several improvements, including updated opacities, 2D effects such as disks, core rotation, outflow cavities, and variable accretion. With all the additions, although predominantly due to the inclusion of episodic accretion, good statistical agreement between the observed and modeled temperature and luminosity distributions was possible. They modeled variability by halting accretion from the disk onto the protostar until the disk mass reached 20% of the star mass, whereupon the star was allowed to accrete at a constant rate of yr until the disk mass was exhausted. Since other terms contributing to the total luminosity (e.g., accretion from the envelope onto the disk and the stellar luminosity) are small compared to the accretion luminosity, the model assumes that protostars exist in a very low luminosity state for most of their formation time and accrete most of their mass during bursts. Statistically, their models predict that protostars must spend % of the time radiating at . This is only marginally consistent with the observed sample of 112 protostars, which contains no sources more luminous than . The models also have a star formation time Myr, which is less than the star formation time inferred by Evans et al. (2009). Future surveys containing more objects are necessary to test the importance of episodic accretion.
There is debate whether allowing for episodic accretion in evolutionary models can reproduce the dispersion of lowluminosity sources on the HertzsprungRussell diagram and thus account for a portion of the age spread inferred for members in young, lowmass star clusters (Baraffe et al., 2009; Hosokawa et al., 2011). If so, episodic accretion could also explain lithium depletion measurements in young stars (Baraffe & Chabrier, 2010). However, this solution is valid only if nearly all the accretion energy is efficiently radiated away from a small fraction of the stellar surface. Other mechanisms such as varying initial conditions and assumptions about the radiative properties of the accretion flow, which are not well constrained, may also contribute to an apparent age spread (Hosokawa et al., 2011).
1.3. Model Assumptions
In this paper we derive the Protostellar Luminosity Function (PLF) predicted by various star formation models and compare with observations (Enoch et al., 2009). Comparison with the mean observed luminosity determines the mean star formation time; comparison with the shape of the PLF (ratio of median to mean and the standard deviation) tests the validity of the models. Our basic formalism (McKee & Offner 2010a, henceforth Paper I) can be adapted to consider a timevarying rate of star formation, bursty accretion,and an arbitrary stellar initial mass function. However, in this paper we base our comparison on four main assumptions. First, we treat only the cases in which the star formation rate is either constant or smoothly accelerating in time (e.g. Palla & Stahler 2000). Second, we assume the accretion rate is a smooth function of time during the accretion phase; in particular, we assume it can be expressed as a function of the current mass and the final mass of the protostar. Rather than excluding stochastic variability completely, we assume that any timevarying component is statistically rare in the samples with which we compare and therefore unlikely to be included the data (see §5.3 for a detailed discussion). Third, we assume that the accretion rate onto the protostar, which can be inferred from the protostellar luminosity, tracks the accretion rate onto the protostardisk system over the majority of the lifetime of the protostar. Finally, we adopt an individualstar Chabrier (2005) IMF truncated at an upper mass limit .
Through our derived PLF models, we investigate the luminosity problem in the context of different star formation theories and explore variations such as an accelerating star formation rate and accretion rates that taper off over time. In §2 we review the Protostellar Mass Function, which is described in detail in Paper I. In §3 we derive the PLF and relative statistics. We compare with observations of local starforming regions in §4 and discuss our results in §5. Approximate analytic results for the PLF for different star formation models are given in the Appendix.
2. The Protostellar Mass Function (PMF)
In Paper I we determined the mass function of protostars (the PMF) in terms of the IMF and the accretion history of the protostars. We assumed that the accretion history could be expressed in terms of the current protostellar mass, , and its final, stellar mass, . Henceforth, and are measured in units of solar masses and is in units of solar masses per year. We expressed the IMF as , where is the fraction of stars born with final masses in the range . We then defined the bivariate PMF, , such that is the fraction of protostars in a region of star formation with masses in the range and final masses in the range . We showed that for steady (i.e., nonaccelerating) star formation
(3) 
where is the formation time for a star of mass and is the average of over the IMF. The PMF, , is just the integral of the bivariate PMF,
(4) 
where the lower limit of integration is
(5) 
and the most and least massive stars formed by the cluster are and , respectively.
To determine the PMF, we required the accretion histories of the protostars. As noted above, we assume that these accretion histories, apply to the gas reaching the protostellar surface as well as to the gas that accretes from the ambient medium onto the protostardisk system. We considered several different models:

Insideout collapse of an isothermal sphere (IS, Shu 1977),
(6) 
The Turbulent Core model (TC, McKee & Tan 2002, 2003), in which the initial core is presumed to be supported by turbulent motions instead of thermal pressure,
(7) (8) where is in yr, is the surface density of the clump in g cm in which the stars are forming, and the parameter is determined by the density profile of the core; we follow McKee & Tan (2003) in setting . Tan & McKee (2004) suggested that the coefficient in the accretion rate in equation (8) be increased by a factor in order to allow for infall at the beginning of the insideout collapse. We retain the value given by McKee & Tan (2003), but note that the overall normalization of the accretion rate remains uncertain.
All these forms for the accretion history can be summarized in the expression
(10) 
where the model exponents are given by:
(11)  
(12)  
(13) 
We also considered a twocomponent turbulent core model (2CTC), a blend of the IS and TC models:
(14) 
where ; we adopt as in Paper I. This model is similar to the TNT model of Myers & Fuller (1992). There is some evidence that a twocomponent accretion model may be more physical in the CA case as well (i.e., a 2CCA model). Smith et al. (2009) show that, at least initially, competitively accreting protostars are surrounded by a small bound envelope. This suggests that at early times the protostars may accrete via a Shulike constant accretion rate. However, the authors find that the envelope mass is not well correlated with the final mass of the star, which they suggest indicates that the CA phase dominates. For comparison, we define a 2CCA model with:
(15) 
where
(16) 
According to Smith et al. (2009), accretion of the core envelope should contribute less than half of the final mass. Therefore, we adopt , which is determined by assuming that a star of average mass accretes half its mass from its envelope.
In the fiducial cases above, the accretion rates increase monotonically over the protostellar lifetime. For the singlecomponent models, the time to form a star, is a powerlaw function of the final mass,
(17) 
where
(18) 
is the time to form a star with . In reality, we expect accretion rates to gradually taper off (e.g., Myers et al., 1998; Myers, 2010), at least for the IS and TC models. Tapered accretion is supported observationally since Class I sources are not significantly brighter than Class 0 sources, and in some cases, the most luminous young source is actually Class 0 (e.g., Evans et al. 2009). Class I sources are also twice as numerous, suggesting that high accretion rates, which depend upon the presence of dense infalling gas, cannot be sustained throughout the majority of the formation time (Enoch et al., 2009). In competitive accretion models, the decline in accretion rates is assumed to be abrupt, so untapered models are likely to be the best approximation for the CA case. If we assume that the decline is linear in time, the tapered accretion case can be written:
(19) 
where is the untapered accretion rate for . We can then express the accretion rate for both the tapered and untapered cases as
(20) 
where is the untapered accretion rate for and
(21) 
The star formation time for singlecomponent models can then written as:
(22) 
In the case of accelerating star formation we assume a birthrate that increases exponentially with time, i.e.,
(23) 
where is the age of a star with mass and final mass and is the characteristic acceleration time. For an accelerating star formation rate, the mean formation time is no longer equal to the observed star formation time. Instead,
(24)  
(25) 
where is the lifetime of Class II sources. For Myr and Myr, the actual mean star formation time of the accelerating cases, , is a factor smaller than the value inferred from observations of the numbers of protostars and young stellar objects, .
3. The Protostellar Luminosity Function (PLF)
3.1. Definition
Our objective is to determine the protostellar luminosity function, , where is the fraction of protostars that have luminosities in the range . The protostellar mass and accretion rate are related to the accretion luminosity through equation (2), which gives . Unfortunately this relation is complicated by the fact that the protostellar radius is a function of both the mass and the accretion rate, , which does not have a simple analytic representation (Stahler, 1988; Hartmann et al., 1997). We use the routine described by Offner et al. (2009), which agrees with the results of Stahler (1988), to determine the protostellar radius (see §3.3). This model exhibits evolution similar to the model with shock boundary conditions in Hosokawa et al. (2011). Consequently, the subsequent evolution of the stars will have ages that are consistent with previously calculated nonaccreting isochrones.
We proceed by using equation (10) to determine from , and then solving the accretion luminosity equation (2) for . Since there are two independent variables, we use bivariate distribution functions: The fraction of protostars in the luminosity range and final mass range is the same as the fraction protostars with masses in the range and final mass range :
(26) 
where was determined in Paper I (see eq. 3). The PLF is then obtained by integrating over ,
(27)  
(28) 
In the Isothermal Sphere case, the accretion rates are independent of the final mass. Consequently, the PLF reduces to:
(29) 
3.2. The Mean and Standard Deviation of the Luminosity
In general, the mean luminosity can be evaluated as . However, it is more instructive to return to the bivariate luminosity function and define the mean accretion luminosity: (26):
(30)  
(31) 
where equation (26) enabled us to eliminate the dependence on for steady star formation and where
(32) 
is the massaveraged harmonic mean radius. Note that equation (31) is only valid in the case of steady star formation, where the accretion luminosity is proportional to the binding energy of the stars formed. In order to evaluate the mean luminosity, we take advantage of the fact that the stellar radius depends primarily on the instantaneous mass and only weakly on the final mass in order to define a mean radius in terms of the current mass ,
(33) 
This approximation is reasonably accurate for the value of the upper mass limit we consider here (); for larger values of , it would be preferable to include a weighting with the IMF.
For protostars with , the stellar luminosity, i.e., the luminosity due to nuclear burning and KelvinHelmholtz contraction, makes only a small contribution to the total luminosity, so that . However, for more evolved protostars beginning late in the Class I stage, accretion diminishes and no longer dominates the total luminosity (White & Hillenbrand, 2004). For the young, lowmass, embedded sources we consider here, the accretion luminosity is much greater than the nuclear luminosity provided . For such stars with tapered accretion, the accretion luminosity dominates during the majority of the evolution, i.e., while . For stars with , the nuclear luminosity becomes important and the accretion luminosity is not a good approximation of the total. However, these stars comprise only a small part of the total sample. For example, the maximum expected protostellar mass, , for a cluster with and (see Figure 7 in Paper I) in the CA case corresponds to a maximum ratio of . In the fiducial IS case, yields , which reflects both the lower accretion rate of the IS model and the higher maximum mass for a cluster with a fixed number of protostars in this model.
The standard deviation is a useful metric for characterizing the breadth of the protostellar luminosity distribution. Since the luminosity spread may encompass several orders of magnitude, we adopt the standard deviation of the logarithm of the protostellar luminosity:
(34)  
This metric has the advantage that it is dimensionless and, like the mean, it can be computed solely in terms of the protostellar mass function if is provided. The ratio of the median to the mean, also a useful dimensionless number, measures the skewness of the distribution.
3.3. Stellar Evolution Model
In order to determine the PLF predicted by a given theory of star formation, we must first construct a comprehensive luminosity model that takes into account both stellar evolution and the accretion history. We adopt the onezone protostellar evolution model developed by Tan & McKee (2004), which is described in detail in Offner et al. (2009). This model has been calibrated to replicate the more detailed calculations of Palla & Stahler (1991, 1992) and updated to reflect the recent recommendations of Hosokawa & Omukai (2009). It corresponds to the “hot” stellar evolution model described in Hosokawa et al. (2011), where an accretion flow directly hits the stellar surface and forms an accretion shock front. Although the gas may be channeled through a disk, the model assumes that it is sufficiently thick such that the accretion column covers most of the stellar surface. In practice, the initial physical conditions and accretion flow properties are not well known, introducing uncertainties in the stellar radius of a factor of for . (e.g., Hosokawa et al. 2011). Figure 1 shows the stellar radius predicted by this model as a function of stellar mass for different accretion histories. For , the model has a maximum disagreement with Hosokawa & Omukai (2009) of % when and generally agrees within 5%. The model has a maximum disagreement of a factor of 2 at for , but it generally is within a factor of 1.2 Hosokawa & Omukai (2009). Although our subgrid model has been thoroughly described in Offner et al. (2009), we give a brief summary here.
The model treats the protostar as a polytrope and follows the stellar contraction by enforcing energy conservation. The model is characterized by six stages, culminating in the arrival of the protostar on the zero age main sequence. In the“precollapse” stage, the collapsing gas densities are substellar. In the second “no burning stage”, the densities are sufficient for the dissociation of H, but are too low for deuterium burning. In the following “core deuterium stage at fixed ,” deuterium burning ignites in the core. Once the supply of deuterium is depleted, the core temperature rises and the protostar enters the “core deuterium burning at variable .” Eventually high temperatures in the core terminate convection and halt the flow of deuterium into the core, whereafter the protostar enters the “shell deuterium burning” stage. Finally, the central temperature reaches K, and the star moves onto the main sequence. In calculating the total luminous output, the zeroage main sequence luminosity serves as an approximation of the interior luminosity, which results from the combination of deuterium burning and KelvinHelmholz contraction. We evolve this model in combination with the accretion rates specified by the star formation models and generate luminosity tables as a function of the instantaneous and final stellar masses for each case with both untapered and tapered accretion. For an accelerating star formation rate, we assume that the physical parameters are set by the initial conditions of the collapse and hence do not themselves accelerate. Thus, for the untapered accelerating and tapered accelerating star formation cases, we use the untapered and tapered luminosity tables, respectively.
We use these tabulated values of to obtain the mean and standard deviation of the luminosity directly. However, in the model, and undergo discontinuous jumps as the deuterium state changes. This is problematic in calculating the PLF, which requires the derivative of . To circumvent this issue, we use a polynomial fit to equation (33). We renormalize the result using the directly derived value of . This strategy simplifies the integral while preserving the trends in the subgrid model and eliminating the weak dependence of on .
3.4. Episodic Accretion
Observations of highluminosity variable sources, such as the prototype FU Ori, suggest that protostars undergo periods of high accretion (e.g., Hartmann & Kenyon 1996). In the extreme case, protostars may spend most of their life in a lowluminosity, lowaccretion phase and accrete nearly all their mass during short, intense accretion bursts. However, only a total of 18 bursting sources have been identified within 1 kpc of the Sun (Greene et al., 2008). These sources have luminosities in the range 20550 , corresponding to accretion rates of yr. The total mass that can be accreted in such bursts is limited by the known star formation rate. Following McKee & Offner (2010a), let be the fraction of mass accreted during episodic accretion periods:
(35) 
where is the total time spent in the bursting state and is the mean stellar mass for a typical IMF. The time spent in the high accretion state can be expressed as:
(36) 
where is the number of protostars accreting in this state within 1 kpc of the Sun (Greene et al., 2008) and the star formation rate, , is 0.016 stars/yr within 1 kpc (Fuchs et al., 2009). This gives yr and 0.25 for an accretion rate of yr. If the protostellar lifetime is 0.5 Myr, then the duty cycle must be 0.2%.
This estimate suggests that our sample of 112 protostars is unlikely to contain any bursting sources. For our observational comparison, we adopt an episodic accretion factor, , to reflect the amount of mass accreted during outbursts. If no protostars are currently undergoing bursts, then the observed sample will have a mean luminosity that is 75% of the true timeaveraged mean. To account for this, we reduce our model mean luminosity by a factor of . In other words, we adopt an effective value for in equation (2) of
(37) 
Note that the Dunham et al. (2010) episodic accretion model corresponds to for a star with . For , this corresponds to An alternative way of correcting for episodic accretion assumes that the protostars accrete at their normal rate for a time and then accrete via bursts for a brief additional time. In this case, the luminosity of the protostars would be unchanged and instead the formation time would be correspondingly shorter. We adopt the former definition of since it seems more plausible and is consistent with simulations exhibiting episodic accretion (e.g., Vorobyov & Basu 2005).
4. Comparison With Observations
4.1. Overview of the Observational Data
Throughout this section, we compare to the Class 0 and Class I source luminosities of Serpens, Ophiuchus, Perseus, Lupus and Chameleon from Enoch et al. (2009) and Evans et al. (2009). This is the same sample used in Dunham et al. (2010). It has been corrected for localized extinction and carefully analyzed to remove nonprotostars and nonembedded sources (Mike Dunham, private communication). The data are comprised of 39 Class 0 and 73 Class I objects for a total sample size of 112. This includes only two sources from Chameleon II and one source from Lupus, such that the sample essentially reflects the properties of Serpens, Ophiuchus, and Perseus. The Class 0 and Class I sources have a 0.50 likelihood of being drawn from the same population, which suggests that the two groups are quite similar. It is unclear if this is due to similar accretion rates or source mixing resulting from inclination effects (e.g., Robitaille et al. 2006).
As discussed by Enoch et al. (2009) and Dunham et al. (2008), the bolometric luminosities have an uncertainty of 50% due to saturation of the 160 m band. Enoch et al. (2009) estimate that 10% error arises from the SED fitting process and 1525% due to finite sampling errors. Although the data we use have improved 350 m observations and more careful source selection, uncertainty arises from a number of factors and is not well constrained. Following Enoch et al. (2009), we adopt a 50% upper uncertainty in the bolometric luminosities, and we use the nonextinction corrected bolometric luminosities, which are naturally underestimates, to set a lower error limit (Melissa Enoch and Mike Dunham, private communication). The extinction–corrected mean and median bolometric luminosities are then:
(38)  
(39) 
Note that the mean is very sensitive to the most luminous sources. The extinction corrected sample of Enoch et al. (2009), which is derived from a nearly identical source list, has a mean of ,which falls within the error of our adopted mean. The disagreement arises from a combination of improved 350 m data used in the Evans et al. (2009) analysis and different treatments of the IRAS fluxes. In either case, the luminosity of a few of the brightest sources may have been overestimated by as much as a factor of two if the 24m flux was omitted due to saturation (Mike Dunham, private communication). Consequently, the upper uncertainty should be considered an upper limit on the luminosity rather than a onesigma error estimate. Since the standard deviation is also sensitive to the distribution outliers, we use the log of the standard deviation for comparison:
(40) 
A second useful nondimensional quantity is the ratio of the median luminosity to the mean:
(41) 
However, the associated uncertainties of this ratio are somewhat uncertain since the median and mean errors are correlated.
In this sample, the mean luminosities of the individual observed star forming regions and their shapes are similar, even though the number of sources in each is different. The PLFs of Perseus, Serpens, and Ophiuchus are shown in Figure 2, where the regions have mean luminosities 5.2 , 5.0 , and 6.2 , respectively. A KolmogorovSmirnov (KS) test indicates that Perseus and Serpens have a 0.39 likelihood of of being drawn from the same parent distribution, while Perseus and Ophiuchus have a 0.45 likelihood and Ophiuchus and Serpens have a 0.86 likelihood. This high level of agreement, despite the apparent dissimilarity of the distributions shown by Figure 2, is partially due to the small number of protostars ( 20) in Serpens and Ophiuchus. While comparing our models with the individual regions would be ideal, small statistics requires that we use the total protostellar distribution over all the clouds.
To derive the model distribution, we require the mass of the most massive star forming in the cluster, . The brightest Class II source in Evans et al. (2009) suggests an upper mass of . Since the number of Class II sources is several times larger than the sample of Class I and Class 0 sources, this should statistically be a good upper limit. However, statistical variations or changes in the star formation rate may impact the upper mass of the currently forming stars. Inspection of the most luminous Class I source also suggests an upper mass of , assuming that its luminosity is dominated by stellar radiation rather than accretion. This latter estimate is also very uncertain since we can’t discount that the source is undergoing an outburst and its luminosity is dominated by accretion. However, we can be fairly confident that is a good upper limit.
The accretion timescale of the models is constrained by the observed protostellar lifetime, i.e., how long protostars spend in the main accretion phase. Both Stage 0 () and Stage I () protostars experience significant accretion, since a significant fraction of the total gas is contained in the envelope (Crapsi et al., 2008). These stages roughly correspond to Class 0 and Class I provided that the Class I envelope mass exceeds 0.1 . Enoch et al. (2009) estimated 0.1 Myr for the Class 0 lifetime. This is longer than the Class 0 lifetime of yr reported by Andre et al. (2000), because these authors base their lifetime solely on data from Ophiuchus, which is now recognized to have a deficit of Class 0 objects compared to other star forming regions (Enoch et al. (2009), who also include one more Class 0 object than Andre et al. (2000) in their Ophiuchus sample.) For local star forming regions, Evans et al. (2009) report an average combined Class 0 and I lifetime of 0.44 Myr. The uncertainty in this estimate arises in part from statistics and Class I/II source confusion, but it is dominated by the uncertainty in the disk lifetime ( Myr), which is necessary to calibrate the ages. Altogether, the mean protostellar lifetime is Myr in the Evans et al. (2009) sample. Their estimates of the mean protostellar lifetimes in individual clouds, Ophiuchus (0.30 Myr), Serpens (0.46 Myr), and Perseus (0.62 Myr), are within the errors of the overall mean. This lifetime is significantly longer than previous estimates, which have adopted either the Class 0 lifetime (e.g., Kenyon et al. 1990) or the core freefall time (0.1 Myr for a 0.5 star; e.g., Myers et al. 1998; Young & Evans 2005). Shorter lifetimes exacerbate the luminosity problem. As discussed by McKee & Offner (2010a), one possible solution is “slow accretion,” in which accretion rates are reduced by protostellar outflows or lengthened disk lifetimes. Estimated protostellar lifetimes on the order of Myr lend credence to a slow accretion picture.
Some studies of Class I protostars have found that stellar luminosity dominates the total luminosity, suggesting that accretion has already diminished Muzerolle et al. (1998); White & Hillenbrand (2004); Connelley & Greene (2010). However, indepth study often reveals that sources classifed as Class I are actually edgeon or heavily obscured Class II sources (White & Hillenbrand, 2004; van Kempen et al., 2009; Heiderman et al., 2010). The observational sample here uses a minimum gas mass criteria of for inclusion in the sample, which is an indicator that accretion is still underway. However, the authors do not probe for tracers of dense, warm gas, such as emission from higher transtions of HCO and CO, which more conclusively separates embedded protostars from obscured Class IIs van Kempen et al. (2009); Heiderman et al. (2010).
4.2. The Mean Luminosity,
The mean of the luminosity distribution serves as a simple onedimensional statistic with which to compare the models and observations directly. In Figure 3, we plot the mean luminosity as a function of the most massive star forming in the cluster for each of the accretion models. We truncate the IS model at since it is unrealistic for describing highmass star formation. In Table 1, we give the model values calculated assuming completeness of the extinctioncorrected data down to for . We plot these tabulated values and the observed mean () from Evans et al. (2009) in an inset plot in Figure 3. As with the mean protostellar mass derived in Paper I, the mean luminosity increases strongly as a function of the upper stellar mass. Unfortunately, for , the model spread is small, making it difficult to discriminate between models. In contrast, the model means diverge significantly for clusters with large in all cases, suggesting that more complete observations of highmass starforming regions would be useful to distinguish models based solely upon the mean luminosity. Note that the model means depend upon the uncertain starformation timescale, so that it is not possible to make an independent comparison of the models and observations (see discussion in §5). However, better constraints on the star formation time in the future should increase the discriminating value of the mean luminosity in comparing models.
As shown in the plot insets, the mean observed luminosity falls above
the models in all the nonaccelerating cases
cases. The untapered TC and CA means, IS tapered mean, and
all accelerating means are consistent with the observational error.
NonAccelerating  Accelerating 


Model  Untapered  Tapered 
Untapered  Tapered 
Isothermal Sphere (IS)  3.27  3.81  4.50  4.49 
2C Turbulent Core (2CTC) 
2.97  3.17  5.25  5.44 
Turbulent Core (TC)  3.63  3.09  5.70  4.44 
2C Competitive Accretion (2CCA) 
2.91  2.68  4.97  4.14 
Competitive Accretion (CA)  4.26  3.26  7.05  4.73 
The mean luminosities in the accelerating cases are up to 30% lower than in the fiducial nonaccelerating cases for a fixed value of because the accelerating cases have more lowmass protostars. However, for a given observed value of , the mean luminosities are raised since their formation time is a factor smaller, as shown in equation (25). It is this adjustment that accounts for the good agreement of the accelerating models with the observations.
Allowing the accretion rates to taper off during the accretion period has a varying effect on the mean luminosities. In the IS case, the accretion rate is no longer independent of mass, increasing the mean luminosity. In comparison, the other models accrete at higher rates at early times, lowering the mean luminosity. Adding tapering also differentiates the models slightly at . However, when the curves are renormalized using the observational completeness limit, the differences are minimized and in the TC and CA cases the mean increases. Among the tapered, nonaccelerating accretion models, only the IS model falls within error under the constraint that the mean star formation time is 0.44 Myr.
4.3. Median Luminosity
While the mean is sensitive to the maximum luminosity, which observationally may be prone to both overestimation and statistical fluctuations, the median serves as a better proxy for the peak of the distribution. In this case, systematic errors in the spectra fitting and extinction corrections are likely to dominate the error in the observed value. Table 2 gives the medians for each of the cases for and Myr, where the observed median is . For the fiducial case, only the 2CTC and TC models are within error. When tapered accretion is included, all but the IS model agree. For an acclerating star formation rate, the untapered 2CTC, 2CTC, TC and tapered CA and TC models are within error.
NonAccelerating  Accelerating 


Model  Untapered  Tapered 
Untapered  Tapered 
Isothermal Sphere (IS)  3.12  3.19  3.86  3.22 
2C Turbulent Core (2CTC) 
1.76  1.85  1.41  2.44 
Turbulent Core (TC)  1.44  1.52  1.91  2.22 
2C Competitive Accretion (2CCA)  0.95  1.75  1.52  2.76 
Competitive Accretion (CA)  0.88  1.25  1.31  1.69 
4.4. The Star Formation Timescale,
We see that in some cases the mean and median luminosities are quite different from the observed values when we fix the starformation time at Myr. However, this time scale is uncertain by a factor 2. Henceforth, rather than fixing the observed lifetime, we derive the average formation time, , by adjusting the models such that in all cases the average model luminosity agrees with the average observed luminosity, . Figure 4 shows the star formation timescale derived from the observed mean luminosity versus the mean model luminosity obtained above from the observational timescale reported by Evans et al. (2009). The timescales are also listed in Table 3, which includes the dimensionless parameter values for each model. For nonaccelerating models, the mean formation time is the same as the value that would be inferred from observation by comparing the number of protostars with the number of Class II sources. For accelerating models, however, the formation time is significantly less, , for the parameters we have adopted; both values are listed in the table. For models with nonaccelerating accretion, the figure strongly suggests that the actual formation time, during which the majority of a protostar’s mass is accreted, is
(42) 
This is longer than the observed Class 0 lifetime of 0.1 Myr (e.g., Enoch et al. 2009), which implies that significant accretion continues into the Class I phase.
Alternatively, one could use the median rather than the mean to adjust the models. This would result in an inferred distribution of more evenly distributed above and below . For example, the untapered, nonaccelerating CA case would require a time half the length of the observed formation time, while the untapered, nonaccelerating IS case would require a time twice as large. However, adjusting the mean gives better agreement between the overall distributions (see Section 4.5).
The vertical error bars on arise from the uncertainty in the adopted disk lifetime: 2 Myr. Another, smaller source of error arises from the possible misclassification of sources. In comparing the observed formation time with that predicted by the models summarized in Table 3, further uncertainty is introduced by the uncertainty in the mean luminosity of the sample, which has been used to normalize the models. Models without acceleration have an average formation time of about 0.3 Myr, which would require disk lifetimes close to 1 Myr, smaller than most of the results cited by Evans et al. (2009). Models with an acceleration time Myr have an average observed formation time of 0.6 Myr, corresponding to a mean disk lifetime of about 3 Myr, at the high end of the observationally inferred values. Improvements in the accuracy of the measured protostellar lifetimes, acceleration times and luminosities will enable more stringent tests of the models.
4.5. The Protostellar Luminosity Function (PLF)
In this section we present the PLF for each model and define four dimensionless parameters to characterize the PLF shape. In Figures 5 and 6, we plot the predicted PLFs and overlay the PLF of the observed distribution of protostars. We exclude sources with (comparable to the limit of the observations after correction for extinction) from the model distributions and then renormalize the PLF to unity.
In Table 3 we give the standard deviation of log luminosity, the ratio of the median luminosity to the mean, the ratio of the maximum luminosity to the mean, and the fraction of very low luminosity objects (VELLOs; see Section 4.5.4) for each model, renormalized such that the means are equal to the observed mean luminosity.
Model 


(Myr) 


Observed  
Isothermal Sphere (IS)  NonAccelerating  
Untapered  0.44  0.55  0.96  5.81  0.08  0.62  
Tapered 
0.45  0.56  0.83  2.97  0.08  0.72  
Accelerating 

Untapered  0.43  0.54  0.86  2.90  0.07  0.85, 0.37  
Tapered  0.39  0.51  0.71  3.67  0.07  0.85, 0.37  
NonAccelerating  
Untapered  0.55  0.64  0.60  11.80  0.09  0.56  
Tapered  0.41  0.53  0.64  5.23  0.08  0.60  
Accelerating  
Untapered  0.56  0.65  0.26  10.57  0.15  0.99, 0.43  
Tapered  0.42  0.53  0.45  6.74  0.07  1.03, 0.45  
Turbulent Core (TC)  NonAccelerating  
Untapered  0.77  0.84  0.42  9.46  0.18  0.69  
Tapered  0.69  0.76  0.51  4.91  0.21  0.58  
Accelerating  
Untapered  0.76  0.83  0.32  9.78  0.21  1.08, 0.47  
Tapered  0.71  0.78  0.53  6.06  0.12  0.84, 0.37  
NonAccelerating  
Untapered  0.53  0.61  0.30  9.46  0.09  0.55  
Tapered  0.49  0.57  0.66  4.04  0.06  0.51  
Accelerating  
Untapered  0.55  0.63  0.30  12.19  0.11  0.94,0.41  
Tapered  0.48  0.57  0.67  4.02  0.05  0.78,0.34  
Competitive Accretion (CA)  NonAccelerating  
Untapered  0.81  0.94  0.21  12.78  0.25  0.81  
Tapered  0.75  0.89  0.44  5.99  0.22  0.62  
Accelerating  
Untapered  0.80  0.93  0.18  13.41  0.27  1.33, 0.58  
Tapered  0.80  0.93  0.43  6.79  0.17  0.89, 0.39 
Standard Deviation of Log Luminosity,
We find that the IS model has the smallest in all cases and is too narrow with respect to the data. In contrast, the TC and CA models have larger values of and encompass the extent of the data. The 2CTC and 2CCA models present a promising compromise and their are slightly too low. Allowing the accretion rate of the protostars to taper off decreases the standard deviation of the distributions in all but the nonaccelerating IS case. This exception occurs because tapering introduces a spread in the accretion rate.
Allowing for time variations of a factor of two, i.e., nonFUOri fluctuations in the accretion rates, would also increase the standard deviation but have little effect on or . To take such lowlevel variability into account, we add 0.3 dex in quadrature to the measured standard deviations. As shown in Table 3, this permits agreement between observations and both the 2CTC and 2CCA models. The TC and CA models remain within error, while the IS models continue to be inconsistent.
In Figure 8, we plot as a function of the most massive forming star. We plot the tabulated values for the case of and the observed standard deviation ( dex) in inset plots. As illustrated by both Figures 8 and 5, the standard deviations of the four models are very distinct. The TC and CA models grow closer together for larger , but remain fairly well separated at . The standard deviations are also relatively insensitive to the upper stellar mass in the cluster. Tapering reduces the width of the distribution for all for these models.
Median to Mean Ratio
The ratio of the median to the mean, ,
virtually eliminates the
dependence of the results on , , and .
Table 3 gives the values of the ratios for each
of the accretion models, where the observed ratio is .
In all cases, the IS model can be ruled
out.
Maximum Luminosity,
To characterize the maximum luminosity of the distributions, we define:
(43) 
where is the PLF derived using the accretion luminosity. For parity between the observations and models, we compare with the ratio of to . The most luminous observed source has a bolometric luminosity of , which is likely an upper limit and may be over estimated by a factor of two. Consequently, we adopt in comparing with the models. Importantly, this source is a Class 0 object, which suggests that the luminosity chiefly arises from accretion. Thus, we derive from , which assumes that the interior luminosity is small compared to the accretion luminosity. We note that the values of in Table 3 are based on the actual model dispersion, , not the enhanced value that allows for factor of 2 variability; allowance for such variability would increase somewhat.
According to Table 3, the untapered CA models have the highest , followed by the . 2CCA and 2CTC cases and the TC ones. The observed maximum luminosity is times higher than for the tapered PLFs, but it is within error of the untapered 2CTC, 2CCA, TC and CA cases. Without a correction for the model mean luminosities, the difference between the observed and model would be a factor of 56. It is unclear whether this discrepancy is an artifact of the larger error of the higher luminosity measurements, variable accretion, or fluctuations in the star formation rate.
The untapered cases, even though directly consistent with , are discrepant with the observations in a different way. For untapered accretion, in which the accretion rates increase monotonically, the maximum luminosities are achieved only when reaches its final value, . Consequently, in order for the untapered cases to be consistent with observations either the Class I phase must have higher observed luminosities than the Class 0 phase or the Class 0 phase must last much longer is currently inferred from observations. The tapered IS, 2CTC, and TC models, though lower as shown by Table 3, have peak luminosities that occur midway through the formation time. In contrast, the untapered 2CCA and CA models, which agree better with the observations, achieve those luminosities only towards the end of accretion, and thus do not appear to be consistent with observation.
Very LowLuminosity Object Fraction,
There is observational evidence in support of a significant population of low luminosity protostars. For example, Dunham et al. (2008) found that 30% of protostars with have luminosities below 0.1 , a sample commonly referred to as very lowluminosity objects (VELLOs). In our sample, which covers the same regions as Dunham et al. (2008) but has been corrected for dust extinction and more carefully pruned to eliminate nonprotostars (e.g., background galaxies), only % of protostars can be considered VELLOs. Note that applying an exinction correction to the data increases the median luminosity by 40% so we define the VELLO fraction as
(44) 
The observational sample contains one source with an extinctioncorrected luminosity of , but it is likely incomplete at luminosities below ; we therefore set . Without correcting for extinction, 20% of the observed sources would have .
In Table 3 we give for each of the models. The IS models, the tapered 2CTC, and the tapered 2CCA models are inconsistent with the data, whereas the TC and CA models are consistent with the observed values in all cases. The CA, and to a lesser extent, the TC models also predict a substantial number of VELLOs below the luminosity completeness limit. Future observations extending below this limit are necessary to confirm or rule out these models. It should be noted that in all prescriptions, VELLOs are produced not by quiescent periods preceeding episodic accretion events, but by the small accretion rates associated with protostars of very low mass. Current observations cannot discriminate between lowluminosity sources in a quiescent phase and those that are simply lowmass objects with low accretion, although it is likely that at least some of the observed protostars fall in the former category. It is therefore possible that the low VELLO fractions of some of the models may be consistent with the actual fraction of the subset of very lowluminosity protostars that are not undergoing episodic accretion.
Constant Radius PLF
In the Appendix we derive the PLF for each model (except 2CCA) assuming only accretion luminosity and adopting a constant protostellar radius, . Figure 7 shows these constantradius curves in the untapered, nonaccelerating case together with the PLFs for which is a polynomial fit to the subgrid protostellar evolution model (as in Figure 5). Figure 7 illustrates that the curve width and maximum luminosity are reduced when the radius is allowed to depend upon mass. The constant radius curves have standard deviations that are 50% larger than when the radius is allowed to vary, a change that, for example, makes the shapes of the constant radius IS curve and varying radius 2CTC curves similar. This indicates that the stellar evolution model does impact the PLF shape and the details of the comparison. In particular, stellar evolutionary models with a narrower range of will have broader luminosity distributions. Since our stellar evolution model depends upon the instantaneous accretion rate and not just the accretion timescale, a star of the same final mass will have a different amount of deuterium remaining at for the different accretion histories. The models fall in the order IS, 2CTC, 2CCA, TC, and CA from most evolved to least evolved, based upon the amount of deuterium burned at a given final mass. However, we note that the model PLF shapes are mainly distinguished by their characteristic accretion models rather than differences in their stellar evolution.
5. Discussion
A number of uncertainties enter into our comparisons. Foremost, the observational data are difficult to obtain and corrections due to obscuration along the line of sight contribute significant error. There is also uncertainty in several of the parameters we have adopted for the models and in the way these parameters are implemented.
The parameter remains one of the most uncertain in our estimation. It actually contains two very different effects: nonradiative loss of accretion energy from the disk (e.g., Ostriker & Shu, 1995) and advection of accretion energy into the stellar interior (Hartmann et al., 1997). The latter effect reduces the accretion luminosity by a factor , where is the fraction of the accretion energy advected into the stellar interior. Most authors agree that this effect is small (Hartmann et al., 1997; Baraffe et al., 2009), and it was not included in the work of Stahler (1988). In a thorough study of the structure of protostellar accretion shocks, Commerçon et al. (2011) have shown that is indeed extremely small provided the accretion flow is optically thin to optical radiation, so we neglect it here. The accretion shock does heat the surface layers of the protostar, however, and if one assumes that this heating is negligible because the accretion is localized onto a small fraction of the protostellar surface, the resulting protostars can be significantly more compact than when the surface of the protostar is heated by the accretion shock (Hartmann et al., 1997). Such models appear to be inconsistent with observation (Baraffe et al., 2009; Hosokawa et al., 2011).
The fraction of mass accreted during bursts, , is also uncertain. Most bursts are observed to occur in the Class I stage, although Class II objects also experience FUOri type events (e.g., Miller et al. 2010). There is some suggestion that bursts increase with age and that Class 0 objects experience smoother accretion (Vorobyov & Basu, 2005; Zhu et al., 2009). There is also evidence that isolated stars experience more frequent outbursts than those in clusters (Greene et al., 2008), a puzzle that highlights a deficit in our understanding of protostellar accretion processes.
If protostars undergo periodic outbursts, then they must exist in a quiescent, lowluminosity stage between bursts. Indeed, some of these quiescent sources may appear as VELLOs. The median luminosity of the FU Ori sources cataloged by Sandell & Weintraub (2001) is , so protostars that increase in luminosity by more than 8 magnitudes would have been in a VELLO state prior to outburst. It is not known how many of the known FU Ori sources had preburst luminosities below , but several of the best known ones did have preburst luminosities above this limit (Hartmann & Kenyon, 1996) and would not have been classified as VELLOs. If the bursting sources were in a lowluminosity state between bursts, they would populate the lower luminosity region of the observed PLF and thus increase the standard deviation, lower the ratio of the median to the mean, and, for sufficiently faint sources, increase of the observed population relative to our models, which do not explicitly include outbursts. If the interburst accretion luminosity is negligible and the nuclear luminosity is approximated by the ZAMS value, then nonaccreting protostars with masses up to about could potentially be classified as VELLOs. The observed value of is thus an upper limit on the number of sources in that luminosity range in our models. As a result, models like the tapered 2CTC and 2CCA models, which fall below the observed , are promising candidate models, while the untapered TC and CA accretion models may well overpredict and are not as promising.
The ratio of the median and mean luminosities can be significantly reduced by the degree of burstiness in the accretion rate of the individual protostars. The same effect can be achieved in models such as the CA and TC models, which achieve a small ratio by virtue of a relatively long phase of lowluminosity accretion. However, as noted above, these models are somewhat artificial, the former because protostars initially accrete some mass from a local reservoir (which leads to the 2CCA model) and the latter because turbulent motions do not exceed thermal motions in such lowmass cores (which leads to the 2CTC model). Nonetheless, this illustrates that it is possible to fit the observed ratio without episodic accretion, and thus, that it is observationally difficult to disentangle the influence of variability from an underlying mean accretion trend shaping the luminosity distribution.
Although only a small number of sources have been observed to undergo large, FU Ori outbursts, many protostellar objects have been observed to undergo variability in luminosity over a factor of two on timescales of months to years (Pech et al., 2010; Covey et al., 2010). Our correction to the dispersion in log discussed in §4.5.1 reflects this, but is very approximate. In fact, we find that some lowlevel timevariability is needed in order for the 2CTC and 2CCA cases to be consistent with observations.
Another source of dispersion in the PLF could stem from variation in the initial conditions for star formation. For example, temperatures in the lowmass starforming regions we study here are in the range K (McKee & Offner, 2010a), which leads to a factor 3 range of accretion luminosity in the IS model. This corresponds to a dispersion of dex, which has a negligible effect when added in quadrature to the dispersion of 0.3 dex for the assumed temporal variability. When samples with a larger range of initial conditions are considered, it is possible that this additional source of dispersion would have to be included.
Our models do not directly take into account other modes of star formation such as fragmentation within accretion disks (Bate, 2009a; Stamatellos & Whitworth, 2009; Kratter et al., 2010). Division of gas accreting onto a shared disk between close companions is a complicated process, which could potentially alter the the accretion dependence on and that we assume. However, the formulation of the models does not necessarily exclude disk fragmentation since we specify the accretion prescription rather than the protostellar origin. We also do not expect disk fragmentation to be common in this observational sample, since heating due to radiative feedback significantly stabilizes lowmass disks (Offner et al., 2009; Bate, 2009b). In the TC and CA senarios, core and filament fragmentation may produce wide binary companions (Offner et al., 2010), but accretion rates for this mode of origin should follow the expected accretion trends. We have not taken binarity into account in our analysis, but its effects are small compared to the large differences among the different accretion models; for example, the difference between Chabrier’s (2005) IMF for individual stars and that for stellar systems is only a factor of 1.25 in the peak mass.
In addition to the accretion history and to uncertainties in parameters, inclination effects may also contribute to the standard deviation of the observed PLF. For example, the extinction correction and therefore the bolometric luminosity of a protostar with a disk that is observed edgeon will be underestimated. Dunham et al. (2010) found that adding inclination effects to their models broadened the bolometric luminosity distribution by %, so that orientation effects alone were not able to fit the data. Moreover, including inclination effects had only a small effect on the high luminosities, at most increasing the maximum by 25%. Consequently, we expect geometry to have a minor effect on the shape of the PLF.
Adopting an accelerating star formation rate is one possible resolution of the apparent discrepancy between the observed star formation timescale of 0.44 Myr, which is based on the assumption of a constant starformation rate (Evans et al., 2009), and the mean luminosity, which suggests a starformation timescale of about 0.3 Myr (Sec. 4.4). We have adopted Myr as the acceleration time, which is the inferred acceleration time for Ophiuchus and is likely to be a lower bound on the acceleration time for the other regions (Palla & Stahler, 2000; McKee & Offner, 2010b). (More recent Ophiuchus data suggest that the region has a deficit of Class 0 objects and therefore a decelerating rate of star formation (Evans et al., 2009).) Insofar as 1 Myr is a lower bound on the acceleration time, it gives the maximum difference between and , and thus the maximum effect of acceleration on protostellar luminosities. IC348, a subregion of Perseus, has an inferred accleration time of 2 Myr (Palla & Stahler, 2000; McKee & Offner, 2010b), but the accelerations for the entire Perseus region and for Serpens are unknown. Given the small statistical sample size, the appreciable uncertainties in the timescale estimates, and our simplified acceleration model, these results should be interpreted as being consistent with, but not demonstrating that, acceleration resolves the apparent discrepancy between the observed starformation time and the mean luminosity by accelerating star formation.
How do the different accretion models compare with the data? Independent of parameter uncertainties, the IS models appear to be inconsistent with observations. Of the five metrics we have adopted—, , , , and —only the last is within the estimated errors. The ratio of the median to mean luminosity is more than 2 standard deviations from the observed value. (However, the error estimates here are quite uncertain, and future data should substantially improve them.) Because protostars accreting according to the TC and CA prescriptions spend significant time at low masses, where ISlike accretion is likely to occur, we believe that the 2CTC and 2CCA models are more realistic physically, even though in some cases the TC and CA models agree quite well with the data. For example, in the CA case lowmass stars spend of their accretion time achieving a mass of 0.1 . It is unlikely that protostars spend this much time at accretion rates significantly less than the IS value. The 2CTC and 2CCA models are thus likely to be a better representation of the TC and CA models, since the constant accretion component dominates the initial protostellar accretion rate.
Likewise, with the exception of the CA models, which are assumed to have accretion terminated over a short time by protostellar feedback, tapered models are likely to be more realistic since accretion from a core declines gradually after the expansion wave reaches the surface of the core (e.g., McLaughlin & Pudritz, 1997). Furthermore, untapered models achieve their maximum luminosity at the end of the protostellar stage, whereas observations show that Class I sources are not noticeably more luminous than Class 0 ones; in addition, the accretion rate for the Class 0 sources is believed to be larger than for the Class I sources, which is consistent with tapering. In our comparison we find that in all cases the tapered models tend to underestimate the maximum luminosity and the standard deviation of the distribution. This may suggest that our parametrization of tapering or our specific choice of is not correct rather than ruling out tapered protostellar accretion rates. It must also be borne in mind that we have not included the effect of variability on , so the tabulated values are lower bounds on the true values.
As a check on our tapering model, we derive tapered PLFs for the IS and TC cases assuming an exponential functional form for the accretion rate:
(45) 
For the IS case, this is similar to the model of Myers et al. (1998), except that in their model approaches as . We modify their model so that at , at which point the accretion rate has declined exponentially by a factor of . For the IS case, we find that , , and change by only a few percent relative to the linearly tapered case. The parameter decreases by 50%, mainly because the exponentially tapered accretion rate does not go to as the linearly tapered case does. Thus, adopting a different tapering model does not alter the fundamental disagreement of the IS case with the observations. For the exponentially tapered TC case, and change by a few percent, while and increase by and %, respectively, and decreases by 40%. Although these changes are more significant, the dimensionless metrics continue to agree with the observations as before, except that becomes slightly too high. While it may be possible to formulate accretion tapering that improves the agreement with the observations, it seems unlikely that a different functional form would modify our conclusion that constant accretion time models agree better with the data than constant accretion rate models.
The comparison with observations could be strengthened by deriving PLFs for the Class 0 and Class I populations separately. In theory, the division between the two classes occurs when the protostellar mass exceeds the envelope mass (Andre et al., 2000; Crapsi et al., 2008). For the IS and TC models, it is straightforward to define Class 0 and Class I PLFs as the distributions of luminosities for which and , respectively. Although such distributions can also be constructed for the CA case using our approximate CA model, one characteristic of competitive accretion is the lack of correspondance between envelope mass and final stellar mass (Smith et al., 2009). Even using this simple definition, comparing the Class PLFs quantitatively with the observations is challenging. Since the protostellar mass is not measurable during the embedded phase and since inclination effects can significantly confuse the classification, the two populations are difficult to distinguish observationally. Because the Class I lifetime is about three times the Class 0 lifetime (Evans et al., 2009), models that do not have a higher rate of accretion as Class 0 sources could be excluded. If the fraction of a protostellar lifetime spent as a Class 0 source were to be less than that found by Evans et al. (2009), the constraint on accretion and early luminosities would become more stringent. Qualitatively, we expect all the untapered models to fail such a test, since accretion, and hence luminosity, in such models rises monotonically with age.
6. Conclusions
We have analyzed the Protostellar Luminosity Function (PLF), which is the presentday luminosity distribution of a cluster of protostars, for several different theories of star formation and compared the results with observation. In our derivation we have assumed that the protostellar masses evolve smoothly onto a truncated Chabrier (2005) IMF, that the accretion rates are a continuous function of the instantaneous and final protostellar masses, and that the star formation rate is either constant or accelerating in time. We also assume that over most of the formation time, the accretion rate onto the protostar tracks the accretion rate from the ambient molecular core onto the protostardisk system. Episodic accretion, such as that occurring in FU Ori outbursts, violates this assumption, but the available observational data indicate that this is not a dominant effect.
The PLF depends explicitly upon the mean formation time of the stars, , and the maximum stellar mass produced, . For the lowmass starforming regions we have compared with, appears to be about . We consider three main accretion prescriptions corresponding to standard models of star formation: Isothermal Sphere (IS, constant accretion rate), Turbulent Core (TC), and Competitive Accretion (CA, constant accretion time). We note that prior to the development of either the CA or TC models, Kenyon et al. (1990) considered both constant accretion rate and constant accretion time in the paper that introduced the luminosity problem. We also consider two hybrid models: the TwoComponent Turbulent Core model (2CTC, a compromise between the IS and TC models) and the TwoComponent Competitive Accretion model (2CCA, a combination of the IS and CA accretion prescriptions; this model was not considered in Paper I). We explore two variations on these models: a case in which the accretion rate smoothly tapers to zero and a case in which the star formation rate accelerates with a characteristic timescale, .
The CA model used here is an approximate analytic representation of the competitive accretion model developed by Bonnell et al. (1997), which begins with protostellar seeds that are produced by a process similar to that in the IS case. As a result, we believe the 2CCA model is a better approximation to their work than the CA model. The TC model was specifically formulated for highmass stars, which we do not focus on here. For the case of lowmass stars, McKee & Tan (2003) proposed the inclusion of an IS stage, which suggests that the 2CTC model is the best representation of their model in this context. We note that this model has some similarities to the TNT model of Fuller & Myers (1993).
We compare our models to protostellar luminosities observed in local lowmass star forming regions (Evans et al. 2009, Enoch et al. 2009). The classical luminosity problem is that observed protostars appear to have luminosities significantly lower than expected theoretically (Kenyon et al., 1990). The extinctioncorrected sample that we adopt from Evans et al. (2009) has a mean luminosity of , which is more than a factor of two larger than the earlier Enoch et al. (2009) results that did not account for extinction. This alone signficantly ameliorates the luminosity problem. In comparing the models and the data, we first used the mean luminosity and the median luminosity. In all permutations of the parameters, the models require that the average protostellar lifetime be Myr, somewhat less than Myr measured by Evans et al. (2009). That is, the model luminosities are too low if the mean lifetime inferred by Evans et al. (2009) is assumed. Thus, with a star formation time of Myr, allowance for nonradiative energy loss in winds () and a modest amount of episodic accretion () is sufficient to lower the mean protostellar luminosity so that there is no longer a “luminosity problem,” in lowmass star formation. We note that this resolution of the luminosity problem is quite consistent with the suggestions of Kenyon et al. (1990), who first pointed out the existence of the problem: among the solutions they proposed were that the formation time was longer than the yr indicated by their data and that some of the accretion was episodic. With a starformation time of 0.3 Myr, the mean accretion rate for a star of mass (the mean mass of the Chabrier, 2005 IMF) is yr.
The discrepancy between our estimate of the mean star formation time and that of Evans et al. (2009), while not large, could be due to a number of factors: the disk lifetime could be shorter than they assumed, the number of Class I sources could less than they assumed, or could be larger than we assumed, and/or the star formation rate could be accelerating as suggested by Followup studies of dense gas tracers, such as HCO, have found that as many as half of previously identified Class I objects are not actually embedded van Kempen et al. (2009); Heiderman et al. (2010). Although Evans et al. (2009) exclude sources embedded in less than 0.1 of gas mass when deriving the lifetime, some sources may nonetheless be misclassified older, heavily obscured objects. Also, while our value of is calculated using all the known bursting sources, deeply embedded FUOri type sources may be missing from the sample. By increasing the ratio of protostars to Class II sources, accelerating star formation produces a longer observationally inferred star formation time than the actual star formation time of the sample. Consequently, the accelerating models have better consistency between the observed star formation time and mean and median luminosities. Although we do expect the star formation rate to vary with time, Palla & Stahler (2000) found evidence for an acceleration time as short as our adopted value, Myr, for only one the observed regions in the protostellar sample (Ophiuchus), and this evidence is contradicted by an apparent deficit of Class O sources in the region.
We then compared four dimensionless quantities that characterize the shape of the PLF: (1) the standard deviation of log (including an allowance for source variability at the factor of two level); (2) the ratio of the median to the mean luminosity; (3) the ratio of the maximum to the mean luminosity; and (4) the fraction of very lowluminosity objects, , defined as the ratio of the number of sources with extinctioncorrected luminosities between and to the number between and . The first three of these are the most strongly discriminating since they are independent of , and the mean protostellar lifetime, ; the fourth quantity, , is only weakly dependent on these factors. We also compared the value of the starformation time required by the models to get the observed mean protostellar luminosity with the starformation time of Myr inferred by Evans et al. (2009) from the ratio of the number of protostars to the number of Class II sources, which were assumed to have a 2 Myr lifetime. Although protostars with different accretion histories have slightly different stellar evolutionary states at the end of accretion, we find that differences between the PLF shapes are driven by the different accretion histories, not the different evolutionary states. We find that the IS model is a poor fit to the data in all cases, mainly due to the strongly peaked nature of its PLF profile. The model results for the four dimensionless parameters are outside the error bars of the data regardless of whether the model is tapered or untapered, accelerating or nonaccelerating. The one parameter the IS model agrees with is the observed starformation time, after renormalizing the accretion rate so that the mean luminosity agrees with observation.
The CA model is in best agreement with the data; only the observed starformation time, , lies outside the error bars, and this is only for the nonaccelerating, untapered case. The CA model predicts a relatively large number of VELLOs below the observational limit of , which will provide a strong test of the model in the future. However, as discussed above, the 2CCA model, which has an initial phase in which the star accretes more rapidly, is closer to the actual competitive accretion model. Furthermore, one of the assumptions of the competitive accretion model is that the gas is cleared out relatively quickly toward the end of the accretion phase, so the untapered version of the 2CCA model is closest to the actual competitive accretion model. However, this model (as well as the untapered CA model) predicts that the maximum luminosity is achieved at late times, not during the Class 0 stage. In general, the length of the Class 0 lifetime and the similarity of Class 0 and Class I luminosities suggests that models that do not accrete a significant fraction of mass during the earliest times are inconsistent. This includes all the untapered models, which acheive their maximum accretion rate, and hence maximum luminosity, at the end of the protostellar lifetime, in what would be the late Class I stage. Otherwise, both the accelerating and nonaccelerating untapered versions of the 2CCA model agree well with the data, although the dispersion is very slightly below the observed value for the nonaccelerating case. The tapered version of the 2CCA model does not compare well with the data.
The TC models are in good agreement with the data, with the exception of the nonaccelerating, tapered case, which has a maximumtomean luminosity ratio that is marginally too low. The observed starformation time, , is marginally too high for the nonaccelerating, untapered case. As for the CA model, however, the 2CTC model, which is similar to the IS model at low masses, is more realistic. The untapered version of this model agrees well with the data, except that the median luminosity is somewhat high for the nonaccelerating case. However, the tapered version of the 2CTC model is the best representation of the model, since the model is essentially a turbulent version of the IS model. The nonaccelerating version of this model marginally agrees with all the data except for ; the accelerating version underpredicts the number of VELLOs (although as remarked above, this may not be a problem if some of the observed VELLOs are episodic sources in a quiescent phase).
We conclude that models that tend towards a constant accretion time and thus produce a greater spread in luminosities (like the CA and TC models), rather than models that have a constant accretion rate (such as the IS model) are in better agreement with the data on the PLF. Ultimately, agreement between a model and the observed PLF is necessary but not sufficient. Models must also reproduce a number of other observed features of protostars and the regions from which they form, including core properties (e.g. Offner et al. 2008; Kirk et al. 2009), the approximate agreement between the luminosities of Class 0 and Class I sources, and the existence of large disks around Class I sources. The CA model, which exhibits good agreement with the PLF, does not do well with any of these additional features; in particular, it does not produce welldefined cores for individual protostars (Enoch et al., 2008). The IS and TC models naturally have cores, and the tapered models can give comparable luminosities for the Class 0 and Class I stages. Nonmagnetic IS and TC models are predicted to have large disks, but the existence of such disks in the presence of magnetic fields is a topic of active investigation (e.g., Mellon & Li, 2009; Ciardi & Hennebelle, 2010).
In this paper, we have developed the PLF as a tool for confronting star formation theories with observation. Ongoing observational efforts should permit a significant improvement in the comparison between theory and observation by providing larger samples, which would reduce statistical fluctuations, and more accurate extinctioncorrected luminosities, which would reduce the current factor of 2 uncertainties. The sample of protostars we have analyzed has an estimated maximum mass ; a larger sample would presumably include more massive stars and enable a stronger test of the theories. If the sample were sufficiently large, one would be able to directly determine the role of large, FU Ori type outbursts on the growth of protostars. Study of the very faint protostars, the VELLOs, should determine the relative proportion of those that are in a quiescent state between outbursts and those that are faint because they have very low mass (as we have assumed). A monitoring program would permit one to characterize the variabliity of the protostars, which we have simply taken as increasing the dispersion in the PLF by a factor of two. Finally, more accurate measurements of the physical conditions in the starforming clumps would enable more accurate theoretical predictions of the accretion histories of the protostars in the sample.
Appendix A PLF Alternative Derivation
The PLF may be obtained either by integrating over or over , where in either case the resulting PLFs are identical. Above we exclusively use the former formalism. For completeness, we give the latter PLF definition here:
(A1)  
(A2) 
where the lower limit of integration is given by equation (5) with . Note that this formulation does not work when the luminosity is independent of the final mass, as in the case of Isothermal Accretion.
Appendix B The PLF for Constant Radius
Accurate evaluation of the PLF requires allowing for the dependence of the protostellar radius on the mass and accretion rate. However, the radius is almost always within a factor 2 of , and if we take to be constant the analysis is simplified considerably.
Combining equations (2) and (20), we express the accretion luminosity as
(B1) 
where
(B2)  
(B3) 
is the luminosity for and in the untapered case. If we let
(B4) 
(where is in units of ), then the relation of and to is
(B5) 
b.1. Isothermal Sphere Accretion
In this case, we have , so that equation (29) gives