The Protostellar Luminosity Function

# The Protostellar Luminosity Function

## Abstract

The protostellar luminosity function (PLF) is the present-day 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 low-mass star formation suggest a mean star formation time of 0.30.1 Myr. Such a timescale, together with some accretion that occurs non-radiatively and some that occurs in high-accretion, episodic bursts, resolves the classical “luminosity problem” in low-mass 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.

stars: formation stars: luminosity function, mass function

## 1. 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 star-forming 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). High-resolution millimeter emission maps tracing dusty envelopes, provide a measure of core masses (Enoch et al., 2008). Combined with mid to far-infrared 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 sub-millimeter 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 pre-main sequence star has accreted or expelled most of the initial envelope mass and produces little sub-millimeter 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. Edge-on 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 non-embedded, Class II objects (Enoch et al., 2008), although objects with thick disks along the line-of-sight 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

 Lacc = facc(Gm˙mr), (2) = 7.8facc(m0.25M⊙)× (˙m2.5×10−6M⊙yr−1)(2.5 R⊙r) L⊙

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 in-falling 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 time-dependent luminosity function for a cluster of evolving protostars. Their model included the luminosity contribution from Class I to main-sequence stars. For steady-state star formation, our constant accretion-rate model produces a protostellar protostellar luminosity function similar to that derived by (Fletcher & Stahler, 1994b) when pre-main sequence stars are excluded (see §3).

Myers et al. (1998) used a semi-analytic 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 zero-age 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 low-luminosity sources were not addressed.

Monte-Carlo 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 line-of-sight. 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 1-D radiative transfer observations of the inside-out 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 low-mass star formation.

