FEEDBACK EFFECTS ON LOW-MASS STAR FORMATION

All simulations have the same initial conditions, also used in Offner et al. (2009b). The initial gas temperature is T_g = 10 K, the box length is L = 0.65 pc, and the average density is $\bar{\rho } = 4.46 \times 10^{-20} \,{\rm g \,cm}^{-3}$, corresponding to n_H = 1.91 × 10⁴ cm⁻³. The total box mass is 185 M_☉. For radiative simulations, the radiation temperature, T_r is initialized to 10 K. The radiation energy density is thus E = aT⁴_r = 7.56 × 10⁻¹¹ erg cm⁻³.

To obtain the turbulent initial conditions, we begin our simulations without self-gravity and apply velocity perturbations to an initially constant density field using the method described in Mac Low (1999). These perturbations correspond to a Gaussian random field with a flat power spectrum in the range 1 ⩽ k ⩽ 2. The application of these perturbations continues for three cloud crossing times and then stops. At this point the turbulence follows a Burgers power spectrum, P(k)∝k⁻², characteristic of supersonic hydrodynamic turbulence. The three-dimensional turbulent Mach number is $\mathcal {M} = 6.6$, which gives a three-dimensional rms velocity dispersion, σ_v = 1.2 km s⁻¹. With this Mach number the cloud is approximately virialized:

$$\begin{equation} \alpha _{{\rm vir}} = \frac{5 \sigma ^2}{G M / R} \simeq 1. \end{equation} \tag{ 1 }$$

This is slightly above the linewidth–size relation σ ≃ 0.7(R/1 pc)^1/2 km s⁻¹ (Solomon et al. 1987; Heyer & Brunt 2004) and is equivalent to σ = 1.2(R/1 pc)^1/2 km s⁻¹, which is well within the observed range (e.g., Falgarone et al. 2009)

After driving for three cloud crossing times, we then turn off driving, turn on gravity, and follow the subsequent gravitational collapse for approximately one global free-fall time:

$$\begin{equation} t_{{\rm ff}} = \sqrt{\frac{3 \pi }{32 G \bar{\rho }}} = 0.315 \,{\rm Myr}, \end{equation} \tag{ 2 }$$

where $\bar{\rho }$ is the mean density of the box. The simulations with radiation become prohibitively computationally expensive at late times and are stopped at t = 0.83t_ff with a total stellar mass of 30 M_☉ for simulation R. The barotropic simulations are continued to t = 1.05t_ff before they are stopped. At this time the total stellar mass in simulation B is 50 M_☉ compared to the total simulation mass of 185 M_☉. There is still gas bound to protostars totaling 11 M_☉ when the simulations end. Our stellar mass estimates may therefore be too low by 20%.

Given our temperature of 10 K, the Jeans length at $\bar{\rho }$ is

$$\begin{equation} \lambda _J = \left(\frac{\pi c_s^2}{G \bar{\rho }}\right)^{1/2} = 0.20 \,{\rm pc}, \end{equation} \tag{ 3 }$$

and the Jeans mass is

$$\begin{equation} M_J = \frac{4 \pi }{3} \left(\frac{\lambda _J}{2}\right)^3 \bar{\rho } = 2.7 \,M_{\odot }. \end{equation} \tag{ 4 }$$

The turbulent Jeans mass, at density $\rho = \mathcal {M}^2 \bar{\rho }$, is 0.4 M_☉.

The calculations have a 256³ base grid with four levels of refinement by factors of two, giving an effective resolution of 4096³. This resolution corresponds to Δx₄ = 32 AU at the finest refinement level.

We use the parallel adaptive mesh refinement (AMR) code ORION for our simulations. The numerical method is nearly identical to what we have used in previous papers (Krumholz et al. 2007a; Offner et al. 2009b; Cunningham et al. 2011; Krumholz et al. 2011; Myers et al. 2011). ORION solves the equations of compressible gas dynamics including self-gravity, radiative transfer, protostellar outflows, and radiating star particles, all on an adaptive grid. Every cell in the grid has four conserved quantities: mass density, ρ, vector momentum density, $\rho \mathbf {v}$, gas energy density, ρe, and radiation energy density, E. These conserved quantities can be used to calculate derived quantities such as velocity, $\mathbf {v}$, and pressure, P. In addition to the gas quantities, we evolve point-mass star particles, each with a position $\mathbf {x}_i$, mass M_i, momentum $\mathbf {p}_i$, angular momentum $\mathbf {j}_i$, and luminosity L_i. The subscript i refers to the star particle number. The particle method is explained in Krumholz et al. (2004, hereafter KKM04), with the addition of radiation (Krumholz et al. 2007a) and outflows (Cunningham et al. 2011). The full set of evolution equations for gas and particles is

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \mathbf {v}) + \sum _i [\dot{M}_{{\rm KKM04}} W(\mathbf {r}_i) - \dot{M}_{w,i} W_w(\mathbf {r}_i) \xi (\theta _i)]= 0, \end{equation} \tag{ 5 }$$

$$\begin{eqnarray} \frac{\partial \rho \mathbf {v}}{\partial t} + \nabla \cdot (\rho \mathbf {vv}) &=& - \nabla P - \rho \nabla \phi - \sum _i (\dot{\mathbf {p}}W(\mathbf {r}_i)\nonumber\\ && - \dot{M}_{w,i} v_{w,i} W_w(\mathbf {r}_i) \xi (\theta _i) \cdot \hat{\mathbf {r}}_i), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\frac{\partial (\rho e)}{\partial t} + \nabla \cdot [(\rho e + P)\mathbf {v}] = \rho \mathbf {v} \nabla \phi - \kappa _R \rho (4 \pi B - c E) \nonumber\\ & &\quad - \left(\frac{\rho }{\mu m_{{\rm H}}} \right)^2 \Lambda (T_g) -\sum _i \bigg[\dot{\varepsilon }_{{\rm KKM04}} W(\mathbf {r}_i)\nonumber\\ &&\quad - \dot{M}_{w,i} W_w(\mathbf {r}_i) \xi (\theta _i) \frac{k_B T_w \,{\rm K}}{\mu (\gamma - 1)}\bigg], \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t} E - \nabla \cdot \left(\frac{c \lambda }{\kappa _R \rho } \nabla E \right) &=& \kappa _P \rho \left(4 \pi B - c E \right) + \left(\frac{\rho }{\mu m_{{\rm H}}} \right)^2 \Lambda (T_g)\nonumber\\ && + \sum _i L_i W(\mathbf {r}_i), \end{eqnarray} \tag{ 8 }$$

$$\begin{equation} \nabla ^2 \phi = - 4 \pi G[\rho + \sum _i M_i \delta (\mathbf {r}_i)], \end{equation} \tag{ 9 }$$

$$\begin{equation} \dot{M}_i = \frac{1}{1+f_w} \dot{M}_{{\rm KKM04}}, \end{equation} \tag{ 10 }$$

$$\begin{equation} \dot{M}_{w,i} = f_w \dot{M}_i = \frac{f_w}{1+f_w} \dot{M}_{{\rm KKM04}}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \dot{\mathbf {p}}_i = \dot{\mathbf {p}}_{{\rm KKM04}}, \end{equation} \tag{ 12 }$$

$$\begin{equation} \mathbf {r}_i = \mathbf {x} - \mathbf {x}_i, \end{equation} \tag{ 13 }$$

