RESOLVING THE LUMINOSITY PROBLEM IN LOW-MASS STAR FORMATION

2.1.1. Basic Processes and Equations

We make use of numerical hydrodynamics simulations in the thin-disk approximation to compute the gravitational collapse of rotating, gravitationally bound, pre-stellar cores. This approximation is an excellent means to calculate the evolution for many orbital periods and many model realizations. It is valid as long as the aspect ratio A of the disk scale height h to radius r is well below unity, which is usually fulfilled in the inner 1000 AU (see Figure 11 in Vorobyov & Basu 2010). Protostellar disks rarely exceed 1000 AU in radius (Vicente & Alves 2005; Vorobyov 2011), making our approach well justified for the purpose of collecting a wide statistical sample of model disks. At larger radial distances, A may approach unity but those regions are not important dynamically and serve as a reservoir for material falling onto the disk during the embedded phase of star formation.

We start our numerical integration in the pre-stellar phase characterized by a collapsing starless core, continue into the embedded phase of star formation during which a star and disk are formed, and terminate our simulations in the T Tauri phase when most of the core has accreted onto the forming star+disk system. The thin-disk approximation allows us to consider spatial scales of order 10,000 AU and compute the evolution of both the core and the star+disk system altogether. An important consequence is that the mass accretion rate onto the disk ($\dot{M}_d$) is not a free parameter of our model but is self-consistently determined by the gas dynamics in the collapsing core.

We introduce a "sink cell" at R_sc = 6 AU and allow for the matter in the computational domain to freely flow into the sink cell but not out of it. We monitor the gas density in the sink cell and when its value exceeds a critical value for the transition from isothermal to adiabatic evolution (∼10¹¹ cm⁻³), we introduce a central gravitating point-mass star. In the subsequent evolution, most of the gas that crosses the inner boundary is assumed to land onto the central star while a small fraction is assumed to be carried away with protostellar jets, with the exact partition being a free parameter (usually 10% is assumed to be ejected).

The main physical processes that are taken into account in our modeling include stellar irradiation, background irradiation with temperature T_bg = 10 K, viscous and shock heating, radiative cooling from the disk surface, and also disk self-gravity. The corresponding equations of mass, momentum, and energy transport are

$$\begin{equation} \frac{{\partial \Sigma }}{{\partial t}} = - \nabla _p \cdot (\Sigma \mbox{\boldmath $ v $}_p ), \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial }{\partial t} (\Sigma \mbox{\boldmath $ v $}_p ) + [ \nabla \cdot (\Sigma \mbox{\boldmath $ v_p $} \otimes \mbox{\boldmath $ v $}_p ) ]_p = - \nabla _p {\cal P} + \Sigma \, \mbox{\boldmath $ g $}_p + (\nabla \cdot \mathbf {\Pi })_p, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial e}{\partial t} +\nabla _p \cdot (e \mbox{\boldmath $ v $}_p ) = -{\cal P} (\nabla _p \cdot \mbox{\boldmath $ v $}_{p}) -\Lambda +\Gamma + \left(\nabla \mbox{\boldmath $ v $}\right)_{pp^\prime }:\Pi _{pp^\prime }, \end{equation} \tag{ 3 }$$

where subscripts p and p' refer to the planar components (r, ϕ) in polar coordinates, Σ is the mass surface density, e is the internal energy per surface area, ${\cal P}=\int ^{h}_{-h} P dh$ is the vertically integrated form of the gas pressure P, h is the radially and azimuthally varying vertical scale height determined in each computational cell using an assumption of local hydrostatic equilibrium, $\mbox{\boldmath $ v $}_{p}=v_r \hat{\mbox{\boldmath $ r $}}+ v_\phi \hat{\mbox{\boldmath $ \phi $}}$ is the velocity in the disk plane, $\mbox{\boldmath $ g $}_{p}=g_r \hat{\mbox{\boldmath $ r $}} +g_\phi \hat{\mbox{\boldmath $ \phi $}}$ is the gravitational acceleration in the disk plane, and $\nabla _p=\hat{\mbox{\boldmath $ r $}} \partial / \partial r + \hat{\mbox{\boldmath $ \phi $}} r^{-1} \partial / \partial \phi$ is the gradient along the planar coordinates of the disk. Viscosity enters the basic equations via the viscous stress tensor $\mathbf {\Pi }$ and the magnitude of kinematic viscosity ν is parameterized using the usual α-prescription of Shakura & Sunyaev (1973). In our models, we use a spatially and temporally uniform α, with its value set to 5 × 10⁻³.

Apart from viscous heating determined by the viscous stress tensor $\mbox{\boldmath $ \Pi $}$, the thermal balance of the disk is also controlled by heating due to stellar and background irradiation and radiative cooling from the disk surface. The latter is calculated using the diffusion approximation of the vertical radiation transport in a one-zone model of the vertical disk structure (Johnson & Gammie 2003):

$$\begin{equation} \Lambda ={\cal F}_{\rm c}\sigma \, T^4 \frac{\tau }{1+\tau ^2}, \end{equation} \tag{ 4 }$$

where σ is the Stefan–Boltzmann constant, T is the midplane temperature of gas, and ${\cal F}_{\rm c}=2\,+\,20\tan ^{-1}(\tau)/(3\pi)$ is a function that secures a correct transition between the optically thick and thin regimes. We use frequency-integrated opacities τ of Bell & Lin (1994).

The heating function is expressed as

$$\begin{equation} \Gamma ={\cal F}_{\rm c}\sigma \, T_{\rm irr}^4 \frac{\tau }{1+\tau ^2}, \end{equation} \tag{ 5 }$$

where T_irr is the irradiation temperature at the disk surface determined by the stellar and background blackbody irradiation as

$$\begin{equation} T_{\rm irr}^4=T_{\rm bg}^4+\frac{F_{\rm irr}(r)}{\sigma }, \end{equation} \tag{ 6 }$$

where T_bg is the uniform background temperature (in our model set to the initial temperature of the natal cloud core) and F_irr(r) is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central star. The latter quantity is calculated as

$$\begin{equation} F_{\rm irr}(r)= \frac{L_\ast }{4\pi r^2} \cos {\gamma _{\rm irr}}, \end{equation} \tag{ 7 }$$

where the incidence angle of radiation γ_irr arriving at the disk surface at radial distance r is calculated using the model's known radial profile of the disk scale height h, and L_* is the sum of the accretion luminosity L_{*, accr} arising from the gravitational energy of accreted gas and the photospheric luminosity L_{*, ph} due to gravitational compression and deuterium burning in the star interior. While the former quantity is calculated from the model's known stellar mass and accretion rate onto the star, the latter is taken from the pre-main-sequence tracks of D'Antona & Mazzitelli (1994). A more detailed explanation of the model can be found in Vorobyov & Basu (2010).

2.1.2. Initial Conditions

The form of the initial gas density, temperature and angular velocity profiles in gravitationally bound, dense pre-stellar cores may vary depending on the environment and physical processes that contribute to the formation of such cores. In this work, we consider two possible gas surface density Σ and angular velocity Ω profiles and assume that the initial gas temperature in the core is equal to T_init ≡ T_bg = 10 K. Most models' cores have Σ and Ω typical for magnetically supercritical cores formed by slow gravitational contraction (Basu 1997):

