RADIATION-HYDRODYNAMIC SIMULATIONS OF MASSIVE STAR FORMATION WITH PROTOSTELLAR OUTFLOWS

Our simulations are initialized to the state of a prestellar core of mass M, each with a power-law density profile given by

$$\begin{equation} \rho (r) \propto r^{-k\rho }, \end{equation} \tag{ 1 }$$

with k_ρ = 3/2, consistent with models (McKee & Tan 2002, 2003) and observation (Beuther et al. 2007), and an initial temperature T_c = 20 K. The average density, initial radius, and free-fall time of the initial core are set by the initial volume average core surface density

$$\begin{equation} \Sigma = \frac{M}{\pi r_c^2} \end{equation} \tag{ 2 }$$

as

$$\begin{equation} \rho _c = \sqrt{\frac{9 \pi \Sigma ^3}{16 M}} \end{equation} \tag{ 3 }$$

$$\begin{equation} r_c = \sqrt{\frac{M}{\pi \Sigma }} \end{equation} \tag{ 4 }$$

and

$$\begin{equation} t_{\rm ff} = \left[\frac{\pi M}{64 G^2 \Sigma ^3}\right]^{1/4}. \end{equation} \tag{ 5 }$$

The initial core is placed in the center of a cubical simulation domain spanning four times the core radius, i.e., L_domain = 4r_c on each side, so that no part of the initial cloud except gas that is entrained into protostellar winds ever approaches the boundary. The initial core is immersed into a uniform ambient environment with density that is 0.01 times that at the edge of the initial core. Pressure balance between the core and environment is maintained by setting the temperature of the ambient gas to 100 times that at the edge of the initial core, T_amb = 2000 K. The Planck mean opacity of the ambient gas is set to zero to ensure that it does not cool or radiatively heat the core. The cores are initialized with a turbulent velocity field chosen to put them in approximate balance between gravity and turbulent motions. Three Gaussian random fields are generated with power spectrum P(k)∝k⁻² for the three velocity components, each normalized to have an integrated norm of unity over the full spectral range sampled. We set the initial velocity in every cell equal to the components of the Gaussian random field times the one-dimensional velocity dispersion,

$$\begin{equation} \sigma _v = \sqrt{\frac{GM}{2(k_\rho -1)r_c}} = \left[\frac{G^2 M \pi \Sigma }{4(k_\rho -1)^2}\right]^{1/4}, \end{equation} \tag{ 6 }$$

corresponding to the velocity at the surface of a singular polytropic sphere (McKee & Tan 2003). The virial parameter α_vir = 5σ²_vGM/r_c (Bertoldi & McKee 1992) is 5 for k_ρ = 3/2, somewhat larger than the value of 15/4 in hydrostatic equilibrium (McKee & Tan 2003). The kinetic energy is therefore initially larger than the gravitational energy, but α_vir decreases with time due to the decay of the turbulence. The initial radiation energy density is set to the value for a blackbody radiation field with radiation temperature T_r = 20 K everywhere.

The earliest stages of high-mass star formation occur in infrared dark clouds (IRDCs; Rathborne et al. 2007) which are detected in absorption against the mid-infrared galactic background (Perault et al. 1996; Egan et al. 1998). Observations of IRDCs indicate a range of surface density in star-forming regions from more tenuous sources with Σ ∼ 0.1 g cm⁻², to more typical galactic star formation conditions with Σ ∼ 1 g cm⁻², to Σ ∼ 10 g cm⁻² (Beuther et al. 2002b; Rathborne et al. 2006; López-Sepulcre et al. 2010) or more in extragalactic superclusters (Turner et al. 2000; McCrady & Graham 2007). Table 1 summarizes the parameters for each of the four computational models presented in this work. The initial conditions for these models have been chosen to study the collapse of galactic IRDCs with high but not atypical mass and surface density. Each of the initial simulation core states are rescaled versions of one another with identical density structure, virial ratio, velocity field, and comparable peak resolution in every run. We have compared each of the simulations at equivalent time in units of the free-fall time of the initial cores. The homology between the runs is broken only by the presence or absence of outflows and by radiative effects. This choice of model parameters therefore probes the surface density dependence of radiative feedback effects, and isolates the effects of protostellar outflows by holding all other parameters constant as they are turned on and off.

Table 1. Simulation Parameters

Σ(g cm−2)	1.0	2.0	2.0	10.0
Wind	on	on	off	on
$M\, (M_{\mathord \odot })$	300	300	300	300
rc (pc)	0.141	0.100	0.100	0.0447
$\bar{n}_H\, ({\rm cm}^{-3})$	7.3 × 105	2.1 × 106	2.1 × 106	2.3 × 107
σv/cs	8.80	10.5	10.5	15.6
tff (kyr)	50.7	30.2	30.2	9.02
Ldomain (pc)	0.565	0.400	0.400	0.179
N0	128	192	192	192
Max level	5	4	4	3
ΔxL (AU)	28.4	26.8	26.8	24.0

Notes. Row 1: initial core surface density; Row 2: initial core mass; Row 3: initial core radius; Row 4: initial core velocity dispersion; Row 5: core free-fall time; Row 6: linear size of the computational domain; Row 7: number of cells per linear dimension on the coarsest level; Row 8: maximum refinement level; Row 9: computational resolution on the finest AMR level.

Download table as:  ASCIITypeset image

We have largely followed the approach set forth by Krumholz et al. (2010) in choosing the parameters for the numerical simulations considered here. However, the simulations in this work differ from the earlier work in several ways. First, these simulations use an initial core mass of $300 \,{M}_{\mathord \odot }$ instead of $100 \,{M}_{\mathord \odot }$. This choice was motivated by our desire to study the evolution of high-mass star systems and our expectation that protostellar outflows, which were not considered in the earlier models, will eject a significant fraction of the initial core. Second, the lowest initial surface density is Σ = 1.0 g cm⁻² instead of Σ = 0.1 g cm⁻² as used in the earlier work. This choice is largely motivated by computational constraints. The high flow speeds present in runs with outflows necessitate smaller numerical time steps than for non-outflow runs and thus increase the computational cost. A simulation with Σ = 0.1 g cm⁻² would be particularly expensive due to its long free-fall time and the need to advance to these later times with numerical time steps that are limited by the cell crossing time of outflow-ejected gas.