$$\begin{equation} \theta _i = {\rm acos} (\hat{\mathbf {r}}_i \cdot \hat{\mathbf {j}}_i). \end{equation} \tag{ 14 }$$

The quantities entering these equations are defined below. Equations (5) and (6) are the fluid equations for mass and momentum, modified to include particles. Equations (7) and (8) are the energy equations for gas and radiation, respectively. The Poisson equation for the gravitational potential, ϕ, is given by Equation (9). The particle evolution is given by Equations (10), (11), and (12). We use periodic boundary conditions for all gas and particle quantities.

For the radiative runs, we adopt Marshak boundary conditions for the radiation field. This allows radiation to escape from the box as it would from a molecular cloud. The equation of state for the gas is given by

$$\begin{equation} P = \frac{\rho k_B T_g}{\mu m_{{\rm H}}} = (\gamma - 1) \rho \left(e - \frac{v^2}{2}\right), \end{equation} \tag{ 15 }$$

where μ = 2.33 is the mean molecular weight for molecular gas of solar composition and γ is the ratio of specific heats. Since most of the simulation domain is too cold to rotationally excite molecular hydrogen, we adopt γ = 5/3, representing a monatomic ideal gas. The term κ_Pρ(4πB − cE) in Equations (7) and (8) represents energy exchanged between the radiation field and the dust in our gas, with B = caT⁴_g/4π representing the Planck emission function integrated over all frequencies. The opacities κ_P and κ_R are Planck and Rosseland means given by the dust opacities from the iron normal, composite aggregates dust model of Semenov et al. (2003). We assume that the gas and the dust are thermally coupled. When the gas temperature exceeds the dust destruction temperature, the energy exchange term goes to zero and the gas and radiation unrealistically decouple. To address cooling from gas above the dust destruction temperature, we use the line cooling function Λ(T_g) from Cunningham et al. (2006). This removes energy from the gas and adds that energy to the radiation field (see Cunningham et al. 2011 for further details). The radiation flux limiter is given by $\lambda = ({1}/{R}) ({\rm coth} \mathit {R} - {1}/{\mathit {R}})$, where R = |∇E|/κ_RρE (Levermore & Pomraning 1981). It should be noted that we have excluded the radiation pressure and radiation enthalpy advection terms from Equations (6), (7), and (8) that appear in the analogous equation in Krumholz et al. (2007a). This approximation is justified in the formation of low-mass stars, as shown in Offner et al. (2009b).

When radiation is neglected, the energy exchange term from Equation (7) disappears, and we close the system of equations with a barotropic equation of state for the gas pressure:

$$\begin{equation} P = \rho c_{s0}^2 \left[1 + \left(\frac{\rho }{\rho _c}\right)^{\gamma -1}\right], \end{equation} \tag{ 16 }$$

where $c_{s0} = \sqrt{k_B T_0 / \mu m_{{\rm H}}}$ is the isothermal sound speed at temperature T₀ = 10 K, γ = 5/3, and ρ_c is the critical density. The critical density determines the transition from isothermal to adiabatic regimes and we adopt ρ_c = 2 × 10⁻¹³ g cm⁻³ to agree with the collapse solution from Masunaga et al. (1998) prior to H₂ dissociation. Simulations that use the barotropic equation of state achieve maximum densities of 5 × 10⁻¹⁵ g cm⁻³, with an effective temperature at ρ_max of 10.8 K. Most of the gas in any given simulations is effectively isothermal.

The particle quantities $\dot{M}_{{\rm KKM04}}$, $\dot{\mathbf {p}}_{{\rm KKM04}}$, and $\dot{\varepsilon }_{{\rm KKM04}}$ in Equations (5–7) represent the sink particle accretion rates of mass, momentum, and energy from the surrounding gas in the absence of winds as given by KKM04. The function W represents a window function over the four-cell accretion zone of the particle and W_w represents the outflow window function. Outflows are implemented as in Cunningham et al. (2011) and summarized here.

Our outflow model is specified by the dimensionless parameter, f_w, which sets the mass flux of outflow as a fraction of the accretion onto a star, and v_{w, i}, the wind launch speed. The mass fraction in our simulations is f_w = 0.3. The wind speed is set by the Keplerian speed at the surface of the star, $v_{k,i} = \sqrt{G M_i / r_{*,i}}$ where r_{*, i} is the protostellar radius, but is capped at 60 km s⁻¹ for computational speed. Specifically, v_{w, i} = min (v_{k, i}, 60 km s⁻¹). The velocity cap has a similar effect to the choice of Cunningham et al. (2011) to use v_{w, i} = v_{k, i}/3 for the most massive stars in the calculation. Wind gas is injected at temperature T_w = 10⁴ K.

The wind is injected over a window function W_w,

$$\begin{equation} W_w = \frac{1}{C_1}\left\lbrace \begin{array}{@{}ll@{}}\displaystyle r^{-2} & \mbox{if } 4 \Delta x < r \le 8 \Delta x \\ 0 & {\rm otherwise} \end{array} \right., \end{equation} \tag{ 17 }$$

which represents a shell just outside of the accretion region. The normalization constant C₁ is computed numerically to avoid numerical aliasing effects that occur from injecting a spherical wind into a Cartesian grid.

The exact angular distribution of the wind is described in the function ξ. The functional form is taken from Matzner & McKee (1999),

$$\begin{equation} \xi (\theta) = \left[ \ln \left(\frac{2}{\theta _0} \right)\big(\sin ^2 \theta + \theta _0^2\big)\right]^{-1}, \end{equation} \tag{ 18 }$$

where θ₀ is a flattening parameter that sets the opening angle of the wind. We use the fiducial value of θ₀ = 0.01. Equation (18) is averaged over the solid angle subtended by a grid cell at the outer radius of W_w. This averaging is particularly important near θ = 0. In addition, ξ is set to zero for θ ∼ π/2, so that winds are not injected directly into the plane of any equatorial disks.

We update the luminosity of each star, L_i, using the protostellar evolution model described in Offner et al. (2009b). In this model, 75% of the accretion energy is radiated away while 25% is nominally used to power winds. The energy of winds in our simulations is already determined by f_w and v_{w, i}, which is typically 15% of the accretion energy. The remaining 10% of the accretion energy is effectively lost. When winds are not present, we still only radiate 75% of the accretion energy for consistency across simulations.

The evolution equations can be described as the fluid and radiation equations from Offner et al. (2009b) combined with the particle equations and line cooling of Cunningham et al. (2011), but there is one important modification. In the KKM04 sink particle methodology, all particles with overlapping accretion zones are merged together. This gives an effective merger radius of eight cells, or 256 AU at a grid resolution of 32 AU. To limit this effect, we changed the merger radius to four cells, representing the point when a particle is in the accretion zone of another particle. This gives an effective merger radius of 128 AU at our resolution. Even with this improvement, our particle algorithm will unrealistically merge stars that pass within 128 AU. To address this, we have implemented a mass limit of m_merge = 0.05 M_☉, above which stars do not merge. This limit is chosen to correspond to the mass of second collapse in the formation of a star with final mass 1 M_☉ (Masunaga et al. 1998; Masunaga & Inutsuka 2000). Second collapse occurs when the protostar's core temperature becomes high enough to dissociate molecular hydrogen. Before second collapse, protostars are extended balls of gas with radii of a few AU and have a much higher collisional cross-section than main-sequence stars. After second collapse, stellar mergers should be extremely rare and we do not allow them. This approach is also used in Myers et al. (2011). The mass of second collapse depends on the accretion history of the protostar and is necessarily less than 0.05 M_☉ for brown dwarfs (Stamatellos & Whitworth 2009; Bate 2012). The effects of numerical merger suppression is explored in Section 4.9.