Recent extensive surveys of five local star-forming 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, 2-D 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 low-luminosity sources on the Hertzsprung-Russell diagram and thus account for a portion of the age spread inferred for members in young, low-mass 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 time-varying 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 time-varying 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 protostar-disk system over the majority of the lifetime of the protostar. Finally, we adopt an individual-star 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 star-forming 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 bi-variate 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., non-accelerating) star formation

 ψp2(m,mf)=mψ(mf)˙m⟨tf⟩, (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 bi-variate PMF,

 ψp(m)=∫mumf,ℓψp2(m,mf)dlnmf, (4)

where the lower limit of integration is

 mf,ℓ=max(mℓ,m) (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 protostar-disk system. We considered several different models:

• Inside-out collapse of an isothermal sphere (IS, Shu 1977),

 ˙m=˙mIS=1.54×10−6(T/10K)3/2. (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,

 ˙m = 3.6×10−5Σ3/4cl(mmf)jmf3/4, (7) ≡ ˙mTC(mmf)jmf3/4, (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 inside-out collapse. We retain the value given by McKee & Tan (2003), but note that the overall normalization of the accretion rate remains uncertain.

• A simplified model to represent competitive accretion (CA, Bonnell et al. 1997, 2001),

 ˙m=˙m1(mmf)2/3mf, (9)

where is the accretion rate for . Note that this accretion rate has the property that stars of all masses form in the same amount of time.

All these forms for the accretion history can be summarized in the expression

 ˙m=˙m1(mmf)jmfjf, (10)

where the model exponents are given by:

 IS: j=jf=0 (11) TC: j=12,jf=34 (12) CA: j=23,jf=1. (13)

We also considered a two-component turbulent core model (2CTC), a blend of the IS and TC models:

 ˙m=˙mIS[1+R2˙m(mmf)2jmf3/2]1/2, (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 two-component 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 Shu-like 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:

 ˙m=˙mIS[1+R2˙m,CA(mmf)4/3mf2]1/2, (15)

where

 R˙m,CA≡˙mCA˙mIS. (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 single-component models, the time to form a star, is a power-law function of the final mass,

 tf=tf1mf1−jf, (17)

where

 tf1=1(1−j)˙m1 (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:

 ˙m=˙m0(m,mf)[1−(ttf)], (19)

where is the untapered accretion rate for . We can then express the accretion rate for both the tapered and untapered cases as

 ˙m=˙m0,1(mmf)jmfjf[1−δn1(mmf)1−j]1/2 (20)

where is the untapered accretion rate for and

 δn1={0untapered,n=0,1tapered,n=1, (21)

The star formation time for single-component models can then written as:

 tf=tf1mf1−jf(1+δn1). (22)

In the case of accelerating star formation we assume a birthrate that increases exponentially with time, i.e.,

 ˙N∗∝e(t−tm)/τ, (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,

 ⟨tf⟩obs = ⟨tII⟩⟨1−e−tf/τ⟩⟨etf/τ⟩(1−e−tII/τ) (24) ≃ [⟨tII⟩τ(1−e−⟨tII⟩/τ)]⟨tf⟩, (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 non-accreting 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 bi-variate 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 :

 Ψp2(L,mf)dlnLdlnmf=ψp2(m,mf)dlnmdlnmf, (26)

where was determined in Paper I (see eq. 3). The PLF is then obtained by integrating over ,

 Ψp(L) = ∫dlnmfΨp2(L,mf), (27) = ∫mumf,ℓ(L)dlnmfψp2[m(L),mf]∣∣∣∂lnL∂lnm∣∣∣. (28)

In the Isothermal Sphere case, the accretion rates are independent of the final mass. Consequently, the PLF reduces to:

 Ψp(L)→m(L)⟨mf⟩∫mumf,ℓdlnmfψ(mf)∣∣1−∂lnr∂lnm∣∣. (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 bi-variate luminosity function and define the mean accretion luminosity: (26):

 ⟨Lacc⟩ = ∫mumℓdlnmf∫mf0dlnmLaccψp2, (30) = 1⟨tf⟩∫mumℓdlnmfψ(mf)[faccGmf22¯r(mf)], (31)

where equation (26) enabled us to eliminate the dependence on for steady star formation and where

 1¯r(mf)≡2mf2∫mf0mdmr(m,mf) (32)

is the mass-averaged 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 ,

 1¯r(m)=2mu2−m2∫mumf,lmfdmfr(m,mf). (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 Kelvin-Helmholtz 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, low-mass, 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:

 σ(logL) = [∫(logL)2Ψp(L)dlnL (34) Missing or unrecognized delimiter for \right

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 one-zone 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 sub-grid 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“pre-collapse” stage, the collapsing gas densities are sub-stellar. 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 zero-age main sequence luminosity serves as an approximation of the interior luminosity, which results from the combination of deuterium burning and Kelvin-Helmholz 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 re-normalize the result using the directly derived value of . This strategy simplifies the integral while preserving the trends in the sub-grid model and eliminating the weak dependence of on .

### 3.4. Episodic Accretion

Observations of high-luminosity 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 low-luminosity, low-accretion 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 20-550 , 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:

 fepi=˙mepiΔtepi⟨mf⟩, (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:

 Δtepi=Np,epi˙N∗, (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 time-averaged 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

 facc,eff=facc(1−fepi)=0.56. (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 non-protostars and non-embedded 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 15-25% 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 non-extinction 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:

 ⟨Lobs⟩ = 5.3+2.6−1.9 L⊙, (38) Lobs,med = 1.5+0.7−0.4 L⊙. (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 one-sigma error estimate. Since the standard deviation is also sensitive to the distribution outliers, we use the log of the standard deviation for comparison:

 σ(logLobs)=0.7+0.2−0.1 dex. (40)

A second useful non-dimensional quantity is the ratio of the median luminosity to the mean:

 Lmed/⟨L⟩=0.3+0.2−0.1. (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 Kolmogorov-Smirnov (K-S) 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.

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, in-depth study often reveals that sources classifed as Class I are actually edge-on 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, ⟨L⟩

The mean of the luminosity distribution serves as a simple one-dimensional 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 high-mass star formation. In Table 1, we give the model values calculated assuming completeness of the extinction-corrected 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 high-mass star-forming regions would be useful to distinguish models based solely upon the mean luminosity. Note that the model means depend upon the uncertain star-formation 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 non-accelerating cases cases. The untapered TC and CA means, IS tapered mean, and all accelerating means are consistent with the observational error.1 This clearly indicates that there is no luminosity problem in the traditional sense. The resolution is a result of the longer protostellar lifetime adopted from Evans et al. (2009) and an effective accretion efficiency, due to a radiative efficiency of 75% and allowance for episodic accretion at the level of 25%. Altogether this reduces the predicted luminosities for the non-accelerating cases by a factor of .

The mean luminosities in the accelerating cases are up to 30% lower than in the fiducial non-accelerating cases for a fixed value of because the accelerating cases have more low-mass 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 re-normalized using the observational completeness limit, the differences are minimized and in the TC and CA cases the mean increases. Among the tapered, non-accelerating 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 over-estimation 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.

### 4.4. The Star Formation Timescale, ⟨tf⟩

We see that in some cases the mean and median luminosities are quite different from the observed values when we fix the star-formation 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 non-accelerating 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 non-accelerating accretion, the figure strongly suggests that the actual formation time, during which the majority of a protostar’s mass is accreted, is

 ⟨tf⟩=0.3±0.1   Myr. (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, re-normalized such that the means are equal to the observed mean luminosity.

#### Standard Deviation of Log Luminosity, σ(logL)

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 non-accelerating 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., non-FU-Ori fluctuations in the accretion rates, would also increase the standard deviation but have little effect on or . To take such low-level 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.18 The non-accelerating 2CTC and the tapered 2CCA models are outside the observational error bars.

#### Maximum Luminosity, Lmax

To characterize the maximum luminosity of the distributions, we define:

 Lmax=∫LuL′Ψ(Lacc) dlnL=1/N∗, (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 5-6. 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 Low-Luminosity Object Fraction, fVELLO

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 low-luminosity 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 non-protostars (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

 fVELLO=N∗(Lmin≤L≤0.14L⊙)N(Lmin≤L≤1.4L⊙). (44)

The observational sample contains one source with an extinction-corrected 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 low-luminosity sources in a quiescent phase and those that are simply low-mass 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 low-luminosity 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 constant-radius curves in the untapered, non-accelerating case together with the PLFs for which is a polynomial fit to the sub-grid 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: non-radiative 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 FU-Ori 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, low-luminosity 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 pre-burst luminosities below , but several of the best known ones did have pre-burst luminosities above this limit (Hartmann & Kenyon, 1996) and would not have been classified as VELLOs. If the bursting sources were in a low-luminosity 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 inter-burst accretion luminosity is negligible and the nuclear luminosity is approximated by the ZAMS value, then non-accreting 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 over-predict 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 low-luminosity 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 low-mass 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 low-level time-variability 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 low-mass star-forming 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 low-mass 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 edge-on 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 star-formation rate (Evans et al., 2009), and the mean luminosity, which suggests a star-formation 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 sub-region 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 star-formation 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 IS-like 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 low-mass 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:

 ˙m=˙m0(m,mf)exp(−2t/tf)    (t≤tf). (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 present-day 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 protostar-disk 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 low-mass star-forming 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 Two-Component Turbulent Core model (2CTC, a compromise between the IS and TC models) and the Two-Component 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 high-mass stars, which we do not focus on here. For the case of low-mass 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 low-mass 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 extinction-corrected 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 non-radiative 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 low-mass 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 star-formation 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 Follow-up 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 FU-Ori 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 low-luminosity objects, , defined as the ratio of the number of sources with extinction-corrected 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 star-formation time required by the models to get the observed mean protostellar luminosity with the star-formation 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 non-accelerating. The one parameter the IS model agrees with is the observed star-formation 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 star-formation time, , lies outside the error bars, and this is only for the non-accelerating, 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 non-accelerating untapered versions of the 2CCA model agree well with the data, although the dispersion is very slightly below the observed value for the non-accelerating 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 non-accelerating, tapered case, which has a maximum-to-mean luminosity ratio that is marginally too low. The observed star-formation time, , is marginally too high for the non-accelerating, 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 non-accelerating 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 non-accelerating 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 well-defined 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. Non-magnetic 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 extinction-corrected 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 star-forming clumps would enable more accurate theoretical predictions of the accretion histories of the protostars in the sample.

We thank Steven Stahler, Jean Francoise Lestrad, Mark Krumholz, Lee Hartmann, Scott Kenyon, Joan Najita, Phil Myers and an anonymous referee for helpful comments. We thank Melissa Enoch, Neal Evans, and Mike Dunham for useful discussions and clarifications of the observational data. This research has been supported by the NSF through grants AST-0908553 (CFM) and AST-0901055 (SSRO).

## 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:

 Ψp(L) = ∫dlnmΨp2(L,m), (A1) = ∫mmaxmmindlnmψp2[m,mf(L,m)]∣∣∣∂lnL∂lnmf∣∣∣, (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

 Missing or unrecognized delimiter for \left (B1)

where

 L(1) = faccG˙m0,1(M2⊙R⊙yr), (B2) = 31.3facc[˙m0,11×10−6M⊙yr−1]   L⊙, (B3)

is the luminosity for and in the untapered case. If we let

 ℓ≡rLL(1) (B4)

(where is in units of ), then the relation of and to is

 ℓ=m1+jmfjf−j[1−δn1(mmf)1−j]1/2. (B5)

### b.1. Isothermal Sphere Accretion

In this case, we have , so that equation (29) gives

 Ψp(L)=ψp(m=ℓ)=ℓ⟨mf⟩∫mumax(mℓ,ℓ)dlnmfψ(m