$$\begin{equation} \Sigma ={r_0 \Sigma _0 \over \sqrt{r^2+r_0^2}}\:, \end{equation} \tag{ 8 }$$

$$\begin{equation} \Omega =2\Omega _0 \left({r_0\over r}\right)^2 \left[\sqrt{1+\left({r\over r_0}\right)^2 } -1\right], \end{equation} \tag{ 9 }$$

where Ω₀ is the central angular velocity and r₀ is the radius of central near-constant-density plateau defined as $r_0 = \sqrt{A} c_{\rm s}^2 /(\pi G\Sigma _0)$, where c_s is the sound speed. These initial profiles are characterized by the important dimensionless free parameter (Ω₀r₀/c_s)² and have the property that the ratio of rotational to gravitational energy β ≈ 0.91(Ω₀r₀/c_s)² (see Basu 1997). We note that the above form of the column density is very similar to the integrated column density of a Bonnor–Ebert sphere with a positive density enhancement A (Dapp & Basu 2009). Throughout the paper, we use A = 1.2. In addition, every core is truncated so that the ratio of the core radius to the radius of the central flat region r_core/r₀ is constant and equal to 6.0. The truncation mechanisms could be erosion of the core outer regions by UV radiation and/or tidal stripping.

The actual procedure for generating a specific core is as follows. First, we fix β between approximately 10⁻³ and 10⁻² based on observational results that find that β ranges between approximately 10⁻⁴ and 10⁻¹ (Goodman et al. 1993; Caselli et al. 2002).⁵ Then we fix the core radius r_core and find r₀ from the condition r_core/r₀ = 6. The central surface density Σ₀ is found from the relation $r_0=\sqrt{A}c_{\rm s}^2/(\pi G \Sigma _0)$ and the resulting core mass M_core is determined from Equation (8). Finally, the central angular velocity Ω₀ is found from the condition β = 0.9(Ω₀r₀/c_s)².

As an alternative to the Σ∝r⁻¹ and Ω∝r⁻¹ profiles, we also consider several models with radially constant Σ and Ω distributions, corresponding to self-gravitating, sheet-like cores with volume density depending only on distance from the midplane ρ(z) = ρ(0) sech²(z/h) (Boss & Hartmann 2001).

In total, we have considered 23 models with initial core masses (M_core) ranging from 0.105 M_☉ to 1.84 M_☉ and initial ratios of rotational to gravitational energy β = (0.275–1.26) × 10⁻². The parameters of these 23 models are summarized in Table 1. The left panel of Figure 1 plots the mass-weighted initial angular velocity versus initial core mass for each model, where we have chosen to plot the mass-weighted rather than central angular velocity for direct comparison to observations. Our models range from about 10⁻¹⁴ rad s⁻¹ to a few times 10⁻¹³ rad s⁻¹, consistent with the observed range found by Goodman et al. (1993). The right panel of Figure 1 plots the distribution of initial core masses for each model. The absence of models with M_core > 2.0 M_☉ is caused by numerical difficulties associated with modeling the collapse of massive cores. We plan to extend our parameter space in a future study. This distribution clearly does not follow the stellar initial mass function (IMF). However, all results in this paper are presented after weighting each individual model by the IMF in order to properly compare to observations (see Section 4 for details on this weighting). Our results are thus robust to the exact number of models at different masses as long as we adequately sample all relevant mass ranges.

Zoom In Zoom Out Reset image size

Table 1. Model Parameters

Model	Ω0	r0	Σ0	Core Outer Radius	Initial Core Mass
	(rad s−1)	(AU)	(g cm−2)	(pc)	(M☉)	β
1	1.6 × 10−13	685	0.18	0.02	0.305	8.75 × 10−3
2	1.95 × 10−13	685	0.18	0.02	0.305	1.26 × 10−2
3	3.0 × 10−14	3770	0.033	0.11	1.684	8.8 × 10−3
4a	2.5 × 10−14	⋅⋅⋅	0.026	0.04	0.612	8.75 × 10−3
5a	0.9 × 10−14	⋅⋅⋅	9.3 × 10−3	0.11	1.686	8.8 × 10−3
6a	1.2 × 10−14	⋅⋅⋅	0.013	0.08	1.22	8.8 × 10−3
7	2.0 × 10−13	445	0.28	0.013	0.194	5.6 × 10−3
8	2.3 × 10−14	3770	0.033	0.11	1.689	5.6 × 10−3
9	2.6 × 10−13	514	0.24	0.015	0.229	1.26 × 10−2
10	2.9 × 10−14	3085	0.04	0.09	1.378	5.6 × 10−3
11	3.2 × 10−14	4115	0.03	0.12	1.84	1.26 × 10−2
12	3.7 × 10−14	2400	0.05	0.07	1.0726	5.6 × 10−3
13	3.8 × 10−14	1645	0.075	0.048	0.7337	2.75 × 10−3
14	3.9 × 10−13	342	0.36	0.01	0.1515	1.26 × 10−2
15	3.9 × 10−14	3430	0.036	0.1	1.5312	1.26 × 10−2
16	3.25 × 10−13	274	0.45	0.008	0.1174	5.6 × 10−3
17	4.7 × 10−14	1885	0.066	0.055	0.8434	5.6 × 10−3
18	4.8 × 10−14	2745	0.045	0.08	1.2422	1.26 × 10−2
19	5.56 × 10−13	240	0.52	0.007	0.1052	1.26 × 10−2
20	1.1 × 10−13	1200	0.1	0.035	0.5349	1.26 × 10−2
21	6.0 × 10−14	1370	0.09	0.04	0.6084	5.6 × 10−3
22	2.8 × 10−14	2230	0.056	0.065	0.999	2.75 × 10−3
23	2.0 × 10−14	2915	0.043	0.085	1.304	2.75 × 10−3

Note. ^aConstant surface density profile model.

Download table as:  ASCIITypeset image

Figure 2 plots the time evolution of the accretion rate onto the protostar ($\dot{M}_s$) for each of the 23 models considered in this paper and listed in Table 1. The accretion is highly variable with episodic bursts of rapid accretion, with the exact amplitude and frequency of the variability differing from one model to the next due to the different initial conditions (in particular initial core mass; see Vorobyov & Basu 2010 for details). The origin of the accretion variability seen in Figure 2 is discussed in detail below in Section 2.1.3. One particular feature worth noting in Figure 2 is that the average accretion rates show a gradual decline with time. On timescales of ∼ 1.0 Myr the decline in $\dot{M}$ may reach three orders of magnitude. This behavior is likely caused by a gradual decline of the disk mass with time, which is particularly prominent starting from the late Class I phase when mass loading from the parent core diminishes (Vorobyov 2009a, 2011).

Zoom In Zoom Out Reset image size

2.1.3. Disk Gravitational Instability and Variable Accretion