The AMR capabilities of the code track the collapsing cores in three dimensions to grid scales of Δx_L ∼ 25 AU on the finest AMR level. This resolution is achieved by discretizing the physical domain on the coarse onto a base grid of N³₀ cells. The placement of finer level grids up to the finest level L was determined by the refinement criteria that any gas denser than one-half the density at the edge of the initial core be refined by at least one level. Further refinement is also triggered wherever the local Jeans number, $J=\sqrt{G \rho \Delta x^2 / (\pi c_s^2)}$ (Truelove et al. 1998), exceeds 0.125, where Δx is the computational cell width on the coarser level, or wherever the local gradient of the radiation energy density |∇E_rad|Δx/E_rad exceeds 0.1.

The simulations use a zero velocity gradient outflow boundary condition for the hydrodynamics and Marshak boundary conditions for the radiation energy density with radiation flux at the edge of the computational domain for a 20 K background radiation field. The gravitational potential is specified at the edge of the domain as the sum of the multipole moments of the mass distribution in the computational volume as a function of time up to the quadrupole term. We adopt an equation of state with γ = 5/3, appropriate for gas too cool for molecular hydrogen to be rotationally excited, but this choice is essentially irrelevant because the gas temperature is set almost purely by radiative effects.

The radiation transport is handled by a frequency-integrated flux limited diffusion approximation. We use the Planck and Rosseland mean dust opacities, κ_P and κ_R, respectively, of the Semenov et al. (2003) iron normal, composite aggregates model plotted in Figure 1. The simulations with protostellar winds introduce strong heating behind wind-driven shocks. When the thermal energy of the gas exceeds that of a gas with molecular weight 0.6m_p and temperature 10⁴ K, we treat the gas as fully ionized with Rosseland opacity κ_R = 0.32 cm² g⁻¹, the value for Thompson scattering at solar metallicity. Our gray flux limited diffusion approximation cannot adequately represent the collisionally excited line cooling processes that dominates at temperatures above the dust destruction temperature. However, it is critical to include them in order to ensure that shocked gas is able to cool. We therefore leave κ_P = 10⁻² cm² g⁻¹ for this gas, ensuring that it does not interact strongly with the ambient radiation field, but we also implement an approximate line cooling function λ(T) to remove energy from gas above the dust destruction temperature and transfer it to the radiation field. We take λ(T) from the function shown in Figure 1 of Cunningham et al. (2006). In each time step, before we perform our ordinary flux limited diffusion radiation solve, we update the gas and radiation energy densities by implicitly solving the operator split-system

$$\begin{eqnarray} \frac{d \rho e}{d t} &=& {-}(\rho /\mu)^2 \lambda (T) \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \frac{d E}{d t} &=& (\rho /\mu)^2 \lambda (T) \end{eqnarray} \tag{ 8 }$$

using the LSODE (Radhakrishnan & Hindmarsh 1993) Gear-type solver. In the above system, e is the specific thermal energy density of the gas, E is the radiation energy density, μ is the mean molecular weight, and T is the gas temperature appropriate for a solar metallicity ionized gas mixture.

Zoom In Zoom Out Reset image size

The ORION code includes a "star particle" algorithm to handle the formation of protostars (Krumholz et al. 2004, 2007a). This algorithm provides for the creation of sub-grid star particles in those cells of the computation that become poised for gravitational collapse to spatial scales smaller than those that can be captured on the computational grid without spurious fragmentation (Truelove et al. 1997). The luminosity, radius and burning state of the star particle is advanced with the simulation according to the protostellar evolution model of McKee & Tan (2003) as updated by Offner et al. (2009). The protostellar evolution model takes as input the mass and accretion history of the star as determined by the simulation, and as output predicts a protostar's radius and luminosity at any given time. The protostellar luminosity prescribed by this model enters the simulation as a source term in the radiation energy density equation, and the protostellar radius is used to compute the Keplerian velocity at the stellar surface, which affects outflows as described below. The protostellar luminosity prescribed by this model enters the simulation as a source term in the radiation energy density equation.

The ORION star particle algorithm has been enhanced for this work to include the driving of bipolar outflows. Our outflow model is specified by the dimensionless parameters f_w and f_v, which set a wind launch speed as a fraction f_v of the Kepler speed at the stellar surface and a mass flux that is a fraction f_w of the rate of accretion onto the star or, equivalently, a fraction f_w/(1 + f_w) of the total mass that is either accreted onto the star or ejected in the wind. Since we are interested in the large-scale impact of the protostellar winds, we assume that the wind is injected over a range of radii determined by a function χ_w(|r|) with an angular dependence given by $\bar{\xi }(\theta _i)$; explicit forms for these functions are given below. The wind driving is imposed by operator split source terms in the gas density, momentum density, and energy density equations with

$$\begin{eqnarray} \left. \frac{d \rho }{d t} \right|_s &=& {-}\dot{M}_{w,i} \chi _w(|{\bf r}_i|) \bar{\xi }(\theta _i), \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \left. \frac{d \rho {\bf v}}{d t} \right|_s &=& {-}f_v v_{k,i} \dot{M}_{w,i} \chi _w(|{\bf r}_i|) \bar{\xi }(\theta _i) \cdot {\bf \hat{r}}_i, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \left. \frac{d \rho e}{d t} \right|_s &=& {-}\dot{M}_{w,i} \chi _w(|{\bf r}_i|) \bar{\xi }(\theta _i) \frac{k_B T_w}{\mu (\gamma -1)}, \end{eqnarray} \tag{ 11 }$$

where

$$\begin{equation} v_{k,i}=\sqrt{\frac{GM_i}{r_{*,i}}} \end{equation} \tag{ 12 }$$

is the Keplerian speed at the surface of the star and r_{*, i} is the protostellar radius; as remarked above, we use the value of r_{*, i} given by the model of McKee & Tan (2003) as updated by Offner et al. (2009). For the simulations presented here we have set the wind-launched gas temperature as T_w = 10⁴ K, appropriate for an ionized wind. The corresponding rate of particle mass growth, particle wind mass ejection rate, acceleration, radial distance, and spatial inclination of the ith star are

$$\begin{eqnarray} \dot{M}_i &=& \frac{1}{1+f_w} \dot{M}_{{\rm KKM04}}, \end{eqnarray} \tag{ 13 }$$

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

$$\begin{eqnarray} {\bf \dot{v}}_{i} &=& {\bf \dot{v}}_{{\rm KKM04}}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} {\bf r}_i &=& {\bf x}-{\bf x}_i, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \theta _i &=& {\rm acos}(\hat{{\bf r}}_i \cdot \hat{{\bf j}}_i), \end{eqnarray} \tag{ 17 }$$

where j_i and x_i are the velocity and position of the ith particle, j_i, $\dot{M}_{{\rm KKM04}}$, and ${\bf \dot{v}}_{{\rm KKM04}}$ are the sink particle angular momentum, accretion rate, and acceleration for the case in which winds are absent as given by the algorithm of Krumholz et al. (2004).