ORION utilizes a second-order Godunov scheme to solve the equations of compressible gas dynamics (Truelove et al. 1998; Klein 1999). These are Equations (5)–(7), excluding terms from stars and radiation. The Poisson equation (9) is solved using a multi-grid iteration scheme (Truelove et al. 1998; Klein 1999; Fisher 2002). The flux-limited diffusion radiation equation (8) and the radiation terms in Equation (7) are solved using the conservative update scheme from Krumholz et al. (2007b) modified to include the pseudo-transient continuation of Shestakov & Offner (2008).

We use the Truelove criterion (Truelove et al. 1997) to determine the addition of new AMR grids so that the gas density in the calculation always satisfies

$$\begin{equation} \rho < \frac{J^2 \pi c_s^2}{G(\Delta x_l)^2}, \end{equation} \tag{ 19 }$$

where Δx_l is the cell size on level l. We adopt a Jeans number of 0.125. In the simulations with radiative transfer, it is necessary to resolve the spatial gradients in the radiation field. Areas of high radiation gradients are near accreting stars, which tend to already be refined under the Truelove criterion. This is not always true for more evolved stars, which have higher luminosities and have accreted the dense gas that would trigger refinement. We find that the radiation gradients are adequately resolved by refining whenever |∇E|∇x_l/E > 0.25.

The evolution of the barotropic simulations with and without winds is depicted in Figure 1. Figures 2 and 3 depict the evolution of the radiative simulations without and with winds, respectively. In all simulations, for t ≲ 0.4t_ff, there are cloud-scale filaments that slowly contract, allowing three turbulent cores of width ∼20, 000 AU and density ∼10⁻¹⁹ g cm⁻³ to form. This length is half the Jeans length at the average cloud density. The first core to form is at the center of each panel in Figures 1–3, the second core is left of the central core, near the left edge, and the third core is near the right edge, at the end of a simulation-wide filament. At this point, the cores begin to fragment, while new cores form, eventually forming six fully developed cores. These cores each have a central stellar system. In simulation R, this central system is all that forms in each core. In all other simulations, the central system represents ∼75% of the stellar mass in the core and an additional group of low-mass stars forms, totaling ∼10 stars per core. These cores with multiple stars generally resemble observed high-stellar-density cores (Teixeira et al. 2007; Chen & Arce 2010). There are an additional 20 stars in cores that are still forming at the end of the simulation, giving a global total of 80 stars. Three of the cores coalesce by the end of the simulation to form a single group of 30 stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the three-dimensional rms velocity dispersion, σ_v, is shown in Figure 4. The global turbulence decays until star formation ramps up at t ∼ 0.5t_ff. There are two main mechanisms for star formation to increase σ_v. First, as stars accrete mass and deepen their gravitational potential, the surrounding gas can convert gravitational energy into kinetic energy as it falls into the stars. This is shown in the gradual increase in Mach number for t > 0.6t_ff in the B simulation. This effect is strong enough to return σ_v to near its original virialized value by itself. In rare cases, a many-body close encounter between stars will eject some gas at high velocities. There is not much momentum injected this way and the energy quickly dissipates, but it causes spikes in σ_v for the barotropic simulations, which have more small-scale fragmentation and therefore more many-body close encounters. The second mechanism occurs when protostellar winds are included. Some mass accreted onto stars is directly injected around the stars at high velocities. This causes the smooth increase in σ_v for simulations BW and RW as well as spikes from events with particularly high accretion rates that lead to bursts in wind momentum.

Zoom In Zoom Out Reset image size

The total momentum injected by winds for model BW is shown in Figure 5. For comparison, a characteristic value of the magnitude of the momentum associated with internal motions in the cloud, M_cloudσ_v is also plotted. For t > 0.8t_ff, the total momentum that has been injected by the winds is greater than the characteristic cloud momentum. At this point, the total amount of turbulent momentum that has been dissipated (including dissipation of wind momentum) is roughly the total amount of momentum that has been injected by the winds. By the end of the BW and RW simulations, the total wind momentum injected into the cloud is over twice the characteristic cloud momentum. The kinetic energy injected from the winds dissipates over time, suggesting a steady-state solution where the velocity dispersion of the cloud is constant with time as the winds replenish energy as quickly as it can dissipate.

Zoom In Zoom Out Reset image size

The total mass in stars and the total number of stars as functions of time are shown are Figures 6 and 7. The realization of the initial turbulence is slightly different between the radiative and barotropic simulations, so star formation begins at different times. This is the result of a slightly higher maximum density due to turbulence and it is unrelated to the choice of radiative or barotropic thermodynamics. The first stars form at t = 0.2t_ff for the radiative simulations and t = 0.35t_ff for the barotropic simulations. The turbulence overlaps enough between the radiative and barotropic cases, however, that at later times the total mass in stars at any given time is similar for the two cases when winds are not included. Winds reduce the mass in stars by about a factor of three in both the radiative and barotropic cases. The number of stars does not change between BW and B, implying winds do not cause or suppress fragmentation by themselves. The number of stars in RW is significantly greater than in R, however, because protostellar luminosity inhibits fragmentation and the winds reduce that luminosity. The three simulations other than R show a dramatic increase in the number of stars at 0.6t_ff < t < 0.8t_ff.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of median stellar mass is shown in Figure 8. This is a rough proxy for the characteristic mass of the protostellar mass function (PMF). Note that it will always be lower than the median mass of the IMF, because not all of the stars have finished accreting. Both BW and RW maintain a median around 0.05 M_☉ (similar to that of the PMFs in McKee & Offner 2010) throughout the simulation. The median does increase for simulation BW around t > t_ff as the formation rate of new stars decreases. The median of B fluctuates more, but is around 0.2 M_☉. Lastly, R maintains a median around 0.5 M_☉. This general behavior should be expected. The median mass is lowest when winds are included and highest when radiation is allowed to suppress fragmentation. The case with both winds and radiation ends up similar to BW because winds reduce protostellar luminosities.

Zoom In Zoom Out Reset image size

The global luminosity evolution for the radiative simulations is plotted in Figure 9. The winds reduce the total luminosity by a factor of up to 10 at any given time. This is expected since the global accretion luminosity is

$$\begin{equation} L=\frac{3}{4}\ \sum _i \frac{G M_{\star i} \dot{M}_{\star i} }{ R_{\star i}} \sim \frac{G M_{\star,\rm tot} \dot{M}_{\star,\rm tot} }{ \langle R_*\rangle }, \end{equation} \tag{ 20 }$$

given that the total mass in stars M_{⋆, tot} is reduced by a factor of three when winds are included, and the total accretion rate of stars $\dot{M}_{\star, \rm tot}$ is also reduced by three, the total luminosity is therefore reduced by a factor of nine assuming the characteristic stellar radius, 〈R_⋆〉, does not change. Main-sequence stars typically have a positive correlation between mass and radius, suggesting the factor of nine is an upper limit. However, at any given point in time, many stars in our simulations are in the degenerate regime where mass and radius are negatively correlated (Chabrier et al. 2009), which counteracts the positive correlation from the higher mass stars and keeps the total luminosity ratio near 9.

