PROTOSTAR MASS FUNCTIONS IN YOUNG CLUSTERS

In this paper "accretion" refers to protostar mass gain from its environment. For most of the accretion history, this mass gain is dominated by accretion through the circumstellar disk. It is assumed that the rate of mass gain from the disk to the protostar is equal to the rate of mass gain from the core-clump environment to the disk, as expected in steady state, and for simplicity both processes are here called accretion.

The accreting protostar of mass m is assumed to gain mass from its surrounding core and clump at a rate $\dot m$, according to

$$\begin{equation} \dot m = \dot m_0 [1 + (m/m_0)^2 ]^{1/2} \end{equation} \tag{ 1 }$$

where the initial mass accretion rate $\dot m_0$ and the mass scale m₀ are parameters independent of mass, time, and position in the cluster. This independence of cluster position neglects the likely increase of accretion rate near the clump center of mass, due to global infall of clump gas (Smith et al. 2009; Wang et al. 2010).

Equation (1) is based on the assumption that all stars begin their accretion from centrally condensed, thermally supported "core" gas. The surrounding turbulent clump gas provides mass for accretion and exerts pressure on the core. Such a two-component form of the accretion rate is needed to account for both low-mass and massive star formation. In the version adopted here, denoted 2CA, the constant term represents low-mass star formation from an isothermal sphere (Shu 1977). The mass-dependent term represents massive star formation, as in competitive accretion models (Bonnell et al 2001). Equation (1) resembles expressions in McKee & Tan (2003), McKee & Offner (2010, hereafter MO), Offner & McKee (2011, hereafter OM), and Paper I. However, here the accretion rate depends on current protostar mass and not on final mass as in MO and OM.

In the high-mass limit of Equation (1), the accretion rate of turbulent clump gas increases as a power of the protostar mass. The value of this power is set by the stopping probability and by the requirement that the mass function approach a declining power law at high mass, as in the IMF. When the stopping probability is equal in each time interval, the high-mass accretion rate varies as $\dot m \sim m$. This accretion rate dependence on mass is weaker than the rate $\dot m \sim m^2$ for stationary accretion from a uniform medium (Bondi 1952).

In this model thermal core gas sets the initial accretion rate of a protostar but not its final mass. The final mass m_f of a protostar may originate from a subregion of the parent core, as when a member of a small group is ejected (Reipurth & Clarke 2001), or it may originate from a region larger than the core, as when a massive star draws from core gas and then from clump gas (Smith et al. 2009; Wang et al. 2010). The final mass m_f is assumed to come from an "original mass" M, which is independent of the core mass.

The mass flow from M to m_f is assumed to be inefficient because of the processes that limit accretion. These include protostar ejection (Reipurth & Clarke 2001), gravitational competition (Bonnell et al. 1997), gas dispersal by stellar feedback due to outflows (Velusamy & Langer 1998), radiative heating and ionization (Offner et al. 2009; Klassen et al. 2012), and exhaustion of the original gas supply. Studies of the CMF and the IMF indicate that the typical star mass is less than the typical core mass by a factor of ∼3, which applies on average but not for each individual star (Rathborne et al. 2009; Michel et al. 2011). It is assumed here that m_f = εM with ε = 1/3 (Alves et al. 2007; Matzner & McKee 2000).

The protostar mass accretion rate is similarly inefficient because of dispersal of some of the original mass before it accretes onto the disk or before it accretes from the disk onto the protostar. An analysis of the TMC-1 system indicates an accretion rate efficiency of 0.25 with respect to spherically symmetric SIS collapse (Terebey et al. 2006). It is assumed here that the efficiencies of accretion and of accretion rate are equal,

$$\begin{equation} \dot m = \varepsilon \dot M. \end{equation} \tag{ 2 }$$

Here $\dot M$ is the accretion rate from the original mass M onto the protostar–disk system in the absence of any dispersal.

The assumption of constant efficiency of mass and mass accretion rate is a highly simplified average over accretion processes that are more complex in space and time, according to high-resolution observations (Lee et al. 2012) and simulations (Girichidis et al. 2011). Also, this assumption neglects neglect temporal inefficiency due to episodic or nonsteady accretion (Dunham & Vorobyov 2012), which is discussed in Section 4.1.

The original mass M₀ of the protostar mass scale m₀ is assumed to equal the mass of a Bonnor–Ebert sphere (Bonnor 1956; Ebert 1955), also known as a CIS, discussed in Section 2. The protostar mass scale m₀ reflects the physics of a critically stable core, which sets M₀, and it reflects the physics of gas dispersal during infall and disk accretion, which sets ε.

The mass scale is chosen to be M_CIS rather than M_SIS, the mass of a truncated singular isothermal sphere (SIS) in the limit of infinite central density (Chandrasekhar 1939). The CIS is more appropriate than the SIS as a mass scale for thermal star formation. The CIS is critically stable and is therefore on the verge of star-forming collapse. In contrast, the SIS is infinitely unstable and is less likely to be achieved as an initial state (Whitworth et al. 1996).

The SIS is nonetheless useful for this paper because of a simple relation between a CIS and a truncated SIS having the same temperature and mean density,

$$\begin{equation} M_{\rm CIS} = \Gamma _{\rm IS} M_{\rm SIS}, \end{equation} \tag{ 3 }$$

where the constant factor Γ_IS is defined by

$$\begin{equation} \Gamma _{\rm IS} \equiv \left({\frac{{\nu _c }}{2}} \right)^{3/2}. \end{equation} \tag{ 4 }$$

Here ν_c is the negative log–log slope of the density profile of an isothermal sphere at its radius of critical stability (Chandrasekhar 1939). At this radius, ν_c = 2.434 (Chandrasekhar & Wares 1949), so that Equation (4) gives Γ_IS = 1.342. This factor is close to the IMF slope Γ_S = 1.35 (Salpeter 1955), which is discussed in Section 3.6.

The protostar mass is obtained as a function of its accretion duration a, or the time since the start of its accretion, from Equations (1) and (2) as

$$\begin{equation} m = m_0 \sinh \left({\frac{{\dot m_0 a}}{{m_0 }}} \right). \end{equation} \tag{ 5 }$$