Values of the parameters f_w and f_v can be both estimated from theory and constrained by observations. Theoretically, the X-wind (Shu et al. 1988) and disk wind (Pelletier & Pudritz 1992) models predict f_w ∼ 0.3, f_v ∼ 1 and f_w ∼ 0.1, f_v ∼ 3, respectively. Both therefore suggest f_wf_v ∼ 0.3. The total momentum P_w carried by an observed outflow from a star of mass M_* is related to f_w and f_v by

$$\begin{equation} f_w f_v = \frac{P_w}{v_k M_*}, \end{equation} \tag{ 18 }$$

where v_k is the Keplerian velocity at the stellar surface. The peak of the stellar initial mass function is at $M_* \approx 0.2 \,{M}_{\mathord \odot }$ (Chabrier 2005), and such stars typically have radii ${\sim} 3 \,{R}_{\mathord \odot }$ during their main accretion phases (e.g., model mC5H of Hosokawa et al. 2011), so typical values of v_k are ∼100 km s⁻¹. This value of v_k together with f_wf_v ∼ 0.3 implies a net wind momentum flux of ∼30 km s⁻¹ per $\,{M}_{\mathord \odot }$ of stars formed.

A number of surveys have used measurements of P_w and estimates of M_* and v_s, together with the relation given in Equation (18), to constrain f_wf_v. Surveying the literature available as of 2000, Richer et al. (2000) estimate f_wf_v ∼ 0.3. More recent observational surveys of several nearby low-mass star-forming regions indicate typical outflow momenta of ∼0.2– $\sim 3.0 \,{M}_{\mathord \odot\ }{\rm km\, s}^{-1}$ (Maury et al. 2009; Arce et al. 2010; Curtis et al. 2010; Ginsburg et al. 2011). The physical properties of the driving sources of most of the surveyed outflows are not very well constrained. However, if we assume that the typical source has accreted half of its final mass $M_* \sim 0.1 \,{M}_{\mathord \odot }$ and has radius $r_* \sim 3 \,{R}_{\mathord \odot }$, then Equations (12) and (18) can be used to extract a range of outflow momentum parameters based on the results of these surveys of 0.025 ≲ f_wf_v ≲ 0.38.

Observationally, f_w and f_v can be better constrained from sources where observational measurements exist for both net outflow momentum and the net mass accreted onto the protostar M_*. Curtis et al. (2010) have surveyed the outflow momentum in the young cluster NGC 1333 and Sandell & Knee (2001) have estimated the mass of warm dusty gas in the collapsing envelope around the deeply embedded protostars that drive several of these outflows. In Table 2 we list the intersection of those sources that have both dust mass measurements from Sandell & Knee (2001) and outflow momentum measurements by Curtis et al. (2010), excluding a few sources that Sandell & Knee (2001) indicated as near the edge of their field of view with unreliable flux densities, and excluding one tight binary source that could not be separately resolved in that work (IRAS 4A). We expect that the envelope masses of early class 0 sources should be somewhat greater than the mass of the embedded protostar and we expect that the mass of later class I sources should exceed that of their envelope. In the absence of better mass constraints for the collection of class I and 0 protostars in the present sample we adopt the assumption M_* ∼ M_Dust as a rough approximation. We assume fiducial stellar radius of $r_* = 3 \,{R}_{\mathord \odot }$ following model mC5H of Hosokawa et al. (2011). From these assumptions, we can constrain f_wf_v from Equations (18) and (12). The datum indicates a range of outflow launch parameters of 0.01 ≲ f_wf_v ≲ 0.15, with no strong statistical correlation of f_wf_v with the source spectral class.

Wind launch speeds represent the Courant time-step constraint in a typical calculation, so large values of f_v impose a particularly onerous computational burden. We therefore choose a wind mass to stellar mass fraction on the high end of the theoretical guidance, f_w = 27% and the wind velocity parameter of, f_v = 1/3. This yields a momentum flux injected by our wind model characterized by f_wf_v = 9%, which is toward the higher end of the observed range of rates of momentum injection by outflows from the low-mass sources tabulated in Table 2.

The function

$$\begin{equation} \chi _w(r) = \frac{1}{C_1}\left\lbrace \begin{array}{@{}ll}r^{-2} & {\rm if }\; 4 \Delta x < r \le 8 \Delta x \\\\ 0 & {\rm otherwise} \end{array} \right. \end{equation} \tag{ 19 }$$

is a normalized weighting kernel that determines the spatial scale of the wind injection where C₁ is a normalization constant to the weighting kernel. C₁ is computed numerically in the ORION code so that numerical aliasing effects of the spherical wind injection region into the Cartesian grid are accounted for exactly.

The remaining function $\bar{\xi }$ describes the angular distribution of the wind mass and momentum flux at the point where it is injected into the computational grid. We take this function from Matzner & McKee (1999), who find that the momentum distribution of prestellar outflow injection asymptotically far away from the protostellar surface as a function of the polar angle θ from the direction of the protostar's rotation is given by

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

where θ₀ is the so-called "flattening parameter" that sets the opening angle of the wind. In the case of low-mass stars, Matzner & McKee (1999) suggest a fiducial value of θ₀ = 0.01 and we adopt the same value here. While stars of type B or later direct momentum in a very well collimated beam (Beuther et al. 2002a, 2003, 2004), O star winds are collimated somewhat more weakly, possibly due to the effect of ionization (Beuther & Shepherd 2005). Since we do not include ionizing radiation in these simulations, we do not attempt to model O star wind broadening.

The large value of ξ near θ = 0 in Equation (20) requires particular care in implementing this model in a numerical code. We implement the driving function by averaging ξ over the polar angle subtended by a grid cell Δθ = atan(1/8) at the outer radius of the weighting kernel (Equation (19)) as

$$\begin{equation} \bar{\xi }(\theta,\theta _0) = \frac{1}{C_2} \left\lbrace \begin{array}{@{}ll}\dsty \frac{1}{\Delta \theta }\int _{\theta -\Delta \theta /2}^{\theta +\Delta \theta /2} \frac{d\theta }{\sin ^2 \theta +\theta _0^2} & {\rm if }\, \big|\sin (\frac{\pi }{2}-\theta)\big| {\,\ge\,} \frac{\Delta x}{r}\\\\ 0 & {\rm otherwise}. \end{array} \right. \end{equation} \tag{ 21 }$$