Zoom In Zoom Out Reset image size

The average stellar luminosity is shown in Figure 10. As was seen in the plot of total luminosity, the mean and median values of protostellar luminosity are much lower when winds are included. The disparity in average luminosity is even greater than the disparity in total luminosity because there are fewer stars when winds are excluded. The average luminosity in simulation R is heavily influenced by a single 6.6 M_☉ star that accounts for over half the total luminosity in the simulation. Unlike the low-mass stars, most of this luminosity is powered by nuclear fusion rather than accretion. Protostellar luminosities will be discussed further in Section 4.3.

Zoom In Zoom Out Reset image size

All simulations start at a background temperature of 10 K and are bathed in 10 K radiation. Stellar radiation and mechanical energy from protostellar outflows can raise this temperature. We have identified gas heated above 12 K as thermally affected by stellar feedback. Turbulent dissipation alone heats almost no gas above 12 K, leaving stellar feedback as the only explanation for this heating. The total mass of this gas is shown in Figure 11.

Zoom In Zoom Out Reset image size

The simulation with winds has significantly less heated gas than the simulation without winds. This is due to the reduced luminosity caused by the winds shown in Figure 9. In each simulation, the mass in heated gas roughly follows the mass in stars. When winds are included, the mass in stars drops by a factor of three, and the mass of heated gas also falls due to the reduced stellar luminosity. Some of the lost radiative heating is replaced by mechanical heating from outflows. This is evident in the ratio of heated gas mass to stellar mass. The mean value of this ratio is 0.9 without winds, but rises to 1.5 with them. Wind gas is injected at a temperature of 10⁴ K, but cools quickly. The mass in gas with temperatures above 1000 K is less than 2% of the mass in stars at any given time.

To further explore the heated gas, temperature–density phase plots with and without winds are shown in Figure 12. The phase plots with and without winds are notably different in two areas. First, the wind gas fills the high-temperature, low-density domain, while the same domain is empty without winds. Second, high-density gas with ρ > 10⁻¹⁶ g cm⁻³ has a higher temperature range without winds than with winds because the extra stellar luminosity heats that gas. When winds are included, that dense gas is less common in addition to being colder; there is more fragmentation, which turns dense gas into stars. In addition, some of the gas is also blown away by the winds themselves.

Zoom In Zoom Out Reset image size

The turbulent evolution shown in Figures 4 and 5 roughly agrees with previous simulations of molecular clouds with outflows, such as those in Nakamura & Li (2007) and Wang et al. (2010). Turbulent energy decays initially, only to be replaced by kinetic energy from winds and from gravity. While these sources can increase the total kinetic energy of gas, the new turbulence is fundamentally different from the isotropic, homogeneous hydrodynamic turbulence it replaces. This result is also seen in Nakamura & Li (2007), who find that the late time turbulent statistics do not match expected isotropic hydrodynamic results. One key difference is the energy from outflows is highly anisotropic. Outflow cavities are marked by long channels with high velocity shear between the fast outflow gas and the slow ambient gas. This shear is detectable as solenoidal energy. There is some compressive energy at the head of the outflow cavity, but most of the surface area is the side walls of the cavity and not the head. To measure the relative importance of solenoidal and compressive motions, we use the ratio 〈|∇ × v|²〉/〈(∇ · v)²〉, shown in Figure 13. In the case of isotropic turbulence, this ratio is 2, which is also the ratio of solenoidal to compressive energies (Elmegreen & Scalo 2004). Wind injection greatly increases the solenoidal velocities, steadily increasing the solenoidal to compressive ratio over the course of the simulation. Bursts in outflows injection around t = 0.5t_ff and t = 0.65t_ff are also visible in the solenoidal velocity. Anisotropic turbulence behaves differently than isotropic turbulence; in particular, it takes longer to decay since it decays on the crossing time calculated from the smallest velocity dispersion (Hansen et al. 2011). It is difficult to measure this increase in the decay time in our simulations, however, because the winds drive turbulence over a wide range of scales.

Zoom In Zoom Out Reset image size

The other major difference between the initial turbulence and wind-driven turbulence is seen in the rms velocity, σ_dense, of gas with $\rho > \bar{\rho }$. Even in isotropic, homogenous, hydrodynamic turbulence, there is a negative correlation between density and velocity, causing σ_dense < 7σ_v (Offner et al. 2009a). The winds themselves are collimated, very low density gas and have difficulty transmitting energy into high-density gas. This means that while σ_v is much greater with winds than without, σ_dense does not change much when winds are included. The evolution of σ_dense compared to σ_v is shown in Figure 14.

Zoom In Zoom Out Reset image size

Because the dense gas is relatively unaffected by the outflows, if our cloud had been centrally concentrated like that of Nakamura & Li (2007) or Wang et al. (2010), the dense part of the cloud would have most likely collapsed on itself even with the support of protostellar outflows. Magnetic fields may have an effect, as shown by Wang et al. (2010). The magnetic fields help transmit outflow energy to a much larger solid angle, so that even a centrally concentrated cloud can achieve a quasi-static balance between outflows and gravity.

Our simulations are marginally able to capture binaries, so most star particles should represent a single stellar object instead of a stellar system. The typical binary separation for main-sequence stars is 50 AU (Mathieu 1994), just slightly larger than our resolution of 32 AU. We cannot form stars within 128 AU of another star due to merging star particles, but stars that form beyond 128 AU can spiral in to 32 AU through interaction with gas. At any given time, about 2/3 of our star particles are in stellar multiples. The multiple properties are dynamic due to unstable high-order multiples and we cannot compare to the observed system properties. We can, however, compare to observed stellar properties. The mass functions of the four main simulations are shown in Figure 15 and compared to the stellar IMF in Chabrier (2005) as well as PMFs from McKee & Offner (2010). The mass functions in Figure 15 are shown at the latest time available for each simulation. The barotropic simulations could be evolved to later times due to the computational expense of flux-limited diffusion with many stellar sources. The mass functions in Figure 15 are not exactly comparable to an IMF because some of the stars are still accreting. The barotropic mass functions have evolved to a later time and should be similar to the IMF. The radiative simulations are still actively accreting and are closer to the PMFs in McKee & Offner (2010). Among the PMFs, the turbulent core and competitive accretion (not shown in Figure 15) models both roughly match RW. The isothermal sphere PMF is the best fit to mass functions of the simulations without winds, but this is merely a function of both mass functions being too top-heavy. The actual accretion is not described by an isothermal sphere and is discussed further in Section 4.3.

Zoom In Zoom Out Reset image size

As expected, the simulations without winds have mass functions skewed to higher masses than those with winds. The mass functions of both simulations without winds have too much mass at the high end compared to the IMF. The best fit to the IMF is from the BW simulation. It should be noted that the normalized mass functions for BW and RW look nearly identical when compared at the same time (explored in the next section). It is a good assumption that RW would also eventually match Chabrier at t ∼ t_ff.