As was discussed in Section 1, episodic accretion is a plausible solution to the luminosity problem and various physical processes have been invoked to produce episodic and repetitive bursts of mass accretion onto the star. Among them, disk gravitational instability and fragmentation in the embedded phase of star formation has been shown to trigger intense bursts of accretion luminosity as the forming fragments are torqued into the disk inner regions and onto the star (Vorobyov & Basu 2005b, 2006, 2010; Machida et al. 2011). Apart from these FU-Ori-like bursts, gravitationally unstable disks are generally characterized by highly variable accretion onto the star caused by the nonlinear interaction of dominant spiral modes (Vorobyov 2009b) or by the presence of a massive planet (Machida et al. 2011). In this section, we take model 12 as a prototype model and use this model to illustrate the role of disk gravitational instability and fragmentation in driving variable accretion onto the star.

First, we must discuss the method by which we calculate basic disk properties. The disk mass and radius are calculated at each time step by disentangling the disk and infalling core on the computational mesh. We do this in practice by adopting a surface density threshold of 0.5 g cm⁻² between the disk and core and also using the radial gas velocity profile. This method is described in detail by Vorobyov (2011), although we note that Vorobyov (2011) adopted a lower surface density threshold of 0.1 g cm⁻². The values of the disk radius and mass both depend on the adopted threshold. Figure 3 plots a histogram showing the fraction of total time the 23 models considered in this paper spend at various disk radii (weighted by the IMF; see Section 4). Our choice of a threshold of 0.5 g cm⁻² is motivated by the fact that the resulting distribution of disk radii peaks between 10² and 10³ AU; adopting the original 0.1 g cm⁻² threshold would shift the entire distribution up by about a factor of 2–3. Sizes of embedded disks are very poorly constrained by observations and are typically assumed to be on the order of 100 AU or less based on simple centrifugal radius arguments. In reality, however, angular momentum transport will cause disks to spread to sizes greater than the centrifugal radii. Furthermore, there are some limited observations supporting the existence of large protostellar disks (i.e., the 1 M_☉, 300 AU disk surrounding Serpens FIRS 1; Enoch et al. 2009a). Ultimately, the correct threshold to adopt in order to disentangle the disk and core in the simulations will remain uncertain until the masses and sizes of protostellar disks are better constrained.

Zoom In Zoom Out Reset image size