Our choice to set $\bar{\xi }$ to zero for angles close to π/2 is driven by numerical considerations. If we allow the outflow to be injected into 4π sr around the star, its mass and momentum are sufficient to disrupt the early development of an equatorial disk. This behavior is an artifact of the necessarily poor resolution inside the wind launching region. We do not resolve the disk scale height, and this artificially puffs up the disk and reduces its mass and momentum density, rendering it far easier for the outflow to disrupt than it would be if its true scale height were resolved. We avoid this problem by reducing the outflow mass and momentum flux to zero in an equatorial belt that is at least one cell thick, ensuring that accreting material always has an uninterrupted path to the star. We note that Schönke & Tscharnuter (2011) have considered the effect of radiative and protostellar outflow feedback on the dynamics of the accretion disk in two dimensions at much higher resolution than the simulations considered here. Their simulations indicate that feedback effects can alter the accretion rate onto the star on shorter timescales and smaller length scales than have been resolved in this study. However, the emphasis of the present study is on the large-scale radiative and feedback effects on the ambient core and we do not attempt to model this small-scale behavior.

The definite integral in Equation (21) evaluates to

$$\begin{eqnarray} &&\int _{\theta -\Delta \theta /2}^{\theta +\Delta \theta /2} \frac{d\theta }{\sin ^2 \theta +\theta _0^2}\nonumber\\ && \quad\,=\, \frac{1}{\theta _0 \sqrt{\theta _0^2+1}} \left[ {\rm atan}\left(\frac{\sqrt{\theta _0^2+1} \tan \left(\theta +\frac{\Delta \theta }{2}\right)}{\theta _0}\right)\right.\nonumber\\ && \quad\, \quad\, \left.-\, {\rm atan}\left(\frac{\sqrt{\theta _0^2+1} \tan \left(\theta -\frac{\Delta \theta }{2}\right)}{\theta _0}\right) \right]. \end{eqnarray} \tag{ 22 }$$

The normalization constant C₂(θ₀) = ∫ξ(θ, θ₀)χ_w(r)d³x is also computed numerically to exactly account for grid aliasing effects. Neglecting the grid aliasing effect, we find C₂ = 8.165 by numerical integration. A second subtlety that arises in the numerical implementation of the wind driving is that care must be taken so that the momentum source terms impart exactly zero net momentum onto the star particles and the gas in the computational domain. If the position of the particle is allowed to vary continuously within its host cell the $\frac{1}{r^2}$ term in Equation (19) may lead to an asymmetric driving. To overcome this problem, we bring the wind driving into symmetry with the numerical grid by rounding the particle position to the nearest half-integer multiple of the grid width

$$\begin{equation} {\bf x_i} \leftarrow 2 \Delta {\bf x}\, {\rm nint}\left[ \frac{{\bf x_i}}{2 \delta {\bf x}} \right] \end{equation} \tag{ 23 }$$

before computing the source terms (Equations (13)–(10)) where nint is the "nearest integer" function.

The wind momentum injection $\bar{\xi }$ is significantly broadened at polar latitudes relative to the analytic prescription for ξ. Consequently less momentum is injected near the equator to satisfy the normalization constraint. To facilitate comparison of our numerical models with the analytic predictions of Matzner & McKee (2000) in Section 3.5 it is convenient to define the effective numerical flattening parameter θ_{0, eff} over a range of angles that separate the outflowing gas from ambient gas:

$$\begin{equation} \xi (\theta,\theta _0=\theta _{0,{\rm eff}})= \bar{\xi }(\theta,\theta _0=0.01). \end{equation} \tag{ 24 }$$

We find that the above expression holds to with 10% for θ > 10° with θ_{0, eff} = 5.75 × 10⁻⁴.

The left column of Figures 2 through 5 show the large-scale evolution of each simulation from t = 0.2t_ff to t = 0.8t_ff in increments of 0.2t_ff. These plots show slices of density with the color-mapping scaled by Σ^3/2 following the scaling given in Equation (3). The spatial scale shown in each plot scaled by the initial core radius with each showing an area (2.5r_c)². The slices are centered on the position of the primary protostar and oriented so that the angular momentum vector of the protostar is upwardly oriented on the page. As expected, once the scaling of cloud radius and surface density is accounted for, the regions that have not been penetrated by the outflow bow shocks collapse in a homologous manner with little dependence on the initial surface density. The protostellar outflows evacuate a shock-bounded cavity through the initial cores and the outflow cavities are the prominent features in the density slices in the left column by t ∼ 0.3t_ff in the lower surface density cases and t ∼ 0.4t_ff in the high surface density case. The propagation speed of the outflow bow shocks through the core and the width of the outflow cavities show a strong dependence on initial surface density. The lower surface density cores show significantly greater disruption due to the protostellar outflow feedback. This effect is due to (1) greater mechanical luminosity of the protostellar winds at lower surface density (an effect that we discuss in detail in Section 3.3) and (2) lower turbulent (σ²_v∝Σ^1/2) and thermal pressures (P_c∝ρ_cT_c ∼ Σ^3/2) in the ambient cores which act to confine the propagation of the outflows. We note that the outflow-evacuated cavities in the cases with outflows emerge from the initial core toward the left side of the density slices in Figure 5. The reason for this is that the primary star particle retains the momentum of the material that it accretes and therefore the star drifts away from the center of mass of the initial core with time (see Section 3.3). The outflowing material therefore emerges first from the thinnest side of the initial cloud, relative to the position of the primary star.

Zoom In Zoom Out Reset image size

The third columns of Figures 2 through 5 show the small-scale evolution of the surface density. Each of the plots are scaled by the initial surface density Σ and the initial core radius, with each showing an area (0.1r_c)² centered on the position of the primary star. Deviations from the homologous scaling with surface density exhibited on larger scales appear by t = 0.2t_ff. By t = 0.4t_ff in Figure 3 and t = 0.6t_ff in Figure 4, notable differences in the disk around the primary protostar emerge. Progressing from the third row, showing the case with highest surface density, to the top row, showing the case of lowest surface density we note an increasing tendency of the disk around the primary star to fragment and generate spiral arms characteristic of a Toomre-unstable disk (Toomre 1964). In the simulations with lower surface densities in the top two rows of Figure 4) we note a prominent increase in disk fragmentation in terms of the multiplicity of star particles and the presence of distinctly separate accretion disks around the primary and secondary protostars in comparison to the highest surface density case in the third row.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The trend toward reduced disk stability and enhanced fragmentation with lower initial surface density is consistent with earlier analytic (Krumholz & McKee 2008) and numerical (Krumholz et al. 2010) works that predict a threshold core surface density for sufficient radiative heating to inhibit disk fragmentation of Σ ∼ 1.0 g cm⁻², as well as with observational data from infrared star-forming clouds that are consistent with this prediction (López-Sepulcre et al. 2010). We will show in Section 3.4 that protostellar outflows provide a mechanism for focusing radiative feedback in the poleward directions, away from the infalling disk in the midplane, as predicted by Krumholz et al. (2005). Protostellar outflows should therefore raise the surface density threshold for high-mass star formation. The simulations presented here do not survey sufficiently low surface densities to quantify this effect, due to the computational cost of simulating protostellar winds inside low surface density cores with commensurately longer free-fall times. Furthermore, we expect that the relative importance of this effect may depend on magnetic fields, which our simulations neglect, and their role in confining the outflow cavity (Hennebelle et al. 2011). We therefore defer more precise quantification of the effect of outflows on the minimum surface density threshold for high-mass star formation to future work that will include magnetic fields and survey lower surface density cores.