Equation (5) indicates that the protostar mass increases linearly at early times, as expected from the collapse of a centrally condensed isothermal core. At late times, the mass increases exponentially as it accretes clump gas, reaching massive-star values if its accretion does not stop first.

Protostar accretion can stop because of the same processes that limit its efficiency. As noted above in Section 3.1, these processes include protostar ejection, gravitational competition, gas dispersal by stellar feedback due to outflows, radiative heating and ionization, and exhaustion of the original gas supply. To date, no detailed study of the combination of these processes is available. Here their combined effect is described in a simple statistical way, as in Fletcher & Stahler (1994) and Paper I.

In this model a cluster forms protostars at constant birthrate b from time t = 0 until the present time t. The probability density that a protostar accretes for duration a and stops accreting between a and a + da is

$$\begin{equation} p(a) = \frac{{\exp (- a/t_s)}}{{t_s [1 - \exp (- t/t_s)]}} \end{equation} \tag{ 6 }$$

where the durations a span the range 0 ⩽ a ⩽ t and t_s is the accretion stopping timescale (Paper I).

When t ≫ t_s, the protostar population approaches a steady state, where the rate of accretion stopping equals the birthrate and the number of protostars approaches bt_s (Fletcher & Stahler 1994). Then the probability density p(a) approaches the condition of being equally like to stop, and the accretion stopping timescale t_s approaches the mean duration of accretion (Basu & Jones 2004; Paper I). This steady state condition t ≫ t_s is satisfied in most embedded clusters, whose oldest stars have an age of ∼1 Myr (Gutermuth et al. 2009), which is much greater than the accretion stopping timescale, ∼0.1 Myr (Dunham & Vorobyov 2012).

In such a steady state cluster, the mass function of still-accreting protostars has the same shape as the final mass function of the cluster but with a smaller amplitude and a smaller maximum mass (Paper I). Assuming that the final mass function matches the IMF, this property allows the IMF to constrain model parameters, as in the next section.

The PMF dN/d log m in steady state is obtained by assuming that the duration of accretion is the most important factor in setting the protostar mass, as found in simulations of competitive accretion (Bate & Bonnell 2005; Bate 2012). Then the probability densities for protostar mass p(m) and accretion duration p(a) are related by p(m)dm = p(a)da, and Equations (5) and (6) give p(m) as

$$\begin{equation} p(m) = \frac{{\Gamma (\mu + \nu)^{ - \Gamma } }}{{m_0 \nu }} \end{equation} \tag{ 7 }$$

where μ ≡ m/m₀ is the protostar mass normalized by the mass scale, ν ≡ (1+μ²)^1/2 is the quadrature sum of 1 and μ, and the exponent is defined by $\Gamma \equiv m_0 /(\dot m_0 t_s) = M_0 /(\dot M_0 t_s)$.

The mass function for N protostars, dN/d log m = ln(10)Nmp(m), can be written from Equation (7) as

$$\begin{equation} \frac{{dN}}{{d\log m}} = \ln (10)N\Gamma \left({\frac{\mu }{\nu }} \right)(\mu + \nu)^{ - \Gamma }. \end{equation} \tag{ 8 }$$

When μ ≪ 1, Equation (8) shows that the mass function increases linearly with μ, because at low mass the probability density p(m) in Equation (7) approaches a constant value. In contrast, when μ ≫ 1, the mass function decreases as μ^−Γ. Thus, Γ equals the negative of the log–log slope of the mass function at high mass.

The power law decline of the mass function at high mass is a consequence of the equally likely to stop probability in Equation (6) and of the linear mass dependence of the accretion rate at high mass in Equation (1). If the mass-dependent clump accretion term was omitted from Equation (1), the resulting mass function would decline at high mass as an exponential and not as a power law.

The constancy of the slope Γ can be understood in terms of the constancy of its defining parameters. The slope $\Gamma \equiv (m_0 /\dot m_0)/t_s$ is the ratio of two timescales—the accretion timescale $m_0 /\dot m_0$ and the accretion stopping timescale t_s. This ratio is expected to be constant when both timescales are controlled by the dynamical time of the accreting medium (Basu & Jones 2004). Observational estimates indicate that these timescales are similar, as discussed in Section 3.6.

The 2CA mass function in Equation (8) is similar to recently published functions that fit observed estimates of the IMF (De Marchi et al. 2011; Parravano et al. 2011; Myers 2011; Paper I). It has the same form as a model of the CMF based on competition between gravitational contraction of a filament and being equally likely to stop because of gas dispersal (Myers 2013). On the other hand, the 2CA mass function is more physical than the foregoing expressions, since it is based on gas properties of a thermal core in a turbulent medium.

The 2CA model resembles IMF models of turbulent fragmentation (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008, 2009; Hopkins 2012) in its reliance on thermal and turbulent gas properties. However, the 2CA model describes accretion onto protostars and specifies when accretion stops. In contrast, the above turbulent fragmentation models predict the CMF and assume that the CMF is a sufficient initial condition to predict the IMF. The models also differ because the 2CA model is simpler and has fewer parameters.

The initial mass accretion rate is related to the high-mass slope of the mass function. From the above definition of Γ and the definition of Γ_IS in Equation (3), the initial mass accretion rate can be written

$$\begin{equation} \dot M_0 = \left({\frac{{\Gamma _{\rm IS} }}{{\Gamma _{\rm S} }}} \right)\left({\frac{{\Gamma _{\rm S} }}{\Gamma }} \right)\left({\frac{{M_{\rm SIS} }}{{t_s }}} \right). \end{equation} \tag{ 9 }$$

In Equation (9), the factor Γ_IS/Γ_S lies within 1% of unity according to Equation (4). Thus, the high-mass slope Γ closely matches the Salpeter value Γ_S if the initial mass accretion rate matches M_SIS/t_s and conversely.

The initial mass accretion rate $\dot M_0$ = M_SIS/t_s needed to obtain Γ ≈ Γ_S resembles the accretion rate of a collapsing SIS. This conclusion follows because M_SIS/t_s can be written as the accretion rate for pressure-free collapse of a SIS, M_SIS/t_f, times the ratio of timescales t_f/t_s. This ratio appears close to unity, since estimates of the typical duration of accretion t_s and of the free-fall time of star-forming gas t_f are each of the order of 0.1 Myr in nearby star-forming regions (Evans et al. 2009; Dunham & Vorobyov 2012; Foster et al. 2009; Lada et al. 2010). The typical departure of t_f/t_s from unity is uncertain, but it is estimated to be less than a factor of three.