Figure 4 presents gas surface density maps for model 12 at several times after the formation of the central star. Only the inner 1200 AU are shown, while the whole computational region amounts to 28,000 AU. The evolutionary times are chosen so as to emphasize three distinct stages of the disk evolution, each depicted by a separate row of images. The upper row shows the initial stage of the disk formation and growth. The disk begins to form at t ≈ 0.015 Myr and several distinct fragments are present as early as at t = 0.1 Myr. The number of fragments varies with time due to the continuing process of disk fragmentation and migration of the fragments onto the star. At around t = 0.16 Myr the second stage of the disk evolution begins when a massive fragment (≈0.15 M_☉) forms in the disk. The fragment is prominent until t = 0.3 Myr when it is ejected from the disk via many-body interaction with other newly born fragments. This ejection event is very transient and is not captured in Figure 4 but is evident in Figure 5 showing the time evolution of the total mass in our computational domain. This is a universal mechanism that seems to occur in massive disks and is reported in other numerical studies of disk fragmentation (Stamatellos & Whitworth 2009; Bate 2009; Basu & Vorobyov 2011). After the ejection event, the last stage of disk evolution ensues, which is shown in the bottom row of Figure 4 and is characterized by a dramatic change in the disk appearance—the disk is no longer prone to fragmentation and shows only a weak spiral structure that diminishes with time. We note that these three stages are only indicative and reflect the overall tendency of a protostellar disk to go through the initially vigorously unstable phase toward a marginally stable configuration. All three stages are often observed in relatively massive disks. In low-mass disks, the first and second stages may be very short or even absent and the ejection event may not be present even in massive enough disks (Basu & Vorobyov 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The three stages of the disk evolution are reflected in the mass accretion rate onto the star, $\dot{M_{\rm s}} = - 2\pi R_{\rm sc} v_r \Sigma$, which we calculate as the mass passing through the sink cell per one time step of numerical integration (which in physical units is usually equal to 10–20 days and the total integration time can go beyond 1 Myr). The latter value is also corrected for the jet efficiency to account for the mass that passes through the sink cell but is later evacuated with the jets. We note that the adopted size of the sink cell R_sc = 6 AU is larger than the stellar radius, except for the very early stage when the forming star is represented by the first hydrostatic core (FHSC). Decreasing R_sc entails a significant increase in the calculation times because the physical size of the computational zones in the ϕ-direction decreases as one approaches the singularity at r = 0, causing a corresponding decrease in the time step. Simultaneously, the number of grid zones in the r-direction needs to be increased to sustain an adequate resolution at r ∼ 100 AU on the log-spaced grid. This effectively limits our choice of R_sc to a few AU at best considering the available computational resources. We acknowledge that our results may be sensitive to this limitation. For instance, additional accretion bursts may be triggered in the inner disk at r < 6 AU due to the thermal ionization instability (Bell & Lin 1994), MRI (Zhu et al. 2009a, 2009b), or disk–planet interaction (Nayakshin & Lodato 2011). We have run a few models with R_sc = 2 AU and found that the general behavior of the accretion rates (variability, burst magnitudes) remains qualitatively similar but further work is needed to extend the computational region still closer to the stellar surface.

The time evolution of the mass accretion rate is shown in Figure 6. $\dot{M}_s$ is negligible in the pre-stellar phase and rises to a maximum value of 1.3 × 10⁻⁵ M_☉ yr⁻¹ when the FHSC forms at t = 0 Myr (the time is counted from this moment). The subsequent short period of evolution is characterized by a gradually declining $\dot{M}_s$ when the material from the core lands directly onto the forming star. A sharp drop in the accretion rate at t = 0.015 Myr manifests the beginning of the disk formation phase when the infalling core material hits the centrifugal barrier near the sink cell and the accretion rates drops to zero. However, the process of mass loading from the core continues, the disk grows in mass and size, and gravitational instability is soon ignited in the disk at t = 0.04 Myr. From this moment and until t = 0.16 Myr, i.e., during the first stage of disk evolution, accretion exhibits variability with periods in the 400–1000 yr range and rates between a few × 10⁻⁷ M_☉ yr⁻¹ and a few × 10⁻⁵ M_☉ yr⁻¹ and also several strong bursts with $\dot{M}_s$ ≳ 2 × 10⁻⁵ M_☉ yr⁻¹. The latter are caused by fragments migrating onto the star due to the loss of angular momentum via gravitational interaction with spiral arms in the disk (Vorobyov & Basu 2005b, 2006, 2010). The typical period of accretion variability implies that its driving force lies in the disk intermediate regions at r = 30–60 AU (for an average stellar mass 0.3 M_☉). Many theoretical and numerical studies, including our own, indicate that this is the minimum distance from the star at which gravitational instability and fragmentation can occur (e.g., Stamatellos & Whitworth 2008; Clarke 2009; Vorobyov & Basu 2010; Meru & Bate 2010), suggesting a causal link between the disk gravitational instability and accretion variability in our model.

Zoom In Zoom Out Reset image size

During t = 0.16–0.3 Myr, i.e., during the second stage of disk evolution, the accretion pattern changes notably and starts to show a higher-amplitude and longer-period variability, with accretion rates lying between 10⁻¹⁰ M_☉ yr⁻¹ and a few × 10⁻⁴ M_☉ yr⁻¹. This change is caused by the formation of a massive fragment clearly visible in the middle row in Figure 4. The typical period of accretion variability is now ≈5 × 10³ yr, implying that the driving force is located at r ≈ 250 AU. This value is an average orbital distance of the fragment, suggesting a causal link between accretion variations and orbital motion of the fragment. A similar correlation between the orbital period of a massive object in the disk and the accretion rate variability was found by Machida et al. (2011). These results demonstrate that the disk gravitational instability and fragmentation are the driving forces for the accretion variability in our models.

During t = 0.3–0.5 Myr, i.e., during the third stage of disk evolution, accretion onto the star again undergoes a dramatic change. There are no more bursts or high-amplitude variations. $\dot{M}_s$ shows a gradual decline with time from a few × 10⁻⁷ M_☉ yr⁻¹ to a few × 10⁻⁸ M_☉ yr⁻¹ with only low-amplitude flickering. This change is caused by the ejection of the massive fragment from the disk and the associated loss of a significant fraction of the disk material. As a result, the disk settles into a state with only marginal gravitational instability characterized by a weak and diffuse spiral structure. At this stage, the accretion rate is mostly determined by viscous torques rather than by gravitational ones (Vorobyov & Basu 2009a).

We caution that accretion rates may be sensitive to our choice of the α-parameter, which determines the magnitude of viscous torques in our model. In the present study, we use a spatially and temporally uniform α = 5 × 10⁻³, a value chosen for consistency with the work of Vorobyov & Basu (2009b, 2010). These authors studied numerically the secular evolution of viscous and self-gravitating disks, with particular emphasis on accretion rates, and found that α < 10⁻² is needed to reproduce both the FU Ori accretion bursts and the range and magnitude of accretion rates found for brown dwarfs and T Tauri stars.

However, the α-parameter may vary in space and time. In particular, if the gas surface density exceeds some critical value (a few × 100 g cm⁻²; Armitage 2011), then the ionization fraction due to cosmic rays and X-rays becomes insufficient to sustain the MRI and the corresponding α-parameter is expected to decrease significantly. In this case, the outburst phenomenon, in particular, and the accretion variability, in general, will be more pronounced as demonstrated by Vorobyov & Basu (2009a, 2010). In our models, surface densities in excess of the critical value are found in the inner 10–20 AU for relatively massive disks (Vorobyov 2010). This means that our choice of a fixed α = 5 × 10⁻³ may cause the outburst magnitude to damp somewhat in the inner parts of the disk.

Figure 7, which plots the time evolution of the global Fourier modes, substantiates our analysis. We calculate the first five modes using the following equation:

$$\begin{equation} C_{\rm m} (t) = {1 \over M_{\rm disk}} \left| \int _0^{2 \pi }\!\! \int _{R_{\rm sc}}^{R_{\rm disc}} \Sigma (r,\phi,t) \, e^{im\phi } r \, dr\, d\phi \right|, \end{equation} \tag{ 10 }$$

where M_disk is the disk mass and R_disk is the disk's physical outer radius.

Zoom In Zoom Out Reset image size

Three distinct patterns of modal behavior are seen in Figure 7, each corresponding to a particular stage of disk evolution. In the early evolution at t = 0.04–0.16 Myr, the nonlinear interaction between competing spiral modes in the gravitationally unstable disk gives rise to the mode variability resembling in magnitude and frequency that of the mass accretion rate.⁶ At t = 0.16–0.3 Myr, global modes are maximal (with lower order modes being marginally dominant) but show little variability, while the corresponding accretion rates are highly variable. The matter that is the dominant source of perturbation is now a massive compact fragment, and the accretion variability is driven by the dynamics of the fragment rather than by the nonlinear interaction of the spiral modes. After t = 0.3 Myr, the global modes show a fast decline by two orders of magnitude caused by a sharp decrease in the disk mass after the ejection event. Although the modes return to a highly variable pattern, their amplitudes are very low, ∼10⁻², meaning that the non-axisymmetric density perturbations in the disk at this stage are of the order of 1% and the mass accretion is now entirely controlled by viscous torques. The latter drive the disk toward an axisymmetric state (Vorobyov & Basu 2009a) and the corresponding accretion rates show only low-amplitude flickering, gradually declining as the disk is drained of its mass reservoir.

The time evolution of the disk, core, and stellar mass as well as the total star+disk+core mass in the computational domain were shown in Figure 5. It is worth noting that the stellar mass is always greater than that of the disk. This is a direct consequence of the regulating nature of the accretion burst phenomenon. The interplay between mass loading from the core and migration of the fragments onto the star helps to sustain the system near the fragmentation boundary, with the disk periodically going above and below this boundary until the infalling core is depleted of matter (Vorobyov & Basu 2006, 2010). The ejection event is well visible in the time evolution of the disk, core, and total mass. Both the disk and core exhibit a sharp change in mass at t ≈ 0.3 Myr, with the disk losing mass and the core gaining mass as the fragment leaves the disk and propagates radially outward through the core. At t ≈ 0.34 Myr, the fragment passes through the outer computational boundary at r = 14,000 AU, which is manifested by a sharp drop in the total mass in the computational domain. From the amplitude of this drop, we can estimate the upper limit on the mass of the fragment as 0.15 M_☉. This estimate includes the fragment itself as well as a mini-disk that surrounds it.

An ideal coupling of hydrodynamic simulations with radiative transfer models involves taking the detailed disk and core structure from the former and using them as inputs to the latter. Our hydro simulations, however, provide the radial structure and the disk vertical scale height but lack the detailed vertical density and temperature distributions (due to the adopted thin-disk approximation) needed to run radiative transfer models. One way to circumvent this difficulty would be to reconstruct the disk and core vertical structure solving for the combined equations of vertical hydrostatic equilibrium and vertical radiation transfer after updating the flow variables in the plane of the disk. This approach, which would render our hydro simulations essentially two plus one dimensional, is currently under development and will be presented in a future paper.

In this work, we have employed a simpler approach by taking analytic profiles for the core and disk structure and re-scaling them according to simulation's known parameters. More specifically, for each of the models listed in Table 1, the hydro simulations described above provide the initial core mass (M_core) and radius⁷ (r_core), along with the time evolution of the core mass, disk mass (M_disk), protostellar mass (M_star), disk outer radius (R_disk), accretion rate onto the protostar ($\dot{M}_s$), and accretion rate onto the disk ($\dot{M}_d$). These parameters are used to construct the core and disk structure according to analytic profiles described in more detail below. Incorporating the exact physical structure at each time step from the simulations will not fundamentally alter our results or conclusions since the mass and radius of the protostar and disk and accretion rates onto both objects are unchanged, but it will affect the detailed comparison to observations since the radiative transfer (and thus resulting core temperature profile, spectral energy distributions (SEDs), and observational signatures) depends on the two-dimensional physical structure.

As in Papers I and II, the core inner radius is held fixed at a value such that the initial optical depth at 100 μm is set equal to 10 (see Paper I, in particular Equation (4)) until the disk outer radius exceeds this value; once this occurs the core inner radius is set equal to the disk outer radius. Following both Papers I and II and Vorobyov & Basu (2010), the protostellar radius (R_star) is calculated from M_star following Palla & Stahler (1991; see Equation (11) of Vorobyov & Basu 2010), except at early times where it is modified to take into account the FHSC (Larson 1969) stage. Following Vorobyov & Basu, R_star is fixed at 5 AU for the duration of the FHSC phase and then smoothly joined to the Palla & Stahler values over a period of 1000 yr when the second core forms. We use a simple polytropic model for the forming star and assume that the transition from the FHSC to the second core (protostar) occurs when the central temperature exceeds 2 × 10³ K. This is an improved implementation of R_star as compared to that in Papers I and II but it yields similar lifetimes of the FHSC, (1–2) × 10⁴ yr.

As in Paper II, we adopt the core density profile given by the TSC84 solution for the collapse of a slowly rotating core to include the effects of rotational flattening. The TSC84 solution results in a core that is initially a spherically symmetric, singular isothermal sphere with a density distribution n∝r⁻²_core, identical to the Shu (1977) solution. As collapse proceeds, the solution takes on two forms: an outer solution that is similar to the non-rotating, spherically symmetric solution and an inner solution that exhibits flattening of the density profile. Since material falling into the central regions originates from larger radii and thus carries more angular momentum as time progresses, the radius where the inner solution must be used, and thus the radius at which flattening becomes significant, increases with time (r_flat∝Ω²₀t³; TSC84).

The TSC84 solution is parameterized by the initial angular velocity of the core and the time since the formation of the protostar. We truncate the solution at the given r_core for each model and then renormalize the density profile so that the core mass matches the current M_core given by the hydro simulations. As in Paper II, we use the velocity profiles given by the TSC84 solution to allow r_core to decrease once the infall radius⁸ exceeds the initial outer radius. By renormalizing the density profile to match the current M_core given by the hydro simulations our model is not completely physically self-consistent. However, this choice preserves both the exact evolution of core mass with time given by the hydro simulations and the qualitative feature of the TSC84 solution that the effects of rotational flattening on the core density profile increase with both time and initial core rotation.

Following Paper II, the assumed disk structure follows a power law in the radial coordinate and a Gaussian in the vertical coordinate, with the density profile given by

$$\begin{equation} \rho _{{\rm disk}}(s,z) = \rho _0 \left(\frac{s}{s_o} \right)^{-\alpha } \exp \left[ -\frac{1}{2} \left(\frac{z}{H_s}\right)^2 \right], \end{equation} \tag{ 11 }$$

where z is the distance above the midplane (z = rcos θ, with r and θ the usual radial and zenith angle spherical coordinates, respectively), s is the distance in the midplane from the origin ($s= \sqrt{r^2 - z^2}$), H_s is the disk scale height, and ρ₀ is the density at the reference midplane distance s₀. All parameters describing the dependence of H_s on s (and thus the flaring of the disk) are the same as in Paper II. We note that this adopted disk flaring is very similar to that found in the hydrodynamical simulations (see Figure 11 in Vorobyov & Basu 2010). On the other hand, the adopted gas surface density profile Σ∝r⁻¹ is somewhat shallower than the Σ∝r^−1.5 usually seen in hydrodynamical simulations of gravitationally unstable disks in the embedded phase of star formation. The density profile is truncated at R_disk and normalized so that the total mass matches M_disk given by the simulations. The inner radius of the disk is set equal to the dust destruction radius, calculated (assuming spherical, blackbody dust grains) as

$$\begin{equation} R_{{\rm disk}}^{{\rm in}} = \sqrt{\frac{L_*}{16 \pi \sigma T_{{\rm dust}}^4}}, \end{equation} \tag{ 12 }$$

where L_* is the protostellar luminosity (see below) and T_dust is the dust destruction temperature (assumed to be 1500 K; e.g., Cieza et al. 2005).

As in both Papers I and II, the total internal luminosity of the protostar and disk at each point in the collapse from core to star contains six components: (1) luminosity arising from accretion from the core directly onto the protostar, (2) luminosity arising from accretion from the core onto the disk, (3) luminosity arising from accretion from the disk onto the protostar, (4) disk "mixing luminosity" arising from luminosity released when newly accreted material mixes with existing disk material, (5) luminosity arising from the release of energy stored in differential rotation of the protostar, and (6) photosphere luminosity arising from gravitational contraction and deuterium burning. The first five components are calculated following Adams & Shu (1986); further details can be found in Papers I and II. In this work, however, the second and third components are calculated using direct input from the simulations, i.e., $\dot{M}_{\rm d}$ and $\dot{M}_{\rm s}$, respectively.

The sixth component, the photosphere luminosity arising from gravitational contraction and deuterium burning, follows the pre-main-sequence tracks of D'Antona & Mazzitelli (1994) with opacities from Alexander et al. (1989). In Papers I and II, an offset of 10⁵ years was assumed between the onset of collapse and the zero time of these tracks (see also Myers et al. 1998). In this work, we take an improved approach by noting that the zero time of the D'Antona & Mazzitelli tracks corresponds to the time when deuterium burning begins in the protostellar interior. Hence, we make use of a simple stellar polytropic model and assume that the zero time of the pre-main-sequence tracks corresponds to the moment when the temperature in the stellar interior reaches 5 × 10⁵ K. In order words, the time offset in our models is now the sum of two quantities: the time needed to form the first hydrostatic core from a pre-stellar cloud core, which may vary depending on the initial core mass, angular momentum, and density distributions, and the time needed to ignite deuterium burning in the forming central protostar, which is usually 2 × 10⁴ yr after the first hydrostatic core formation.

Neither the time evolution of the protostellar radius nor that of the photosphere luminosity is included in our models in a self-consistent manner. R_star is calculated from M_star following Palla & Stahler (1991), who assume a constant accretion rate of 10⁻⁵ M_☉ yr⁻¹. The photosphere luminosity is included by adopting the pre-main-sequence tracks of D'Antona & Mazzitelli (1994). These tracks do not include accretion and suffer from the uncertainty of how to define the zero time relative to the onset of collapse and formation of the protostar. In reality, the accretion luminosity depends on the protostellar radius (L_acc ∝ 1/R_star), and both the magnitude and zero time of the photosphere luminosity depend on the accretion history (Baraffe et al. 2009; Hosokawa et al. 2011; Hartmann et al. 2011). Thus, the accuracy of both components of the total model luminosity can be improved by including self-consistent calculations of the protostellar radius and photosphere luminosity. Such calculations can now be performed for any arbitrary accretion history (e.g., Baraffe et al. 2009) and will be explored in a future paper.

Finally, there is also external luminosity arising from heating of the core by the interstellar radiation field (ISRF). We adopt the same ISRF as Papers I and II: the Black (1994) ISRF, modified in the ultraviolet to reproduce the Draine (1978) ISRF, and then extincted by A_V = 0.5 of dust with properties given by Draine & Lee (1984) to simulate extinction by the surrounding lower density environment. We input the mean intensity of this ISRF (J_ext) into the radiative transfer code as an additional source of heating. The luminosity added to L_bol from this external heating, L_ext, is determined after each radiative transfer model is run by subtracting the total internal luminosity (the sum of the above six components) from the total model luminosity.

For each model, we use the two-dimensional, axisymmetric, Monte Carlo dust radiative transfer package RADMC (Dullemond & Turolla 2000; Dullemond & Dominik 2004) to calculate the two-dimensional dust temperature profile of the core throughout the duration of each model. Each model begins when the FHSC forms and terminates when 10% of the initial core mass remains, which we have taken to be the Class I/II boundary (Vorobyov 2009b; Vorobyov & Basu 2010). Some models with initially low core masses are extended until 1% of the initial core mass remains. The simulation output is resampled onto a time step grid appropriate for the radiative transfer models; a detailed description of this resampling and its effects on the results are given in the Appendix. For the dust properties, we adopt the dust opacities of Ossenkopf & Henning (1994) appropriate for thin ice mantles after 10⁵ yr of coagulation at a gas density of 10⁶ cm⁻³ (OH5 dust), which have been shown to be appropriate for cold, dense cores (e.g., Evans et al. 2001; Shirley et al. 2005). Isotropic scattering off dust grains is included in the model as described in Paper II. SEDs at each time step are then calculated at nine different inclinations ranging from i = 5°–85° in steps of 10°. An inclination of i = 0° corresponds to a pole-on (face-on) system, while an inclination of i = 90° corresponds to an edge-on system.

We use the model SEDs to calculate observational signatures of the models at each time step for each inclination. We calculate the bolometric luminosity (L_bol), the ratio of bolometric to submillimeter luminosity (L_bol/L_smm), and the bolometric temperature (T_bol). L_bol is calculated by integrating over the full SED,

$$\begin{equation} \mbox{$L_{{\rm bol}}$}= 4\pi d^2 \int _0^{\infty }\! S_{\nu }d\nu \,, \end{equation} \tag{ 13 }$$

while the submillimeter luminosity is calculated by integrating over the SED for λ ⩾ 350 μm,

$$\begin{equation} \mbox{$L_{{\rm smm}}$}= 4\pi d^2 \int _0^{\nu =c/350\, \mu {\rm m}} S_{\nu }d\nu \,. \end{equation} \tag{ 14 }$$

The bolometric temperature is defined to be the temperature of a blackbody with the same flux-weighted mean frequency as the source (Myers & Ladd 1993). Following Myers & Ladd, T_bol is calculated as

$$\begin{equation} \mbox{$T_{{\rm bol}}$}= 1.25 \times 10^{-11} \, \frac{\int _0^{\infty } \nu S_{\nu } d\nu }{\int _0^{\infty } S_{\nu } d\nu } \ \rm {K} . \end{equation} \tag{ 15 }$$

T_bol can be thought of as a protostellar equivalent of T_eff; T_bol starts at very low values (∼10 K) for cold, starless cores and eventually increases to T_eff once the core and disk have fully dissipated. The integrals defined in Equations (13)–(15) are calculated using the trapezoid rule to integrate the finitely sampled model SEDs.

We use the 1024 young stellar objects (YSOs) in the five large, nearby molecular clouds surveyed by the Spitzer Space Telescope Legacy Project "From Molecular Cores to Planet Forming Disks" (Evans et al. 2003) as our observational data set. Evans et al. (2009) compiled photometry and calculated L_bol and T_bol for all 1024 YSOs in the same manner as described above after correcting the photometry for foreground extinction (see Evans et al. 2009 for details). They concluded that 112 of the 1024 YSOs are embedded protostars based on association with a millimeter continuum emission source tracing a core. Since the models presented here are those of cores collapsing to form protostars, it is to these 112 objects that we compare. Core masses for all objects in Perseus, Ophiuchus, and Serpens are taken from Enoch et al. (2009b).

As shown above in Figures 8 and 9, the models considered here reproduce the full range of observations in L_bol–T_bol space and provide a reasonable match to the observed protostellar luminosity distribution. Thus we conclude that the accretion process predicted by the Vorobyov & Basu (2005b, 2006, 2010) simulations of collapsing cores resolves the luminosity problem, although we caution that these results must be revisited in the future once simulations that feature smaller sink cells closer to the stellar surface and fully capture both the physics in the inner disk and the actual protostellar accretion rates are possible (see Section 2.1.3 for further discussion). In Paper II, we argued that episodic accretion is both necessary and sufficient to resolve the luminosity problem, but this conclusion is subject to uncertainty since the actual prescription for episodic accretion included in that paper was quite simple and did not fully capture the accretion process in the simulations. Here we have coupled the exact evolution of the collapsing cores with radiative transfer calculations and demonstrated that the Vorobyov & Basu simulations resolve the luminosity problem, although we defer the question of whether or not episodic accretion itself is necessary to Section 5.4 below. We also note here that these models predict a smooth distribution in L_bol–T_bol space, whereas the models in Paper II featured white "excluded zones" where the models spent no time but sources were observed to exist. The fact that only three initial mass cores were considered in Paper II (0.3, 1, and 3 M_☉) was argued to artificially create these zones; the increased sampling of core masses in this paper confirms this argument.

Most of the remaining discrepancy between the model predictions and observations of protostellar luminosities actually appears in the form of a "reverse luminosity problem," evident in both Figures 8 and 9 as an overabundance of time spent at L_bol ≲ 0.1 L_☉ compared to observations. Indeed, a K-S test on the observed and model luminosity distributions that only compares the distributions above 0.1 L_☉ returns a value of 0.87, slightly higher than the value of 0.85 obtained for the full distributions. What causes the disagreement at low luminosities? Only 7% of the observed sources have L_bol ⩽0.1 L_☉ whereas the models spent 21% of the total time at such luminosities. Furthermore, the majority of this time (19% out of 21%) is spent at T_bol ⩽100 K, where there are no observed sources. At least some, and possibly all, of this discrepancy can be explained by observational completeness limits. As noted above, the c2d survey is generally complete to protostars with L_bol ≳ 0.05 L_☉, although the exact completeness limit varies depending on distance and the detailed shape of the SED of each source. Indeed, several extremely low luminosity sources undetected in the Spitzer c2d survey, with both internal and bolometric luminosities significantly below 0.1 L_☉, have recently been discovered, most through detections of outflows driven by cores with no associated Spitzer c2d sources (Chen et al. 2010; Enoch et al. 2010; Dunham et al. 2011; Pineda et al. 2011). With such low luminosities at least some of these sources may be first hydrostatic cores. Sensitive interferometer outflow surveys and very deep Herschel and James Webb Space Telescope infrared surveys directed toward cores currently classified as starless are needed to fully identify and characterize the population of such extremely low luminosity protostars and/or first cores before an accurate comparison to the models can be made for luminosities below 0.1 L_☉.

As noted above in Section 4, the models do not provide a good match to the observed T_bol distribution, with most of the discrepancy in a population of embedded objects at high T_bol (≳ 1000 K) predicted by the models but lacking in the observations. In the models, most of this time spent at T_bol ≳ 1000 K arises when R_disk is larger than a few hundred AU and the line of sight does not pass through the disk. As described in detail in Section 2.2, we adopt analytic profiles for the disk and core density profiles since the hydro simulations do not provide the full vertical density structure. In order to do this, we define the core inner radius to be equal to the disk outer radius so that there is no overlap between where the disk and core density profiles are defined. However, as a consequence of this method, large cavities devoid of material exist above the surface of the disk but within the core inner radius, and these cavities increase as the disk sizes increase. Lines of sight that pass through these cavities have reduced optical depths, allowing more short-wavelength emission to escape and thus increasing T_bol.

In reality, such large cavities are unlikely to exist; instead, the disk and core density profiles should smoothly join together. As we have mentioned in Section 2.2, an improved methodology that reconstructs the disk and core vertical structure and incorporates this exact structure into the radiative transfer calculations (rather than adopting simple analytic profiles) is currently under development and will be presented in a future paper. This method will result in a more accurate distribution of material above the disk surface and will likely remove much of the discrepancy between observed and model values of T_bol. Since the distribution of luminosities is set mainly by the accretion rates and protostellar masses, we argue that including a more accurate physical structure should not significantly alter our results on the resolution of the luminosity problem, although this will be explicitly tested in a future paper. We also note here the possibility that these models contain too much rotation and angular momentum, since this would cause both the rotational flattening of the cores and sizes of the disks to be overestimated, resulting in the models overpredicting T_bol compared to observations. As discussed in some detail in Section 2.1.2 above, we do not consider this to be very likely, but we acknowledge that it is a possibility since β, the ratio of rotational to gravitational energy, is not a directly observable quantity and observed ranges of β could be overestimated if infall and/or outflow motion is mistakenly attributed to rotation.

Even after including a more accurate physical structure, some discrepancy between observed and model values of T_bol may remain, particularly once the effects of outflow cavities are included. Figure 11, which plots an example SED for a model with high T_bol, shows that such models feature a double-peaked SED. This figure also demonstrates that, from optical wavelengths to either the 24 or 70 μm bands probed by Spitzer observations, such models resemble transition disk Class II sources (sometimes also referred to as cold disk sources; see Merín et al. 2010 and references therein). Given that our models predict that 40% of embedded sources are classified as Class II by T_bol, and that the observed fraction of Class 0+I to Class II sources is 0.19 when classifying via extinction-corrected values of T_bol (Evans et al. 2009), our models predict that 12% of Class II YSOs are actually embedded objects with SEDs like that shown in Figure 11. While this is consistent with the upper range of the observed fraction of Class II sources with transition disk SEDs (3%–12%; Merín et al. 2010; Furlan et al. 2011), the models clearly overpredict this fraction, as described above. Nevertheless, we caution that a small fraction of Class II sources with transition disk SEDs could in fact be embedded sources.

Zoom In Zoom Out Reset image size

If the above statement is true, why are they missing from the observed sample of embedded sources? As discussed in detail in Paper II, whether or not a population of embedded objects with high T_bol exists remains an open question. The Evans et al. (2009) sample is based primarily on the association of Spitzer sources featuring rising SEDs and red colors with millimeter continuum sources (Enoch et al. 2009b; Dunham et al. 2008) and is likely biased against such objects since they would often not meet these criteria and would instead be assumed to be chance alignments between millimeter continuum sources and background sources and/or later-stage YSOs. Furthermore, the extinction corrections applied by Evans et al. (2009) to the embedded sources are average extinctions over each individual cloud rather than true line-of-sight extinctions, and thus may underestimate the true extinction since current, active star formation (and thus the position of the youngest, embedded sources) is associated with the densest parts of molecular clouds (e.g., Heiderman et al. 2010; Lada et al. 2010). As shown in Paper II, such underestimates could, in the worst cases, artificially lower T_bol from several thousand K to several hundred K. Future work must revisit the observational samples and carefully evaluate whether or not a population of embedded sources with high enough T_bol to be classified as Class II or Class III exists.

For the 23 models listed in Table 1 and considered in this paper, the duration of the embedded phase (t_emb) ranges from 0.052 (model 19) to 0.951 Myr (model 11). For each model, Figure 12 plots both t_emb (left panel) and the average accretion rate over the embedded phase duration (right panel, calculated as final stellar mass divided by t_emb) versus the final stellar mass produced by the model. The duration of the embedded phase increases approximately linearly with the final stellar mass produced, and thus the average mass accretion rate (not to be confused with the instantaneous mass accretion rate, which is highly variable; see Figure 2) is approximately constant and does not depend on final stellar mass. Offner & McKee (2011) argued that models that tend toward constant accretion time rather than constant accretion rate are necessary to match the observed protostellar luminosity distribution, consistent with the conclusion by Myers (2010) that accretion rates that increase with mass can at least partially resolve the luminosity problem. We disagree that such models are necessary. Our models tend toward a constant average accretion rate of ∼1–3 × 10⁻⁶ M_☉ yr⁻¹ rather than a constant accretion time yet still provide an excellent match to the observed protostellar luminosity distribution.

Zoom In Zoom Out Reset image size

The average duration of the embedded phase, weighted by final stellar mass as described in Section 4, is 0.12 Myr. In contrast, Evans et al. (2009) derived an embedded phase lifetime of t_emb = 0.44 Myr based on the ratio of Class 0+I sources to Class II sources in the c2d sample and the assumption of a Class II lifetime of 2 Myr. Thus, these models predict a significantly shorter t_emb than suggested by recent observations. However, a number of caveats apply to this comparison: (1) we have taken the Class I/II boundary to be the point at which 10% of the initial core mass remains and terminate the models at this point. In reality the exact point at which to set this class boundary is uncertain and could easily shift the duration of the embedded phase by factors of ∼ 2 in either direction; (2) the observationally determined t_emb is pinned to a Class II lifetime of 2 Myr, but the uncertainty in this lifetime is about 1 Myr (50%; see discussion in Evans et al. 2009); and (3) the number of Class 0+I sources in the Evans et al. sample may be overestimated.

The third point above is emphasized by three recent studies. First, van Kempen et al. (2009) observed 22 Class I sources in Ophiuchus in HCO⁺ J = 4–3, a tracer of warm and dense gas, and showed that 11 (50%) are not detected and thus show no evidence of being surrounded by a dense core. Second, McClure et al. (2010) used their revised extinction law and classification method to show that greater than 50% of the Class I sources in Ophiuchus are highly extincted disk sources no longer embedded within cores. Third, Heiderman et al. (2010) observed 53 Class I sources in a variety of nearby star-forming regions in HCO⁺ J = 3–2, a tracer of dense gas, and showed that 31 (58%) are not detected and thus not associated with a dense core. All three studies likely represent upper limits to the true fraction of "fake" Class I sources. In the cases of van Kempen et al. (2009) and McClure et al. (2010) this is because there is substantial evidence that Ophiuchus is located behind an extinction screen that, if not properly accounted for, will redden source SEDs and artificially increase the number of Class I sources (see, e.g., Figure 12 of Evans et al. 2009). In the case of Heiderman et al. (2010) this is because their study specifically targeted suspicious Class I sources located in low extinction regions of clouds. We thus consider 50% as an upper limit to the fraction of Class I sources in the Evans et al. (2009) sample that are actually misclassified Class II sources. Shifting 50% of the Class I sources to Class II would decrease the Evans et al. (2009) observationally determined value of t_emb from 0.44 Myr to 0.2 Myr. Combined with the other caveats mentioned above, the duration of the embedded phase predicted by our models is marginally consistent with observations but should be revisited in the future as uncertainties in the observations are improved.

Finally, the Stage 0 ([M_star + M_disk]/[M_star + M_disk + M_core] ⩽0.5; André et al. 1993) duration for the 23 models considered in this paper ranges from 0.009 Myr (model 19) to 0.256 Myr (model 15), with an average (again weighted by final stellar mass as described in Section 5) of 0.027 Myr. Compared to the total embedded duration of 0.12 Myr, our models predict that the Stage 0 phase is only 23% of the total embedded duration. This is a natural consequence of the fact that these models feature average mass accretion rates that decrease with time (Vorobyov & Basu 2010; see Figure 2), thus the first 50% of the mass will accrete from the core faster than the second 50% of the mass. Our results are generally consistent with recent observationally determined estimates of the lifetime of Stage 0 relative to the total embedded phase. For example, Enoch et al. (2009b) found 39 Class 0 and 89 Class I sources in Perseus, Ophiuchus, and Serpens, giving a relative Stage 0 lifetime of 30% of the total embedded phase duration. Additionally, Maury et al. (2011) found that between 9 and 12 of the 57 protostars they identified in the Serpens South cluster (Gutermuth et al. 2008) were Class 0 sources, giving a relative Stage 0 lifetime of 16%–21% of the total embedded phase duration. While the Enoch et al. and Maury et al. results are difficult to compare quantitatively since they use different classification methods that trace the underlying physical stage to different degrees of accuracy (T_bol in the case of Enoch et al. and position in L_bol–M_core space in the case of Maury et al.), the general agreement between our models and these observational results is encouraging.

The 23 models considered in this paper feature between 0 (model 16) and 97 (model 15) accretion bursts, where the exact criteria for defining bursts are given in the Appendix. The percentage of total time spent in bursts ranges from 0% (model 16) to 11.8% (model 15), and the percentage of total mass accreted in bursts ranges from 0% (model 16) to 35.5% (model 4). On average (where the average is weighted by final stellar mass as described in Section 4), 1.3% of the total time is spent in bursts and 5.3% of the total mass is accreted in these bursts. These values represent the statistical average values of the fraction of total time spent and mass accreted in bursts assuming a standard Kroupa IMF (see Section 4 for details). We caution that the exact values depend on the criteria used for defining bursts and would increase if a lower accretion rate floor were used (see the Appendix).

The simple, parameterized models presented in Paper II spend between 1.5% and 2% of their total time in bursts, consistent with the 1.3% of total time featured by these models. However, in the Paper II models between 50% and 91% of the final stellar mass accretes in bursts, in stark contrast to the 5.3% of total mass that accretes in bursts in the models considered here. The explanation for this large difference lies in the detailed implementation of accretion bursts in Paper II. As described in Section 1, a burst was triggered each time the ratio of M_disk to M_star exceeded 0.2. At this point the accretion rate onto the star was increased from 0 M_☉ yr⁻¹ to 10⁻⁴ M_☉ yr⁻¹ until the disk was fully drained of mass, allowing the cycle to begin anew. However, in reality only the most extreme bursts reach accretion rates of 10⁻⁴M_☉ yr⁻¹, as evident from Figure 2. Furthermore, only a very small amount of the mass in the disk accretes in a given burst (about 10⁻² M_☉on average), whereas in Paper II the assumption that the disk was fully drained in each burst led to situations where up to 0.1–0.2 M_☉ were accreting in single bursts. The models considered here, based on actual hydrodynamic simulations rather than simple parameterizations, are significantly more realistic.

Comparison to the results presented in Papers I and II demonstrates that the accretion process predicted by the Vorobyov & Basu (2005b, 2006, 2010) simulations essentially resolves the luminosity problem inherent in models with constant mass accretion. As first noted by Kenyon et al. (1990), there are two types of non-steady mass accretion that could potentially resolve the luminosity problem: (1) accretion rates that start high and then decrease with time, and (2) generally low accretion rates punctuated by short, episodic bursts of high accretion. Figure 2 clearly illustrates that the Vorobyov & Basu simulations feature both declining accretion rates with time and short, episodic accretion bursts. Which of these is responsible for resolving the luminosity problem?

Given that, on average, only 5.3% of the total mass accretes in bursts, one might suspect that it is the declining accretion rates with time rather than the bursts that lower model luminosities and improve the match to observations. However, this 5.3% excludes all of the mass that accretes in lower-amplitude accretion rate increases that do not meet the criteria for a burst as defined in the Appendix. Thus, to properly evaluate whether the bursts are required in order to resolve the luminosity problem, we have time-averaged the accretion rates to filter out the effects of the bursts and variability and re-run all 23 models. We averaged all models over 20,000 yr durations unless the total model duration was less than 0.2 Myr, in which case we decreased the averaging window to either 10,000 or 5000 yr in order to preserve at least 10 time steps between the formation of the protostar and the end of the embedded phase.

Figure 13 plots histograms comparing the fraction of total time the time-averaged models spend at various L_bol and T_bol to the observed distributions, similar to Figure 9. While the time-averaged models provide a much better match to observations than models with constant mass accretion (see Papers I and II for such models), comparing Figures 9 and 13 shows that the time-averaged models feature a small shift to higher luminosities and do not do quite as good of a job resolving the luminosity problem. This is confirmed by a K-S test on the observed and time-averaged model luminosity distributions, which returns a value of 0.59, lower than the value of 0.85 returned for the original models.

Zoom In Zoom Out Reset image size

These results indicate that the declining accretion rates are not solely responsible for resolving the luminosity problem, a finding that is only reinforced by the fact that, even with the large windows over which we have time-averaged, the variability and bursts are not fully filtered out¹⁰ (see Figure 14 for an example). We thus conclude that it is a combination of both the accretion rates that decline with time and the variability and episodic bursts that resolve the luminosity problem. We consider this to be a plausible result given that the Vorobyov & Basu simulations self-consistently predict both, and we argue that the fact that the bursts are required is in general agreement with the other indirect evidence for accretion variability and bursts described in Section 1.

Zoom In Zoom Out Reset image size