We can isolate the effect of protostellar outflows on the small-scale evolution by comparing the second and fourth rows of Figures 2 through 5 which show the results of our numerical experiments with and without protostellar outflow ejection, both with the same initial surface density of Σ = 2.0 g cm⁻². By t = 0.4t_ff enhanced radiation trapping in the case with outflows has lead to substantially warmer circumstellar gas in comparison to the case without outflows. The temperature structure in the Σ = 2.0 g cm⁻² case without protostellar outflows in the fourth row is more similar to the Σ = 10.0 g cm⁻² case with outflows in the third row than to the case at the same surface density with winds in the second row. As we will show in Section 3.4, protostellar outflow cavities carve a path of reduced optical depth through the initial core that channel radiative flux away from the center of the core. This escape of radiative energy reduces the efficacy of radiative heating in the central regions of the collapsing core.

Zoom In Zoom Out Reset image size

Enhanced disk fragmentation associated with outflows and decreasing surface density is due to the reduced effectiveness of radiative heating. The right columns of Figures 2 through 5 show the mass-weighted, column-averaged radiation temperature T_r over the same regions as the column density projections in the center column. We note that the dense infalling gas is strongly coupled to the radiation field, so that the gas temperature T ≈ T_r. Only the tenuous outflow-evacuated regions achieve sufficient post-shock temperatures to break the radiation-gas coupling by increasing the gas temperature beyond the dust destruction temperature. The column-averaged radiation temperature plots therefore allow us to probe the temperature structure of the dense infalling gas without confusion from the projection of shock heated layers along the surfaces of the outflow cavity. The temperature structure of the gas on small scales shows even stronger dependence on the initial surface density than the column density. We note enhanced temperature with increasing surface density as early as t = 0.2t_ff when no stars are producing significant power via nuclear burning. The dependence of temperature on surface density at these times is solely due to the factors pointed out by Krumholz & McKee (2008): (1) higher surface density cores have higher accretion rates, thereby generating higher accretion luminosity and (2) higher surface density cores have higher optical depth to more effectively trap radiation. The trend toward increasing gas temperature with increasing surface density is also present at all later times, due to additional radiative output from nuclear burning in the primary star in addition to enhanced accretion luminosity and enhanced radiative trapping.

Stars with masses ${\gtrsim} 15\,{M}_{\mathord \odot }$ generate sufficient radiation pressure to exceed their gravitational acceleration. The second columns of Figures 3–5 show that by then the spherically averaged radiation force from the central protostar exceeds the inward gravitation attraction acting on the dusty envelope of infalling gas in the case of Σ = 10.0 g cm⁻² and in the case of Σ = 2.0 g cm⁻² without outflows. In the case without outflows this strong radiative force drives the expansion of a bubble of circumstellar gas away from the central source. The early development of this bubble can be seen in the lower left panel of Figure 3 at t = 0.4t_ff and the radial extent of the bubble grows to a size scale comparable to that of the initial core by t = 0.8t_ff as shown in Figure 5. The radiation bubble emerges from the initial core on the left side of the density slices in Figure 5. This is due to the drift of the primary star away from the center of mass of the initial core with time. The radiation bubble emerges first from the thinnest side of the initial cloud, relative to the position of the primary star. Accretion onto the primary star continues through the radiatively supported bubble via Rayleigh–Taylor unstable modes (Krumholz et al. 2009; Jacquet & Krumholz 2011) that develop dense, radiatively self-shielding spikes of infalling gas. The evolution of radiative bubbles in similar simulations without winds and without initial turbulence are discussed in detail by Krumholz et al. (2009). The sole difference between the radiation bubbles presented here and those in the earlier work is that the bubbles presented here are considerably less symmetrical about the central source owing to the turbulent ambient environment.

In the case with winds, the regions where the net force is dominated by the outward radiation force lie within the outflow cavity. These regions are dominated by outflow irrespective of the radiation force. With protostellar outflow, in no case does radiation force exceed that of gravity acting on the infalling core gas. Consequently, no such radiation supported bubbles form in any of the models with protostellar winds. As predicted in Krumholz et al. (2005), the cavities evacuated by protostellar outflows provide sufficient focusing of the radiative flux in the poleward directions that accretion continues through the regions of the infalling envelope onto the disk that are not disrupted by the protostellar wind shocks, and the infalling motion of this gas is not interrupted by the effects of radiation pressure.

The upper left panel of Figure 6 shows the time dependence of mass accretion onto star particles for each simulation, and the middle left panel shows the time dependence of mass accretion onto the primary protostar. Abrupt jumps in the primary protostellar mass occur due to the merger of other particles that fall toward the primary star. The code has been constructed to merge star particles that cannot be resolved on the resolution scale of four computational zones on the finest level (see Krumholz et al. 2004). Because these mergers are resolution-dependent effects, we have taken care to assure that the mass contribution to stellar sources due to mergers is small. As we have discussed earlier, the ambient core in the case of Σ = 1.0 g cm⁻² is subject to the least efficient radiative heating and is the most susceptible to gravitational fragmentation and therefore the most difficult to resolve. The particle mergers account for ∼15% of the total primary protostellar mass by the end of the simulation in the Σ = 1.0 g cm⁻² case and about ∼10% in the higher surface density cases. Most of the mass accumulated by merger events in the case of Σ = 1.0 g cm⁻² occurs due to the merger of a secondary star particle of $3.0 \,{M}_{\mathord \odot }$ at t = 0.52t_ff. We do note, however, that Myers et al. (2011) have found that imposing a maximum mass threshold for mergers enhanced fragmentation and limited the rate of mass growth of the primary star in more highly resolved but otherwise similar simulations. We therefore regard the stellar mass predictions in the models presented here as an upper limit. In comparing the cases with protostellar outflows we note that the total system mass in stars, M, exhibits a weak trend toward more rapid accretion with higher surface density even when the time is scaled by the free-fall time. In other words, $\dot{M} t_{{\rm ff}}$ increases weakly with the surface density of the initial core, Σ. This is because outflows entrain and unbind less gas from the ambient core in cases with higher surface density. Therefore, the luminosity output from the more massive protostars is enhanced in the higher surface density cases. This contributes to more effective heating and decreased fragmentation of the ambient core in the higher surface density cases as noted in Section 3.2. Consequently, the simulations at lower surface density fragment into low multiple systems earlier and are characterized by significantly reduced accretion to the primary star in comparison to the models at higher surface density.