Thus, the initial mass accretion rate is constrained to approximate that of SIS collapse. However, this constraint is not specific enough to discriminate among SIS collapse models with and without thermal pressure, with varying initial overdensity, or with varying initial velocity (Shu 1977; Fatuzzo et al. 2004). A detailed model of CIS initial state and SIS accretion rate was presented by Lee et al. (2004).

To make specific predictions, it is assumed that the timescales t_f and t_s are equal, i.e. the mean duration of accretion equals the free fall time of the mass scale gas. Their common value is denoted as the star formation timescale τ = t_f = t_s. Then the mass function slope Γ equals the ratio of timescales $\Gamma = (m_0 /\dot m_0)/\tau$. This slope matches the Salpeter slope when the initial mass accretion rate onto the protostar is that of inefficient pressure-free collapse of a SIS, $\dot m_0 = {{8\varepsilon \sigma ^3 }/{(\pi G)}}$, where σ is the thermal velocity dispersion and G is Newton's gravitational constant.

On the other hand, if the timescales t_f and t_s are not equal, the foregoing points indicate that a CIS mass scale and an initial mass accretion rate due to pressure-free SIS collapse give a mass function slope Γ = Γ_IS(t_f/t_s). The slope can thus be "top-heavy" or "bottom-heavy" compared with the IMF depending on the ratio of t_f, the free-fall time of mass scale gas, and t_s, the mean accretion duration.

The analytic form of the mass function model allows simple expressions for the mean, median, and modal masses. They are similar to each other, and they reflect the peak but not the width of the broad mass distribution.

The mean mass is obtained from the probability density p(m) in Equation (7) by integrating mp(m) over all masses. The median mass is obtained by setting the cumulative probability to 1/2, and the modal mass is obtained by setting the derivative of the mass function in Equation (8) to zero. The mean and median mass become independent of their mass limits when the limits are sufficiently far from the mass scale, i.e. when m_min ≪ m₀ and m_max ≫ m₀. Equations (10) and (11) give expressions for these limiting cases. The modal mass, in Equation (12), is independent of the limiting values. In each case the mass is given in its dimensionless form, μ = m/m₀:

$$\begin{equation} \bar \mu = \frac{{\Gamma}} {{\Gamma^{2} - 1}}, \end{equation} \tag{ 10 }$$

$$\begin{equation} \mu _{{\rm med}} = \sinh (\ln 2/\Gamma), \end{equation} \tag{ 11 }$$

$$\begin{equation} \mu _{\rm mod } = \left[ {\frac{{(1 + 4/\Gamma ^2)^{1/2} - 1}}{2}} \right]^{1/2}. \end{equation} \tag{ 12 }$$

For Γ = 1.342 as in Section 3.1, Equations (10)–(12) give $\bar \mu$ = 1.68, μ_med = 0.540, and μ_mod = 0.630. These values indicate that the mass scale lies between the median and mean mass, and it is greater than the mode of the mass function by the factor 1/0.630 = 1.59. Also, the modal mass is significantly less than the mean mass, m_mod/$\bar m$ = 0.375, since the asymmetric mass function is steep at low mass and shallow at high mass.

For mass limits m_min = 0.01 M_☉ and m_max = 100 M_☉, and for mass scale m₀ = 0.22 M_☉ as adopted below, the mean mass exceeds the limiting value in Equation (10) by 1%, and the median mass is less than the limiting value in Equation (11) by 9%.

The mass scale, initial accretion rate, and stopping probability in Equations (1)–(6) require three independent parameters, which can be specified as the thermal gas velocity dispersion σ; the mass efficiency ε of star-forming gas, and the star-forming timescale τ. This timescale denotes the mean accretion duration and the free fall time of the mass scale gas, which are assumed equal to each other. In combination with the model equations, these three parameters completely specify the mass function model.

Parameter values that match observational estimates and that yield a mass function matching standard representations of the IMF are σ = 0.19 km s⁻¹, ε = 1/3, and τ = 0.12 Myr. These are based on observations of nearby regions of star formation. There NH₃ line observations of starless cores in Perseus give a median thermal velocity dispersion 0.19 km s⁻¹ corresponding to temperature 11 K (Foster et al. 2009), the mass function of starless cores in the Pipe Nebula appears shifted from the IMF by a factor three (Alves et al. 2007), and modelling of the evolution of the spectral energy distributions of protostars indicates that the embedded phase has a duration of 0.12 Myr (Dunham & Vorobyov 2012).

These parameter values yield a mass scale of m₀ = 8εσ³τΓ/(πG) = 0.22 M_☉. Applying Equations (10)–(12) gives a median, modal, and mean mass of 0.12, 0.14, and 0.37 M_☉, respectively. The mean mass agrees closely with the mean mass estimated from observations, 0.36 M_☉ (Weidner & Kroupa 2006).

Variation of parameter values from the above choices affect the mass function mass scale but not its shape. In this model, the mass function shape is constant because its log–log slopes, unity at low mass and −Γ at high mass, are independent of parameter values. The mass scale is also unchanged if its parameters vary so that the product εσ³τ is unchanged. Neglecting variation in ε, this condition is equivalent to varying the mean gas density and temperature as $\bar n \sim T^3$, i.e., where star-forming gas is hotter, it must also be denser. Observational evidence for small variation in mass scale is discussed in Section 5.1.

The 2CA mass function based on the above equations and parameters agrees well with the IMF estimates of Kroupa (2002) and Chabrier (2005). Each of these segmented IMF estimates is represented in Figure 1 by a smooth function with no discontinuities (Myers 2010). The 2CA mass function lies between the two IMF estimates and has essentially the same high-mass slope and the same mass scale. The low-mass limit is a power law of slope unity, due to the combination of ELS and constant initial mass accretion rate. At low mass the Kroupa and Chabrier estimates diverge, while the 2CA model approximates the geometric mean of the two.