Theoretical predictions of protostellar luminosities are often too high compared to observations of regions of low-mass star formation (Kenyon et al. 1990; Young & Evans 2005; Enoch et al. 2009), so it is important to compare our own simulations with observations. The mean and median luminosities of protostars observed in nearby clusters are 〈L_obs〉 = 5.3^+2.6_−1.9 L_☉ and L_{obs, med} = 1.5^+0.7_−0.4 L_☉, respectively (Enoch et al. 2009; Evans et al. 2009; Offner & McKee 2011). The typical mean and median luminosity in our simulation with winds are 〈L〉 = 6.9 L_☉ and L_med = 1.4 L_☉, in agreement with the observations. An additional useful quantity is the standard deviation of the luminosity, σ(L). This is sensitive to outliers in the luminosity distribution, but this can be mitigated by using the log of the luminosity. We find σ(log L) = 0.77 dex. This matches the observed value σ(log L_obs) = 0.7^+0.2_−0.1. For these comparisons, we have discarded stars with L < 0.05 L_☉, below the detection limit of the observations.

Theoretical frameworks that do match the observed protostellar luminosities are presented in McKee & Offner (2010) and Offner & McKee (2011). These models differ from the straightforward, isothermal sphere models (Fletcher & Stahler 1994a, 1994b) in three important ways. First, Offner & McKee (2011) assumed that 1/4 of the energy of the gas that accreted onto the star was removed non-radiatively. This effect is captured in our simulations. Second, they considered accretion rates that rise over time, such as predicted in core accretion and competitive accretion. These accretion rates all take the form,

$$\begin{equation} \dot{m} = \dot{m_1} \left(\frac{m}{m_f}\right)^j m_f^{j_f}, \end{equation} \tag{ 21 }$$

where m is the instantaneous mass of a protostar, m_f is the final mass of a protostar, and $\dot{m_1}$ is a constant throughout a cloud. More realistic accretion rates will rise at early times and slowly decline over time, in what McKee & Offner (2010) call "tapered accretion." There are multiple approaches to this; the one taken by McKee & Offner (2010) is to assume a linear decrease in the accretion rate with time, $\dot{m} \propto 1-t/t_f$, which implies that

$$\begin{equation} \dot{m} = \dot{m_1} \left(\frac{m}{m_f}\right)^j m_f^{j_f}\left[1 - \left(\frac{m}{m_f}\right)^{1-j}\right]^{1/2}. \end{equation} \tag{ 22 }$$

The values of j and j_f for several different accretion theories are shown in Table 2 as well as the values measured in our simulations with winds. Radiative and barotropic simulations produce the same j and j_f. The fits shown are for Equation (22). Our data do not agree with the functional form for untapered accretion (Equation (21)), but if this were used, j_f would remain the same while j ∼ −0.1. Our measured values of j = 0.3 ± 0.2 and j_f = 0.6 ± 0.2 marginally agree with the turbulent core accretion model, but do not match any of the other theories. Our j and j_f actually lie in between isothermal sphere and turbulent core accretion. This suggests that the most appropriate theories may be the two-component turbulent core model, which is approximately isothermal sphere accretion at early times and turbulent core accretion at later times. Other candidates include the TNT model (Myers & Fuller 1992), which includes both thermal and nonthermal support, and the core-clump accretion model from Myers (2011b) which is approximately isothermal sphere accretion at early times and reduced Bondi accretion at late times. Our low j does not agree with competitive accretion, which accelerates with mass.

While we agree with tapered accretion from McKee & Offner (2010), we do not rule out other tapered models, such as the exponential tapering in Schmeja & Klessen (2004). Our average accretion rates are lower than Schmeja & Klessen (2004) by an order of magnitude, partially due to our inclusion of protostellar outflows, but we see a similar rise and fall of accretion, which can be fit by many functional forms.

The last effect considered by McKee & Offner (2010) and Offner & McKee (2011) is episodic accretion from FU Ori type protostars (e.g., Hartmann & Kenyon 1996). If protostars spend most of their life in a low-accretion, low-luminosity phase and accrete nearly all their mass during short intense accretion bursts, the median stellar luminosity can be greatly reduced. This episodic accretion is thought to arise from disk instabilities (Vorobyov & Basu 2006). Offner & McKee (2011) estimated that 25% of the mass was accreted in this manner. We do not resolve disks and therefore do not see episodic accretion in our simulations. Given that our simulations match the observed luminosity function without episodic accretion, models that rely on substantial episodic accretion (Dunham et al. 2010; Stamatellos et al. 2011) may no longer be necessary.

Fragmentation in the simulations occurs in two distinct phases. First, the cloud as a whole forms cores of size ∼20, 000 AU (0.1 pc). This size scale is half the Jeans length at density $\bar{\rho }$. These large-scale cores each have a major filament within them that is above the critical line density for stability (Larson 1985; Inutsuka & Miyama 1992, 1997),

$$\begin{equation} \lambda _{{\rm crit}} = \frac{2 k T}{\mu m_{{\rm H}} G}. \end{equation} \tag{ 23 }$$

At our temperature, T = 10 K, the critical line density is λ_crit = 1.0 × 10¹⁶ g cm⁻¹. The line densities of the filaments in our cores range from 1.7 × 10¹⁶ to 4.0 × 10¹⁶ g cm⁻¹, similar to filament line densities seen in Serpens (André et al. 2010). The general morphology is similar to the hub-filament structure in Myers (2009) and Myers (2011a) with the hub at the center of each core.

The large-scale fragmentation is not affected by radiation. In contrast to high-mass star formation (Krumholz et al. 2007a; Cunningham et al. 2011; Krumholz et al. 2011), there is simply not enough protostellar luminosity to affect the large scales except possibly at late times in the R simulation. The winds do travel through the entire simulation domain and could theoretically affect the fragmentation, but this does not happen in practice since the winds do not couple well to the cores.

The second stage of fragmentation occurs as the filament in each core contracts under self-gravity. The filaments can then fragment into many stars around the central stellar system. Unlike the large-scale core fragmentation, this small-scale fragmentation can be significantly suppressed by radiation (Krumholz et al. 2007a; Offner et al. 2009b). This is best demonstrated in Figure 7. Simulation R stagnates at 30 stars, while B and BW finish with 80 stars. At the end of R, the total protostellar luminosity is 2500 L_☉ and this is currently heating 30 M_☉ of gas. Winds by themselves do not affect small-scale fragmentation. The total number of fragments, and therefore the total number of stars, are the same in the barotropic simulations with and without winds, as shown in Figure 7. It should be noted that our fragments occur in the core and not in accretion disks. This agrees with higher resolution core evolution simulations (Offner et al. 2010).

When winds and radiation are combined, the winds have a significant indirect effect. The winds lower the mass and accretion rate of the protostars, which lowers the luminosity as seen in Figure 9. This means that the fragmentation suppression seen in comparing R to B should be reduced when comparing RW to BW. To investigate this in more detail, we have shown the mass functions of RW and BW, both at t = 0.83t_ff, in Figure 16. For purposes of comparison, we have excluded stars with mass M < 0.05 M_☉, since these stars are usually short-lived and possibly subject to details of numerical sink particle formation or merging. The two mass functions are nearly the same. A two sample Kolmogorov–Smirnov test gives a 50% chance that both samples are drawn from the same underlying function. The primary difference between the two mass functions lies in the heaviest star in the simulation. The core that forms the heaviest star fragments early in BW, turning a star that is 2.8 M_☉ in RW into a 2.0 M_☉ star in BW with two extra M ∼ 0.4 M_☉ stars. The reduced fragmentation in RW for that star might be expected from radiative feedback, since that is the most luminous star in the simulation.