Zoom In Zoom Out Reset image size

The upper right and middle right panels of Figure 6 show the stellar luminosity and the speed of the protostellar wind driving from the primary protostar as a function of time. Nuclear burning dominates the radiative output from stars with mass exceeding $5 \,{M}_{\mathord \odot }$ and the primary protostar dominates the stellar feedback into the core. The rapid increase in mechanical feedback from winds in all of the models is driven by the Kelvin–Helmholtz contraction of protostar (Shu et al. 1987). Stellar contraction causes the escape speed at the protostellar surface to increase which in turn increases the wind ejection speed with v_w = v_esc/3 (see Section 2.4). By t = 0.2t_ff all of the models have undergone sufficient contraction to ignite deuterium burning in the stellar cores, leading to a rapid increase in luminosity. Rapid contraction of the stellar surface continues as deuterium is exhausted in the prestellar core, giving rise to a commensurate increase in wind speed from t = 0.2t_ff to t = 0.3t_ff.

Contrary to the increase in the effects of radiative feedback with surface density, we find that the effects of mechanical feedback on the ambient cores decrease with the initial surface density. For t > 0.4t_ff the primary wind speeds show a noticeable decrease with increasing surface density. This occurs due to a mismatch between the Kelvin time for the contraction of the stellar surface and the free-fall time of the ambient molecular cloud core. The former depends mostly on the protostellar mass whereas the latter scales as t_ff∝Σ^−3/4 (Equation (5)). This means that lower surface density simulations form stars that undergo more rapid Kelvin contraction per unit free-fall time and thereby eject more powerful winds at equivalent stages of collapse. We note that this result is driven by the assumption built into our numerical model that the wind ejection speed is proportional to the escape speed at the protostellar surface. As discussed in Section 3.2, the more powerful winds contribute to the enhanced disruption of the ambient core in the lower surface density simulations. However, the overarching trend toward enhanced disruption of the ambient core with lower surface density would remain even if the wind speed as a function of core free-fall time were independent of surface density. First, the outflow break-out time from the core, t_outflow = r_c/v_w would scale, relative to the free-fall time, as t_outflow/t_ff ∼ Σ^1/4, resulting in later outflow emergence from higher surface density cores. Second, higher surface density cores enhance the confinement of outflow cavities due to higher ambient pressures.

The cascade of turbulent motions through the collapsing cloud strongly affects the motion of the primary star as a function of time. The lower left panel of Figure 6 plots the distance between the primary star, at position r_p, and the center of mass of the system, at position r_COM. The drift of the primary star away from the center of the system is a result of the random velocity field. Because most of the turbulent energy is in large wavelength modes, the denser center of the cloud will tend to have a different velocity than the lower density edges. The star initially forms in a compressional mode that is much larger than the local reservoir of gas that starts the protostar, but smaller than the diameter of the core. This initial reservoir of gas that starts the star is not co-moving with the center of mass of the core. Therefore, while the star does initially form at the bottom of the gravitational potential well it has a finite kinetic energy and can oscillate within the well. As a result the location of the star drifts from where it first forms in our simulations.

The cascade of turbulent motions through the collapsing cloud also strongly affect the orientation of the primary star as a function of time. The lower right panel of Figure 6 shows the angle of inclination of the angular momentum accreted by the primary protostar particle relative to the $\hat{{\bf z}}$ direction. We note that the total angular momentum of the initial cloud has an orientation that is inclined 27° to the $\hat{{\bf z}}$ direction. Because no radiative heating feedback is present in the initial state of the cloud, the early evolution of the simulations are characterized by fragmentation of the densest portion of the initial cloud into 2–4 gravitationally bound particles, as shown in Figure 2. By t = 0.3t_ff, the radiative feedback from the primary star is sufficient to suppress local fragmentation and the initial fragments merge. The angular momentum of the primary particle at early time varies over ∼90° as a consequence of variation in the angular momentum of the surrounding gas and as a consequence of the coalescence of the initial fragments. After t = 0.3t_ff the primary star has built up a sufficient moment of inertia relative to the rate of angular momentum deposition by accretion that the orientation of the star changes less rapidly. Subsequent evolution (t = 0.3t_ff to t = 0.8t_ff) of the orientation of the primary star is characterized by a rotation of 20°–40°. We note that the change in angular momentum in the lower surface density cases (Σ = 1.0 g cm⁻² and Σ = 2.0 g cm⁻²) occurs in abrupt jumps, whereas the higher surface density case is less susceptible to numerical fragmentation and therefore are characterized by relatively smoother change.

The cumulative distribution of stellar masses at t = 0.5t_ff for each of the runs is shown in Figure 7. Both the highest surface density case with Σ = 10.0 g cm⁻² and the case without outflows at surface density Σ = 2 g cm⁻² collapse to a single star of 70% of the total mass accreted, consistent with the qualitative similarity in the small-scale temperature structure noted in Section 3.2. The lower surface density cases on the other hand fragment into binary systems with ${>} 1 \,{M}_{\mathord \odot }$ secondaries present by t = 0.4t_ff. These results demonstrate that the absence of outflows and/or higher initial surface densities result in less fragmentation and the production of fewer, more massive stars.

Zoom In Zoom Out Reset image size

The radiative feedback from stellar sources into the ambient cloud is not isotropic. Four possible causes of anisotropic stellar output include (1) the clearing of an optically thin outflow-swept cavity along the poles of the star, (2) shielding in the midplane due to the presence of an optically thick accretion disk, (3) non-uniform optical thickness of the ambient cloud due to the displacement of the star relative to the center of mass of the ambient cloud, and (4) non-uniform optical thickness of the ambient cloud due to the turbulent density structure. In this section we consider the distribution of radiative output from the primary star in each simulation as a function of solid angle to elucidate the importance of each of these effects in shaping the radiative feedback into the ambient cloud.