Zoom In Zoom Out Reset image size

The good agreement between the model PMF and the IMF curves in Figure 1 indicates that the model can reproduce the low-mass shape and the mass scale of the IMF, with parameter choices based on observed properties of star-forming gas. The high-mass slope matches that of the IMF by construction, since it was assumed in Section 3.6 to constrain the initial mass accretion rate. On the other hand, the similarity between the PMF and the IMF as in Paper I is a special case. In general, the PMF may vary significantly from the IMF (MO; OM).

The luminosity of an accreting protostar is dominated by its accretion luminosity for a broad range of protostar masses. This property allows the observable accretion luminosity to be used as a probe of the protostar mass, which is not generally observable by other means. The accretion luminosity can be written as

$$\begin{equation} L = \frac{{\gamma Gm\dot m}}{{R_H }} \end{equation} \tag{ 13 }$$

where the luminosity efficiency γ is the ratio of the accretion luminosity from the star–disk system to that due to spherical accretion onto a spherical surface of radius R_★.

Most models assume that accretion occurs through a circumstellar disk. For steady accretion of infalling material through an optically thick disk, γ = 1/2 (Hartmann 1998). It is less clear how much kinetic energy of accretion reaches the protostar surface to be radiated as accretion luminosity. Assuming disk accretion and that half the kinetic energy is radiated away before it reaches the protostar, γ = 3/4 (Offner et al. 2009). Assuming disk accretion and that none of the kinetic energy is radiated away, γ = 1 (Klassen et al. 2012). These estimates neglect the possibility of nonsteady or episodic accretion, which can further reduce γ (Dunham & Vorobyov 2012).

The radius of the accreting protostar also depends on how much energy the accreting mass brings to the protostellar surface. For spherical accretion, the maximum energy is deposited, and the radius increases with both stellar mass and mass accretion rate (Stahler et al. 1986; Hosokawa & Omukai 2009). In contrast, it is expected that accretion of the same mass through a circumstellar disk, or through magnetospheric accretion columns, allows more energy to be radiated away, bringing less energy to the protostar surface.

In the limit of "cold disk accretion," the protostar radius is significantly smaller than for spherical accretion, and it varies less with protostar mass and accretion rate (Palla & Stahler 1992; Hartmann et al. 1997; Hosokawa et al. 2010). These properties of spherical and cold disk accretion apply during the "convection" phase, until the protostar has acquired a mass of several M_☉, depending on the accretion rate. Then the "swelling" phase begins, as the surface layers of the protostar rapidly expand (Stahler 1988), in response to the outward transport of entropy (Hosokawa et al. 2010).

Here simplified descriptions of these limiting cases are used to bracket the likely protostar radius and its corresponding accretion luminosity. It is expected that the most accurate description lies between these limiting cases. In the following expressions, subscripts s and c refer respectively to spherical and cold-disk accretion.

In the spherically symmetric limit, the radius is approximated by

$$\begin{equation} R_{\bigstar s} = R_{0s} \left({\frac{m}{{m_1 }}} \right)^\alpha \left({\frac{{\dot m}}{{\dot m_1 }}} \right)^\beta \end{equation} \tag{ 14 }$$

where R_0s = 26 R_☉, m₁ = 1 M_☉, $\dot m_1$ = 10⁻³ M_☉ yr⁻¹, α = 0.27, and β = 0.41 as given by Stahler et al. (1986). This relation agrees within a factor ∼2 with the stellar models of Hosokawa & Omukai (2009) for masses greater than 0.1 M_☉. For the accretion rate in Equation (1), Equation (14) applies up to a maximum mass ∼5 M_☉, when the radius begins to swell.

In the cold-disk limit, the radius is assumed to have the constant value R_★c = R_0c = 2.5 R_☉. This radius value is typical for low-mass spherical models, and it lies below most higher-mass models of cold-disk accretion (Stahler 1988; Hartmann et al. 1997; Hosokawa et al. 2010). It is assumed to apply for masses of 0.05–5 M_☉. However, unlike Equation (14) it is not an analytic fit to cold-disk accretion models. It should be viewed as a constant-radius model rather than a true cold-disk accretion model.

The constant radius approximation reflects the relatively small radius variation with mass at fixed accretion rate in the models of Hosokawa et al. (2010), but it does not match the typical radius value at each accretion rate. This approximation should be considered a simple basis for comparison rather than an accurate representation of cold-disk accretion models.

For spherically symmetric accretion, the accretion luminosity can be written using Equations (13) and (14),

$$\begin{equation} L_s = L_{0s} \mu ^{1 - \alpha } \nu ^{1 - \beta } \end{equation} \tag{ 15 }$$

where the luminosity scale is

$$\begin{equation} L_{0s} \equiv \frac{{\gamma _s Gm_1^\alpha \dot m_1^\beta m_0^{2 - \alpha - \beta } }}{{R_{0s} (\Gamma _s \tau)^{1 - \beta } }}. \end{equation} \tag{ 16 }$$

In the cold-disk approximation, the luminosity is

$$\begin{equation} L_c = L_{0c} \mu \nu \end{equation} \tag{ 17 }$$

where the luminosity scale is

$$\begin{equation} L_{0c} \equiv \frac{{\gamma _c Gm_0^2 }}{{R_{0c} \Gamma _c \tau }}. \end{equation} \tag{ 18 }$$

Equations (15)–(18) show that at low mass compared to the mass scale, the luminosity increases with mass to the power of 0.71–1.0, while at high mass, the luminosity increases more rapidly with mass to the power of 1.3–2.0. For the same mass, the cold disk luminosity exceeds the spherical luminosity because the radius is generally greater in the spherically symmetric case.

To evaluate typical luminosities, the parameter values σ = 0.19 km s⁻¹, ε = 1/3, τ = 0.12 Myr, and Γ_c = Γ_s = 1.35 are adopted, which match the 2CA model to the IMF. The luminosity efficiencies are γ_c = 1/2 for cold disk accretion and γ_s = 3/4 for spherical accretion, following the above discussion. Then the luminosity scales are L_0s = 6.03 L_☉ and L_0c = 1.88 L_☉. For protostar masses 0.05, 1, and 5 M_☉, Equations (15) and (17) give the accretion luminosity as 0.44–2.1 L_☉, 30–60 L_☉, and 250–1500 L_☉. These accretion luminosities exceed the zero-age main sequence luminosity due to nuclear burning for all masses considered.