Zoom In Zoom Out Reset image size

To understand why radiation does not significantly alter star formation for regions of low-mass star formation such as we are considering, we expand on a line of reasoning developed by Bate (2009b). He introduced an effective Jeans length, which we write as λ_eff, defined by the condition

$$\begin{equation} \lambda _{\rm eff}\equiv \frac{GM}{c_s^2}, \end{equation} \tag{ 24 }$$

so that λ_eff is the radius at which the escape velocity equals the isothermal sound speed, c_s. Since the mass inside λ_eff is M = (4π/3)ρλ³_eff, one finds that Bates' effective Jeans length is smaller than the standard Jeans length, λ_J, which is given in Equation (3), by a factor π(4/3)^1/2 = 3.63. On the other hand, λ_eff is comparable to the radius of the fragments we find in our simulations, which have a diameter of λ_J/2, corresponding to a radius of λ_J/4. Bates' effective Jeans mass, which is the mass within a radius λ_eff, is smaller than the standard value in Equation (4), by a factor π³/√27 ≃ 6.0.

To estimate the dust temperature at a distance λ_eff from a star of luminosity L, we follow Bate (2009b) and ignore the frequency dependence of the dust absorption coefficient; λ_eff is then given by the condition

$$\begin{equation} L=4\pi \lambda _{\rm eff}^2\sigma _{\rm SB}T^4, \end{equation} \tag{ 25 }$$

where σ_SB is the Stefan–Boltzmann constant. Since c_s = (kT/μm_H)^1/2, we have

$$\begin{equation} T=\left(\frac{G\mu m_{\rm H}\rho L}{3k\sigma _{\rm SB}}\right)^{1/5}=5.3\left(\frac{\rho }{10^{-19}\mbox{ g cm$^{-3}$}} \cdot \frac{L}{L_\odot }\right)^{1/5} \mbox{K}. \end{equation} \tag{ 26 }$$

Bate (2009b) adopted a fiducial luminosity of 150 L_☉, which gives T ≃ 15 K for ρ ≃ 10⁻¹⁹ g cm⁻³. This is larger than the background temperature of 10 K that we have assumed; the Jeans mass is therefore increased and fragmentation is suppressed. However, the observed median luminosity in local star-forming regions is far smaller: Enoch et al. (2009) find a median luminosity of only 1.5 L_☉, which gives a temperature of only 5.8 K at λ_eff. Since this is less than the background temperature, we conclude that protostars typically do not suppress fragmentation at densities of order 10⁻¹⁹ g cm⁻³, which are characteristic of low-mass star-forming regions. At higher densities, the accretion luminosity can raise the temperature enough to begin to suppress fragmentation in clouds with T ≃ 10 K.

Our cloud fragments at about half the Jeans length (i.e., at ∼0.1 pc), but then continues to fragment below this point. Fragmentation at the Jeans length is commonly observed (Blitz & Williams 1997; Enoch et al. 2008). In instances where observers have the resolution and sensitivity to resolve fragmentation at scales below the Jeans length, however, even more fragmentation is found (Motte et al. 1998; Johnstone et al. 2000; Teixeira et al. 2007; Chen & Arce 2010; Bontemps et al. 2010). Fragmentation can be quantified as a function of scale, r. Given a set of clump locations in a cloud, one can calculate a set of clump pair separations. Let the differential number of pairs separated by distance r be dN_pair = H(r)dln r. The number of clump pair separations for randomly distributed clumps is H_ran(r). The two-point correlation function, w, can be calculated from these quantities (Johnstone et al. 2000),

$$\begin{equation} w(r) = \frac{H(r)}{H_{{\rm ran}}(r)} - 1. \end{equation} \tag{ 27 }$$

The two-point correlation function has been measured for the central parsec of ρ Ophiuchus by Johnstone et al. (2000). In this measurement, excess fragmentation (w > 0) is found below r ∼ 3 × 10⁴ AU, similar to the Jeans length of the cloud. There is a power-law fit, w(r)∝r^−0.75, in this regime. Larson (1995) also measured clustering of stars in Taurus and found a power-law fit with a break at 8000 AU, but attributed the break to binary stars. The separation between stars has had time to evolve since the initial fragmentation, so we narrow our focus to comparisons with Johnstone et al. (2000).

To compare our simulations to the observed w, we first created optically thin column density maps of our simulations and convolved them with a Gaussian with an FWHM of 1600 AU. The resolution was chosen to be similar to that from Johnstone et al. (2000). We used the Clumpfind algorithm from Williams et al. (1994) on the convolved column density map from simulation RW to obtain a list of clumps and their positions. To investigate the possibility of time evolution of w, this is performed at an early time in the simulation and then again at a late time, t ∼ 0.4t_ff and t ∼ 0.75t_ff, respectively. The results are shown in Figure 17. As expected, the correlation function drops off above r ∼ 4 × 10⁴ AU, about 2/3 the Jeans length at $\bar{\rho }$ for our simulation. More remarkably, the correlation function in our simulation matches that measured in ρ Ophiuchus quite well at all scales below this drop off. The early and late time simulations also generally agree with each other, suggesting there is not much time evolution in w. There is a discrepancy between the two times at r ∼ 5 × 10³ AU. At these small scales, fragmentation does increase in time, as high-density regions have more time to form and fragment.

Zoom In Zoom Out Reset image size

It should be noted that at t ∼ 0.4t_ff, our stars have not provided very much feedback, and the simulation can be considered solely gravito-hydrodynamic. Simulations with just hydrodynamics and gravity are scale-free. The exact match with ρ Ophiuchus is partially due to our choice of cloud parameters to mimic ρ Ophiuchus. The scale-free results are the break in w(r) at the Jeans length the w(r)∝r^−0.75 functional form.

There is a wealth of observations cataloguing masses of cores in star-forming regions (Motte et al. 1998; Testi & Sargent 1998; Johnstone et al. 2000, 2001; Motte et al. 2001; Beuther & Schilke 2004; Stanke et al. 2006; Alves et al. 2007; Enoch et al. 2008; Sadavoy et al. 2010). Most of these observations are unable to resolve the small-scale fragmentation seen in ρ Ophiuchus, but still provide valuable information. To recreate the observations, we produced optically thin column density maps of our simulations in all three directions and convolved them with Gaussian beams chosen to match the observations. We then applied Clumpfind to the convolved column density, similar to our comparison to ρ Ophiuchus.