Figure 8 shows the radiative flux from the primary protostar as a function of spherical angle in each model at times t = 0.3, 0.4, and 0.6t_ff, normalized to the radiative flux that would be expected in an isotropic environment ${\bf F}_{\rm isotropic} \cdot \hat{{\bf r}}= L_p/(4 \pi r^2)$. The plots shown in Figure 8 have been constructed from slices of the radial radiation flux at a distance of 1500 AU from the primary star, rotated into the coordinate system where the angular momentum of the primary star points upward. The angular distribution of the radiative flux is relatively insensitive to the position of the spherical slice, provided that the slice has radius greater than that of the disk around the primary protostar. In varying the radius of the spherical slice from 1500 AU to 10⁴ AU, we note a decreased contrast between regions of low and high radial flux by about ∼25%, which we attribute, in part, to the flux limited diffusion approximation used for the radiation transport. However, the location and extent in solid angle of the large-scale features is independent of the position of the slice. The azimuthal coordinate facing the thinnest edge of the cloud due to the motion of the star relative to the center of mass of the system is centered in each of the plots. Regions of peak outward radiative flux correlate very well with outflow ejection, and regions of low outward radiative flux correlate very well with the presence of the accretion disk near the midplane. The clearing of low optical depth paths of escape by outflow ejection and shielding by the dense accretion disk in the midplane are therefore the dominant effects in focusing the radiative feedback from the star. As gas falls in toward the primary protostar, it carries its angular momentum with it. In the highest surface density gas the accretion flow is fairly smooth because radiative heating raises the pressure near the primary star and prevents gas around it from clumping up. At lower surface densities, however, the accreting gas may be partially collapsed under its own gravity, or may even have collapsed completely to form stars that then merge with the primary. As a result, angular momentum tends to be accreted in distinct lumps, leading to rapid reorientation of the accreting star over short timescales. The radiative flux is far more focused as bipolar in the case of Σ = 10.0 g cm⁻¹, consistent with the narrower geometry of the outflow cavity in this case.

Zoom In Zoom Out Reset image size

Figure 9 shows the azimuthally averaged radiative flux from the primary protostar as a function of polar angle in each model at times t = 0.3, 0.4, and 0.6t_ff. At each of these times, the effect of the presence of protostellar outflow cavities is clearly evident showing polar radiation flux 1.7–15 times that at the midplane. The models with protostellar outflow show enhancement of the polar flux relative to the case without outflows by comparable factors. The degree of poleward focusing of the radiation flux diminishes with time due to broadening of the outflow evacuated cavities that focus the radiation. We note that polar flux in the Σ = 1.0 g cm⁻² and Σ = 2.0 g cm⁻² cases at t = 0.4t_ff underrepresents the flux focusing due to outflow cavities because the outflow cavities are tilted by ∼20° relative to the orientation of the primary protostar shown in Figure 8. Similar radiation focusing effects were also shown in Krumholz et al. (2005), where the authors considered the effects of the presence of outflow cavities on radiation escape from the infalling envelope around massive protostars. Using static radiative transfer models they showed that focusing of the radiative force from the central star throughout the outflow cavity results in a reduction of the equatorial radiative flux relative to a control model without outflows by factors of 1.7–14, depending on the width and shape of the region of low optical depth in the outflow cavity. Therefore, the models presented here support the Krumholz et al. (2005) prediction that outflows reduce the Eddington radiation pressure barrier to high-mass star formation by reducing the radiation force exerted in the infalling cloud gas. However, it has also been shown that fully three-dimensional Rayleigh–Taylor modes can remove the Eddington barrier even when protostellar outflows are neglected (Krumholz et al. 2009).

Zoom In Zoom Out Reset image size

Matzner & McKee (2000) give an analytic model for the core to star formation efficiency

$$\begin{equation} \epsilon _{\rm core}=1-M_{\rm ej}/M, \end{equation} \tag{ 25 }$$

where M_ej is the net mass ejected from the core by entrainment into the protostellar outflow and M is the initial core mass. For the case of an unmagnetized core the model predicts

$$\begin{equation} \epsilon _{\rm core} = \frac{2 X(c_g)}{1+\sqrt{1+4(1+f_w)^2 X^2(c_g)}}, \end{equation} \tag{ 26 }$$

where the dimensionless function X is defined by

$$\begin{equation} X = c_g \ln \left(\frac{2}{\theta _0}\right) \frac{v_{\rm esc}}{f_w \bar{v}_w}, \end{equation} \tag{ 27 }$$

v_esc is the escape speed from the edge of the core, and the parameter c_g depends on the core density profile and free-fall time relative to the effective timescale for the wind driving. For an unmagnetized core with profile $\rho \propto r^{k_\rho }$, Equation (A19) of Matzner & McKee (2000) provides an estimate of

$$\begin{equation} c_{g,\gg 1}=\frac{\pi (9-3k_\rho)(4-k_\rho)}{8-3k_\rho }\left(\frac{t_w}{t_{\rm ff}}\right) \end{equation} \tag{ 28 }$$

for steady winds. The estimate depends on the age of the steady wind, t_w, and is valid for c_g ≫ 1. In the limit of an impulsively driven wind, t_w → 0, Matzner & McKee (2000) provide an estimate of

$$\begin{equation} c_{g,0}=\sqrt{\frac{9-3k_\rho }{8-3k_\rho }}. \end{equation} \tag{ 29 }$$

To compare our simulations with the analytic model, we choose c_g by interpolating between the two limits as

$$\begin{equation} c_g=\frac{\pi (9-3k_\rho)(4-k_\rho)}{8-3k_\rho }\left(\frac{t_w}{t_{\rm ff}}\right)+\sqrt{\frac{9-3k_\rho }{8-3k_\rho }}. \end{equation} \tag{ 30 }$$

For k_ρ = 3/2 this expression becomes

$$\begin{equation} c_g=2.52\frac{t_w}{t_{\rm ff}}+1.13. \end{equation} \tag{ 31 }$$

The mass-weighted average wind speed that characterizes the wind momentum injection into the core is

$$\begin{equation} \bar{v}_w = \sum _{\rm stars} \frac{1}{f_w M_i} \int \dot{M}_{w,i} v_{w,i} dt. \end{equation} \tag{ 32 }$$

Values of $\bar{v}_w$ and v_esc are given in Table 3. Because our numerical wind injection approach is based on volume-averaged quantities inside of an eight-cell radius wind injection sphere as described in Section 2.4, the effective numerical flattening parameter is θ_{0, eff} = 5.75 × 10⁻⁴ for winds with an opening angle >32° and we shall use this as the flattening parameter for the purposes of comparing the simulations to the analytic model. We note that X depends logarithmically on the flattening parameter (Equation (26)) and therefore the model prediction is not very sensitive to the estimate of the effective numerical flattening angle.

Table 3. Outflow Ejection