For each observed protostar luminosity function (PLF), the cold-disk and spherical accretion models are used to estimate the underlying PMF, from

$$\begin{equation} \frac{{dN}}{{d\log m}} = \frac{{dN}}{{d\log L}}\frac{{d\log L}}{{d\log m}}, \end{equation} \tag{ 19 }$$

provided the luminosity increases monotonically with mass. The term d log L/d log m is evaluated for each of the limiting luminosity–mass relations (15) and (17).

The relations between accretion luminosity and mass in Equations (13)–(19) are based on an accretion history that is one of several that can fit observed luminosity distributions (OM). These relations assume that the mass accretion is a smooth function of time, neglecting fluctuations due to episodic accretion (Dunham & Vorobyov 2012). Thus, the present inference of protostar mass from its accretion luminosity is specific to the 2CA model, and it should be considered more useful for a statistical sample than for individual protostars.

The 2CA model predicts typical luminosities that match those in nearby cluster-forming regions. When Equation (19) is used with the model parameters in Section 3.8, the resulting luminosity functions have modal values 0.78 L_☉ for cold disk accretion and 3.8 L_☉ for spherical accretion. These estimates are similar to modal luminosities in nearby star-forming regions, typically 1–3 L_☉ (Evans et al. 2009; Dunham & Vorobyov 2012; Kryukova et al. 2012).

The 2CA model matches typical observed luminosities because it assumes that accretion is inefficient. The 2CA model specifies inefficiency ε in the mass transfer, ε in the rate of mass transfer from core to protostar, and γ in the energy transfer from the disk to the protostar, as described in Section 3.1. Thus, 2CA model luminosities do not exceed typical observed luminosities by an order of magnitude, in contrast to the simple spherical collapse models used to describe the "luminosity problem" (Kenyon & Hartmann 1995, hereafter KH; MO; OM).

Similarly inefficient departures from simple accretion flow were detailed in a study of the protostar and outflow in the dense core TMC-1 (Terebey et al. 2006). The efficiency of mass transfer from the core to the protostar was estimated to be 0.6, the efficiency of the mass accretion rate from the core to the protostar was 0.3, and the efficiency of energy transfer from the disk to the protostar was 0.6, for a combined luminosity efficiency of 0.09 with respect to spherical isothermal collapse onto a spherical protostar. This paper adopts the corresponding factors of 0.3, 0.3, and 0.50–0.75 for a combined efficiency of γε² = 0.06–0.09.

The luminosity expression of KH (their Equation (2)) assumes perfectly efficient spherical flow onto a spherical protostar surface. Application of the above efficiency estimates to KH Equation (2) yields accretion luminosities of 0.83–2.1 L_☉, which is in good agreement with typical luminosities in the KH sample.

A better understanding of the physical basis of the inefficiencies assumed in this paper is still needed, in terms of the structure and duration of accretion flows and dispersing flows. Furthermore, there is significant evidence for temporal inefficiency of mass transfer to the protostar in the form of episodic or intermittent accretion, due to disk instabilities (KH; Hartmann et al. 1997; Evans et al. 2009; OM; Dunham & Vorobyov 2012). The present model assumes only steady accretion. Thus, a complete treatment of the luminosity problem requires a better understanding of both the spatial and temporal inefficiencies with respect to idealized models of steady spherical accretion.

Recent observations at mid-infrared and far-infrared wavelengths with the Spitzer and Herschel telescopes have greatly improved our knowledge of protostar luminosity distributions in nearby star-forming regions (Evans et al. 2009; Dunham et al. 2010; Kryukova et al. 2012). These distributions are referred to here as PLFs (OM). In several regions the PLF has a single peak near 1 L_☉ and a high-luminosity tail extending to 10²–10³ L_☉ (Evans et al. 2009; Dunham et al. 2010; Kryukova et al. 2012; Kryukova 2011).

Observed PLFs have been compared to theoretical predictions based on various accretion models of star formation and on the requirement that the model reproduce the IMF (MO; OM; Paper I). Here PLFs are fit in order to compare the underlying mass functions with each other and with the IMF. It is assumed as before that the underlying mass distribution has the form of Equation (8) based on the 2CA accretion model and on random accretion stopping. However, the mass scale and high-mass slope are now allowed to depart from those that match the IMF. The corresponding model PLF is obtained from Equation (19) for each version of the protostar radius. For each cluster, each of these two model PLFs is fit to the observed luminosity distribution. Then each fit PLF is used to infer the underlying mass function.