The observations have a wide range of beam sizes due to the range in distances to star-forming regions, so we also used a range of beam sizes for comparison. We find that the CMF derived from our simulated observations is highly sensitive to the beam size used. As the beam size increases, the smallest cores are no longer detectable and drop out of the CMF. In addition, tight clusters of cores become unresolved and look like new, much larger cores. Both effects increase the median clump mass. The effect of overlapping cores has been explored further in Kainulainen et al. (2009b). This sensitivity of the CMF to resolution is also seen in the observations and shown in Figure 18. To gather the median mass of a range of CMFs, we used the data tabulated in Reid & Wilson (2006). All cores from Reid & Wilson (2006) are detected using Clumpfind. To supplement that data, we also used the three clouds (Serpens, Perseus, and ρ Ophiuchus) from Enoch et al. (2008). Enoch et al. (2008) do not rely solely on Clumpfind, but use a method that returns similar results. Clumpfind has many limitations (Pineda et al. 2009; Goodman et al. 2009), but is still useful for the purposes of comparison. When compared to these observational data, our simulated CMFs match quite well. The simulated CMF with 3200 AU and 6400 AU beams are interesting in particular, as this is the usual range of telescope beams for nearby clouds. The extent of our simulation domain (130,000 AU) prevents useful comparisons to more distant observations. In addition, our base grid resolution, 512 AU, causes discrepancies with observations of cores at small beam sizes. Bound cores will trigger adaptive refinement and can go to smaller scales, but smaller, ephemeral, cores that would show up in observations are not always captured in the simulations. This means the simulated median masses at 800 AU and 1600 AU are too high compared to observations which include all cores, bound or not. This is most apparent when comparing to the observations of ρ Ophiuchus with a beam size of 1600 AU. The net effect is a flattening in the clump mass versus beam size relation at small beam sizes for the simulated clump observations.

Zoom In Zoom Out Reset image size

Figure 18 does not include the CMF from the Pipe Nebula (Alves et al. 2007). The resolution is comparable to the best observations of ρ Ophiuchus, but the median mass is near 1 M_☉. This makes it stand apart from the other observations. This region has a low density and is not actively forming stars. In addition, the observed cores are largely unbound (Lada et al. 2008). The Pipe Nebula should perhaps be considered a pre-star-forming region, unlike Serpens or ρ Ophiuchus. These pre-star-forming regions have much higher Jeans lengths and masses and have fundamentally different column density distribution functions (Kainulainen et al. 2009a).

Given that we have identified clumps at the beginning of star formation and have stars at the end of the simulation, a natural question is how the initial clumps relate to the final stars. We will first consider the cores identified by clumpfind discussed in the previous section. These cores represent what an observer might find in a similar cloud. For this comparison, we need to use simulations that have run to completion, so cannot use the radiative simulations. As a first pass at the correlation between the CMF and IMF, we will start with simulation B, where winds are absent and one might expect the CMF and IMF to overlap. The initial CMF and the final IMF for B are shown in Figure 19. For proper comparison, the mass functions are not normalized, and total counts at each mass are shown. To avoid triple counting cores in the CMF, the synthetic observations are taken in only one line-of-sight direction instead of all three. Otherwise these CMFs are the same as those in previous sections. When focusing on the CMF found with the 1600 AU beam size, the initial CMF and final IMF are well correlated. The typical mass and total number of objects match well between cores and stars. Unfortunately, this comparison only holds for the 1600 AU beam size. The CMF is highly sensitive to beam size, as shown in Figure 18. When a slightly smaller, 800 AU, beam size is used, the clumps are too small and too numerous to all correlate with stars. It should be noted that this analysis is performed using a single clump identification method (Clumpfind) on a crowded field. It is possible that different methods of clump identification or clouds with well-separated clumps do not have as much sensitivity to telescope resolution in the resulting CMF. Even without resolution effects, we are considering all cores, bound and unbound. Low-mass cores are less likely to be bound or to form stars (Myers 2009; Padoan & Nordlund 2011), which will skew the CMF to lower masses. CMFs that only include bound cores should be more correlated with the IMF.

Zoom In Zoom Out Reset image size

We do not need to limit ourselves to synthetic observations of our simulations and can identify cores from the full six-dimensional phase space. This is done using the "find_clumps" routine in the yt analysis toolkit (Turk et al. 2011). The algorithm is described in more detail in Smith et al. (2009). It uses density contours to return a hierarchy of clumps, where each clump can contain smaller child clumps. Our simulations have hundreds of local maxima in density large enough to be considered clumps. At t = 0.3t_ff, only 40 of these are bound, defined as |potential energy| > (thermal energy + kinetic energy). We will only consider the bound clumps in this analysis.

The bound clumps cover a large range of sizes with the largest clump filling almost the entire box (181 M_☉). This clump contains five child clumps. Each of these clumps contain many generations of children and grandchildren. When comparing to the IMF, one should ignore clumps with multiple child clumps and count the children instead. Throwing out clumps with masses less than 0.05 M_☉, there are only five bound childless clumps at t = 0.3t_ff, for a total mass of 2.5 M_☉. At t = 0.4t_ff, some of the parent clumps split into more children, resulting in 16 bound childless clumps, though the total mass is largely unchanged, at 2.6 M_☉. These clumps are notably smaller than the primary turbulent cores described in Section 3.1 and will be called "sub-cores." These sub-cores lead to the burst of star formation at t = 0.5t_ff and their mass function is shown in Figure 20. The median mass of these sub-cores is 0.13 M_☉. Even if each sub-core forms exactly one star, they do not explain the stellar mass function of the simulation, which eventually has a median mass of 0.3 M_☉ when winds and radiation are not included. This discrepancy can be explained by the 20, 000 AU turbulent cores. These objects are bound, but they cannot be detected with density contours due to their supersonic turbulence. Each core has multiple pockets of high-density gas. When using density contours, one sees many smaller unrelated, unbound clumps instead of a few larger bound clumps that correspond to the physical cores. The cores are visible by eye and should be identifiable with a more advanced density search. We identify them by stellar clustering. The central stellar systems (usually binary systems) in these cores are all much more massive than the sub-cores. If these central stellar systems are ignored, the median mass of stars moves from 0.3 M_☉ to 0.16 M_☉, much closer to the median sub-core mass. The median mass of non-central stars for the case with winds is 0.07 M_☉. The turbulent core and central star properties are summarized in Table 3.

Zoom In Zoom Out Reset image size

The cores are bound at late times even when only gravity due to the stars is considered, with the kinetic energy in stars approximately half the potential energy of stars in each core. In addition, the core-to-core velocity dispersion is typically 0.4–0.5 km s⁻¹. This is notably lower than the cloud velocity dispersion, which starts at 1.2 km s⁻¹. There is a strong anti-correlation between velocity and density, as the densest gas occurs at stagnation points in a turbulent flow. This means that the core-to-core velocity dispersion will naturally be much lower than that of the gas (Offner et al. 2009a).

Once the cores have formed, each core is carved out by the outflows of its own protostellar system. This yields the core to star efficiency factor, 0.2 < _core < 1.0. The amount of mass lost from a spheroidal core can be calculated from the total momentum output and opening angle of the winds (Matzner & McKee 2000), but the cores in our simulations are more complicated.

The best way to calculate _core is to compare the mass of stars in simulations with and without winds. In the simulations with the barotropic equation of state (B and BW), the total, the mean, and the median mass of stars are all approximately 3:1 comparing the non-wind simulation to the wind simulation at any point in time. This means _core ≃ 1/3. The 2/3 of the core mass that is lost in the presence of winds is mostly core gas that is entrained in the winds. This was not precisely quantified, and it is possible that other mechanisms account for some mass loss, such as unbinding outer core gas when the mass from the interior core is lost. In the radiative simulations (R and RW), the total mass of stars is also 3:1 comparing the two simulations. The mean and median masses are 3:1 for simulation RW, but closer to 10:1 in simulation R due to unphysically strong radiative suppression of fragmentation.

In the event that cores are accreting mass, _core is a function of time, and may be larger than 1 if the instantaneous core mass is less than the final stellar system mass (Padoan & Nordlund 2011). Given that the sub-cores in our simulations are massive enough at early times to create their final stellar systems, the sub-cores probably do not accrete much mass. It is extremely difficult to track these cores over time, however, so we cannot quantify the amount of sub-core accretion.