Σ (g cm−2)	1.0	2.0	10.0
tend (tff)	0.6	0.8	0.8
vesc (km s−1)	4.27	5.08	7.60
$\bar{v}_w|_{t=t_{\rm end}}\, ({\rm km\, s}^{-1})$	87.7	72.2	71.0
core	0.70	0.73	...
wind, simulation	1.06	0.342	0.0563
wind	0.42	0.370	...

Notes. Simulation results (rows 1–4 and 6) and analytic model predictions (rows 5 and 7). The columns indicate the cases of Σ = 1.0 g cm⁻², Σ = 2.0 g cm⁻², and Σ = 10.0 g cm⁻² from left to right. As discussed in the text, the Σ = 10.0 g cm⁻² simulation was not evolved sufficiently far in time to compare with the analytic model.

Download table as:  ASCIITypeset image

Due to constraints on computational time, we have not run numerical simulations sufficiently long to determine the final star formation efficiency. To facilitate comparison between the numerical and analytic model, we focus our attention to the ratio of the mass ejected by winds to the total stellar mass,

$$\begin{equation} \epsilon _{\rm wind} = \frac{M_{\rm ej}}{\sum _{\rm stars} M_i} = \frac{1-\epsilon _{\rm core}}{\epsilon _{\rm core}}, \end{equation} \tag{ 33 }$$

where M_ej is the total mass ejected from the system. This quantity can be computed as a function of time throughout the simulation. For the purpose of comparing our numerical simulations to the analytic prediction, we heuristically define the ejected mass as any mass that has been either ejected from the simulation domain or that is propagating with a sufficient radial component of velocity away from the center of mass of the system to overcome its gravitational binding to the system. The left panel of Figure 10 shows a plot of the total wind mass ejected in each simulation as a function of the total mass in stars, and the right panel shows the simulation result for the wind ejection efficiency _wind = M_ej/∑_starsM_i as a function of time for each simulation. We note that wind-launched gas in the highest surface density simulation with Σ = 10.0 g cm⁻² emerges from the initial core relatively late in the simulation at t ∼ 0.6t_ff. At the end of the simulation, much of the wind-launched gas is still entrained in the limbs of the outflow cavity. This is a transient effect that is not included in the analytic model and we therefore will focus our attention in comparing the analytic prediction to only the Σ = 1.0 g cm⁻² and Σ = 2.0 g cm⁻² models. By inspection of Figure 10, we adopt the values of t_w = 0.5t_ff and t_w = 0.3t_ff as characteristic of the age of the winds in the Σ = 1.0 g cm⁻² and Σ = 2.0 g cm⁻² simulations, respectively.

Zoom In Zoom Out Reset image size

The analytic predictions given by Equations (26) and (33) for the outflow mass-weighted wind speed $\bar{v}_w$ for each simulation are given in Table 3. The wind ejection efficiency at the end of each simulation (t = t_end) are listed in the table as _{wind, simulation} for comparison to the model predictions.

In drawing conclusions from comparing the analytic model result for the wind ejection ratio _wind to the simulation result _{wind, simulation} it is important to bear in mind that the analytic model predicts the fraction of the initial core that is ejected by a wind over the entire course of evolution of the initial core, whereas the numerical model is only averaged over the duration of the core evolution up to the end of the simulation. Furthermore, the analytic model cannot account for the time dependence associated with the propagation of protostellar wind shocks through the highly inhomogeneous ambient core.

For the lowest surface density case Σ = 1.0 g cm⁻², the simulation has ejected more mass with the wind per unit of mass accretion than predicted by the analytic model by a factor of 2.3. This is likely due to efficient entrainment of core gas into the wind via the interaction of the protostellar winds as they propagate through the network of dense filaments in the ambient core at early time. By t = 0.4t_ff, a wide wind cavity has been cleared through the initial core in the simulation (see the top row of Figure 3). The time-integrated outflow ejection behavior in the simulation is biased toward this early-time behavior because the simulations are only advanced to t = 0.6t_ff in the Σ = 1.0 g cm⁻² case and t = 0.8t_ff in the other models. In Figure 10, we note that evolution of the system past 0.5t_ff carries forward with less efficient entrainment of core gas into the outflows as the outflowing gas at later time propagates unimpeded through the wind channel. We expect that with continued evolution the system would asymptote to a steady state value of _core that is closer to the model prediction. However, constraints of computational cost have prevented us from testing this expectation.

The intermediate surface density case Σ = 2.0 g cm⁻² exhibits low ejection efficiency due to during the initial core-crossing of the outflow from t = 0 to t = 0.3t_ff. The ejection efficiency plateaus at later time with _core varying between 0.34 and 0.42 at later time. We note that by the end of the simulation at t = 0.8t_ff, the outflow ejection efficiency in the simulation is 92% of the analytic prediction.

We note that the star formation rate per free-fall time in our simulations, _ff = _coret_ff/t_end, is much greater than the observationally estimated value of a few percent found by Krumholz & Tan (2007) and Evans et al. (2009). This apparent discrepancy is quite easy to understand, both observationally and theoretically.

On the observational side, the Krumholz & Tan (2007) and Evans et al. (2009) estimates were for gas clumps at densities of at most ∼10⁵ cm⁻³, and the typical objects at this density are ∼10⁴  M_☉ pc-sized clumps that are forming entire star clusters. In comparison, the prestellar cores that we have simulated are much smaller and denser: n ∼ 10⁶ cm⁻³, r ∼ 0.1 pc, M ∼ 10² M_☉. In the terminology of McKee & Ostriker (2007), they are "cores" rather than "clumps." There are no observational measurements for the value of _ff in such structures. Indeed, Krumholz & Tan (2007) commented that _ff must reach values ∼1 rather than ∼0.01 at some density higher than what then current observations probed. Furthermore, it is clear that 100 M_☉ cores forming massive stars must be rare exceptions, since massive stars are rare. Since measurements of _ff only provide statistical averages, it is possible that a few n ∼ 10⁶ cm⁻³ cores like the ones we have studied undergo rapid collapse, but that there are more numerous structures at similar density that are not undergoing rapid monolithic collapse, so that the average value of _ff is much lower than in the core we have simulated.

On the theoretical side, values of _ff ∼ 0.01 are expected only in regions where there is turbulence at roughly virial levels (Krumholz et al. 2005). In our simulations, while we start out with such turbulence, this decays rapidly. Since these simulations contain a single massive star with a dominant outflow, there is nothing to drive turbulence in the core, and this allows _ff to rise rapidly as the turbulence becomes sub-virial. One can see this effect in Figure 6: the total mass in stars rises very slowly at first, and accelerates as time goes on and the turbulence decays.