Observed PLFs and model fits are presented for all currently available regions that have observed luminosity estimates for at least ∼100 protostars: Cygnus (Figure 2; Kryukova et al. 2012), Orion A (Figure 3), Mon R2 (Figure 4), and Cep OB3 (Figure 5; Kryukova et al. 2012). A PLF comparison with the "c2d" ensemble drawn from several nearby clouds was given in Dunham et al. (2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PLF fit requires three parameters—the number of protostars, the luminosity scale, and the high-mass slope. The luminosity scale and high-mass slope were chosen so that the model PLF passes through all or nearly all of the histogram bins. This goal was achieved on the high side of the luminosity peak, and one or two bins are missed on the low side, where the data are less complete (Kryukova et al. 2012). The best-fit PLFs for cold-disk and spherical accretion agree well with each other for each cluster. The PLF shape for spherical accretion is slightly more asymmetrical than for cold disk accretion.

Over all four clusters, the luminosity scales L_0c and L_0s span the narrow range of 0.6–1.8 L_☉. The underlying mass functions have negative log–log slope Γ_c = 1.0–1.8 for cold disk accretion, bracketing the Salpeter slope of the IMF, and significantly shallower slopes of Γ_s = 0.6–0.9 for spherical accretion. Table 2 summarizes these fit parameters.

Table 2. Properties of Cluster PLFs and PMFs

Cluster	NPS	Γc	L0c	γc	τc	m0c	$\dot m_{0c}$	Tc	Γs	L0s	γs	τs	m0s	$\dot m_{0s}$	Ts
			(L☉)		(Myr)	(M☉)	(M☉ Myr−1)	(K)		(L☉)		(Myr)	(M☉)	(M☉ Myr−1)	(K)
Cyg	1695	1.05	1.5	0.5	0.1	0.14	1.3	10	0.60	0.78	0.75	0.1	0.030	0.49	3.6
Ori	229	1.20	1.6	0.5	0.1	0.18	1.5	12	0.65	0.90	0.75	0.1	0.034	0.52	4.0
Mon	96	1.80	2.1	0.5	0.1	0.25	1.4	15	0.85	0.85	0.75	0.1	0.037	0.44	4.2
Cep	100	1.00	1.2	0.5	0.1	0.14	1.4	10	0.55	0.55	0.75	0.1	0.022	0.44	3.0

Note. In each cluster, N_PS is the number of protostars; Γ is the high mass exponent in Equation (8) which gives best model fit to the observed PLF; L₀ is the best-fit luminosity scale; γ is the assumed luminosity efficiency with respect to spherical accretion; τ is the assumed star formation timescale; m₀ is the derived mass function mass scale; $\dot m_0$ is the derived initial mass accretion rate; and T is the derived kinetic temperature of mass scale gas. Quantities associated with the protostar radius models for cold disk accretion or spherical accretion are subscripted with c or s, respectively.

Download table as:  ASCIITypeset image

The fits of model PLFs to observed PLFs in Section 4.2 suggest that each underlying PMF can be inferred by fitting the model luminosity function to the observed PLF. The PMF is useful because it reflects the birth history of cluster stars and can test different theories of star formation (MO; OM). For constant birthrate, 2CA accretion, and random stopping, it has the same shape as the final stellar mass function (Paper I).

However, it is not clear whether the typical embedded cluster satisfies the assumptions of constant birthrate, a steady state between birthrate and accretion stopping, and other conditions of the 2CA model with random stopping. In this section, mass functions are inferred from the PLFs fit in Section 4.2 to compare model mass functions with each other and with the IMF.

In estimating a PMF from an observed PLF, it is not required here that the PMF match the IMF. The accretion law and stopping probability each have the same form as in Section 3, so the PMF again tends to a power law at high mass. However, the temperature of the mass scale CIS is allowed to depart from the initially assumed 11 K, and the initial mass accretion rate is allowed to depart from the rate of SIS collapse. Thus, the mass scales m_0s and m_0c and the power law exponents Γ_s and Γ_c are not constrained to the standard IMF values.

To obtain the PMF, the exponents Γ_s and Γ_c are determined directly from the fits to the PLF. For a PLF with negative log–log slope Γ_L at high luminosity, Γ_s = (2-a-β)Γ_L and Γ_c = 2Γ_L, based on Equations (8) and (19). To obtain the PMF mass scales, the PLF fit values of L_0s and Γ_s for spherical accretion are combined with an assumed value of γ_s/τ^1−β in Equation (16) to obtain the mass scale m_0s. The fit values of L_0c and Γ_c for cold disk accretion are combined with an assumed value of γ_c/τ in Equation (18) to obtain the mass scale m_0c. The two limiting PMFs are then obtained for each cluster from Equation (8).

In each case, the star-forming timescale τ is assumed to be 0.1 Myr as in Section 3.8, and the luminosity efficiencies are γ_c = 0.50 and γ_s = 0.75 as discussed above. The mass range is 0.1–5 M_☉ for spherical accretion, limited by the range where Equation (14) is valid. The mass range is 0.05–5 M_☉ for cold disk accretion, where the lower minimum value is chosen to reveal the turnover mass of the mass function.

For cold-disk and spherical accretion, the mass functions are shown for each cluster in Figures 6–9. For each cluster and accretion model, Table 2 gives the assumed star-forming timescale and luminosity efficiency and the inferred mass scale, initial mass accretion rate, and initial mass scale gas temperature. The mass scale is obtained from Equations (15)–(18). The initial accretion rate is obtained from $\dot m_0 = m_0 /(\Gamma \tau)$, based on Equation (7). The kinetic temperature of the mass scale gas is obtained from T = (m_p/k)[πm₀G/(8Γετ)]^2/3, based on Equation (3), where m_p = 2.33m_H is the mean particle mass and k is Boltzmann's constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PMFs in Figures 6–9 and their properties in Table 2 show that the mass functions for spherical accretion have a mass scale mass lower by a factor of ∼3 and a high-mass slope shallower by a factor of ∼2 compared with the mass functions for cold disk accretion. The mass scale for spherical accretion, typically 0.04 M_☉, causes the modal mass to lie below 0.1 M_☉, the minimum mass where the radius approximation in Equation (14) is verified by the stellar structure calculations of Hosokawa & Omukai (2009). In contrast, the modal mass for cold-disk accretion lies above its minimum mass. Thus, in Figures 6–9 the mass functions for cold-disk accretion show both a high-mass power law and a low-mass turnover, while the mass functions for spherical accretion show only a high-mass power law. The low-mass turnover for spherical accretion is not shown because its location is too uncertain. The mass functions for spherical accretion would require relatively more stars with masses near the peak of the IMF in order to resemble the mass functions for cold-disk accretion.

In contrast to the mass functions for spherical mass accretion, the cold-disk mass functions are relatively similar in mass scale and shape to the IMF. However, evidence for early depletion of Li expected from cold-disk accretion models with high burst accretion rates was not found in a recent search (Sergison et al. 2013).

Table 2 indicates that initial mass accretion rates are 0.4–0.5 M_☉ Myr⁻¹ for spherical accretion and 1.3–1.5 M_☉ Myr⁻¹ for cold disk accretion. Each of these rates is within a factor of ∼2 of standard estimates for low-mass star formation. However, the initial temperatures of mass scale gas are 3–4 K for spherical accretion. These temperatures are significantly less than the 10–15 K values for cold disk accretion, which agree well with NH₃ line observations of dense cores (Rosolowsky et al. 2008).

Such low inferred temperatures could increase to match observed values if the product of the timescale τ and the mass efficiency ε decreased by a factor of ∼3. Similarly, the mass scale and initial temperature would increase toward observed values if the luminosity efficiency γ_s decreased by a similar factor of ∼3. However, it is difficult to justify any of these changes.

Thus, the spherical accretion model of the protostar radius appears to be inconsistent with the models of accretion and accretion stopping when they are used to fit observed PLFs and to compare with initial gas temperatures. In contrast, the constant protostar radius of the cold-disk accretion model is consistent with observed dense core temperatures but not with calculations of massive protostar evolution, which predict a substantial increase in the protostar radius with mass and mass accretion rate (Stahler et al. 1986; Yorke & Bodenheimer 2008; Hosokawa & Omukai 2009).

The simplest resolution of these discrepancies may be that the protostar radius lies between the models of spherical accretion and cold-disk accretion in a way that allows the typical mass scale to be 0.1–0.2 M_☉. This result might arise if low-mass protostars have a nearly constant radius and if the radii of more massive protostars increase less rapidly with mass and mass accretion rate than in Equation (14).

For a PMF inferred from a PLF fit, the slope depends more simply on model parameters than does the mass scale. The slopes inferred for spherical and cold-disk accretion are directly related to the slope of the PLF, as noted in Section 4.3. They have the ratio Γ_s/Γ_c = 1− (α + β)/2 = 0.66 for identical fits to a PLF that declines as a power law at high luminosity. Thus, the ratio of PMF slopes obtained from the same PLF differs from unity, mainly because the protostar radius depends differently on the mass and accretion rate between spherical accretion, where α + β = 0.68, and cold-disk accretion, where α + β = 0. This difference between PMF slopes does not depend on model efficiencies, timescales, or the choice of initial velocity dispersion.

If the likely radius dependence on mass and accretion rate lies between those of the cold-disk and spherical accretion models, Figures 6–9 and Table 2 suggest that the likely PMF is top-heavy with respect to the IMF. The eight best-fit values of Γ in Table 2 have a mean ± standard error of 1.0 ± 0.1. This typical slope lies below nearly all of the slopes in the compilation of Bastian et al. (2010). The only way the typical inferred PMF could be consistent with the IMF is if the protostar radii generally remain close to the constant value assumed here as a limiting case.

Such top-heavy PMFs are inconsistent with the steady state cluster model of Section 3 if the mass scale and the initial mass accretion rate are due to thermal equilibrium and collapse and if the mean accretion duration t_s is equal to t_f, the free-fall time of mass scale gas. The top-heavy PMFs cannot be explained by a cluster that is too young to have reached steady state, since then the cluster would have too few stars to be recognized as a large cluster, and it would also have too few massive protostars to match the IMF (Paper I).

A simple explanation for a top-heavy PMF in a steady state cluster with a thermal mass scale and thermal initial accretion rate as in Section 3 is that t_s exceeds t_f. If accretion typically lasts longer than the free-fall time of the mass scale gas, the slope of the PMF decreases, as noted in Section 3.6. Such a decrease would be accompanied by an increase in the modal mass, even if the mass scale remained constant. Equation (12) indicates that a decrease in slope from 1.35 to 0.67 would be accompanied by an increase in modal mass by a factor 2.0.

Alternately, a top-heavy PMF might arise because the cluster evolution differs from steady state. First-born stars are expected to experience less stellar feedback from neighbor stars than do later-born stars. Then first-born stars could have an average accretion duration that is greater than that of later-born stars and greater than the free-fall time of the typical core. If so, the slope of the mass function could start out shallow and then steepen as more low-mass stars form. These speculations are discussed further in Section 5.5.

The key features of a mass function model are its mass scale and the high-mass slope. In this model the mass scale is the mass of a CIS, which has the same dependence on mean density and kinetic temperature as the Jeans mass but a smaller coefficient. When Equation (3) is expressed in terms of core mean density n and kinetic temperature T, the mass scale mass is

$$\begin{equation} m_0 = \varepsilon \frac{{\Gamma _{IS} }}{{m_p^2 }}\left({\frac{6}{{\pi n}}} \right)^{1/2} \left({\frac{{kT}}{G}} \right)^{3/2} \end{equation} \tag{ 20 }$$

where m_p = 2.33m_H is the mean particle mass, k is Boltzmann's constant, and ε, Γ_IS, and G have the same values as in Section 3.8.

This mass function mass scale should not vary significantly if the model is to match the IMF under a wide variety of star-forming conditions. It is difficult to examine this question thoroughly, because the mean density and temperature inferred from observations can vary significantly because of differences in tracer and in resolution. Here the distribution of mass scales is examined for resolution-matched samples of starless cores in low-mass and massive star-forming regions.

The Perseus molecular cloud complex hosts the embedded low-mass clusters NGC 1333 and IC348 and a number of more isolated young stars (Bally et al. 2008). The maximum bolometric luminosity of its protostars is 9 L_☉ according to near- and mid-infrared observations (Kryukova et al. 2012). In contrast, a sample of 15 regions within 3.5 kpc of the Sun was selected to reveal conditions for formation of stars of intermediate and high mass in clusters (Sánchez-Monge et al. 2013). Their minimum and median bolometric luminosities are 340 and 3100 L_☉. If the minimum luminosity originates from a single protostar, its mass is at least 3–5 M_☉ according to the models of Section 3.9. Each region has significant clustering, but it is expected that the cluster luminosity is dominated by emission from its most massive members.

The regions were mapped in the NH₃ (1, 1) and (2, 2) molecular lines that trace density and kinetic temperature in star-forming gas (Ho & Townes 1983). These lines were observed with the Green Bank Telescope and the Very Large Array, giving essentially the same linear resolution, 0.04 pc, in 39 starless cores in Perseus (Foster et al. 2009) and in 23 starless cores in the massive star-forming regions (Sánchez-Monge et al. 2013).

Despite the differences in the protostar properties of these samples, the mass scales of their star-forming gas lie within a factor of two. The number distributions of m₀ have first-quartile, median, and third-quartile values of 0.15, 0.17, and 0.19 M_☉ for the low-mass sample and 0.19, 0.28, and 0.53 M_☉ for the high-mass sample. The mass scale distributions each span a factor of ∼2 and are marginally distinct. Their widths and separation are much narrower than the mass function itself. Thus, the mass function that allows for their spread is negligibly different from a mass function based on their median mass scale. Such a small range in mass scales may arise from the dependence of thermal dust–gas coupling on density (Elmegreen et al. 2008) or possibly from the reduction of fragmentation due to radiative heating in a cluster environment (Bate 2009; Krumholz 2011).

The consistency of mass scales between low-and high-mass star-forming regions within a factor of ∼2 would not hold if the effective temperature of turbulent motions in the NH₃ line widths were used instead of the kinetic temperature. Since the cores in massive star-forming regions have much broader line widths than in low-mass star-forming regions, the median mass scale masses based on line widths differ by a typical factor of ∼8.

In this model, the high-mass slope is due to both turbulent and thermal gas properties. The requirement that the mass function approach a power law at high mass is met by the linear dependence of the mass accretion rate on mass at high mass. This dependence refers to clump gas, which is more turbulent than core gas. However, the value of the power law slope, $\Gamma = (m_0 /\dot m_0)/\tau$, is due to the mass scale m₀ and initial mass accretion rate $\dot m_0$, which are set by the physics of the thermally supported core.

It is important to understand more physically how properties of turbulent clump gas lead to the accretion rate needed for the power law part of the IMF. The dependence of accretion rate on mass found here, $\dot m$ ∼ m, is significantly weaker than for point-mass accretion of nonturbulent polytropic gas, $\dot m$ ∼ m² (Bondi 1952). Such a reduction in accretion rate may apply to accretion of gas whose turbulence includes vorticity. An analysis of this process suggests an accretion rate weaker than in pure Bondi accretion and indicates that inclusion of magnetic fields may further weaken this dependence (Krumholz et al. 2006).

The association of thermal gas properties with the high-mass slope of the IMF might seem to conflict with the observations of massive protostars in gas whose line widths are supersonic (McKee & Tan 2003) and with observations of starless cores with turbulent motions in regions of massive star formation (Wang et al. 2008; Ragan et al. 2012). However, in this model, massive protostars are born in thermal cores, and they become massive by accreting both thermal core gas and turbulent clump gas. By the time they are identified as being massive because of their high luminosity, they are associated with gas that is turbulent because of its clump origin and also with core gas that has become more turbulent because of massive star heating, winds, and ionization.

It is necessary to specify when accretion stops in order to specify protostar final mass. However, the statistical model of being equally likely to stop used here is justified mainly by its simplicity. Its mean accretion duration is set by observational estimates, but the form of the distribution itself is neither supported nor refuted by observations and simulations. The likely mechanisms of accretion stopping, including ejection, ionization, winds and outflows, and exhaustion of initial gas, have been identified and discussed, but their relative importance over the history of a forming cluster remains unclear.

It is difficult to see how observations can clarify the distribution of accretion durations, but simulations may be helpful, provided they include the relevant physical processes and can be run long enough to achieve mass functions of final masses. Studies of how simulated protostar masses depend on accretion duration in the presence and absence of outflows, radiative heating, photoionization, and magnetic fields may be useful in developing more realistic distributions of accretion duration.

Furthermore, the same mechanisms that cause accretion to stop are invoked to justify the assumption of a constant efficiency factor ε between the original gas mass and protostar mass and between the corresponding mass accretion rates. Although this constant efficiency is a useful simplification, a more serious treatment of dispersal should account for the gradual stopping of accretion rather than assuming a constant efficiency followed by a sudden stop.

This model assumes that final protostar masses depend more on differences in accretion duration from one protostar to the next than on differences in initial core properties. Thus, the initial accretion rate and the mass scale in Equation (1) are each constants rather than distributions due to core-to-core variations in temperature and mass. Consequently, this model is independent of the core mass function, and the resulting mass function is due to accretion and accretion stopping rather than to core-to-core variations.

Nonetheless it is important to understand the origin of the core mass function and to understand why the CMF has a similar shape to the IMF. Several models of turbulent fragmentation reproduce the CMF (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hennebelle & Chabrier 2009; Hopkins 2012), but these models assume that the IMF masses are a fixed fraction of the corresponding CMF masses, without explaining how this fraction arises. Alternatively, it has been suggested that the IMF and the CMF have similar shapes because each is the result of the competition between accretion and dispersal. In one such model, cores grow as filaments contract radially under their self-gravity, and these cores stop growing as outer filament gas is dispersed by stellar feedback, photoionization, and evaporation (Myers 2013).

As discussed in Section 4, the PMFs inferred from the 2CA model and from observed cluster PLFs have a typical slope of Γ ≈ 1.0, which is significantly shallower than the slope of the IMF.

This result is of interest because it is the first of such conclusions based on protostars having a typical accretion age of the order of 0.1 Myr. In contrast, previous studies of young stellar object (YSO) mass functions are based on near-infrared photometry, optical and infrared spectroscopy, and pre-main sequence evolutionary tracks. Such studies of the L1688 cluster in Oph, the Orion Nebula Cluster, the Taurus complex, IC 348, σ Ori, λ Ori, and Cha I indicate typical YSO ages of a few Myr (Da Rio et al. 2012; Bastian et al. 2010; De Marchi et al. 2011).

The top-heavy mass functions in Figures 6–9 are reminiscent of mass functions in the Arches cluster and in NGC 3603. These massive starburst clusters with ages of 2–4 Myr show significant mass segregation and flattening of their inner mass function with respect to their outer mass function (Bastian et al. 2010). However, as noted above, the mass functions described here refer to much younger stars than those used to obtain mass functions in the Arches and NGC 3603.

A top-heavy PMF might arise from the preferentially early birth of massive stars. Such early birth is a common feature of simulations of cluster formation (Bate et al. 2003; Smith et al. 2009; Wang et al. 2010). In these simulations, the first protostars appear near the center of gravity of their parent clump. Their accretion is relatively long-lived because it is not yet limited by the dispersing effects of winds and outflows from later-born stars. These later-born stars contribute more to the mass function at low mass than at high mass, thereby increasing Γ toward the eventual IMF slope. A factor of ∼2 decrease in accretion timescale over the history of the cluster appears sufficient to account for the shallow slopes in Figures 6–9 and Table 2. However, a more detailed formulation of the accretion stopping probability is needed to test these ideas quantitatively.