It is useful to place the stellar accretion in our simulations in the context of existing star formation models. Two current popular models are turbulent fragmentation and competitive accretion. In the turbulent fragmentation model, (Padoan 1995; Padoan & Nordlund 2002; McKee & Tan 2003; Padoan et al. 2007; Hennebelle & Chabrier 2008, 2009), supersonic turbulence in molecular clouds creates many cores. Each bound core then collapses into a single stellar system (if the structure of the core is dominated by internal turbulence, the resulting model for the formation of the star was termed the "turbulent core model" by McKee & Tan 2003). In this scenario, the mass from each star is accreted almost entirely from its natal core. In the competitive accretion model (Zinnecker 1982; Bonnell et al. 1997, 2001; Bate & Bonnell 2005; Bonnell et al. 2006), the bound cores begin with a mass ∼0.1 M_☉. The molecular cloud undergoes a global collapse and all stars accrete from the entire cloud. Protostars exhaust their cores at low-masses and then grow by Bondi–Hoyle accretion. The protostars compete with each other for mass from the host cloud and the dynamics of this competition lead to the IMF. Roughly speaking, the virial parameter of the cloud decides which model is applicable (Krumholz et al. 2005; Bonnell & Bate 2006; Offner et al. 2008). In clouds with sufficiently sub-virial turbulence, global collapse is possible and competitive accretion prevails. In virialized clouds, core accretion dominates. Most simulations of star formation start with virialized clouds, but turbulence quickly dissipates and simulations that do not regenerate turbulence become sub-virial and demonstrate competitive accretion. Turbulence can be regenerated by protostellar outflows, H ii regions, or a cascade from larger scales; in simulations, the cascade can be regenerated by large-scale driving.

Our simulations largely agree with turbulent fragmentation models, while introducing a hierarchical aspect of sub-cores within turbulent cores. We form turbulent cores on the cloud Jeans length and each of those cores forms a central binary or single star with mass roughly equal to the core mass. Even the smaller stars are formed in their own sub-cores and do not accrete from the cloud at large. Our turbulent cores do accrete from the larger cloud (increasing their initial mass by ∼75% over the course of the simulations), which was not originally part of the core accretion theory, but recent turbulence simulations suggest turbulent cores do accrete from their host cloud (Falceta-Gonçalves & Lazarian 2011).

The accretion from the larger cloud onto cores in our simulation is possibly caused by the fact that we have not included turbulent driving from large scales. Turbulence is regenerated to some extent by protostellar outflows, but this is relatively ineffective in denser gas. Our simulations are then similar to the undriven simulations in Offner et al. (2008), which also show accretion onto cores. The simulations in Offner et al. (2008) with external driving produce cores that do not accrete much from the cloud. If magnetic fields were included in our simulations, the outflows would couple to much more of the gas (Wang et al. 2010), which would reduce accretion onto the cores.

The stellar mass functions in our simulations are influenced by the details of our sink particle merger process. Mergers are necessary because all codes will introduce numerical fragmentation once they can no longer resolve the Jeans length on the finest scale (Truelove et al. 1997). More stringent sink particle conditions can reduce the number of unwanted sink particles (Federrath et al. 2010), but numerical fragments are unavoidable. If sink particles are not allowed to merge, these numerical fragments will steal mass from physical fragments and masquerade as real stars, artificially lowering the IMF. On the other hand, allowing sink particles to merge has a similar effect on the IMF as suppressing gas fragmentation. If all sink particles that pass near each other are merged together, the particles will consolidate over time. Eventually, all IMFs look similar to the top heavy IMF from R, where radiation suppressed most fragmentation. Our decision to only merge protostars with masses less than 0.05 M_☉, the mass of the first core of a 1 M_☉ star, is a compromise between the two extremes of no mergers and all mergers.

This suppression of mergers is seen comparing the IMFs of RW and BW. The radiative simulation should suppress some fragmentation while the barotropic simulation fragments all the way down to the Jeans mass at the self-opacity limit, m ∼ 0.004 M_☉. Nevertheless, they have the same IMF. The barotropic simulation does fragment much more than the radiative simulation near the resolution limit. This is seen in the total number of sink particles created, where the barotropic simulation creates seven times more particles than the radiative simulation. These particles are nearly all very small in mass and immediately merged. The net effect of the extra mergers is to suppress fragmentation by combining fragments below the merger radius (128 AU). The two simulations are nearly identical above this scale and therefore produce the same IMF.

Our choice of merger mass (m_merge = 0.05 M_☉) is based on calculations of first collapse of a solar mass star (Masunaga et al. 1998; Masunaga & Inutsuka 2000), but the correct mass is not certain. In addition, smoothed particle hydrodynamics simulations with resolutions of a few AU (compared to our 32 AU) move their numerical fragmentation to much smaller scales and do not show stellar mergers of this type (Stamatellos & Whitworth 2009; Bate 2012). To investigate the effect of our mass choice, we repeated RW with a merger mass of 0.01 M_☉ out to t = 0.55t_ff. At this point in the simulations, the total mass in stars is 2 M_☉. In the simulation with m_merge = 0.05 M_☉, there are 17 total sink particles; whereas in the simulation with m_merge = 0.01 M_☉, there are 30 total sink particles, even though the total mass in particles is the same. Given that the final IMF for m_merge = 0.05 M_☉ is at slightly lower masses than the observed IMF, nearly doubling the number of sink particles with m_merge = 0.01 M_☉ would move the IMF to masses less than half the observed value.

For most isothermal simulations, there is a free parameter which allows rescaling of the mass, M, density, ρ, or temperature T, while maintaining all of the same dimensionless parameters α, $\mathcal {M}$, and M_J/M. When the merger mass is included, m_merge also scales with M. Now if we scale the simulation to a higher mass, we are also increasing m_merge. There is some uncertainty in m_merge, but it would be difficult to justify increasing it much more than our current level. Even increasing m_merge by a factor of two would bring it to uncomfortably close to the characteristic IMF mass of 0.2 M_☉. This means that the mass scales of our simulation are relatively stationary. When protostellar winds are included, they introduce a new fixed dimensionless number, the Mach number of the winds $\mathcal {M}_{{\rm wind}}$. Since the speed of the winds is proportional to the escape velocity from the stellar surface (i.e., ∝M^1/2 McKee & Ostriker 2007), $\mathcal {M}_{{\rm wind}}$ sets the quantity M/T. In practice, this is not a very tight constraint because the wind speed itself is quite uncertain (Downes & Cabrit 2007). When radiation is important, the luminosity of each star is set by complicated stellar models that depend on M as well as the accretion history. The timescale, which goes into the accretion rate, is set by t_ff and therefore ρ. In addition, the resulting radiation hydrodynamics depends on the temperature. This uniquely sets all scales in the simulations. Even when radiation is not dynamically important, we do match the observed protostellar luminosities and cannot change our masses without jeopardizing the agreement. Using the approximation L∝M/t_ff∝M^3/2, our cloud mass is constrained to 160 M_☉ < M < 200 M_☉ before it no longer falls in the error bars luminosities. Cloud masses below the lower limit do not match the observed median luminosity and masses above the upper limit do not match the mean. Our general match to the IMF, which sets M, provides an additional, looser constraint.