# The effects of magnetic fields and protostellar feedback on low-mass cluster formation

The effects of magnetic fields and protostellar feedback on low-mass cluster formation Abstract We present a large suite of simulations of the formation of low-mass star clusters. Our simulations include an extensive set of physical processes – magnetohydrodynamics, radiative transfer, and protostellar outflows – and span a wide range of virial parameters and magnetic field strengths. Comparing the outcomes of our simulations to observations, we find that simulations remaining close to virial balance throughout their history produce star formation efficiencies and initial mass function (IMF) peaks that are stable in time and in reasonable agreement with observations. Our results indicate that small-scale dissipation effects near the protostellar surface provide a feedback loop for stabilizing the star formation efficiency. This is true regardless of whether the balance is maintained by input of energy from large-scale forcing or by strong magnetic fields that inhibit collapse. In contrast, simulations that leave virial balance and undergo runaway collapse form stars too efficiently and produce an IMF that becomes increasingly top heavy with time. In all cases, we find that the competition between magnetic flux advection towards the protostar and outward advection due to magnetic interchange instabilities, and the competition between turbulent amplification and reconnection close to newly formed protostars renders the local magnetic field structure insensitive to the strength of the large-scale field, ensuring that radiation is always more important than magnetic support in setting the fragmentation scale and thus the IMF peak mass. The statistics of multiple stellar systems are similarly insensitive to variations in the initial conditions and generally agree with observations within the range of statistical uncertainty. stars: low-mass, stars: protostars, ISM: jets and outflows, ISM: magnetic fields 1 INTRODUCTION The formation of stars through the collapse of pre-stellar molecular clouds is a rich process, mediated by the coupled effects of gravity, magnetic fields, turbulence, and feedback from protostellar evolution in the form of outflows and radiation. To date, numerical models of star cluster formation have mostly probed the coupling of at most two of these effects. While most contemporary protostellar cluster models include turbulence and self-gravity, the remaining effects have received a more piecemeal treatment. Models that include magnetic effects and protostellar outflow feedback have typically ignored radiative transfer (Li & Nakamura 2006; Wang et al. 2010; Federrath 2015), while simulations that include radiative transfer with protostellar outflows have neglected magnetic fields (Hansen et al. 2012; Krumholz, Klein & McKee 2012). Yet other works have included magnetic fields and radiative transfer but not protostellar outflows (Price & Bate 2008, 2009; Peters et al. 2011). Two exceptions that includes all three effects are the work of Myers et al. (2014) and Li, Klein & McKee (2018), who consider the formation of star clusters consistent with observations of regions of high-mass star formation with mean column densities Σ = 1 g cm−3. In this work, we perform an analogous study of regions of low-mass star formation, characterized by a mean column density Σ = 0.1 g cm−3. Under these conditions the coupling of radiation from young stars to the infalling gas is weaker (Krumholz & McKee 2008; Krumholz et al. 2010; Cunningham et al. 2011), and it is possible that magnetic and protostellar outflow effects are more prominent. A consistent and comprehensive theory of star formation must simultaneously confront a number of observational constraints, as summarized in the review by Krumholz (2014). The relative strengths of the magnetic field and gravity in a clump of molecular gas are measured by the mass-to-flux ratio. Statistical surveys of Zeeman splitting measurements indicate that this ratio is 2–3 times larger than the critical value at which the magnetic field is strong enough to prevent gravitational collapse (Falgarone et al. 2008; Crutcher et al. 2010). Recent modelling that allows for the difference between the theoretical mass-to-flux ratio, which is area weighted, and the observed one, which is density weighted, show that the data summarized by Crutcher et al. (2010) are consistent with a value for this ratio that is three times the critical value (Li, McKee & Klein 2015). Next, the collapse of a molecular cloud must give rise to a distribution of stellar masses consistent with the observed stellar initial mass function (IMF). The observed IMF is well represented by a power law for stars more massive than 1 M⊙ (Salpeter 1955) and a lognormal distribution peaked at 0.2 M⊙ for lower mass sources (Chabrier 2005). Simulations suggest that magnetic fields influence the process by which the IMF is established by adding support to otherwise super-Jeans scale masses, thereby suppressing fragmentation and influencing the distribution of stellar masses (Padoan et al. 2007; Commerçon, Hennebelle & Henning 2011; Hennebelle et al. 2011; Federrath & Klessen 2012; Myers et al. 2013). Radiative heating from protostellar sources similarly acts to increase the local Jeans mass that characterizes the gas in the protostars’ vicinity, further reducing fragmentation and inhibiting the formation of lower mass protostars (Bate 2009; Offner et al. 2009; Myers et al. 2011; Krumholz, Klein & McKee 2011, 2012; Krumholz et al. 2016). Outflows also indirectly play a role in mediating the IMF by reducing the rate of accretion on to stars, which in turn reduces the rate of radiative heating. In spite of the complexity of the coupling among magnetic, radiative and outflow feedback, the shape of the IMF is remarkably insensitive to the star-forming environment, particularly for the majority of stars that are heavier than brown dwarfs but not so massive as to be exceptionally rare, where observational statistics are best constrained (Offner et al. 2014). Not only must the complex physical processes associated with star formation collude to fragment molecular clouds into a distribution of masses that is consistent with the observed IMF, but the star formation rate must also agree with observation. Star formation is remarkably inefficient. Absent turbulent or magnetic support, a molecular cloud should collapse entirely into stars on a free-fall time-scale. In fact, the star formation efficiency per free-fall time εff, defined as the fraction of gas mass that is converted to stars per free-fall time, was found to lie between 0.001 and 0.09 by Krumholz & Tan (2007); a much larger and more recent compilation of data gives εff = 0.015 with a scatter of ∼0.5 dex for both Galactic and extragalactic regions (Krumholz, Dekel & McKee 2012), and subsequent authors have obtained similarly low values (e.g. García-Burillo et al. 2012; Evans, Heiderman & Vutisalchavakul 2014; Salim, Federrath & Kewley 2015; Usero et al. 2015; Heyer et al. 2016).1 Turbulence provides a significant source of support against rapid collapse. Gas velocity dispersions inferred from observed molecular cloud line width–size relations (Larson 1981) are consistent with the inertial range scaling of both supersonic hydrodynamic and magnetohydrodynamical turbulence (Gammie & Ostriker 1996; Ostriker, Gammie & Stone 1999) and the magnitude of the velocity dispersions is consistent with a near virial equipartition between turbulent support and gravitational contraction (McKee, Li & Klein 2010). Numerical studies of magnetized turbulence show that turbulence at virial levels, in conjunction with protostellar outflows, can produce values of εff close to the observed, low values (Li & Nakamura 2006; Wang et al. 2010; Padoan, Haugbølle & Nordlund 2012; Krumholz, Klein & McKee 2012; Myers et al. 2014; Federrath 2015). Magnetohydrodynamical turbulence, however, undergoes exponential decay via dissipation in shocks on a turbulent crossing time (Stone, Ostriker & Gammie 1998; Mac Low 1999), which is comparable to the free-fall time in a virialized cloud. Gravitational collapse can drive turbulence, but only at the cost of increasing the density and thus lowering the free-fall time (Robertson & Goldreich 2012; Murray et al. 2015). As a result, in the absence of an external or internal energy source to drive the turbulence, hydrodynamical models fall into a state of rapid, near free-fall collapse (Hansen et al. 2012; Krumholz et al. 2012). A number of forcing mechanisms applicable to low-mass star clusters have been proposed, including accretion from larger scales (e.g. Klessen & Hennebelle 2010; Goldbaum et al. 2011; Lee & Hennebelle 2016a,b), internal forcing by protostellar outflows (e.g. Li & Nakamura 2006; Matzner 2007; Wang et al. 2010; Cunningham et al. 2011; Peters et al. 2014; Federrath 2015), and forcing by external H ii regions (e.g. Matzner 2002; Krumholz, Matzner & McKee 2006; Goldbaum et al. 2011) and supernova explosions (e.g. Padoan et al. 2016, 2017; Pan et al. 2016). In this paper, we present a series of simulations of the formation of low-mass stellar systems, which we compare to the various observational constraints on the star formation process we have outlined above. Our goal is to search for models that are capable of simultaneously reproducing the observed IMF and star formation efficiency, determine the relationship between these two constraints, and identify the physical mechanisms that are responsible for satisfying them. In order to accomplish this, we carry out a suite of simulations All our simulations use radiative transfer and protostellar radiation feedback, since these are likely key to setting the IMF. In order to explore the influence of outflows, we toggle these on and off. We also vary the initial magnetization of our regions, and consider both models with turbulence that is continuously driven to maintain virial equilibrium and those where no artificial driving is included. We describe our numerical method and simulation setup in Section 2, we describe the results in Section 3, we explore the question of magnetic fields in more detail in Section 4, and we summarize our findings in Section 5. 2 SIMULATION SETUP 2.1 Initial conditions We perform all our simulations using the orion2 code . Our setup shares several initial parameters with Offner et al. (2009) and Hansen et al. (2012). The initial state consists of a solar metallicity molecular gas with a mean molecular weight 2.33mp and a temperature Tg = 10 K, corresponding to an initial sound speed cs = 0.19 km s−1. The length of the cubical domain is L = 0.65 pc and the average density is $$\bar{\rho }=4.46\times 10^{-20}\ {\rm g\ cm}^{-3}$$ ($$n_{\rm H_2} = 1.15\times 10^4$$ cm−3) corresponding to a total mass of M = 185 M⊙. The initial mean state corresponds to a minimum stable Jeans length λJ = 0.20 pc, and a Jeans mass MJ = 2.7 M⊙. Our models are carried forward with periodic boundary conditions for the gravity and magnetohydrodynamics (MHD). We use Marshak boundary conditions with an external radiation field corresponding to a 10 K isotropic blackbody for the radiation field. The ability of a magnetized cloud to resist gravitational collapse is measured by the magnetic critical mass (Mouschovias & Spitzer 1976),   $$M_\Phi = c_\Phi \frac{\Phi }{\sqrt{G}},$$ (1)where Φ is the total magnetic flux threading the cloud and $$c_\Phi$$ is a dimensionless constant that depends weakly on the equilibrium shape of the cloud in the absence of gravity. A cloud of mass less than $$M_\Phi$$ cannot collapse, whereas a cloud with a greater mass cannot be supported against collapse by magnetic fields alone. Because the value of $$c_\Phi$$ is only weakly sensitive to the cloud geometry (McKee et al. 2010), we adopt the value appropriate for a slab threaded by a uniform, perpendicular field, $$c_\Phi =1/2\pi$$, consistent with the convention chosen in several contemporary theoretical (Li et al. 2015) and observational works (Crutcher et al. 2010). We then define the normalized mass-to-flux ratio   $$\mu _\Phi = \frac{M}{M_\Phi } = 2\pi \sqrt{G}\left(\frac{M}{\Phi }\right)=2\pi \sqrt{G}\left(\frac{\Sigma }{\bar{B}}\right),$$ (2)where $$\bar{B}\equiv \Phi /S$$ is the magnitude of the mean magnetic field, Σ is the mass surface density, and S the cloud's cross-sectional area in the plane orthogonal to the magnetic field direction. Since mass and magnetic flux are conserved for the entire simulation box, the value of $$\mu _\Phi$$ for the box is constant and can be used to specify the initial mean magnetic field. We define initial fields large enough to render the system close to magnetically critical, $$\mu _\Phi \approx 1$$, ‘strong,’ initial fields that make the system moderately supercritical, $$\mu _\Phi \approx 2\text{--}3$$ ‘moderate,’ and initial fields corresponding to $$\mu _\Phi \gg 1$$ ‘weak.’ It is convenient to express the strength of the magnetic field relative to the kinetic energy in terms of the Alfvèn Mach number,   $${\cal M}_{\rm A}=\frac{v_{\rm rms}}{v_{\rm A}},$$ (3)where $$v_{\rm A}^2\equiv B_{\rm rms}^2/(4\pi \bar{\rho })$$ is the Alfvèn velocity, Brms is the root-mean-square (rms) volume-weighted field strength, and vrms is the rms mass-weighted, 3D velocity dispersion. Relative to thermal pressure, the strength of the magnetic field can be expressed in terms of the plasma-β parameter   $$\beta = \frac{8\pi \bar{\rho }c^2_s}{B_{\rm rms}^2} = 2\left(\frac{{\cal M}_{\rm A}}{{\cal M}}\right)^2.$$ (4) In this work, we carry out simulations with several different mean magnetic field strengths. Our models begin with an initially uniform magnetic field B0 in the $$\hat{z}$$-direction. We consider cases with mass-to-flux ratios $$\mu _\Phi =1.56\textrm{,}2.17\textrm{,}23.1\textrm{,and} \infty$$. The corresponding initial Alfvèn Mach numbers are $${\cal M}_{\rm A,0}=1.0\textrm{,}1.4\textrm{,}15\textrm{,and} \infty$$, and the initial plasma β's are β0 = 0.046, 0.089, 10, and ∞. We note that the recent models by Federrath (2015) consider a somewhat more weakly magnetized model ($$\mu _\Phi =3.8$$) than our two most strongly magnetized cases, and the model by Wang et al. (2010) has an initial cloud with $$\mu _\Phi =1.4$$, comparable to our strongest field model. Turbulent conditions for starting the cluster simulation are obtained by applying the driving recipe of Mac Low (1999). The procedure consists of applying a force consisting of a purely solenoidal, uniform spectrum of modes spanning wave numbers $$\boldsymbol {k}$$ such that $$1 \le |\boldsymbol {k}|L/2\pi \le 2$$ to a uniform initial gas. This turbulent driving force is intended to mock up the effect of sources of mechanical energy on scales larger than the simulation box, via one of the numerous possible mechanisms discusses in Introduction. The observed velocity dispersion to size relation can be met by models of compressive turbulence that span a range of relative contributions from solenoidal and compressive forcing (Federrath, Klessen & Schmidt 2009). Recent simulations of supernova-driven turbulence (Padoan et al. 2016; Pan et al. 2016) show that even when the initial driving is purely compressive, the turbulence becomes primarily solenoidal in the inhomogeneous interstellar medium (ISM); indeed, these authors conclude that external supernova-driven turbulence is best characterized as 84 per cent solenoidal after mediation through the ISM to adjacent star-forming regions. The work of Federrath et al. (2010) indicates that the statistical properties of turbulent flow are not sensitive to ≲ 20 per cent compressive contributions to compressive driving. This suggests that our purely solenoidal forcing is a plausible approximation for the driving of turbulence in a molecular cloud. The turbulent driving force is scaled to maintain a roughly constant mass-weighted rms thermal Mach number of $${\cal M}= v_{\rm rms}/c_{\rm s} \equiv \sqrt{3}\sigma _{\rm v}/c_{\rm s}= 6.6$$, where σv is the 1D velocity dispersion. The virial parameter,   $$\alpha _{\rm vir}\equiv \frac{5\sigma _{\rm v}^2 R}{GM},$$ (5)measures the balance between kinetic and gravitational energies; a spherical cloud of uniform density has αvir = 1 in virial equilibrium. With R = L/2 and our chosen value of $${\cal M}$$, the computational box has αvir = 1.05. The gas is initially evolved under the action of the driving force without gravity according to the equations of ideal, isothermal MHD for two crossing times,   $$t_{\rm cross}=\frac{L}{v_{\rm rms}}=0.51\ \textrm{Myr},$$ (6)on a uniform mesh with 5123 zones with periodic boundary conditions. This resolution has been shown by Li et al. (2012) to be sufficient to capture a well-resolved inertial cascade down to scales of  L/30 = 0.02 pc, an order of magnitude smaller than the initial Jeans scale of the turbulent cloud. Consequently, the turbulent flows on the scales that give rise to the initial fragmentation of the simulated cluster are well resolved. Henceforth, we denote the state of the turbulent gas after evolving for two crossing times as the time t = 0.2 We take these conditions as the initial state of the turbulent cluster simulation, which includes a more comprehensive set of physics that we describe in Section 2.3. 2.2 Defining the free-fall time and star formation rate Because our simulations take place in a periodic box that represents only a portion of a cloud, we must give some care to how we go about defining the free-fall time and the efficiency of star formation. Absent magnetic or thermal support to oppose self-gravity, the free-fall time of gas structures of characteristic density ρ is   $$t_{\rm ff} = \sqrt{\frac{3\pi }{32G \rho }}.$$ (7)The coefficient $$\sqrt{3\pi /32}$$ applies to spherical geometries, but by convention we retain it for all geometries. The mean cloud density $$\bar{\rho }$$ corresponds to a free-fall time tff, IC = 0.315 Myr, where we use the subscript IC to refer to quantities defined for the uniform-density state before driving the initial turbulence to steady state. We can measure star formation rates relative to this free-fall time via the usual dimensionless star formation rate parameter   $$\epsilon _{\rm ff,IC} = \frac{\dot{M}_* t_{\rm ff,IC}}{M}.$$ (8) While εff, IC provides one measure of the star formation rate, it can be somewhat misleading. As we drive the turbulence to statistical steady state, the gas density distribution becomes non-uniform, and a majority of the mass ends up residing at densities much higher than $$\bar{\rho }$$. This leads to two biases. First, because most of the mass is at higher density, we might expect εff, IC to end up relatively large even if star formation in the dense gas is inefficient, simply because the dense gas has a much shorter free-fall time. Second, because we are using a periodic box to handle gravity, much of the low-density gas is not self-gravitating at all, and thus would not collapse even with efficient star formation. Including this mass in a calculation of εff would lead to an artificially low value. To counteract these biases, we define a new reference density that provides a better description of the gas after the turbulence is fully developed. Let ρmed be the mass-weighted median density – i.e. the density at which half the mass is denser and half less dense – at the end of the turbulent driving phase (t = 0). Then, we define ρ>med as the average density of the regions with ρ > ρmed at that time. In effect, ρ>med gives an estimate of the mean density of the gas that is actually self-gravitating prior to any star formation. We denote the associated free-fall time as tff, >med. The dimensionless star formation rate is defined as   $$\epsilon _{\rm ff,>med} = \frac{\dot{M}_* t_{\rm ff,>med}}{M_{\rm >med}}.$$ (9)Since M = 2M>med, this can also be expressed as   $$\epsilon _{\rm ff,>med} = \frac{2 \dot{M}_* t_{\rm ff,>med}}{M}.$$ (10) While εff, >med is a natural physical measure of star formation efficiency, it cannot readily be compared to observations. This is both because our simulation at t = 0 represents a somewhat artificial state of driven turbulence but no gravitational effects, and because observations cannot usually measure gas densities directly in any event. Instead, volume densities must be inferred indirectly from molecular line critical densities, or, more commonly for Galactic observations, using column densities that can be measured almost directly. Following Krumholz et al. (2012), we define the observationally inferred characteristic density by taking projected column densities and extracting those regions denser than 0.048 g cm−2 = 230  M⊙ pc−2, corresponding to those regions with K-band extinction AK > 0.8 mag, a threshold used by a number of previous authors (e.g. Lada, Lombardi & Alves 2010; Lada et al. 2012, 2013). We compute the projected area of these regions S>0.8, and from this determine the projected mass enclose M>0.8 and the effective radius $$R_{\rm >0.8} = \sqrt{S_{\rm >0.8}/\pi }$$. We then define the characteristic observed density of the projection as $$\rho _{\rm obs, proj} = 3 M_{\rm >0.8} / (4 \pi R_{\rm >0.8}^3)$$, and take the mean value from projections in the three cardinal directions to be the characteristic density ρ>0.8 for each simulation. We compute this quantity at the final time in each simulation (in contrast to ρ>med, which is calculated at t = 0, when gravity is turned on), and we denote the corresponding free-fall time as tff, >0.8. The star formation rate per tff, >0.8 is then   $$\epsilon _{\rm ff,>0.8} = \frac{\dot{M}_* t_{\rm ff,>0.8}}{M_{\rm >0.8}},$$ (11)and this is the quantity that is most natural to compare to observed star formation rates. In Table 1, we summarize the two characteristic densities (ρ>med, and ρ>0.8) that we have defined for all our simulations. The former is typically 7–9 times greater than the mean density $$\bar{\rho }=4.46 \times 10^{-20}\, {\rm g\ cm}^{-3}$$, due to turbulent compression. The observationally defined ρ0.8 ranges from $$2\ \textrm{to}\ 7 \times \bar{\rho }$$ due to the influence of the combined effects of turbulent compression, gravitational collapse, and stellar feedback. We also report the mean values of the Alfvèn Mach number and plasma β at the end of the driving phase and at the final time in the simulation, to make it clear how these evolve with time as well. This evolution is a result of field amplification by turbulence and gravitational collapse. For consistency, we use the same notation for these other quantities as for the density and free-fall time, i.e. quantities measured in the initial, uniform state have a subscript ‘IC’ (for initial condition), quantities averaged over gas denser than the median density at the end of the driving phase have subscript ‘>med’, and those defined based on averages within the 0.8 mag contour are defined by the subscript ‘>0.8’. Table 1. Model parameters and initial conditions. The outflows and driving columns indicate whether we included protostellar outflows, and whether we continued driving the turbulence after gravity was turned on. For the remaining parameters, the subscript ‘IC’ denotes the state in the uniform initial condition, the subscript ‘>med’ denotes quantities averaged over gas denser than the median density at the instant when self-gravity is switched on, and the subscript ‘>0.8’ denotes averages over gas with a K-band extinction AK > 0.8 mag at the termination of the model. See the text for the details of how the averaging regions are defined. The mean density in all cases is 4.46 × 10−20 g cm−3.   Physics included  IC  >med  >0.8  μΦ  Outflows  Driving  βIC  $$\mathcal {M}_{A,\rm IC}$$  ρ>med (g cm−3)  β>med  $$\mathcal {M}_{\rm >med}$$  ρ>0.8 (g cm−3)  β>0.8  $${\cal M}_{\rm A,>0.8}$$  1.56  Y  N  0.046  1.0  2.36 × 10−19  0.041  0.95  3.21 × 10−19  0.039  1.09  1.56  N  N  0.046  1.0  2.36 × 10−19  0.041  0.95  1.45 × 10−19  0.038  1.14  1.56  Y  Y  0.046  1.0  2.36 × 10−19  0.041  0.95  9.76 × 10−20  0.037  1.05  2.17  Y  N  0.089  1.39  2.18 × 10−19  0.057  1.14  1.53 × 10−19  0.074  1.92  2.17  Y  Y  0.089  1.39  2.18 × 10−19  0.057  1.14  2.35 × 10−19  0.067  1.70  23.1  Y  N  10.0  14.8  1.90 × 10−19  0.46  3.14  1.62 × 10−19  0.41  6.67  23.1  Y  Y  10.0  14.8  1.90 × 10−19  0.46  3.14  1.07 × 10−19  0.75  4.45  ∞  Y  N  ∞  ∞  2.66 × 10−19  ∞  ∞  1.33 × 10−19  ∞  ∞  ∞  Y  Y  ∞  ∞  2.66 × 10−19  ∞  ∞  9.24 × 10−20  ∞  ∞    Physics included  IC  >med  >0.8  μΦ  Outflows  Driving  βIC  $$\mathcal {M}_{A,\rm IC}$$  ρ>med (g cm−3)  β>med  $$\mathcal {M}_{\rm >med}$$  ρ>0.8 (g cm−3)  β>0.8  $${\cal M}_{\rm A,>0.8}$$  1.56  Y  N  0.046  1.0  2.36 × 10−19  0.041  0.95  3.21 × 10−19  0.039  1.09  1.56  N  N  0.046  1.0  2.36 × 10−19  0.041  0.95  1.45 × 10−19  0.038  1.14  1.56  Y  Y  0.046  1.0  2.36 × 10−19  0.041  0.95  9.76 × 10−20  0.037  1.05  2.17  Y  N  0.089  1.39  2.18 × 10−19  0.057  1.14  1.53 × 10−19  0.074  1.92  2.17  Y  Y  0.089  1.39  2.18 × 10−19  0.057  1.14  2.35 × 10−19  0.067  1.70  23.1  Y  N  10.0  14.8  1.90 × 10−19  0.46  3.14  1.62 × 10−19  0.41  6.67  23.1  Y  Y  10.0  14.8  1.90 × 10−19  0.46  3.14  1.07 × 10−19  0.75  4.45  ∞  Y  N  ∞  ∞  2.66 × 10−19  ∞  ∞  1.33 × 10−19  ∞  ∞  ∞  Y  Y  ∞  ∞  2.66 × 10−19  ∞  ∞  9.24 × 10−20  ∞  ∞  View Large 2.3 Physics included and numerical method Our protostellar cluster simulations are performed using the orion2 adaptive mesh refinement (AMR) code (Li et al. 2012) beginning from the initial state generated following the procedure described in the previous section. The code solves the equations of ideal MHD using the scheme of Mignone et al. (2012), along with coupled self-gravity (Truelove et al. 1998; Klein et al. 1999) and radiation transfer (Krumholz et al. 2007) in the two-temperature, mixed-frame, grey, flux-limited diffusion approximation. The exact set of equations solved is the same as those given in Myers et al. (2014). orion2 contains many options for the numerical advection scheme, but we have found the best stability with the Harten-Lax-Van Leer-Discontinuities Riemann solver (Miyoshi & Kusano 2005) coupled to the constrained transport magnetic flux advection of Londrillo & del Zanna (2004), particularly for magnetically dominated flows. Our radiative transfer calculations use the frequency-integrated grey dust opacities from the iron-normal, composite aggregates model of Semenov et al. (2003). For each magnetic field strength parameter considered in this study, $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞, we present two models – one with decaying turbulence where the large-scale driving forces described in the previous section are turned off, and one where the turbulent forcing continues with a constant rate of energy injection that balances the rate of turbulent decay as estimated for a given mean magnetic field strength and turbulent Mach number by Mac Low (1999). In the case without turbulent forcing, self-gravity is switched on and turbulent forcing is switched off simultaneously. Our cluster evolution models neglect non-ideal MHD effects. These effects have been shown to impact turbulent density profiles on small scales Li, Myers & McKee (2012) and therefore could influence the fragmentation dynamics of strongly magnetized cores. Recent studies on the impact of non-ideal MHD effects on the evolution and fragmentation of self-gravitating cores indicate that the evolution of supercritical cores is most strongly influenced by their initial conditions and that the impact of non-ideal effects on their subsequent evolution and fragmentation are small by comparison Wurster, Price & Bate (2017). Consequently, our models should reliably capture the gross impacts of varying large-scale magnetic field strength on the cloud fragmentation and the statistics of its collapse. There are two limitations of our radiative transfer method to which we should call attention. First, we assume that the gas and dust temperatures are the same. In the dense regions simulated by Myers et al. (2014), this is a very good assumption, because the dust and gas temperatures become nearly identical at densities above ∼104–105 cm−3 (e.g. Goldsmith 2001). The assumption of strong coupling is still valid for most of the gas in our simulation domain that is actively participating in star formation (cf. ρ>, med in Table 1), but it is not valid for the lower density, non-self-gravitating regions. In these parts of the flow, we are somewhat overestimating the dust cooling rate for the gas. The second approximation we make is to assume that the absorption opacity is the same as the emission opacity, which is equivalent to assuming that the dust and radiation temperatures are the same. This approximation is clearly invalid very close to the star, where dust is exposed to direct stellar radiation and reprocesses the radiation into the infrared. However, the dust opacity to direct stellar radiation is so high that all this reprocessing occurs in a thin zone at the dust destruction front, which for low-mass star formation is confined to ≲1 au from the star (except perhaps over the small range in solid angle evacuated by outflows). Such small structures are unresolved in our simulations. Our approximation would only become problematic for stars larger than ∼20  M⊙, massive enough that radiation pressure can inflate bubbles that push the dust destruction front out to thousands of au (e.g. Rosen et al. 2016). While orion2 does include a hybrid radiative transfer method to handle this situation (Rosen et al. 2017), we do not use it here because our simulations do not form any stars massive enough to require it. Once the radiation field is dominated by reprocessed radiation, the dust and radiation temperatures will remain close so long as the dust is opaque. Outside the dust photosphere, the emission opacity declines while the absorption opacity remains constant, but this affects only the low-frequency part of the emission spectrum, which has only a small fraction of the energy (Chakrabarti & McKee 2005).3 Our models include several source terms to capture the effects of protostellar feedback that originate on smaller length scales than those resolved in the model. We initialize a sink particle in any zone on the finest AMR level, and only on the finest level, that becomes dense enough to reach a local Jeans number   $$J=\sqrt{\frac{G \rho \Delta x^2}{\pi c_{\rm s}^2}} > \frac{1}{4}.$$ (12)Sink particles then evolve and interact with the cluster through gravity according to the methodology of Krumholz, McKee & Klein (2004), updated to include the effects of magnetic fields on the rate of gas accretion on to sink particles (see the appendix of Lee et al. 2014). The key assumption in our treatment is that the point mass accretes mass, but not flux. This assumption is motivated by observations that show that the magnetic flux in young stellar objects is orders of magnitude less than that in the gas that formed these objects, implying that flux accretion is very inefficient (McKee & Ostriker 2007). The sink particles couple to their surrounding gas gravitationally and through accretion and feedback source terms that operate within a radius of 4Δx on the finest AMR mesh level around each sink particle. The sink particles emit radiation according to the protostellar evolution model of Offner et al. (2009). Our models include the atomic line-cooling sources described in Cunningham et al. (2011) to treat strongly heated regions behind outflow shocks that are sufficiently fast ∼30 km s−1 to dissociate molecules. Feedback due to protostellar outflows is also included via momentum sources around each sink particle following the procedure described in Cunningham et al. (2011): the outflow mass ejection rate is set by the fraction of accreted gas that is ejected into the outflow fw and the outflow ejection speed vw. In this work, we make the same wind model parameter choices as Hansen et al. (2012) by setting fw = 0.3 and vKep = min (vKep,  60 km s−1), where vKep is the Kepler speed at the surface of the protostar. We note that these parameters differ from the parameters used in Cunningham et al. (2011) and Li et al. (2018) that set vw = vKep/3. Consequently, the models presented here impart relatively more outflow feedback during the earlier phases of pre-stellar evolution relative to Cunningham et al. (2011) and Li et al. (2018). Computational expediency motivates the cap on the wind velocity so that isolated massive protostars do not dominate the time-step of the cluster model. The AMR hierarchy is initialized on a 2563 base grid that we denote as level $$\mathcal {L}=0$$. Our initial state is taken from turbulent initial conditions that were established on a 5123 grid with no adaptivity. After t = 0, regions of the flow within 2 base-grid zones of an interface with a density jump Δρ/ρ > 0.5 on the base level are refined to $$\mathcal {L}=1$$, so that the 5123 resolution of the initial state is retained in these regions. On all AMR levels, refinement is triggered in regions in which the local Jeans scale is resolved by fewer than eight zones, J > 1/8 (Truelove et al. 1997). Note that this refinement criterion triggers refinement to the finest level, $$\mathcal {L}_{\rm MAX}$$, at a density one-fourth that required to trigger insertion of a sink particle. This ensures that strongly collapsing regions are refined to the finest level. Refinement is also triggered when there is a large jump in radiation energy density, |∇Er|Δx/Er > 0.125. A resolution sensitivity test was conducted for the case with $$\mu _\Phi =2.17$$ by first performing this run with $$\mathcal {L}_{\rm MAX}=4$$ and then repeating it without the density gradient refinement triggering on level $$\mathcal {L}_0$$ and with $$\mathcal {L}_{\rm MAX}=3$$. We found that the total mass accreted on to protostars out to t = 1.35tff in the two models agreed to within 3 per cent and that the mass distribution of protostars had a similar level of agreement. Our cluster models with decaying turbulence begin with $$\mathcal {L}_{\rm MAX}=4$$. The effective resolution of all collapsing and/or heated regions is therefore $$L/(256 \times 2^{\mathcal {L}_{\rm MAX}}) = 32$$ au. Given the good agreement between the runs with $$\mathcal {L}_{\rm max}=3$$ and 4, the $$\mu _\Phi =1.56$$ and 2.17 models with decaying turbulence were switched to $$\mathcal {L}_{\rm MAX}=3$$ with 65 au effective resolution at t = 1.7tff in order to reduce the run time. The models with turbulent forcing were also performed with $$\mathcal {L}_{\rm MAX}=3$$ and 65 au effective resolution, except for the $$\mu _\Phi =23.1$$ case which retained $$\mathcal {L}_{\rm MAX}=4$$ throughout the evolution of the model. We note that the limited resolution of these models results in a gain in computational expediency but comes at some cost to the fidelity of the model in treating the fragmentation of low-mass cores. The density threshold triggering sink formation at our background 10 K temperature is given by inverting equation (12) as 1.1 × 10−15–4.5 × 10−15 g cm−3 for 65 and 32 au effective resolution. The isothermal-to-adiabatic transition where compressional heating balances cooling occurs at a substantially higher density of ∼ 10− 14(T/10 K)6 g cm−3 (see chapter 16 of Krumholz 2015), and thus our resolution would not be sufficient to capture fragmentation in a simulation without radiative transfer. However, once stars form, radiative heating powered by accretion luminosity heats material at much lower densities than one would infer simply by balancing adiabatic compression against cooling. This substantially eases the resolution required to capture fragmentation, and we show below that our resolution is sufficient that in practice we essentially always resolve the region where gas transitions from the background temperature to elevated temperatures around most protostars in the simulations. None the less, we caution that our limited resolution may still cause us to miss fragmentation in rare regions with little radiative heating. 3 RESULTS In this section, we examine the results of the simulations summarized in Table 1. In the discussion and figure legends that follow we denote the various models by their initial mass-to-flux ratio; cases that include turbulent driving are labelled with the designation ‘Driven’. We have one model with outflows turned off, to which we refer as ‘no wind’; this model also is not driven. 3.1 Protostellar mass accretion Fig. 1 shows the temporal evolution of the total mass M* assembled into protostars for each model as a function of time. The left-hand panel shows the simulations with decaying turbulence and the right-hand panel shows the simulations with large-scale turbulent driving, a convention we shall follow throughout this paper. Fig. 2 shows the star formation efficiency per initial free-fall time εff, >med (see Section 2.2 for the precise definition) as a function of time. Fig. 3 shows the number of protostars formed in each model as a function of time. Figure 1. View largeDownload slide Total mass accreted for the models without (left) and with (right) turbulent forcing. Figure 1. View largeDownload slide Total mass accreted for the models without (left) and with (right) turbulent forcing. Figure 2. View largeDownload slide Instantaneous value of the dimensionless star formation rate εff, >med for the models without (left) and with (right) turbulent forcing (right). The data have been boxcar-smoothed over a time-scale of 0.05 tff, >med to remove short-time-scale variability. Figure 2. View largeDownload slide Instantaneous value of the dimensionless star formation rate εff, >med for the models without (left) and with (right) turbulent forcing (right). The data have been boxcar-smoothed over a time-scale of 0.05 tff, >med to remove short-time-scale variability. Figure 3. View largeDownload slide Total number of protostars formed for the models without (left) and with (right) turbulent forcing. Figure 3. View largeDownload slide Total number of protostars formed for the models without (left) and with (right) turbulent forcing. We also report, in Table 2, both εff, >med and the observed value of the star formation efficiency per free-fall time, εff, >0.8 (again see Section 2.1). The quantities reported in the table are time-weighted averages over the last 0.5tff, >med of the simulation. For most models, we run the simulation until the value of εff, >med stabilizes, and thus the values reported in Table 2 are the steady-state ones. There are two exceptions to this statement. First, the two weakest field models with decaying turbulence ($$\mu _\Phi =23.1$$ and ∞) do not appear to have well-converged values of εff, and instead appear to approach near-free-fall collapse. We run these models for less time than the others, as such efficient collapse is generally inconsistent with the observational constraints. For them, the values of εff reported in Table 2 are only lower limits, and we flag them as such in the table. The second exception are the models with $$\mu _\Phi =2.17$$. We are forced to halt these runs due to computational constraints, as the number of sink particles and highly refined zones eventually becomes such that exploring these models further is prohibitively expensive. While we would, absent consideration of computation constraints, prefer to continue the $$\mu _\Phi =2.17$$ runs further, we note that these models have begun to exhibit stabilized εff.0 by the termination time. Thus, we do not report these as upper limits in Table 2, but warn readers here that they are less secure than the values reported for the other runs. Table 2. Star formation efficiency per free-fall time in each run. The various values of εff are computed for different definitions of the free-fall time, as discussed in Section 2.2. The quantities reported are averaged over the last 0.5tff, >med of each simulation. A value preceded by ≳ indicates that εff was still rising strongly and the end of the simulation, and the values listed are therefore lower limits. μΦ  Outflows  Driving  εff, IC(per cent)  εff, >med(per cent)  εff, >0.8(per cent)  $$\mu _\Phi = 1.56$$  Y  N  12  8.2  6.0  $$\mu _\Phi = 1.56$$  N  N  23  15.8  19  $$\mu _\Phi = 1.56$$  Y  Y  4.8  3.2  4.2  $$\mu _\Phi = 2.17$$  Y  N  26  18.8  19  $$\mu _\Phi = 2.17$$  Y  Y  15  11.0  8.7  $$\mu _\Phi = 23.1$$  Y  N  ≳ 30  ≳ 24  ≳ 22  $$\mu _\Phi = 23.1$$  Y  Y  8.2  6.4  6.6  $$\mu _\Phi = \infty$$  Y  N  ≳ 50  ≳ 34  ≳ 37  $$\mu _\Phi = \infty$$  Y  Y  6.0  4.0  5.2  μΦ  Outflows  Driving  εff, IC(per cent)  εff, >med(per cent)  εff, >0.8(per cent)  $$\mu _\Phi = 1.56$$  Y  N  12  8.2  6.0  $$\mu _\Phi = 1.56$$  N  N  23  15.8  19  $$\mu _\Phi = 1.56$$  Y  Y  4.8  3.2  4.2  $$\mu _\Phi = 2.17$$  Y  N  26  18.8  19  $$\mu _\Phi = 2.17$$  Y  Y  15  11.0  8.7  $$\mu _\Phi = 23.1$$  Y  N  ≳ 30  ≳ 24  ≳ 22  $$\mu _\Phi = 23.1$$  Y  Y  8.2  6.4  6.6  $$\mu _\Phi = \infty$$  Y  N  ≳ 50  ≳ 34  ≳ 37  $$\mu _\Phi = \infty$$  Y  Y  6.0  4.0  5.2  View Large For comparison, Krumholz et al. (2012) combine data from galactic clouds, nearby galaxies, and galaxies at high redshift, and find that all are consistent with a star formation efficiency εff, >0.8 = 0.015, with a scatter of ∼0.5 dex. As noted in Introduction, numerous other studies have confirmed this basic result. Comparing to the values shown in Fig. 2 and reported in Table 2, we see that the models with $$\mu _\Phi = 1.56$$ fall within the observational envelope if we include either driving or outflows, as do all the more weakly magnetized models with driven turbulence, although these tend to lie towards the top end of the observationally permitted range. We therefore find that the star formation rate can be kept low enough to be compatible with observations either if there are external mechanical energy inputs to maintain the turbulence or if there are no such influences, but the magnetic field lies close to the critical value. 3.2 Cluster morphology Fig. 4 shows the mean density along lines of sight in the $$\hat{\boldsymbol {x}}$$-direction at the termination time of each model, overlaid with a semi-transparent rendering of regions with |v| > 2 vrms. Recall that the initial magnetic field is in the $$\hat{\boldsymbol {z}}$$-direction, and is therefore upwardly oriented on the page in Fig. 4. We note that the mean density along the line of sight, which is proportional to the surface density, becomes increasingly filamentary as the strength of the magnetic field increases. In the cases where the magnetic field is the most dynamically significant, the outflows are oriented predominately perpendicularly to the orientation of the densest structures. The $$\mu _\Phi =1.56$$ case with decaying turbulence indicates a final cluster geometry that is most dominated by the magnetic field. In this case, the cluster is arranged mostly in a planar geometry normal to the $$\hat{\boldsymbol {z}}$$-direction. The magnetic fields in this case are strong enough to constrain large-scale collapse to be predominantly in the direction of the initial magnetic field, with less harassment from large-scale turbulent forcing than the equivalent model with turbulent forcing. Figure 4. View largeDownload slide Volume-weighted mean densities projected along the $$\hat{\boldsymbol {x}}$$-direction for the models with turbulent forcing (right) and with decaying turbulence (left). The column density can be recovered by multiplying the mean density shown by the computational domain length 0.65 pc; for reference, this means that a density of 10−18 g cm−3 corresponds to a column density of 2.0 g cm−2. The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ runs, from top to bottom. Regions with gas velocity greater than 2vrms are shaded in red with the opacity of the red regions increasing to a maximum for velocities exceeding 5vrms. White circles indicate the projected position of the protostellar sink particles. Figure 4. View largeDownload slide Volume-weighted mean densities projected along the $$\hat{\boldsymbol {x}}$$-direction for the models with turbulent forcing (right) and with decaying turbulence (left). The column density can be recovered by multiplying the mean density shown by the computational domain length 0.65 pc; for reference, this means that a density of 10−18 g cm−3 corresponds to a column density of 2.0 g cm−2. The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ runs, from top to bottom. Regions with gas velocity greater than 2vrms are shaded in red with the opacity of the red regions increasing to a maximum for velocities exceeding 5vrms. White circles indicate the projected position of the protostellar sink particles. 3.3 Protostellar mass functions We next examine the mass functions for the stars formed in our simulations. As a first step to this, we verify that these distributions have come close to being stable in time. To demonstrate this, in Fig. 5 we plot the temporal evolution of the 25th percentile, median and 75th percentile of the distribution of protostellar masses in each model as a function of time. Similar to the temporal evolution of εff, we see that the protostellar mass percentiles have mostly stabilized by the end of the runs, except for the cases of a weak initial field ($$\mu _\Phi =23.1$$ and ∞) with decaying turbulence. Consequently, we can interpret the protostellar mass distributions of the moderate and strong initial field models with decaying turbulence, and all of the models with driven turbulence, as being only weakly time-dependent. We do warn, however, that this stabilization is a result of an equilibrium between new stars forming at low mass, and existing stars growing. As noted in Krumholz et al. (2016), if some process were to halt the formation of new stars but the existing stars were to continue accreting their envelopes, the mass distribution would shift upward. Figure 5. View largeDownload slide The evolution of the stellar mass distribution for models with turbulent forcing (right) and with decaying turbulence (left). In these figures, the central lines show the 50th percentile mass, while the shaded regions indicate the 25th–75th percentile range. Colours correspond to different values of $$\mu _\Phi$$, as indicated in the legend. Figure 5. View largeDownload slide The evolution of the stellar mass distribution for models with turbulent forcing (right) and with decaying turbulence (left). In these figures, the central lines show the 50th percentile mass, while the shaded regions indicate the 25th–75th percentile range. Colours correspond to different values of $$\mu _\Phi$$, as indicated in the legend. Fig. 6 shows the protostellar mass distributions at the termination of each model. The black curve in the figure shows the observed IMF of Chabrier (2005). The mass distributions in our models have stabilized over dynamical time-scales (tff, >med), as shown in Fig. 5. However, because the sources in our models are still accreting, our models are not exactly comparable to the IMF that is representative of the distribution of final stellar masses. For this reason, we also compare our models to the protostellar mass functions (PMFs) from the analytic theory of McKee & Offner (2010); these PMFs are the mass distributions that one expects a series of protostars to possess if their final mass distribution matches the Chabrier (2005) IMF, and if the numbers of stars formed versus time N*(t) and the accretion histories of individual stars follow some specified functional form. We observe that N*(t) is a roughly linear function in all simulations except those that display runaway collapse, and thus we adopt a linear function for N*(t). For the accretion rate versus time of individual protostars, McKee & Offner (2010) consider a range of cases, all of which we plot together with the simulation results in Fig. 6: competitive accretion (CA), turbulent core (TC), 2-component competitive accretion (2CCA), and 2-component turbulent core (2CTC); the latter two both approach the constant accretion rate expected for an isothermal sphere for low-mass sources (IS), and we also plot the pure IS case.4 The only other parameter required by the McKee & Offner (2010) PMF analytic model is the upper limit to the most massive star that will form in the cluster mu. We have generated the theoretical PMFs in Fig. 6 by setting mu to the larger of 3  M⊙ or the value of mu necessary so that the theoretical PMF has exactly one source heavier than the second most massive particle present in the termination of the simulation. The choice of mu largely affects the shape of the high-mass tail of the PMF and the conclusions that follow from comparing our simulation to the theoretical PMFs are insensitive to this choice. The p-value of the Kolmogorov–Smirnov (K-S) statistic comparing each simulation to each mass function is reported in Table 3. The relative disagreement between simulations and mass functions pairs with p-values less than 0.05 is unlikely due to random chance. These distribution pairs are different with high statistical confidence. We denote the p-value of these pairs in bold text in the table. Figure 6. View largeDownload slide Histograms of protostellar masses at the end of each of the models with turbulent forcing (right) and with decaying turbulence (left). The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ from top to bottom. The black line shows the observed IMF (Chabrier 2005). The other curves denoted in the figure legend show the analytic PMFs of McKee & Offner (2010) for the single- and two-component turbulent core models, the single- and two-component competitive accretion models and the isothermal sphere model. Figure 6. View largeDownload slide Histograms of protostellar masses at the end of each of the models with turbulent forcing (right) and with decaying turbulence (left). The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ from top to bottom. The black line shows the observed IMF (Chabrier 2005). The other curves denoted in the figure legend show the analytic PMFs of McKee & Offner (2010) for the single- and two-component turbulent core models, the single- and two-component competitive accretion models and the isothermal sphere model. Table 3. K-Sstatistic p-values for the distribution of protostars in each simulation as compared to the observed IMF (Chabrier 2005) and several analytical PMFs including the TC, 2CTC, CA, 2CCA, and the isothermal sphere models (McKee & Offner 2010). Simulations that include turbulent forcing are denoted as ‘Driven’ in the left-hand column. The other models are the cases of decaying turbulence. Simulation to mass function pairs with p-values less than 0.05 can be confidently interpreted as a sample drawn from a different distribution than the given mass function. The p-value for these ‘rejected’ models is denoted in boldface in the table. Comparison analytic model  IMF  TC  2CTC  CA  2CCA  IS  $$\mu _\Phi =1.56$$  3.2 × 10−4  0.094  1.7 × 10−3  0.69  1.0 × 10−3  1.1 × 10−7  $$\mu _\Phi =1.56$$ No wind  1.1 × 10−3  7.7 × 10−5  1.0 × 10−3  6.0 × 10−6  4.4 × 10−4  0.054  $$\mu _\Phi =2.17$$  0.19  0.92  0.27  0.56  0.25  0.013  $$\mu _\Phi =23.1$$  0.058  4.2 × 10−3  0.014  1.1 × 10−3  9.9 × 10−3  0.35  $$\mu _\Phi =\infty$$  0.50  1.3 × 10−3  0.011  1.4 × 10−4  9.2 × 10−3  0.35  $$\mu _\Phi =1.56$$ Driven  3.0 × 10−3  0.52  0.13  0.91  0.080  2.1 × 10−3  $$\mu _\Phi =2.17$$ Driven  3.5 × 10−3  0.72  0.086  0.94  0.050  1.6 × 10−4  $$\mu _\Phi =23.1$$ Driven  0.033  0.19  0.032  0.63  0.033  2.9 × 10−5  $$\mu _\Phi =\infty$$ Driven  4.4 × 10−3  0.43  0.036  0.60  4.4 × 10−3  1.9 × 10−5  Comparison analytic model  IMF  TC  2CTC  CA  2CCA  IS  $$\mu _\Phi =1.56$$  3.2 × 10−4  0.094  1.7 × 10−3  0.69  1.0 × 10−3  1.1 × 10−7  $$\mu _\Phi =1.56$$ No wind  1.1 × 10−3  7.7 × 10−5  1.0 × 10−3  6.0 × 10−6  4.4 × 10−4  0.054  $$\mu _\Phi =2.17$$  0.19  0.92  0.27  0.56  0.25  0.013  $$\mu _\Phi =23.1$$  0.058  4.2 × 10−3  0.014  1.1 × 10−3  9.9 × 10−3  0.35  $$\mu _\Phi =\infty$$  0.50  1.3 × 10−3  0.011  1.4 × 10−4  9.2 × 10−3  0.35  $$\mu _\Phi =1.56$$ Driven  3.0 × 10−3  0.52  0.13  0.91  0.080  2.1 × 10−3  $$\mu _\Phi =2.17$$ Driven  3.5 × 10−3  0.72  0.086  0.94  0.050  1.6 × 10−4  $$\mu _\Phi =23.1$$ Driven  0.033  0.19  0.032  0.63  0.033  2.9 × 10−5  $$\mu _\Phi =\infty$$ Driven  4.4 × 10−3  0.43  0.036  0.60  4.4 × 10−3  1.9 × 10−5  View Large We find that our runs with weak initial fields ($$\mu _\Phi =23.1$$ and ∞) and no driving have peak masses that exceed the observed mode of the IMF (m* = 0.2  M⊙ – Chabrier 2005), and that are continuing to rise (Fig. 5), at the termination of the models. As discussed in Section 3.1, these models also approach a state of rapid, near global free-fall collapse, contrary to observational constraints. Furthermore, the shape of the mass distributions is closest to that of the IS case. The ever-increasing peak of the PMFs in these simulations appears to be a manifestation of the overheating problem identified by Krumholz et al. (2011), whereby values of εff close to unity produce such high accretion luminosities that fragmentation is strongly suppressed, leading to a mass function that shifts to ever-higher masses as collapse continues but no new stars can be created. Furthermore, these simulations have low K-S statistic p-values when compared against all analytic models except for the IS one. In contrast, the moderate and strong initial field cases all achieve stable or slowly increasing median masses in the range 0.05 M⊙ < m* < 0.15 M⊙, independent of the presence of external turbulent forcing. All of the models with continuously driven turbulence also show stable peak masses independent of the strength of the initial magnetic field. At the termination of these models, the median masses lie between the analytically predicted peak masses for the single CCA and single CTC models of 0.05 and 0.095  M⊙, respectively. Because of the median masses for the CA and TC models are close to those of the numerical simulations, these simulations also have high K-S p-values when compared to the CA and TC theoretical models. The shape of the mass functions are not amenable to reliably differentiating among the moderately magnetized models or the weakly magnetized models with continuous turbulent forcing to within the statistical limitations of the number of sources in these simulations. Fig. 7 shows the protostellar mass distribution at the termination of the strongly magnetized $$\mu _\Phi =1.56$$ without turbulent driving or protostellar outflow ejection. The effect of outflows on the protostellar mass distribution can be seen clearly by comparing this figure to the top left panel of Fig. 6, which shows the undriven $$\mu _\Phi =1.56$$ model with outflows. Clearly omission of outflows leads to a significant increase in the typical stellar mass and very low K-S test p-values. This difference is partly due to direct removal of mass by the outflows, but this effect alone cannot explain the difference in mass distribution, partly because it is too large (a factor of ≳10 in mass), and partly because there is also a dramatic drop in the total number of stars when we turn-off outflows, something that cannot easily be explained by mass ejection. Instead, the shift in the mass function here appears to result from the interaction between outflows and radiative heating. Outflows indirectly reduce the luminosity of individual stars by lowering their accretion rates (Krumholz et al. 2012; Hansen et al. 2012; Myers et al. 2014), and create paths of reduced optical depth, allowing more efficient escape of radiation from the gas immediately around the protostars (Krumholz, McKee & Klein 2005; Cunningham et al. 2011). Due to the omission of these two effects in the no-outflow model, radiative heating is much more effective at suppressing fragmentation, leading to the production of fewer, much more massive stars. We showed in Section 3.1 that the strong initial magnetic support in the no-outflow model was sufficient to limit the star formation efficiency to values consistent with observation, but clearly the shape of the mass distribution of this model is not consistent with observation. Figure 7. View largeDownload slide Histogram of the protostellar masses at the end of the $$\mu _\Phi =1.56$$ model without protostellar feedback. Figure 7. View largeDownload slide Histogram of the protostellar masses at the end of the $$\mu _\Phi =1.56$$ model without protostellar feedback. It is also interesting to compare our results to previous simulations making use of a similar set of physical processes. Our $$\mu _\Phi =1.56$$ model with continuous driving is comparable to the recent models by Li et al. (2018); the simulations have similar initial (βIC and $$\mathcal {M}_{{\rm A},\rm IC}$$) and post-driving (β>med, $$\mathcal {M}_{{\rm A},\rm >med}$$ and ρ>med) conditions. However, while in our case this model produces an IMF peak at slightly less than 0.1  M⊙, the Li et al. (2018) simulations yield at peak at 0.2  M⊙, comparable to the observed IMF peak. Given the similarities of the dimensionless conditions of the two cluster models, we attribute this difference to the outflow ejection parameters discussed in Section 2.3. The models presented here eject three times the momentum flux in the early phases of pre-stellar collapse and therefore have more disruptive effects on their parent cores. While Cunningham et al. (2011) point out that observations do not precisely constrain the outflow model ejection parameters, the fact that the model in Li et al. (2018) achieves better agreement with the well-constrained median IMF suggests that the outflow ejection parameters used in that model are probably more realistic. This conclusion, however, does not affect the relative differences between the various runs that we have identified; indeed, a shift to higher masses as a result of somewhat weaker outflows would worsen the agreement between the observed IMF and those produced in the weakly magnetized, non-driven models. The other obvious prior work to which we can compare are the simulations of Hansen et al. (2012). Indeed, our non-magnetized, decaying turbulence model is identical to theirs apart from the random phases used during the initial turbulent driving phase. Both models ran for comparable times as well. However, the model presented in this work exhibits rapidly increasing accretion rates and commensurately increasing source masses at late time, while the simulations of Hansen et al. do not. While this might at first seem quite surprising, the result can be understood by examining Fig. 1. This figure shows that the stellar mass at first increases relatively slowly, but then rapidly increases as the turbulence decays and free-fall collapse begins. The behaviour in the Hansen et al. simulation is consistent with what we see just before the onset of rapid collapse in our simulation. Since the exact time at which runaway collapse begins is almost certainly sensitive to the random driving, the most likely explanation for the difference is that the Hansen et al. simulation was simply not run quite long enough to reach the runaway collapse phase where the collapse begins to alter the IMF. 3.4 Stellar multiplicity The statistics of multiple systems provides another observational constraint against which our models can be tested. To compute the stellar system multiplicity, we use the algorithm of Bate (2009). The algorithm works recursively: one identifies the most strongly bound pair of stars in the simulation, and replaces it with a single object at the centre of mass position of the pair with the total mass and momentum of its constituent protostars, then repeats until there are no more bound pairs that can be combined. In cases where combining the most bound pair of objects would create a system with a multiplicity greater than four, one skips it and combines the next most bound combination instead. This is motivated by the fact that high-multiple bound systems are dynamically unstable and such systems are likely to break apart if the model were continued further in time. At the end of this procedure, we obtain S single star systems, B binary systems, T triple systems, and Q quadruple systems. For the set of bound systems, the multiplicity fraction is defined as (Bate 2009; Krumholz et al. 2012; Myers et al. 2014)   $$\textrm{MF} = \frac{B+T+Q}{S+B+T+Q}$$ (13) Fig. 8 shows the multiplicity fraction as a function of the primary star mass in each of our simulations. Given the limited number of systems in at least some of our simulations, computing the multiplicity fraction requires some care. We do so in two ways. First, we form running averages, meaning that, for each primary star (i.e. single star or the most massive member of a multiple), we evaluate equation (13) for the set of systems consisting of that primary and the next most massive and next least massive primaries. This quantity is shown as the blue curve in Fig. 8. Second, we can compute multiplicity in bins. To do so, we start with the least massive primary, and then place the primaries in bins of mass, with the bin width taken to be either sufficient to contain four primaries, or 0.2 dex, whichever is greater. For each bin, we regard the systems in it as samples drawn from a binomial distribution of multiples or singles, and from those samples, we compute a central estimate and a 68 per cent confidence interval on the true multiple system fraction in that bin (see Krumholz et al. 2012 for full details on how this computation is carried out). We plot the binned 68 per cent confidence intervals as the shaded regions in Fig. 8. Figure 8. View largeDownload slide Multiple stellar system fraction at the end of each of the models with turbulent driving (right) and with decaying turbulence (left), computed as a function of primary star mass. The four rows show the $$\mu _\Phi =1.56, 2.17, 23.1$$, and ∞, from top to bottom. For the simulations, the blue curve indicates evolution of a running average over three systems as a function of mass, while shaded regions indicate the mean multiplicity and central 68 per cent probability range on the multiplicity in bins – see the main text for details on how the running average and binned values are computed. Black crosses denote the mass range and mean multiplicity found in observational surveys. The two highest mass surveys are statistical lower bounds. The horizontal marker on these observations indicate the lower bound. The observational data, from low to high mass, are compiled from Basri & Reiners (2006) and Allen (2007) (combined into the left most point), Fischer & Marcy (1992), Raghavan et al. (2010), Preibisch et al. (1999), Mason et al. (2009), and these are the same data used in Bate (2012) and Krumholz et al. (2012). Figure 8. View largeDownload slide Multiple stellar system fraction at the end of each of the models with turbulent driving (right) and with decaying turbulence (left), computed as a function of primary star mass. The four rows show the $$\mu _\Phi =1.56, 2.17, 23.1$$, and ∞, from top to bottom. For the simulations, the blue curve indicates evolution of a running average over three systems as a function of mass, while shaded regions indicate the mean multiplicity and central 68 per cent probability range on the multiplicity in bins – see the main text for details on how the running average and binned values are computed. Black crosses denote the mass range and mean multiplicity found in observational surveys. The two highest mass surveys are statistical lower bounds. The horizontal marker on these observations indicate the lower bound. The observational data, from low to high mass, are compiled from Basri & Reiners (2006) and Allen (2007) (combined into the left most point), Fischer & Marcy (1992), Raghavan et al. (2010), Preibisch et al. (1999), Mason et al. (2009), and these are the same data used in Bate (2012) and Krumholz et al. (2012). Comparing either the running averages or the binned results to the observed multiplicity fractions that we also plot in Fig. 8, we find that all of the models generally agree with observations within the range of statistical uncertainties. It is noteworthy that the weakly magnetized models without turbulent forcing are consistent with observed multiplicity statistics even though their PMFs (Section 3.3) are strongly inconsistent with observations. The driven turbulence $$\mu _\Phi =2.17$$ and 23.1 models somewhat overpredict the multiplicity fraction for systems with primary mass Mp ∼ 1  M⊙, but this may be the result of poor statistics in this mass range in those models, and we note that the observational surveys are based on field stars, while there are strong hints that multiplicity fraction is higher for still-embedded systems (Tobin et al. 2016), which is the more relevant comparison for our simulations.5 We also note that, while the bins at low mass are generally consistent with the observations individually, taken together it is clear that we systematically underproduce multiple systems at low mass relative to the observational data. This is a likely a resolution effect. Our models do not resolve gas-particle gravity forces below the scale where sink source terms couple to the gas –4Δx. Low-mass multiple systems must necessarily be close and consequently more difficult to resolve to remain bound. This limitation is also noted in the multiplicity fractions extracted from the models of Bate (2009), Krumholz et al. (2012), and Myers et al. (2014). 3.5 Turbulent decay and protostellar feedback rates Fig. 9 shows the temporal evolution of the volume-weighted rms gas velocity for each model. As intended, the models that include steady turbulent forcing remain within 10 per cent of a constant rms velocity with αvir ∼ 1. For the runs without driving, at early time, t ≲ tff, turbulent decay causes the rms velocity to drop. However, after 1–2 free-fall times (depending on the level of magnetic support), the system begins to undergo a global collapse, and this causes the velocity dispersion to rise again. Figure 9. View largeDownload slide Volume-weighted root-mean-squared gas velocity, normalized by the gas sound speed, as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). Note the difference in scales between the left- and right-hand panels. Figure 9. View largeDownload slide Volume-weighted root-mean-squared gas velocity, normalized by the gas sound speed, as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). Note the difference in scales between the left- and right-hand panels. Fig. 10 shows the temporal evolution of the mass-weighted rms velocity dispersion. This averaging is more analogous to mass-biased observational line width–size measurements. However, it is also heavily biased towards dense, locally collapsing regions, which somewhat mutes the effect of turbulent decay relative to infall. Consequently, the clear drop in velocity dispersion seen in Fig. 9 for the non-driven runs is absent here. We infer that mass-weighted line-width diagnostics are surprisingly insensitive to interruptions in large-scale driving source of durations up to ≲ 1tcross. Figure 10. View largeDownload slide Same as Fig. 9, but now showing mass-weighted rather than volume-weighted velocities. As with Fig. 9, note the difference in scales between the two panels. Figure 10. View largeDownload slide Same as Fig. 9, but now showing mass-weighted rather than volume-weighted velocities. As with Fig. 9, note the difference in scales between the two panels. It is also interesting to compare the $$\mu _\Phi =1.56$$ cases with and without outflows. In principle, the outflows have more than enough mechanical power to alter the turbulence. We illustrate this in Fig. 11, which shows the outflow mechanical luminosity versus time, $$L_{\rm outflow} = (1/2) \dot{M}_{\rm wind} v_{\rm wind}^2$$. This is typically ∼1–10 L⊙ at late times in our simulations, For comparison, the turbulent decay rate is $${\sim } \bar{\rho } v_{\rm rms}^2 L^3/t_{\rm cross} \sim 10^{-3}\ {\rm L}_{{\odot }}$$, roughly three orders of magnitude smaller. However, when we examine the volume-weighted rms velocity, the $$\mu _\Phi =1.56$$ runs with and without outflows are nearly identical, and in the mass-weighted plot they are only slightly different. These results indicate that outflow mechanical power is efficiently dissipated by radiative shocks and does not quickly couple to the volume-filling turbulence of the surrounding molecular cloud, consistent with earlier works exploring the coupling of continuously driven jets through a stratified environments (Banerjee, Klessen & Fendt 2007). What effects the outflows do have are limited to the dense gas surrounding protostellar cores, which is why the mass-weighted average shows an effect from outflows but the volume-weighted one does not. However, we emphasize that even in a mass-weighted sense the effect is modest. In our simulations, outflows do not appear to be efficient drivers of turbulence. Figure 11. View largeDownload slide Total mechanical luminosity of the protostellar outflows as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). The data have been boxcar-smoothed with an interval of 0.05 tff. Figure 11. View largeDownload slide Total mechanical luminosity of the protostellar outflows as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). The data have been boxcar-smoothed with an interval of 0.05 tff. 3.6 Local magnetic support The turbulent driving phase introduces stretching and amplification of the magnetic field, relative to the initial magnetic field strength. The Alfvèn Mach number captures the strength of the turbulence relative to the initial magnetic field. As discussed in Section 2, the magnetized models in this work span the range $$1 < {\cal M}_{\rm A}< 15$$. However, the weakly magnetized models endure more significant stretching of the magnetic field and more small-scale magnetic field amplification in the turbulent driving phase than the more strongly magnetized models. The tangled magnetic field is further amplified due to collapse after self-gravity is switched on. Therefore, it is only the mean, large-scale field that is weak in our ‘weak-field’ model. On the scale of the pre-stellar dense cores, the magnetic field is significantly more tangled and stretched than in stronger field models. Similar field amplification due to turbulence and the early gravitational collapse phase into starless dark clumps was noted in Li et al. (2015, section 8 in particular). In this section, we examine the magnetic field structure in the vicinity of protostars, with the goal of understanding how the combination of turbulent and gravitational amplification maps from the initial, large-scale magnetic field to the final, small-scale one in the dense gas that is bound to and actively feeding accretion on to protostellar sources. We seek to characterize the mean mass-to-flux ratio, averaged over the bound gas around each protostar that is poised for imminent accretion. We take this gas to be that which satisfies the condition   $$\frac{(\boldsymbol {v}-\boldsymbol {v_*})\cdot (\boldsymbol {x}-\boldsymbol {x_*})}{|\boldsymbol {x}-\boldsymbol {x_*}| } + c_{\rm s} < \sqrt{\frac{Gm_*}{|\boldsymbol {x}-\boldsymbol {x_*}|}},$$ (14)where $$\boldsymbol {x}$$, $$\boldsymbol {v}$$, and cs are the position, velocity, and sound speed of the gas in the computational domain and $$\boldsymbol {x_*}$$, $$\boldsymbol {v_*}$$, and m* are the position, velocity, and mass of each protostar. This criterion selects gas with insufficient thermal or kinetic energy to escape eventual accretion to the protostar, but distinguishes between infalling gas (for which $$\boldsymbol {v}-\boldsymbol {v_*}$$ is anti-aligned with $$\boldsymbol {x} - \boldsymbol {x_*}$$, and thus the dot product is negative), which is more bound, and outflowing gas (vectors aligned), which is less so. We assign gas that is bound to more than one star under this criterion only to the star to which it is most strongly bound. For each star i, we have now defined a mass mgas, i of ‘bound gas.’ We next compute the volume-weighted mean magnetic field direction in this gas, $$\hat{\boldsymbol B}_{\rm i}$$, and we define the associated magnetic flux Φi as the net flux through the intersection of the plane defined by $$\hat{\boldsymbol B}_{\rm i}$$ and the volume of assigned gas. Finally, we characterize the local mass-to-flux ratio of the gas assigned to each star as μΦ, i = mgas, i/MΦ, i where MΦ, i is the magnetic critical mass computed from equation (1). It should be noted that the significance of $$\mu _\Phi$$ changes if a protostar is present in the cloud: The protostellar gravity also acts on the gas, so that accretion can occur normal to the field even for $$\mu _\Phi <1$$. Fig. 12 shows the temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The average is weighted by the mass of each protostar. Due to the computational burden associated with post-processing this result from our models, the plots in Fig. 12 are sampled only every 0.2tff, >med in time. In the figure, the curve for each model begins after the formation of the first sink particle. Figure 12. View largeDownload slide Temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The averaging is weighted by protostellar mass. The plot is presented in units of the magnetically critical mass, for the models with turbulent forcing (right) and with decaying turbulence (left). Figure 12. View largeDownload slide Temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The averaging is weighted by protostellar mass. The plot is presented in units of the magnetically critical mass, for the models with turbulent forcing (right) and with decaying turbulence (left). In interpreting Fig. 12, it is important to bear in mind that this diagnostic characterizes the net impact of multiple competing processes. Under our sink particle prescription (Lee et al. 2014), infalling gas is released from its associated magnetic flux the instant that it deposited on to the sink particle. This prescription is intended to mimic resistive effects near the surface of the protostar at length scales below that which can be resolved by our models. Consequently, accretion causes the mass-to-flux ratio associated with locally bound dense gas clumps to decrease. The accreted mass is removed from the gas but the flux remains. However, the magnetic flux may subsequently escape, because it is no longer anchored to the gravitationally bound material – this is a form of magnetic interchange instability (Zhao et al. 2011; Krasnopolsky et al. 2012; Cunningham et al. 2012; Li et al. 2014). If the magnetic flux tubes escape the gravitationally bound region, this drives the mass-to-flux ratio back up. Similarly, turbulent magnetic reconnection acts to diffuse magnetic flux and increases the local mass-to-flux ratio in regions that undergo reconnection. In this case, the gas mass remains but the magnetic flux is displaced outward (Santos-Lima, de Gouveia Dal Pino & Lazarian 2012). Outflow ejection and turbulent forcing also impart mechanical energy to the system, and the resulting flows can further amplify, entangle, and ultimately reconnect magnetic flux tubes (Sur et al. 2012; Li et al. 2015). In this way, accretion, outflow, field tangling, and reconnection all influence the mass-to-flux ratio. Since accretion itself is episodic, so too is the evolution of the mean mass-to-flux ratio of the bound gas in Fig. 12. We anticipate the degree of tangling and episodic reconnection to be most pronounced in the weak-field models as weaker fields are more readily tangled and amplified by the turbulent flow. This too is borne out in Fig. 12. In the undriven models shown in the left-hand panel of Fig. 12, the mass-to-flux ratios in the bound gas systematically decreases with time for each of the models. This indicates that buoyant flux tubes rising away from the protostars as a consequence of accretion act enhance the degree of magnetic support to the bound gas, and this enhancement overwhelms the effect of turbulent reconnection. The $$\mu _\Phi =1.56$$ model without outflows indicates systematically stronger magnetic support in bound gas than the case with outflows. This indicates that outflows contribute to magnetic reconnection and consequent outward transport of magnetic flux, and explains why these models achieve comparable star formation efficiency (εff, IC = 12 per cent in case without outflows and εff, IC = 23 per cent in case with outflows). The enhanced magnetic support in case without outflows offsets, the loss of mechanical feedback support when the outflows are turned off. The $$\mu _\Phi =2.17$$ model is also strongly magnetically supported, with the mass-to-flux ratio in the bound gas approaching $$\mu _\Phi \sim 1$$ at late time, consistent with its accretion rate εff, IC = 14 per cent, slightly higher than the more strongly magnetized $$\mu _\Phi =1.56$$. The model that began with a field characterized globally as $$\mu _\Phi =23.1$$ attains a characteristic mass-to-flux ratio $$\mu _\Phi \approx 5$$ via field amplification and gravitational compression by the time the first protostellar particle forms. While the level of magnetic support in the bound gas increases with time, the degree of magnetic support achieved by the end of the model is insufficient to significantly reduce the net accretion rate from that of the purely hydrodynamical case. The evolution of the mass-to-flux ratio of the bound gas in the models with driven turbulence models is shown in the right-hand panel of Fig. 12. In this case, the model that began with a global $$\mu _\Phi =2.17$$ retains a roughly constant mass-to-flux ratio in the bound gas of $$\mu _\Phi \sim 2$$ throughout its evolution. The model that began with a global $$\mu _\Phi =1.56$$ attains stronger magnetic stabilization in the dense gas with $$\mu _\Phi < 1$$ by the time the first protostars have formed. Subsequent evolution indicates a slowly decreasing level of magnetic support in the bound gas, stabilizing to $$\mu _\Phi \gtrsim 1$$. The model that began with a global $$\mu _\Phi =23.1$$ enhances the degree of magnetic support in the bound gas with time, reaching to $$\mu _\Phi \sim 4$$ by the end of the model. Taken in aggregate, our results indicate that small-scale dissipation effects near the surface of a protostar provide a feedback loop for stabilizing the star formation efficiency. Bound gas cores that are undersupported initially undergo more rapid collapse than cores that are better supported, independent of whether the initial support is predominately turbulent or magnetic. The more rapidly collapsing regions reconnect and expel magnetic flux to the core from the newly formed protostar at a faster rate. The end result of this process over a time-scale >1tff is that cores with continuous mechanical support that remain in near virial equilibrium approach mass-to-flux ratios of $$\mu _\Phi < 5$$ independent of the large-scale mean field strength. Models with decaying or otherwise sub-virial levels of turbulent support will lead to greater local enhancement of the level of magnetic support in the bound gas, approaching $$\mu _\Phi < 2$$ at late time. 4 MAGNETIC AND RADIATIVE EFFECTS IN SHAPING THE IMF In Section 3.3, we demonstrated that our models have protostellar mass distributions that are stable or evolving slowly towards higher mass only slowly compared to cloud-collapse time-scales. The models with weak initial fields or a lack of driving all have median masses that are too large in comparison to the observed IMF, while the strong field and driven cases have smaller median masses that are at least roughly compatible with these models reproducing the observed IMF. It is therefore interesting to investigate in more depth what mechanisms shape this mass distributions. To explore this question, we turn to the techniques described in Krumholz et al. (2016), who analysed the simulations of Myers et al. (2014); these simulations use the same physics we include here, but focus on a much denser region appropriate to sites of high-mass star formation. Comparing the two sets of simulations therefore enables us to understand how this question depends on the star-forming environment. 4.1 Analysis methodology We begin with a brief summary of the Krumholz et al. (2016) analysis, and refer readers interested in the details to that paper for a full description. The central idea of this analysis is to determine what physical process is responsible for halting fragmentation in the vicinity of each protostar, thereby leaving the mass available for accretion and the build-up of more massive stars, rather than the production of additional low-mass stars. To this end, we examine spherical regions of gas within a $$4000\ {\rm {\rm au}}$$ radius around each protostellar sink particle at time intervals of 0.05tff, >med. We divide the volume of interest into 128 spherical shells of radius r concentric with the position of the sink particle. For each shell, we compute the total gas enclosed mgas, excluding the sink particle mass, the mean density ρ of the enclosed material, and its mass-weighted mean sound speed cs.6 From the resulting profiles in ρ and cs, we compute an effective temperature   $$T_{\rm eff} = \frac{\mu m_{\rm H}}{k_{\rm B}} c_{\rm s}^2 ,$$ (15)and, more importantly from the standpoint of determining whether gas will fragment, the Bonnor–Ebert mass (Ebert 1955; Bonnor 1956),   $$m_{\rm BE} = 1.86 \sqrt{\frac{c_{\rm s}^3}{G^3\rho }}.$$ (16)The Bonnor–Ebert mass characterizes the degree of thermal support. Absent other forces, objects less massive than mBE are stable against collapse and objects more massive than mBE are unstable against collapse. To elucidate the importance of heating by protostars in stabilizing the core, we repeat the computation with the sound speed cs = 0.19 cm s−1, the sound speed of molecular gas at 10 K used as the background isothermal temperature in the simulations. We shall denote this quantity as mBE, 10. Comparison of mBE with mBE, 10 enables us to determine the importance of radiative heating in limiting fragmentation. To characterize the importance of magnetic forces, we compute the magnetic flux Φ threading each spherical shell, defined so that $$\Phi = \int |\boldsymbol {B}\cdot \hat{n}|\, d^2x$$ is computed along a circular cross-section that is concentric with the spherical shell and has normal $$\hat{n}$$. The absolute value is taken in the integration on the assumption that oppositely directed fluxes do not reconnect and cancel. We consider 12 possible orientations for $$\hat{n}$$, uniformly distributed on the unit sphere, with one value aligned with the star's angular momentum vector, and use the largest value of Φ, we find for all the $$\hat{n}$$ values considered. From this, we can compute an effective mean magnetic field strength inside a radius r,   $$B_{\rm eff}(<\!r) = \frac{\Phi }{\pi r^2},$$ (17)a magnetic critical mass $$m_\Phi$$ (equation 1), and the magnetically supported mass (Mouschovias & Spitzer 1976; McKee & Ostriker 2007)   $$m_{\rm B}=\frac{m^3_\Phi }{m^2_{\rm gas}}.$$ (18)We prefer mB to $$m_\Phi$$ for this analysis as the former is an intensive quantity that does not depend on the size of the volume considered. It has been well established that magnetic fields influence the spectrum of stars formed and reduce the star formation efficiency (Price & Bate 2008). Here, we elucidate the relative contribution or radiative feedback to magnetic support by comparing the enclosed mass mgas with the thermally and magnetically supported masses mBE and mB. This provides direct insight into the role of thermal support (augmented by radiative heating) and magnetic fields in shaping the IMF. For sufficiently small shells around each protostar, we always find mgas ≪ mBE and mgas ≪ mB, i.e. the enclosed mass is too small to fragment and form another star given the level of thermal and magnetic supports. This gas will almost certainly be accreted on to the star that it surrounds, augmenting this star's mass. At sufficiently large radii, mgas ≫ mBE and mgas ≫ mB, i.e. the mass enclosed is large enough that is able to fragment and produce additional stars, rather than being accreted on to the star it surrounds. Thus, the shell at which mgas ≈ max (mBE, mB) marks a rough dividing line between gas that will end up in the existing star and gas that will end up elsewhere, and m* + mgas provides a rough estimate of the mass to which the star is likely to grow. We define the thermal critical mass mBE, crit as the smallest (and almost always sole) mass where mgas = mBE, and analogously for mB, crit and mBE, 10, crit. Thus these critical masses provide a rough estimate of the future mass supply available for each star. Krumholz et al. (2016) show that, in the simulations of Myers et al. (2014), stars with mass m* ≪ 0.2  M⊙ almost invariably have stabilized gaseous envelopes such that mBE, crit ≫ m*, explaining the relative rarity of brown dwarfs as compared to stars: brown dwarf mass fragments can form, but when they do, they stabilize enough gas around themselves that they usually continuing accreting well out of the brown dwarf mass regime. Moreover, Krumholz et al. (2016) find that their stabilized regions have mBE, crit ≫ mB, crit ≫ mBE, 10, crit, indicating that, while magnetic support is more important than thermal support in the absence of radiative heating, once one considers the effects of stellar radiation, it is the dominant factor in limiting fragmentation and thus establishing the location of the IMF peak. 4.2 Profiles of mean density, temperature, and magnetic field Fig. 13 shows the average density profile around central stars with three different masses; the averages plotted are taken over all protostars at all times. The density profiles always show a central dip at 130 au < r < 260 au, which is a result of our sink particle accretion zones; we therefore ignore this region in our analysis. At larger radii, we find that the moderately and weakly magnetized models ($$\mu _\Phi \ge$$ 2) without driving have density profiles ρ ∝ r−3/2. This profile is consistent with that expected for free-fall collapse on to a point mass. However, for Bondi-type flow, the density at a fixed distance from the central object should increase with mass. The fact that it is not inconsistent with pure Bondi flow, but is consistent with the TC model of McKee & Tan (2003) and with the model of Murray et al. (2015) for the for collapsing, self-gravitating, TCs. In contrast, the strongly magnetized ($$\mu _\Phi = 1.56$$) undriven cases, and all driven cases for a central star mass of m* = 0.15  M⊙, show slightly shallower density profiles. This is consistent with a build-up of magnetic flux resulting in a slower than free-fall collapse. Figure 13. View largeDownload slide Mean density profiles ρ(r) averaged in spherical shells around protostars. The panels show the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. The bottom row depicts results for the model with continued turbulence driving after the onset of self-gravity. Different colours show means for protostars of mass 0.025, 0.05, and 0.15  M⊙. For each mass bin, the central line shows the mean, while the shaded band shows the 1σ dispersion for protostars in that mass bin. The dashed black line shows the density scaling ρ ∝ r−3/2, expected for free-fall collapse. Figure 13. View largeDownload slide Mean density profiles ρ(r) averaged in spherical shells around protostars. The panels show the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. The bottom row depicts results for the model with continued turbulence driving after the onset of self-gravity. Different colours show means for protostars of mass 0.025, 0.05, and 0.15  M⊙. For each mass bin, the central line shows the mean, while the shaded band shows the 1σ dispersion for protostars in that mass bin. The dashed black line shows the density scaling ρ ∝ r−3/2, expected for free-fall collapse. Fig. 14 shows averaged profiles of Teff. Except for the $$\mu _\Phi = 1.56$$ undriven model, these profiles indicate increasing Teff with central star mass. As Krumholz et al. (2016) point out, to maintain a constant density profile, infall velocities and accretion rate must increase as the progenitor mass increases. Since accretion luminosity is the dominant source of radiative heating in low-mass sources, the temperature around them rises with their mass. The weakly magnetized models show the most rapid increase in heating with mass. Krumholz et al. (2016) noted a similar trend in models of denser high-mass star-forming regions, and also found that the temperature profiles were always close to Teff ∝ r−0.3. Here we find substantially shallower profiles, likely because the lower optical depth environment we are simulating is less effective at trapping the radiation (Krumholz & McKee 2008). Figure 14. View largeDownload slide Same as Fig. 13, but showing the mean effective temperature Teff(r) (equation 15). The black dashed line shows a Teff ∝ r−0.3. The initial temperature of the models is T = 10 K. Figure 14. View largeDownload slide Same as Fig. 13, but showing the mean effective temperature Teff(r) (equation 15). The black dashed line shows a Teff ∝ r−0.3. The initial temperature of the models is T = 10 K. In Fig. 15, we show the corresponding average profiles of Beff. Consistent with our findings in Section 3.6, despite the large differences in large-scale magnetic flux, there is very little difference between the small-scale magnetic field strengths around protostars in the different simulations. Nor does the effective field strength increase with star mass, as one might expect from a naive picture in which accretion drags magnetic flux towards a star. Figure 15. View largeDownload slide Same as Fig. 13, but showing the mean effective magnetic field Beff(<r) (equation 17). The black dashed horizontal line indicates the initial magnetic field strength before driving to the initial turbulent state. Figure 15. View largeDownload slide Same as Fig. 13, but showing the mean effective magnetic field Beff(<r) (equation 17). The black dashed horizontal line indicates the initial magnetic field strength before driving to the initial turbulent state. The non-dependence of the small-scale magnetic fields on either the large-scale flux or the amount of gas that has already been accreted strongly suggests that the magnetic field profiles we find are the result of a self-regulation mechanism operating on small scales. We can think of this self-regulation as due to two pairs of competing processes, although given our limited resolution and the complex nature of the flows in our simulations, we cannot easily identify which mechanisms dominate. First, there is a competition between dynamo amplification and reconnection. In a turbulent dynamo, magnetic fields amplify until these two effects balance (Sur et al. 2012; Li et al. 2015). Second, there is a competition between advection and magnetic interchange instabilities. As gas accretes, magnetic flux tubes are advected into the cores around the protostars, which tends to drive up Beff. The countervailing process is escape of flux through magnetic interchange instabilities. In reality, the interchange instabilities are triggered by resistive effects close to the surface of the protostar, which decouple flux tubes from the mass that anchors them and allows them to drift outward (Zhao et al. 2011; Krasnopolsky et al. 2012; Cunningham et al. 2012; Li et al. 2014); we model this by not accreting any magnetic flux into the star particles. If the flux liberated from the protostar diffused out to a radius R, then it would contribute a mean field   $$B_{\rm eff,*}=2G^{1/2}\left(\frac{L}{Z}\right)\frac{m_*}{\mu _\Phi R^2},$$ (19)where we have assumed that the protostar has formed from a cylinder of gas with a height Z less than the box height L. For m* = 0.15  M⊙ and $$\mu _\Phi =1.56$$, the results in Fig. 15 place a lower limit R > 1000 au on the distance that the field has diffused outward. 4.3 Critical mass In Fig. 16, we show the critical mass of gas supported by thermal pressure, mBE, thermal pressure fixed to the initial temperature, 10 K mBE, 10, and magnetic fields, mB, as a function of central stellar mass m* as described in Section 4.2. When mcrit > m*, the mass of the stabilized envelope around a protostar is much greater than the mass of the star itself, and thus is likely to greatly increase its mass. Figure 16. View largeDownload slide Critical mass of gas supported by thermal pressure mBE, thermal pressure fixed to the initial temperature 10 K mBE, 10 and magnetic fields mB as a function of central stellar mass m*. Shaded bands indicate 1σ dispersion over the protostars. The dashed black line indicates where the mass supported against collapse by pressure or magnetic forces is sufficient to double the present stellar mass, m* = mcrit. The panels shows the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. Figure 16. View largeDownload slide Critical mass of gas supported by thermal pressure mBE, thermal pressure fixed to the initial temperature 10 K mBE, 10 and magnetic fields mB as a function of central stellar mass m*. Shaded bands indicate 1σ dispersion over the protostars. The dashed black line indicates where the mass supported against collapse by pressure or magnetic forces is sufficient to double the present stellar mass, m* = mcrit. The panels shows the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. First consider thermal support. Krumholz et al. (2016) found that mBE was typically a factor of ∼10 larger than mBE, 10, indicating that radiative heating increased the amount of mass that could be thermally supported by a factor of ∼10. In our lower mass cluster models, the difference between mBE and mBE, 10, is less pronounced, closer to 0.5 dex than 1 dex. However, the difference is smaller not because mBE is smaller, but rather because mBE, 10 is larger. That is to say, in our lower density simulation, the characteristic density is lower, and thus the Bonnor–Ebert mass in unheated gas is larger. In comparison, the mass that can be supported after heating, mBE, is almost unchanged and independent of environment – the increase in the temperature of the gas has been compensated by an increase in its density. Next consider magnetic support. The results in Fig. 16 are consistent with the findings of Krumholz et al. (2016) that the magnetic fields are relatively unimportant in supporting the cores compared to thermal pressure. Again, this points to the importance of local processes in regulating the magnetic field structure in the vicinity of protostars. The magnetic flux in the $$\mu _\Phi =1.56$$ case is sufficient to support a mass $$m_\Phi =119\ {\rm M}_{{\odot }}$$ against collapse, but the field threading the local region around the protostars is strong enough to support only 0.01–0.1  M⊙. This result explains why the mode of the IMFs (Section 3.3) in our models are relatively insensitive to the large-scale magnetic field. The critical mass supported against fragmentation is determined primarily by thermal effects enhanced by radiation. Magnetic fields thus have a significant effect on the rate of star formation, but not on the IMF. 5 CONCLUSIONS We present a numerical study of low-mass star cluster formation, and how this process is affected by the strength of the large-scale magnetic field, the presence or absence of turbulent energy input, and the effects of protostellar outflows. Our simulations include MHD, radiative transfer and protostellar radiation feedback, and protostellar outflows, a combination of physical processes that has received little attention in the literature to date. Our primary finding is that our models are able to reproduce the observed (low) efficiency of star formation and the observed location of the IMF peak (to within a factor 2) only in models that include outflows and that do not undergo global collapse, either because they are supported by external energy input that drives the turbulence, or because they are close enough to the line of magnetic criticality for magnetic fields to inhibit collapse. In contrast, multiplicity statistics are insensitive to all variations in the initial conditions. Krumholz et al. (2011, 2012) had previously conjectured that there might be a link between the low efficiency of star formation and the location of the IMF peak, but this was based on a single set of non-magnetized simulations. Our much more extensive suite of radiation-MHD simulations puts this conclusion on a much firmer footing, and suggests a general principle: the IMF peak location is directly linked to the inefficiency of star formation. The primary mechanism controlling the IMF peak appears to be accretion-powered radiation from forming protostars, which produces enhanced thermal pressure on small scales in the vicinity of forming stars, as noted earlier by Bate (2009) as well. In comparison, the small-scale magnetic field structure around forming protostars is both insensitive to the large-scale magnetic field strength, and dynamically sub-dominant compared to radiation when it comes to regulating gas fragmentation. These results strongly favour a picture in which the small-scale magnetic field structure around protostars is regulated by local competitions between accretion of magnetized material and dynamo amplification, which tend to enhanced the field strength, and magnetic interchange instability and turbulent reconnection, which tend to suppress it. The competition between these processes tends to force the local mass-to-flux ratio in gas bound to protostars to ∼2 (normalized to the critical value), regardless of the large-scale field. Our models lack resolution to capture the effects of thermal support against gravitational contraction deep within very low-mass cores and we lack non-ideal MHD effects that are most strongly coupled to turbulent density profiles on small scales. However, these effects should influence our modelled protostellar mass distributions only on the low-mass tail m* ≲ 0.05 M⊙ and our conclusions are not sensitive to these effects. We acknowledge that higher resolution models with non-ideal effects and multiple realizations to capture better statistics at protostellar mass extrema could improve the comparisons with observation made here by capturing a wider range of source masses with higher fidelity. Acknowledgements Support for this research was provided by NASA through NASA ATP grants NNX13AB84G (RIK and CFM), NNX15AT06G (MRK), the US Department of Energy at the Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 (AJC and RIK), the NSF through grant AST-1211729 (RIK and CFM), and the Australian Research Council's Discovery Projects funding scheme (project DP160100695, MRK). This research was also supported by grants of high-performance computing resources from the National Center of Supercomputing Application through grant TGMCA00N020, under the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation grant number OCI-1053575, the computing resources provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and resources and services from the National Computational Infrastructure, which is supported by the Australian Government. LLN-JRNL-737724. Footnotes 1 In contrast to these results, Murray (2011) and Lee, Miville-Deschênes & Murray (2016) argue εff ∼ 1 in the most luminous star clusters, based on their high ratios of ionizing luminosity to molecular gas mass. This interpretation has been disputed by Feldmann & Gnedin (2011) and Krumholz (2014), who suggest that these ratios are high not because εff is large, but rather because very massive clusters have strong feedback that destroys their parent clouds efficiently, thereby rendering invalid Murray's and Lee et al.'s assumption that one can identify and measure the gas mass of these clusters’ progenitors. However, this dispute is immaterial for the low-mass, low-density clusters with which we are concerned here, for which there is a strong observational consensus that εff is low. 2 We note that the time required to saturate the magnetic energy depends on the initial field strength. For models initialized with magnetic fields much weaker than those considered in this paper, significantly more than two dynamical times are required to fully saturate the magnetic energy of the system due to small-scale dynamos (Federrath et al. 2011; Tricco, Price & Federrath 2016). However, the models considered in this work are too strongly magnetized for small-scale eddies to generate significant field amplification, and the magnetic power spectra have been shown to saturate within a few crossing times in these regimes (Li et al. 2012). 3 The choice of the temperature at which one should evaluate the mean opacity for a grey method such as ours is in fact very subtle, even in calculations that separately track the dust, gas, and radiation temperatures. This is because tabulated opacity tables, such as the ones from Semenov et al. (2003) that we use, generally assume that all three temperatures are equal, and different temperatures matter for different physical effects – for example, condensation of ice mantles on grains depends on the gas temperature, evaporation of grains depends on the dust temperature, and the radiation spectrum depends on the radiation temperature. Absent tabulated opacities that consider out-of-equilibrium dust, gas, and radiation fields, there is no single choice of temperature at which the opacity can be computed that properly captures all these divergent effects. However, greater accuracy could be obtained by distinguishing the absorption opacity, which is primarily affected by the radiation temperature, and the emission opacity, which is primarily determined by the dust temperature. 4 We do not consider any of McKee & Offner's (2010) ‘tapered accretion’ analytic models, since the median accretion rate in our simulations shows no signs of tapering (Fig. 5). 5 We unfortunately cannot place the data from Tobin et al. (2016) on Fig. 8, because the masses of the individual stars are largely unconstrained. 6 As in Krumholz et al. (2016), we focus on cumulative rather than differential quantities because they are less subject to noise; see appendix A of Krumholz et al. (2016) for a demonstration that this choice does not alter the qualitative results. REFERENCES Allen P. R., 2007, ApJ , 668, 492 CrossRef Search ADS   Banerjee R., Klessen R. S., Fendt C., 2007, ApJ , 668, 1028 CrossRef Search ADS   Basri G., Reiners A., 2006, AJ , 132, 663 CrossRef Search ADS   Bate M. R., 2009, MNRAS , 392, 590 CrossRef Search ADS   Bate M. R., 2012, MNRAS , 419, 3115 CrossRef Search ADS   Bonnor W. B., 1956, MNRAS , 116, 351 CrossRef Search ADS   Chabrier G., 2005, in Corbelli E., Palle F., eds, Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later . Springer, Dordrecht, p. 41 Chakrabarti S., McKee C. F., 2005, ApJ , 631, 792 CrossRef Search ADS   Commerçon B., Hennebelle P., Henning T., 2011, ApJ , 742, L9 CrossRef Search ADS   Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ , 725, 466 CrossRef Search ADS   Cunningham A. J., Klein R. I., Krumholz M. R., McKee C. F., 2011, ApJ , 740, 107 CrossRef Search ADS   Cunningham A. J., McKee C. F., Klein R. I., Krumholz M. R., Teyssier R., 2012, ApJ , 744, 185 CrossRef Search ADS   Ebert R., 1955, Z. Astrophys. , 37, 217 Evans N. J., II, Heiderman A., Vutisalchavakul N., 2014, ApJ , 782, 114 CrossRef Search ADS   Falgarone E., Troland T. H., Crutcher R. M., Paubert G., 2008, A&A , 487, 247 CrossRef Search ADS   Federrath C., 2015, MNRAS , 450, 4035 CrossRef Search ADS   Federrath C., Klessen R. S., 2012, ApJ , 761, 156 CrossRef Search ADS   Federrath C., Klessen R. S., Schmidt W., 2009, ApJ , 692, 364 CrossRef Search ADS   Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010, A&A , 512, A81 CrossRef Search ADS   Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ , 731, 62 CrossRef Search ADS   Feldmann R., Gnedin N. Y., 2011, ApJ , 727, L12 CrossRef Search ADS   Fischer D. A., Marcy G. W., 1992, ApJ , 396, 178 CrossRef Search ADS   Gammie C. F., Ostriker E. C., 1996, ApJ , 466, 814 CrossRef Search ADS   García-Burillo S., Usero A., Alonso-Herrero A., Graciá-Carpio J., Pereira-Santaella M., Colina L., Planesas P., Arribas S., 2012, A&A , 539, A8 CrossRef Search ADS   Goldbaum N. J., Krumholz M. R., Matzner C. D., McKee C. F., 2011, ApJ , 738, 101 CrossRef Search ADS   Goldsmith P. F., 2001, MNRAS , 557, 736 Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, ApJ , 747, 22 CrossRef Search ADS   Hennebelle P., Commerçon B., Joos M. et al.  , 2011, A&A , 528, A72 CrossRef Search ADS   Heyer M., Gutermuth R., Urquhart J. S., Csengeri T., Wienen M., Leurini S., Menten K., Wyrowski F., 2016, A&A , 588, A29 CrossRef Search ADS   Klein R. I., Fisher R. T., McKee C. F., Truelove J. K., 1999, Numer. Astrophys. , 240, 131 CrossRef Search ADS   Klessen R. S., Hennebelle P., 2010, A&A , 520, A17 CrossRef Search ADS   Krasnopolsky R., Li Z.-Y., Shang H., Zhao B., 2012, ApJ , 757, 77 CrossRef Search ADS   Krumholz M. R., 2014, Phys. Rep. , 539, 49 CrossRef Search ADS   Krumholz M. R., 2015, preprint (arXiv:1511.03457) Krumholz M. R., McKee C. F., 2008, Nature , 451, 1082 CrossRef Search ADS PubMed  Krumholz M. R., Tan J. C., 2007, ApJ , 654, 304 CrossRef Search ADS   Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ , 611, 399 CrossRef Search ADS   Krumholz M. R., McKee C. F., Klein R. I., 2005, ApJ , 618, L33 CrossRef Search ADS   Krumholz M. R., Matzner C. D., McKee C. F., 2006, ApJ , 653, 361 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., Bolstad J., 2007, ApJ , 667, 626 CrossRef Search ADS   Krumholz M. R., Cunningham A. J., Klein R. I., McKee C. F., 2010, ApJ , 713, 1120 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ , 740, 74 CrossRef Search ADS   Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ , 745, 69; Erratum 2013, ApJ, 779, 89 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ , 754, 71 CrossRef Search ADS   Krumholz M. R., Myers A. T., Klein R. I., McKee C. F., 2016, MNRAS , 460, 3272 CrossRef Search ADS   Lada C. J., Lombardi M., Alves J. F., 2010, ApJ , 724, 687 CrossRef Search ADS   Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ , 745, 190 CrossRef Search ADS   Lada C. J., Roman-Zuniga C., Lombardi M., Forbrich J., Alves J. F., 2013, ApJ , 77, 133 CrossRef Search ADS   Larson R. B., 1981, MNRAS , 194, 809 CrossRef Search ADS   Lee A. T., Cunningham A. J., McKee C. F., Klein R. I., 2014, ApJ , 783, 50 CrossRef Search ADS   Lee E. J., Miville-Deschênes M.-A., Murray N. W., 2016, ApJ , 833, 229 CrossRef Search ADS   Lee Y.-N., Hennebelle P., 2016a, A&A , 591, A30 CrossRef Search ADS   Lee Y.-N., Hennebelle P., 2016b, A&A , 591, A31 CrossRef Search ADS   Li P. S., Martin D. F., Klein R. I., McKee C. F., 2012, ApJ , 745, 139 CrossRef Search ADS   Li P. S., Myers A., McKee C. F., 2012, ApJ , 760, 33 CrossRef Search ADS   Li P. S., McKee C. F., Klein R. I., 2015, MNRAS , 452, 2500 CrossRef Search ADS   Li P. S., McKee C. F., Klein R. I., 2018, MNRAS , 473, 4220 CrossRef Search ADS   Li Z.-Y., Nakamura F., 2006, ApJ , 640, L187 CrossRef Search ADS   Li Z.-Y., Banerjee R., Pudritz R. E., Jørgensen J. K., Shang H., Krasnopolsky R., Maury A., 2014, in Beuthe H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI . Univ. Arizona Press, Tucson, AZ, p. 173 Londrillo P., del Zanna L., 2004, J. Comput. Phys. , 195, 17 CrossRef Search ADS   Mac Low M.-M., 1999, ApJ , 524, 169 CrossRef Search ADS   Mason B. D., Hartkopf W. I., Gies D. R., Henry T. J., Helsel J. W., 2009, AJ , 137, 3358 CrossRef Search ADS   Matzner C. D., 2002, ApJ , 556, 302 CrossRef Search ADS   Matzner C. D., 2007, ApJ , 659, 1394 CrossRef Search ADS   McKee C. F., Offner S. S. R., 2010, ApJ , 716, 167 CrossRef Search ADS   McKee C. F., Ostriker E. C., 2007, ARA&A , 45, 565 CrossRef Search ADS   McKee C. F., Tan J. C., 2003, ApJ , 585, 850 CrossRef Search ADS   McKee C. F., Li P. S., Klein R. I., 2010, ApJ , 720, 1612 CrossRef Search ADS   Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS , 198, 7 CrossRef Search ADS   Miyoshi T., Kusano K., 2005, J. Comput. Phys. , 208, 315 CrossRef Search ADS   Mouschovias T. C., Spitzer L., Jr, 1976, ApJ , 210, 326 CrossRef Search ADS   Murray N. W., 2011, ApJ , 729, 133 CrossRef Search ADS   Murray D. W., Chang P., Murray N. W., Pittman J., 2015, MNRAS , 465, 1316 CrossRef Search ADS   Myers A. T., Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ , 735, 49 CrossRef Search ADS   Myers A. T., McKee C. F., Cunningham A. J., Klein R. I., Krumholz M. R., 2013, ApJ , 766, 97 CrossRef Search ADS   Myers A. T., Klein R. I., Krumholz M. R., McKee C. F., 2014, MNRAS , 439, 3420 CrossRef Search ADS   Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ , 703, 131 CrossRef Search ADS   Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI . Univ. Arizona Press, Tucson, AZ, p. 53 Ostriker E., Gammie C. F., Stone J. M., 1999, ApJ , 513, 259 CrossRef Search ADS   Padoan P., Nordlund Å., Kritsuk A. G., Norman M. L., Li P. S., 2007, ApJ , 661, 972 CrossRef Search ADS   Padoan P., Haugbølle T., Nordlund, Å., 2012, ApJ , 759, L27 CrossRef Search ADS   Padoan P., Pan L., Haugbølle T., Nordlund Å., 2016, ApJ , 822, 11 CrossRef Search ADS   Padoan P., Haugbølle T., Nordlund Å., Frimann S., 2017, ApJ , 840, 48 CrossRef Search ADS   Pan L., Padoan P., Haugbølle T., Nordlund Å., 2016, ApJ , 825, 30 CrossRef Search ADS   Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 2011, ApJ , 729, 72 CrossRef Search ADS   Peters T., Klaassen P. D., Mac Low M.-M., Schrön M., Federrath C., Smith M. D., Klessen R. S., 2014, ApJ , 788, 14 CrossRef Search ADS   Preibisch T., Balega Y., Hofmann K.-H., Weigelt G., Zinnecker H., 1999, New Astron. , 4, 531 CrossRef Search ADS   Price D. J., Bate M. R., 2008, MNRAS , 385, 1820 CrossRef Search ADS   Price D. J., Bate M. R., 2009, MNRAS , 398, 33 CrossRef Search ADS   Raghavan D. et al.  , 2010, ApJS , 190, 1 CrossRef Search ADS   Robertson B., Goldreich P., 2012, ApJ , 750, L31 CrossRef Search ADS   Rosen A. L., Krumholz M. R., McKee C. F., Klein R. I., 2016, MNRAS , 463, 2553 CrossRef Search ADS   Rosen A. L., Krumholz M. R., Oishi J. S., Lee A. T., Klein R. I., 2017, J. Comp. Phys. , 330, 924 CrossRef Search ADS   Salim D. M., Federrath C., Kewley L., 2015, ApJ , 806, L36 CrossRef Search ADS   Salpeter E. E., 1955, ApJ , 121, 161 CrossRef Search ADS   Santos-Lima R., de Gouveia Dal Pino E. M., Lazarian A., 2012, ApJ , 747, 21 CrossRef Search ADS   Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A , 410, 611 CrossRef Search ADS   Stone J. M., Ostriker E. C., Gammie C. F., 1998, ApJ , 508, L99 CrossRef Search ADS   Sur S., Federrath C., Schleicher D. R. G., Banerjee R., Klessen R. S., 2012, MNRAS , 423, 3148 CrossRef Search ADS   Tobin J. J. et al.  , 2016, ApJ , 818, 73 CrossRef Search ADS   Tricco T. S., Price D. J., Federrath C., 2016, MNRAS , 461, 1260 CrossRef Search ADS   Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., II, Howell L. H., Greenough J. A., 1997, ApJ , 489, L179 CrossRef Search ADS   Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., II, Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ , 495, 821 CrossRef Search ADS   Usero A. et al.  , 2015, AJ , 150, 115 CrossRef Search ADS   Wang P., Li Z.-Y., Abel T., Nakamura F., 2010, ApJ , 709, 27 CrossRef Search ADS   Wurster J., Price D. J., Bate M. R., 2017, MNRAS , 466, 1788 CrossRef Search ADS   Zhao B., Li Z.-Y., Nakamura F., Krasnopolsky R., Shang H., 2011, ApJ , 742, 10 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# The effects of magnetic fields and protostellar feedback on low-mass cluster formation

, Volume 476 (1) – May 1, 2018
22 pages

/lp/ou_press/the-effects-of-magnetic-fields-and-protostellar-feedback-on-low-mass-Vo0qVrL5ja
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty154
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present a large suite of simulations of the formation of low-mass star clusters. Our simulations include an extensive set of physical processes – magnetohydrodynamics, radiative transfer, and protostellar outflows – and span a wide range of virial parameters and magnetic field strengths. Comparing the outcomes of our simulations to observations, we find that simulations remaining close to virial balance throughout their history produce star formation efficiencies and initial mass function (IMF) peaks that are stable in time and in reasonable agreement with observations. Our results indicate that small-scale dissipation effects near the protostellar surface provide a feedback loop for stabilizing the star formation efficiency. This is true regardless of whether the balance is maintained by input of energy from large-scale forcing or by strong magnetic fields that inhibit collapse. In contrast, simulations that leave virial balance and undergo runaway collapse form stars too efficiently and produce an IMF that becomes increasingly top heavy with time. In all cases, we find that the competition between magnetic flux advection towards the protostar and outward advection due to magnetic interchange instabilities, and the competition between turbulent amplification and reconnection close to newly formed protostars renders the local magnetic field structure insensitive to the strength of the large-scale field, ensuring that radiation is always more important than magnetic support in setting the fragmentation scale and thus the IMF peak mass. The statistics of multiple stellar systems are similarly insensitive to variations in the initial conditions and generally agree with observations within the range of statistical uncertainty. stars: low-mass, stars: protostars, ISM: jets and outflows, ISM: magnetic fields 1 INTRODUCTION The formation of stars through the collapse of pre-stellar molecular clouds is a rich process, mediated by the coupled effects of gravity, magnetic fields, turbulence, and feedback from protostellar evolution in the form of outflows and radiation. To date, numerical models of star cluster formation have mostly probed the coupling of at most two of these effects. While most contemporary protostellar cluster models include turbulence and self-gravity, the remaining effects have received a more piecemeal treatment. Models that include magnetic effects and protostellar outflow feedback have typically ignored radiative transfer (Li & Nakamura 2006; Wang et al. 2010; Federrath 2015), while simulations that include radiative transfer with protostellar outflows have neglected magnetic fields (Hansen et al. 2012; Krumholz, Klein & McKee 2012). Yet other works have included magnetic fields and radiative transfer but not protostellar outflows (Price & Bate 2008, 2009; Peters et al. 2011). Two exceptions that includes all three effects are the work of Myers et al. (2014) and Li, Klein & McKee (2018), who consider the formation of star clusters consistent with observations of regions of high-mass star formation with mean column densities Σ = 1 g cm−3. In this work, we perform an analogous study of regions of low-mass star formation, characterized by a mean column density Σ = 0.1 g cm−3. Under these conditions the coupling of radiation from young stars to the infalling gas is weaker (Krumholz & McKee 2008; Krumholz et al. 2010; Cunningham et al. 2011), and it is possible that magnetic and protostellar outflow effects are more prominent. A consistent and comprehensive theory of star formation must simultaneously confront a number of observational constraints, as summarized in the review by Krumholz (2014). The relative strengths of the magnetic field and gravity in a clump of molecular gas are measured by the mass-to-flux ratio. Statistical surveys of Zeeman splitting measurements indicate that this ratio is 2–3 times larger than the critical value at which the magnetic field is strong enough to prevent gravitational collapse (Falgarone et al. 2008; Crutcher et al. 2010). Recent modelling that allows for the difference between the theoretical mass-to-flux ratio, which is area weighted, and the observed one, which is density weighted, show that the data summarized by Crutcher et al. (2010) are consistent with a value for this ratio that is three times the critical value (Li, McKee & Klein 2015). Next, the collapse of a molecular cloud must give rise to a distribution of stellar masses consistent with the observed stellar initial mass function (IMF). The observed IMF is well represented by a power law for stars more massive than 1 M⊙ (Salpeter 1955) and a lognormal distribution peaked at 0.2 M⊙ for lower mass sources (Chabrier 2005). Simulations suggest that magnetic fields influence the process by which the IMF is established by adding support to otherwise super-Jeans scale masses, thereby suppressing fragmentation and influencing the distribution of stellar masses (Padoan et al. 2007; Commerçon, Hennebelle & Henning 2011; Hennebelle et al. 2011; Federrath & Klessen 2012; Myers et al. 2013). Radiative heating from protostellar sources similarly acts to increase the local Jeans mass that characterizes the gas in the protostars’ vicinity, further reducing fragmentation and inhibiting the formation of lower mass protostars (Bate 2009; Offner et al. 2009; Myers et al. 2011; Krumholz, Klein & McKee 2011, 2012; Krumholz et al. 2016). Outflows also indirectly play a role in mediating the IMF by reducing the rate of accretion on to stars, which in turn reduces the rate of radiative heating. In spite of the complexity of the coupling among magnetic, radiative and outflow feedback, the shape of the IMF is remarkably insensitive to the star-forming environment, particularly for the majority of stars that are heavier than brown dwarfs but not so massive as to be exceptionally rare, where observational statistics are best constrained (Offner et al. 2014). Not only must the complex physical processes associated with star formation collude to fragment molecular clouds into a distribution of masses that is consistent with the observed IMF, but the star formation rate must also agree with observation. Star formation is remarkably inefficient. Absent turbulent or magnetic support, a molecular cloud should collapse entirely into stars on a free-fall time-scale. In fact, the star formation efficiency per free-fall time εff, defined as the fraction of gas mass that is converted to stars per free-fall time, was found to lie between 0.001 and 0.09 by Krumholz & Tan (2007); a much larger and more recent compilation of data gives εff = 0.015 with a scatter of ∼0.5 dex for both Galactic and extragalactic regions (Krumholz, Dekel & McKee 2012), and subsequent authors have obtained similarly low values (e.g. García-Burillo et al. 2012; Evans, Heiderman & Vutisalchavakul 2014; Salim, Federrath & Kewley 2015; Usero et al. 2015; Heyer et al. 2016).1 Turbulence provides a significant source of support against rapid collapse. Gas velocity dispersions inferred from observed molecular cloud line width–size relations (Larson 1981) are consistent with the inertial range scaling of both supersonic hydrodynamic and magnetohydrodynamical turbulence (Gammie & Ostriker 1996; Ostriker, Gammie & Stone 1999) and the magnitude of the velocity dispersions is consistent with a near virial equipartition between turbulent support and gravitational contraction (McKee, Li & Klein 2010). Numerical studies of magnetized turbulence show that turbulence at virial levels, in conjunction with protostellar outflows, can produce values of εff close to the observed, low values (Li & Nakamura 2006; Wang et al. 2010; Padoan, Haugbølle & Nordlund 2012; Krumholz, Klein & McKee 2012; Myers et al. 2014; Federrath 2015). Magnetohydrodynamical turbulence, however, undergoes exponential decay via dissipation in shocks on a turbulent crossing time (Stone, Ostriker & Gammie 1998; Mac Low 1999), which is comparable to the free-fall time in a virialized cloud. Gravitational collapse can drive turbulence, but only at the cost of increasing the density and thus lowering the free-fall time (Robertson & Goldreich 2012; Murray et al. 2015). As a result, in the absence of an external or internal energy source to drive the turbulence, hydrodynamical models fall into a state of rapid, near free-fall collapse (Hansen et al. 2012; Krumholz et al. 2012). A number of forcing mechanisms applicable to low-mass star clusters have been proposed, including accretion from larger scales (e.g. Klessen & Hennebelle 2010; Goldbaum et al. 2011; Lee & Hennebelle 2016a,b), internal forcing by protostellar outflows (e.g. Li & Nakamura 2006; Matzner 2007; Wang et al. 2010; Cunningham et al. 2011; Peters et al. 2014; Federrath 2015), and forcing by external H ii regions (e.g. Matzner 2002; Krumholz, Matzner & McKee 2006; Goldbaum et al. 2011) and supernova explosions (e.g. Padoan et al. 2016, 2017; Pan et al. 2016). In this paper, we present a series of simulations of the formation of low-mass stellar systems, which we compare to the various observational constraints on the star formation process we have outlined above. Our goal is to search for models that are capable of simultaneously reproducing the observed IMF and star formation efficiency, determine the relationship between these two constraints, and identify the physical mechanisms that are responsible for satisfying them. In order to accomplish this, we carry out a suite of simulations All our simulations use radiative transfer and protostellar radiation feedback, since these are likely key to setting the IMF. In order to explore the influence of outflows, we toggle these on and off. We also vary the initial magnetization of our regions, and consider both models with turbulence that is continuously driven to maintain virial equilibrium and those where no artificial driving is included. We describe our numerical method and simulation setup in Section 2, we describe the results in Section 3, we explore the question of magnetic fields in more detail in Section 4, and we summarize our findings in Section 5. 2 SIMULATION SETUP 2.1 Initial conditions We perform all our simulations using the orion2 code . Our setup shares several initial parameters with Offner et al. (2009) and Hansen et al. (2012). The initial state consists of a solar metallicity molecular gas with a mean molecular weight 2.33mp and a temperature Tg = 10 K, corresponding to an initial sound speed cs = 0.19 km s−1. The length of the cubical domain is L = 0.65 pc and the average density is $$\bar{\rho }=4.46\times 10^{-20}\ {\rm g\ cm}^{-3}$$ ($$n_{\rm H_2} = 1.15\times 10^4$$ cm−3) corresponding to a total mass of M = 185 M⊙. The initial mean state corresponds to a minimum stable Jeans length λJ = 0.20 pc, and a Jeans mass MJ = 2.7 M⊙. Our models are carried forward with periodic boundary conditions for the gravity and magnetohydrodynamics (MHD). We use Marshak boundary conditions with an external radiation field corresponding to a 10 K isotropic blackbody for the radiation field. The ability of a magnetized cloud to resist gravitational collapse is measured by the magnetic critical mass (Mouschovias & Spitzer 1976),   $$M_\Phi = c_\Phi \frac{\Phi }{\sqrt{G}},$$ (1)where Φ is the total magnetic flux threading the cloud and $$c_\Phi$$ is a dimensionless constant that depends weakly on the equilibrium shape of the cloud in the absence of gravity. A cloud of mass less than $$M_\Phi$$ cannot collapse, whereas a cloud with a greater mass cannot be supported against collapse by magnetic fields alone. Because the value of $$c_\Phi$$ is only weakly sensitive to the cloud geometry (McKee et al. 2010), we adopt the value appropriate for a slab threaded by a uniform, perpendicular field, $$c_\Phi =1/2\pi$$, consistent with the convention chosen in several contemporary theoretical (Li et al. 2015) and observational works (Crutcher et al. 2010). We then define the normalized mass-to-flux ratio   $$\mu _\Phi = \frac{M}{M_\Phi } = 2\pi \sqrt{G}\left(\frac{M}{\Phi }\right)=2\pi \sqrt{G}\left(\frac{\Sigma }{\bar{B}}\right),$$ (2)where $$\bar{B}\equiv \Phi /S$$ is the magnitude of the mean magnetic field, Σ is the mass surface density, and S the cloud's cross-sectional area in the plane orthogonal to the magnetic field direction. Since mass and magnetic flux are conserved for the entire simulation box, the value of $$\mu _\Phi$$ for the box is constant and can be used to specify the initial mean magnetic field. We define initial fields large enough to render the system close to magnetically critical, $$\mu _\Phi \approx 1$$, ‘strong,’ initial fields that make the system moderately supercritical, $$\mu _\Phi \approx 2\text{--}3$$ ‘moderate,’ and initial fields corresponding to $$\mu _\Phi \gg 1$$ ‘weak.’ It is convenient to express the strength of the magnetic field relative to the kinetic energy in terms of the Alfvèn Mach number,   $${\cal M}_{\rm A}=\frac{v_{\rm rms}}{v_{\rm A}},$$ (3)where $$v_{\rm A}^2\equiv B_{\rm rms}^2/(4\pi \bar{\rho })$$ is the Alfvèn velocity, Brms is the root-mean-square (rms) volume-weighted field strength, and vrms is the rms mass-weighted, 3D velocity dispersion. Relative to thermal pressure, the strength of the magnetic field can be expressed in terms of the plasma-β parameter   $$\beta = \frac{8\pi \bar{\rho }c^2_s}{B_{\rm rms}^2} = 2\left(\frac{{\cal M}_{\rm A}}{{\cal M}}\right)^2.$$ (4) In this work, we carry out simulations with several different mean magnetic field strengths. Our models begin with an initially uniform magnetic field B0 in the $$\hat{z}$$-direction. We consider cases with mass-to-flux ratios $$\mu _\Phi =1.56\textrm{,}2.17\textrm{,}23.1\textrm{,and} \infty$$. The corresponding initial Alfvèn Mach numbers are $${\cal M}_{\rm A,0}=1.0\textrm{,}1.4\textrm{,}15\textrm{,and} \infty$$, and the initial plasma β's are β0 = 0.046, 0.089, 10, and ∞. We note that the recent models by Federrath (2015) consider a somewhat more weakly magnetized model ($$\mu _\Phi =3.8$$) than our two most strongly magnetized cases, and the model by Wang et al. (2010) has an initial cloud with $$\mu _\Phi =1.4$$, comparable to our strongest field model. Turbulent conditions for starting the cluster simulation are obtained by applying the driving recipe of Mac Low (1999). The procedure consists of applying a force consisting of a purely solenoidal, uniform spectrum of modes spanning wave numbers $$\boldsymbol {k}$$ such that $$1 \le |\boldsymbol {k}|L/2\pi \le 2$$ to a uniform initial gas. This turbulent driving force is intended to mock up the effect of sources of mechanical energy on scales larger than the simulation box, via one of the numerous possible mechanisms discusses in Introduction. The observed velocity dispersion to size relation can be met by models of compressive turbulence that span a range of relative contributions from solenoidal and compressive forcing (Federrath, Klessen & Schmidt 2009). Recent simulations of supernova-driven turbulence (Padoan et al. 2016; Pan et al. 2016) show that even when the initial driving is purely compressive, the turbulence becomes primarily solenoidal in the inhomogeneous interstellar medium (ISM); indeed, these authors conclude that external supernova-driven turbulence is best characterized as 84 per cent solenoidal after mediation through the ISM to adjacent star-forming regions. The work of Federrath et al. (2010) indicates that the statistical properties of turbulent flow are not sensitive to ≲ 20 per cent compressive contributions to compressive driving. This suggests that our purely solenoidal forcing is a plausible approximation for the driving of turbulence in a molecular cloud. The turbulent driving force is scaled to maintain a roughly constant mass-weighted rms thermal Mach number of $${\cal M}= v_{\rm rms}/c_{\rm s} \equiv \sqrt{3}\sigma _{\rm v}/c_{\rm s}= 6.6$$, where σv is the 1D velocity dispersion. The virial parameter,   $$\alpha _{\rm vir}\equiv \frac{5\sigma _{\rm v}^2 R}{GM},$$ (5)measures the balance between kinetic and gravitational energies; a spherical cloud of uniform density has αvir = 1 in virial equilibrium. With R = L/2 and our chosen value of $${\cal M}$$, the computational box has αvir = 1.05. The gas is initially evolved under the action of the driving force without gravity according to the equations of ideal, isothermal MHD for two crossing times,   $$t_{\rm cross}=\frac{L}{v_{\rm rms}}=0.51\ \textrm{Myr},$$ (6)on a uniform mesh with 5123 zones with periodic boundary conditions. This resolution has been shown by Li et al. (2012) to be sufficient to capture a well-resolved inertial cascade down to scales of  L/30 = 0.02 pc, an order of magnitude smaller than the initial Jeans scale of the turbulent cloud. Consequently, the turbulent flows on the scales that give rise to the initial fragmentation of the simulated cluster are well resolved. Henceforth, we denote the state of the turbulent gas after evolving for two crossing times as the time t = 0.2 We take these conditions as the initial state of the turbulent cluster simulation, which includes a more comprehensive set of physics that we describe in Section 2.3. 2.2 Defining the free-fall time and star formation rate Because our simulations take place in a periodic box that represents only a portion of a cloud, we must give some care to how we go about defining the free-fall time and the efficiency of star formation. Absent magnetic or thermal support to oppose self-gravity, the free-fall time of gas structures of characteristic density ρ is   $$t_{\rm ff} = \sqrt{\frac{3\pi }{32G \rho }}.$$ (7)The coefficient $$\sqrt{3\pi /32}$$ applies to spherical geometries, but by convention we retain it for all geometries. The mean cloud density $$\bar{\rho }$$ corresponds to a free-fall time tff, IC = 0.315 Myr, where we use the subscript IC to refer to quantities defined for the uniform-density state before driving the initial turbulence to steady state. We can measure star formation rates relative to this free-fall time via the usual dimensionless star formation rate parameter   $$\epsilon _{\rm ff,IC} = \frac{\dot{M}_* t_{\rm ff,IC}}{M}.$$ (8) While εff, IC provides one measure of the star formation rate, it can be somewhat misleading. As we drive the turbulence to statistical steady state, the gas density distribution becomes non-uniform, and a majority of the mass ends up residing at densities much higher than $$\bar{\rho }$$. This leads to two biases. First, because most of the mass is at higher density, we might expect εff, IC to end up relatively large even if star formation in the dense gas is inefficient, simply because the dense gas has a much shorter free-fall time. Second, because we are using a periodic box to handle gravity, much of the low-density gas is not self-gravitating at all, and thus would not collapse even with efficient star formation. Including this mass in a calculation of εff would lead to an artificially low value. To counteract these biases, we define a new reference density that provides a better description of the gas after the turbulence is fully developed. Let ρmed be the mass-weighted median density – i.e. the density at which half the mass is denser and half less dense – at the end of the turbulent driving phase (t = 0). Then, we define ρ>med as the average density of the regions with ρ > ρmed at that time. In effect, ρ>med gives an estimate of the mean density of the gas that is actually self-gravitating prior to any star formation. We denote the associated free-fall time as tff, >med. The dimensionless star formation rate is defined as   $$\epsilon _{\rm ff,>med} = \frac{\dot{M}_* t_{\rm ff,>med}}{M_{\rm >med}}.$$ (9)Since M = 2M>med, this can also be expressed as   $$\epsilon _{\rm ff,>med} = \frac{2 \dot{M}_* t_{\rm ff,>med}}{M}.$$ (10) While εff, >med is a natural physical measure of star formation efficiency, it cannot readily be compared to observations. This is both because our simulation at t = 0 represents a somewhat artificial state of driven turbulence but no gravitational effects, and because observations cannot usually measure gas densities directly in any event. Instead, volume densities must be inferred indirectly from molecular line critical densities, or, more commonly for Galactic observations, using column densities that can be measured almost directly. Following Krumholz et al. (2012), we define the observationally inferred characteristic density by taking projected column densities and extracting those regions denser than 0.048 g cm−2 = 230  M⊙ pc−2, corresponding to those regions with K-band extinction AK > 0.8 mag, a threshold used by a number of previous authors (e.g. Lada, Lombardi & Alves 2010; Lada et al. 2012, 2013). We compute the projected area of these regions S>0.8, and from this determine the projected mass enclose M>0.8 and the effective radius $$R_{\rm >0.8} = \sqrt{S_{\rm >0.8}/\pi }$$. We then define the characteristic observed density of the projection as $$\rho _{\rm obs, proj} = 3 M_{\rm >0.8} / (4 \pi R_{\rm >0.8}^3)$$, and take the mean value from projections in the three cardinal directions to be the characteristic density ρ>0.8 for each simulation. We compute this quantity at the final time in each simulation (in contrast to ρ>med, which is calculated at t = 0, when gravity is turned on), and we denote the corresponding free-fall time as tff, >0.8. The star formation rate per tff, >0.8 is then   $$\epsilon _{\rm ff,>0.8} = \frac{\dot{M}_* t_{\rm ff,>0.8}}{M_{\rm >0.8}},$$ (11)and this is the quantity that is most natural to compare to observed star formation rates. In Table 1, we summarize the two characteristic densities (ρ>med, and ρ>0.8) that we have defined for all our simulations. The former is typically 7–9 times greater than the mean density $$\bar{\rho }=4.46 \times 10^{-20}\, {\rm g\ cm}^{-3}$$, due to turbulent compression. The observationally defined ρ0.8 ranges from $$2\ \textrm{to}\ 7 \times \bar{\rho }$$ due to the influence of the combined effects of turbulent compression, gravitational collapse, and stellar feedback. We also report the mean values of the Alfvèn Mach number and plasma β at the end of the driving phase and at the final time in the simulation, to make it clear how these evolve with time as well. This evolution is a result of field amplification by turbulence and gravitational collapse. For consistency, we use the same notation for these other quantities as for the density and free-fall time, i.e. quantities measured in the initial, uniform state have a subscript ‘IC’ (for initial condition), quantities averaged over gas denser than the median density at the end of the driving phase have subscript ‘>med’, and those defined based on averages within the 0.8 mag contour are defined by the subscript ‘>0.8’. Table 1. Model parameters and initial conditions. The outflows and driving columns indicate whether we included protostellar outflows, and whether we continued driving the turbulence after gravity was turned on. For the remaining parameters, the subscript ‘IC’ denotes the state in the uniform initial condition, the subscript ‘>med’ denotes quantities averaged over gas denser than the median density at the instant when self-gravity is switched on, and the subscript ‘>0.8’ denotes averages over gas with a K-band extinction AK > 0.8 mag at the termination of the model. See the text for the details of how the averaging regions are defined. The mean density in all cases is 4.46 × 10−20 g cm−3.   Physics included  IC  >med  >0.8  μΦ  Outflows  Driving  βIC  $$\mathcal {M}_{A,\rm IC}$$  ρ>med (g cm−3)  β>med  $$\mathcal {M}_{\rm >med}$$  ρ>0.8 (g cm−3)  β>0.8  $${\cal M}_{\rm A,>0.8}$$  1.56  Y  N  0.046  1.0  2.36 × 10−19  0.041  0.95  3.21 × 10−19  0.039  1.09  1.56  N  N  0.046  1.0  2.36 × 10−19  0.041  0.95  1.45 × 10−19  0.038  1.14  1.56  Y  Y  0.046  1.0  2.36 × 10−19  0.041  0.95  9.76 × 10−20  0.037  1.05  2.17  Y  N  0.089  1.39  2.18 × 10−19  0.057  1.14  1.53 × 10−19  0.074  1.92  2.17  Y  Y  0.089  1.39  2.18 × 10−19  0.057  1.14  2.35 × 10−19  0.067  1.70  23.1  Y  N  10.0  14.8  1.90 × 10−19  0.46  3.14  1.62 × 10−19  0.41  6.67  23.1  Y  Y  10.0  14.8  1.90 × 10−19  0.46  3.14  1.07 × 10−19  0.75  4.45  ∞  Y  N  ∞  ∞  2.66 × 10−19  ∞  ∞  1.33 × 10−19  ∞  ∞  ∞  Y  Y  ∞  ∞  2.66 × 10−19  ∞  ∞  9.24 × 10−20  ∞  ∞    Physics included  IC  >med  >0.8  μΦ  Outflows  Driving  βIC  $$\mathcal {M}_{A,\rm IC}$$  ρ>med (g cm−3)  β>med  $$\mathcal {M}_{\rm >med}$$  ρ>0.8 (g cm−3)  β>0.8  $${\cal M}_{\rm A,>0.8}$$  1.56  Y  N  0.046  1.0  2.36 × 10−19  0.041  0.95  3.21 × 10−19  0.039  1.09  1.56  N  N  0.046  1.0  2.36 × 10−19  0.041  0.95  1.45 × 10−19  0.038  1.14  1.56  Y  Y  0.046  1.0  2.36 × 10−19  0.041  0.95  9.76 × 10−20  0.037  1.05  2.17  Y  N  0.089  1.39  2.18 × 10−19  0.057  1.14  1.53 × 10−19  0.074  1.92  2.17  Y  Y  0.089  1.39  2.18 × 10−19  0.057  1.14  2.35 × 10−19  0.067  1.70  23.1  Y  N  10.0  14.8  1.90 × 10−19  0.46  3.14  1.62 × 10−19  0.41  6.67  23.1  Y  Y  10.0  14.8  1.90 × 10−19  0.46  3.14  1.07 × 10−19  0.75  4.45  ∞  Y  N  ∞  ∞  2.66 × 10−19  ∞  ∞  1.33 × 10−19  ∞  ∞  ∞  Y  Y  ∞  ∞  2.66 × 10−19  ∞  ∞  9.24 × 10−20  ∞  ∞  View Large 2.3 Physics included and numerical method Our protostellar cluster simulations are performed using the orion2 adaptive mesh refinement (AMR) code (Li et al. 2012) beginning from the initial state generated following the procedure described in the previous section. The code solves the equations of ideal MHD using the scheme of Mignone et al. (2012), along with coupled self-gravity (Truelove et al. 1998; Klein et al. 1999) and radiation transfer (Krumholz et al. 2007) in the two-temperature, mixed-frame, grey, flux-limited diffusion approximation. The exact set of equations solved is the same as those given in Myers et al. (2014). orion2 contains many options for the numerical advection scheme, but we have found the best stability with the Harten-Lax-Van Leer-Discontinuities Riemann solver (Miyoshi & Kusano 2005) coupled to the constrained transport magnetic flux advection of Londrillo & del Zanna (2004), particularly for magnetically dominated flows. Our radiative transfer calculations use the frequency-integrated grey dust opacities from the iron-normal, composite aggregates model of Semenov et al. (2003). For each magnetic field strength parameter considered in this study, $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞, we present two models – one with decaying turbulence where the large-scale driving forces described in the previous section are turned off, and one where the turbulent forcing continues with a constant rate of energy injection that balances the rate of turbulent decay as estimated for a given mean magnetic field strength and turbulent Mach number by Mac Low (1999). In the case without turbulent forcing, self-gravity is switched on and turbulent forcing is switched off simultaneously. Our cluster evolution models neglect non-ideal MHD effects. These effects have been shown to impact turbulent density profiles on small scales Li, Myers & McKee (2012) and therefore could influence the fragmentation dynamics of strongly magnetized cores. Recent studies on the impact of non-ideal MHD effects on the evolution and fragmentation of self-gravitating cores indicate that the evolution of supercritical cores is most strongly influenced by their initial conditions and that the impact of non-ideal effects on their subsequent evolution and fragmentation are small by comparison Wurster, Price & Bate (2017). Consequently, our models should reliably capture the gross impacts of varying large-scale magnetic field strength on the cloud fragmentation and the statistics of its collapse. There are two limitations of our radiative transfer method to which we should call attention. First, we assume that the gas and dust temperatures are the same. In the dense regions simulated by Myers et al. (2014), this is a very good assumption, because the dust and gas temperatures become nearly identical at densities above ∼104–105 cm−3 (e.g. Goldsmith 2001). The assumption of strong coupling is still valid for most of the gas in our simulation domain that is actively participating in star formation (cf. ρ>, med in Table 1), but it is not valid for the lower density, non-self-gravitating regions. In these parts of the flow, we are somewhat overestimating the dust cooling rate for the gas. The second approximation we make is to assume that the absorption opacity is the same as the emission opacity, which is equivalent to assuming that the dust and radiation temperatures are the same. This approximation is clearly invalid very close to the star, where dust is exposed to direct stellar radiation and reprocesses the radiation into the infrared. However, the dust opacity to direct stellar radiation is so high that all this reprocessing occurs in a thin zone at the dust destruction front, which for low-mass star formation is confined to ≲1 au from the star (except perhaps over the small range in solid angle evacuated by outflows). Such small structures are unresolved in our simulations. Our approximation would only become problematic for stars larger than ∼20  M⊙, massive enough that radiation pressure can inflate bubbles that push the dust destruction front out to thousands of au (e.g. Rosen et al. 2016). While orion2 does include a hybrid radiative transfer method to handle this situation (Rosen et al. 2017), we do not use it here because our simulations do not form any stars massive enough to require it. Once the radiation field is dominated by reprocessed radiation, the dust and radiation temperatures will remain close so long as the dust is opaque. Outside the dust photosphere, the emission opacity declines while the absorption opacity remains constant, but this affects only the low-frequency part of the emission spectrum, which has only a small fraction of the energy (Chakrabarti & McKee 2005).3 Our models include several source terms to capture the effects of protostellar feedback that originate on smaller length scales than those resolved in the model. We initialize a sink particle in any zone on the finest AMR level, and only on the finest level, that becomes dense enough to reach a local Jeans number   $$J=\sqrt{\frac{G \rho \Delta x^2}{\pi c_{\rm s}^2}} > \frac{1}{4}.$$ (12)Sink particles then evolve and interact with the cluster through gravity according to the methodology of Krumholz, McKee & Klein (2004), updated to include the effects of magnetic fields on the rate of gas accretion on to sink particles (see the appendix of Lee et al. 2014). The key assumption in our treatment is that the point mass accretes mass, but not flux. This assumption is motivated by observations that show that the magnetic flux in young stellar objects is orders of magnitude less than that in the gas that formed these objects, implying that flux accretion is very inefficient (McKee & Ostriker 2007). The sink particles couple to their surrounding gas gravitationally and through accretion and feedback source terms that operate within a radius of 4Δx on the finest AMR mesh level around each sink particle. The sink particles emit radiation according to the protostellar evolution model of Offner et al. (2009). Our models include the atomic line-cooling sources described in Cunningham et al. (2011) to treat strongly heated regions behind outflow shocks that are sufficiently fast ∼30 km s−1 to dissociate molecules. Feedback due to protostellar outflows is also included via momentum sources around each sink particle following the procedure described in Cunningham et al. (2011): the outflow mass ejection rate is set by the fraction of accreted gas that is ejected into the outflow fw and the outflow ejection speed vw. In this work, we make the same wind model parameter choices as Hansen et al. (2012) by setting fw = 0.3 and vKep = min (vKep,  60 km s−1), where vKep is the Kepler speed at the surface of the protostar. We note that these parameters differ from the parameters used in Cunningham et al. (2011) and Li et al. (2018) that set vw = vKep/3. Consequently, the models presented here impart relatively more outflow feedback during the earlier phases of pre-stellar evolution relative to Cunningham et al. (2011) and Li et al. (2018). Computational expediency motivates the cap on the wind velocity so that isolated massive protostars do not dominate the time-step of the cluster model. The AMR hierarchy is initialized on a 2563 base grid that we denote as level $$\mathcal {L}=0$$. Our initial state is taken from turbulent initial conditions that were established on a 5123 grid with no adaptivity. After t = 0, regions of the flow within 2 base-grid zones of an interface with a density jump Δρ/ρ > 0.5 on the base level are refined to $$\mathcal {L}=1$$, so that the 5123 resolution of the initial state is retained in these regions. On all AMR levels, refinement is triggered in regions in which the local Jeans scale is resolved by fewer than eight zones, J > 1/8 (Truelove et al. 1997). Note that this refinement criterion triggers refinement to the finest level, $$\mathcal {L}_{\rm MAX}$$, at a density one-fourth that required to trigger insertion of a sink particle. This ensures that strongly collapsing regions are refined to the finest level. Refinement is also triggered when there is a large jump in radiation energy density, |∇Er|Δx/Er > 0.125. A resolution sensitivity test was conducted for the case with $$\mu _\Phi =2.17$$ by first performing this run with $$\mathcal {L}_{\rm MAX}=4$$ and then repeating it without the density gradient refinement triggering on level $$\mathcal {L}_0$$ and with $$\mathcal {L}_{\rm MAX}=3$$. We found that the total mass accreted on to protostars out to t = 1.35tff in the two models agreed to within 3 per cent and that the mass distribution of protostars had a similar level of agreement. Our cluster models with decaying turbulence begin with $$\mathcal {L}_{\rm MAX}=4$$. The effective resolution of all collapsing and/or heated regions is therefore $$L/(256 \times 2^{\mathcal {L}_{\rm MAX}}) = 32$$ au. Given the good agreement between the runs with $$\mathcal {L}_{\rm max}=3$$ and 4, the $$\mu _\Phi =1.56$$ and 2.17 models with decaying turbulence were switched to $$\mathcal {L}_{\rm MAX}=3$$ with 65 au effective resolution at t = 1.7tff in order to reduce the run time. The models with turbulent forcing were also performed with $$\mathcal {L}_{\rm MAX}=3$$ and 65 au effective resolution, except for the $$\mu _\Phi =23.1$$ case which retained $$\mathcal {L}_{\rm MAX}=4$$ throughout the evolution of the model. We note that the limited resolution of these models results in a gain in computational expediency but comes at some cost to the fidelity of the model in treating the fragmentation of low-mass cores. The density threshold triggering sink formation at our background 10 K temperature is given by inverting equation (12) as 1.1 × 10−15–4.5 × 10−15 g cm−3 for 65 and 32 au effective resolution. The isothermal-to-adiabatic transition where compressional heating balances cooling occurs at a substantially higher density of ∼ 10− 14(T/10 K)6 g cm−3 (see chapter 16 of Krumholz 2015), and thus our resolution would not be sufficient to capture fragmentation in a simulation without radiative transfer. However, once stars form, radiative heating powered by accretion luminosity heats material at much lower densities than one would infer simply by balancing adiabatic compression against cooling. This substantially eases the resolution required to capture fragmentation, and we show below that our resolution is sufficient that in practice we essentially always resolve the region where gas transitions from the background temperature to elevated temperatures around most protostars in the simulations. None the less, we caution that our limited resolution may still cause us to miss fragmentation in rare regions with little radiative heating. 3 RESULTS In this section, we examine the results of the simulations summarized in Table 1. In the discussion and figure legends that follow we denote the various models by their initial mass-to-flux ratio; cases that include turbulent driving are labelled with the designation ‘Driven’. We have one model with outflows turned off, to which we refer as ‘no wind’; this model also is not driven. 3.1 Protostellar mass accretion Fig. 1 shows the temporal evolution of the total mass M* assembled into protostars for each model as a function of time. The left-hand panel shows the simulations with decaying turbulence and the right-hand panel shows the simulations with large-scale turbulent driving, a convention we shall follow throughout this paper. Fig. 2 shows the star formation efficiency per initial free-fall time εff, >med (see Section 2.2 for the precise definition) as a function of time. Fig. 3 shows the number of protostars formed in each model as a function of time. Figure 1. View largeDownload slide Total mass accreted for the models without (left) and with (right) turbulent forcing. Figure 1. View largeDownload slide Total mass accreted for the models without (left) and with (right) turbulent forcing. Figure 2. View largeDownload slide Instantaneous value of the dimensionless star formation rate εff, >med for the models without (left) and with (right) turbulent forcing (right). The data have been boxcar-smoothed over a time-scale of 0.05 tff, >med to remove short-time-scale variability. Figure 2. View largeDownload slide Instantaneous value of the dimensionless star formation rate εff, >med for the models without (left) and with (right) turbulent forcing (right). The data have been boxcar-smoothed over a time-scale of 0.05 tff, >med to remove short-time-scale variability. Figure 3. View largeDownload slide Total number of protostars formed for the models without (left) and with (right) turbulent forcing. Figure 3. View largeDownload slide Total number of protostars formed for the models without (left) and with (right) turbulent forcing. We also report, in Table 2, both εff, >med and the observed value of the star formation efficiency per free-fall time, εff, >0.8 (again see Section 2.1). The quantities reported in the table are time-weighted averages over the last 0.5tff, >med of the simulation. For most models, we run the simulation until the value of εff, >med stabilizes, and thus the values reported in Table 2 are the steady-state ones. There are two exceptions to this statement. First, the two weakest field models with decaying turbulence ($$\mu _\Phi =23.1$$ and ∞) do not appear to have well-converged values of εff, and instead appear to approach near-free-fall collapse. We run these models for less time than the others, as such efficient collapse is generally inconsistent with the observational constraints. For them, the values of εff reported in Table 2 are only lower limits, and we flag them as such in the table. The second exception are the models with $$\mu _\Phi =2.17$$. We are forced to halt these runs due to computational constraints, as the number of sink particles and highly refined zones eventually becomes such that exploring these models further is prohibitively expensive. While we would, absent consideration of computation constraints, prefer to continue the $$\mu _\Phi =2.17$$ runs further, we note that these models have begun to exhibit stabilized εff.0 by the termination time. Thus, we do not report these as upper limits in Table 2, but warn readers here that they are less secure than the values reported for the other runs. Table 2. Star formation efficiency per free-fall time in each run. The various values of εff are computed for different definitions of the free-fall time, as discussed in Section 2.2. The quantities reported are averaged over the last 0.5tff, >med of each simulation. A value preceded by ≳ indicates that εff was still rising strongly and the end of the simulation, and the values listed are therefore lower limits. μΦ  Outflows  Driving  εff, IC(per cent)  εff, >med(per cent)  εff, >0.8(per cent)  $$\mu _\Phi = 1.56$$  Y  N  12  8.2  6.0  $$\mu _\Phi = 1.56$$  N  N  23  15.8  19  $$\mu _\Phi = 1.56$$  Y  Y  4.8  3.2  4.2  $$\mu _\Phi = 2.17$$  Y  N  26  18.8  19  $$\mu _\Phi = 2.17$$  Y  Y  15  11.0  8.7  $$\mu _\Phi = 23.1$$  Y  N  ≳ 30  ≳ 24  ≳ 22  $$\mu _\Phi = 23.1$$  Y  Y  8.2  6.4  6.6  $$\mu _\Phi = \infty$$  Y  N  ≳ 50  ≳ 34  ≳ 37  $$\mu _\Phi = \infty$$  Y  Y  6.0  4.0  5.2  μΦ  Outflows  Driving  εff, IC(per cent)  εff, >med(per cent)  εff, >0.8(per cent)  $$\mu _\Phi = 1.56$$  Y  N  12  8.2  6.0  $$\mu _\Phi = 1.56$$  N  N  23  15.8  19  $$\mu _\Phi = 1.56$$  Y  Y  4.8  3.2  4.2  $$\mu _\Phi = 2.17$$  Y  N  26  18.8  19  $$\mu _\Phi = 2.17$$  Y  Y  15  11.0  8.7  $$\mu _\Phi = 23.1$$  Y  N  ≳ 30  ≳ 24  ≳ 22  $$\mu _\Phi = 23.1$$  Y  Y  8.2  6.4  6.6  $$\mu _\Phi = \infty$$  Y  N  ≳ 50  ≳ 34  ≳ 37  $$\mu _\Phi = \infty$$  Y  Y  6.0  4.0  5.2  View Large For comparison, Krumholz et al. (2012) combine data from galactic clouds, nearby galaxies, and galaxies at high redshift, and find that all are consistent with a star formation efficiency εff, >0.8 = 0.015, with a scatter of ∼0.5 dex. As noted in Introduction, numerous other studies have confirmed this basic result. Comparing to the values shown in Fig. 2 and reported in Table 2, we see that the models with $$\mu _\Phi = 1.56$$ fall within the observational envelope if we include either driving or outflows, as do all the more weakly magnetized models with driven turbulence, although these tend to lie towards the top end of the observationally permitted range. We therefore find that the star formation rate can be kept low enough to be compatible with observations either if there are external mechanical energy inputs to maintain the turbulence or if there are no such influences, but the magnetic field lies close to the critical value. 3.2 Cluster morphology Fig. 4 shows the mean density along lines of sight in the $$\hat{\boldsymbol {x}}$$-direction at the termination time of each model, overlaid with a semi-transparent rendering of regions with |v| > 2 vrms. Recall that the initial magnetic field is in the $$\hat{\boldsymbol {z}}$$-direction, and is therefore upwardly oriented on the page in Fig. 4. We note that the mean density along the line of sight, which is proportional to the surface density, becomes increasingly filamentary as the strength of the magnetic field increases. In the cases where the magnetic field is the most dynamically significant, the outflows are oriented predominately perpendicularly to the orientation of the densest structures. The $$\mu _\Phi =1.56$$ case with decaying turbulence indicates a final cluster geometry that is most dominated by the magnetic field. In this case, the cluster is arranged mostly in a planar geometry normal to the $$\hat{\boldsymbol {z}}$$-direction. The magnetic fields in this case are strong enough to constrain large-scale collapse to be predominantly in the direction of the initial magnetic field, with less harassment from large-scale turbulent forcing than the equivalent model with turbulent forcing. Figure 4. View largeDownload slide Volume-weighted mean densities projected along the $$\hat{\boldsymbol {x}}$$-direction for the models with turbulent forcing (right) and with decaying turbulence (left). The column density can be recovered by multiplying the mean density shown by the computational domain length 0.65 pc; for reference, this means that a density of 10−18 g cm−3 corresponds to a column density of 2.0 g cm−2. The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ runs, from top to bottom. Regions with gas velocity greater than 2vrms are shaded in red with the opacity of the red regions increasing to a maximum for velocities exceeding 5vrms. White circles indicate the projected position of the protostellar sink particles. Figure 4. View largeDownload slide Volume-weighted mean densities projected along the $$\hat{\boldsymbol {x}}$$-direction for the models with turbulent forcing (right) and with decaying turbulence (left). The column density can be recovered by multiplying the mean density shown by the computational domain length 0.65 pc; for reference, this means that a density of 10−18 g cm−3 corresponds to a column density of 2.0 g cm−2. The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ runs, from top to bottom. Regions with gas velocity greater than 2vrms are shaded in red with the opacity of the red regions increasing to a maximum for velocities exceeding 5vrms. White circles indicate the projected position of the protostellar sink particles. 3.3 Protostellar mass functions We next examine the mass functions for the stars formed in our simulations. As a first step to this, we verify that these distributions have come close to being stable in time. To demonstrate this, in Fig. 5 we plot the temporal evolution of the 25th percentile, median and 75th percentile of the distribution of protostellar masses in each model as a function of time. Similar to the temporal evolution of εff, we see that the protostellar mass percentiles have mostly stabilized by the end of the runs, except for the cases of a weak initial field ($$\mu _\Phi =23.1$$ and ∞) with decaying turbulence. Consequently, we can interpret the protostellar mass distributions of the moderate and strong initial field models with decaying turbulence, and all of the models with driven turbulence, as being only weakly time-dependent. We do warn, however, that this stabilization is a result of an equilibrium between new stars forming at low mass, and existing stars growing. As noted in Krumholz et al. (2016), if some process were to halt the formation of new stars but the existing stars were to continue accreting their envelopes, the mass distribution would shift upward. Figure 5. View largeDownload slide The evolution of the stellar mass distribution for models with turbulent forcing (right) and with decaying turbulence (left). In these figures, the central lines show the 50th percentile mass, while the shaded regions indicate the 25th–75th percentile range. Colours correspond to different values of $$\mu _\Phi$$, as indicated in the legend. Figure 5. View largeDownload slide The evolution of the stellar mass distribution for models with turbulent forcing (right) and with decaying turbulence (left). In these figures, the central lines show the 50th percentile mass, while the shaded regions indicate the 25th–75th percentile range. Colours correspond to different values of $$\mu _\Phi$$, as indicated in the legend. Fig. 6 shows the protostellar mass distributions at the termination of each model. The black curve in the figure shows the observed IMF of Chabrier (2005). The mass distributions in our models have stabilized over dynamical time-scales (tff, >med), as shown in Fig. 5. However, because the sources in our models are still accreting, our models are not exactly comparable to the IMF that is representative of the distribution of final stellar masses. For this reason, we also compare our models to the protostellar mass functions (PMFs) from the analytic theory of McKee & Offner (2010); these PMFs are the mass distributions that one expects a series of protostars to possess if their final mass distribution matches the Chabrier (2005) IMF, and if the numbers of stars formed versus time N*(t) and the accretion histories of individual stars follow some specified functional form. We observe that N*(t) is a roughly linear function in all simulations except those that display runaway collapse, and thus we adopt a linear function for N*(t). For the accretion rate versus time of individual protostars, McKee & Offner (2010) consider a range of cases, all of which we plot together with the simulation results in Fig. 6: competitive accretion (CA), turbulent core (TC), 2-component competitive accretion (2CCA), and 2-component turbulent core (2CTC); the latter two both approach the constant accretion rate expected for an isothermal sphere for low-mass sources (IS), and we also plot the pure IS case.4 The only other parameter required by the McKee & Offner (2010) PMF analytic model is the upper limit to the most massive star that will form in the cluster mu. We have generated the theoretical PMFs in Fig. 6 by setting mu to the larger of 3  M⊙ or the value of mu necessary so that the theoretical PMF has exactly one source heavier than the second most massive particle present in the termination of the simulation. The choice of mu largely affects the shape of the high-mass tail of the PMF and the conclusions that follow from comparing our simulation to the theoretical PMFs are insensitive to this choice. The p-value of the Kolmogorov–Smirnov (K-S) statistic comparing each simulation to each mass function is reported in Table 3. The relative disagreement between simulations and mass functions pairs with p-values less than 0.05 is unlikely due to random chance. These distribution pairs are different with high statistical confidence. We denote the p-value of these pairs in bold text in the table. Figure 6. View largeDownload slide Histograms of protostellar masses at the end of each of the models with turbulent forcing (right) and with decaying turbulence (left). The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ from top to bottom. The black line shows the observed IMF (Chabrier 2005). The other curves denoted in the figure legend show the analytic PMFs of McKee & Offner (2010) for the single- and two-component turbulent core models, the single- and two-component competitive accretion models and the isothermal sphere model. Figure 6. View largeDownload slide Histograms of protostellar masses at the end of each of the models with turbulent forcing (right) and with decaying turbulence (left). The four rows show the $$\mu _\Phi =1.56$$, 2.17, 23.1, and ∞ from top to bottom. The black line shows the observed IMF (Chabrier 2005). The other curves denoted in the figure legend show the analytic PMFs of McKee & Offner (2010) for the single- and two-component turbulent core models, the single- and two-component competitive accretion models and the isothermal sphere model. Table 3. K-Sstatistic p-values for the distribution of protostars in each simulation as compared to the observed IMF (Chabrier 2005) and several analytical PMFs including the TC, 2CTC, CA, 2CCA, and the isothermal sphere models (McKee & Offner 2010). Simulations that include turbulent forcing are denoted as ‘Driven’ in the left-hand column. The other models are the cases of decaying turbulence. Simulation to mass function pairs with p-values less than 0.05 can be confidently interpreted as a sample drawn from a different distribution than the given mass function. The p-value for these ‘rejected’ models is denoted in boldface in the table. Comparison analytic model  IMF  TC  2CTC  CA  2CCA  IS  $$\mu _\Phi =1.56$$  3.2 × 10−4  0.094  1.7 × 10−3  0.69  1.0 × 10−3  1.1 × 10−7  $$\mu _\Phi =1.56$$ No wind  1.1 × 10−3  7.7 × 10−5  1.0 × 10−3  6.0 × 10−6  4.4 × 10−4  0.054  $$\mu _\Phi =2.17$$  0.19  0.92  0.27  0.56  0.25  0.013  $$\mu _\Phi =23.1$$  0.058  4.2 × 10−3  0.014  1.1 × 10−3  9.9 × 10−3  0.35  $$\mu _\Phi =\infty$$  0.50  1.3 × 10−3  0.011  1.4 × 10−4  9.2 × 10−3  0.35  $$\mu _\Phi =1.56$$ Driven  3.0 × 10−3  0.52  0.13  0.91  0.080  2.1 × 10−3  $$\mu _\Phi =2.17$$ Driven  3.5 × 10−3  0.72  0.086  0.94  0.050  1.6 × 10−4  $$\mu _\Phi =23.1$$ Driven  0.033  0.19  0.032  0.63  0.033  2.9 × 10−5  $$\mu _\Phi =\infty$$ Driven  4.4 × 10−3  0.43  0.036  0.60  4.4 × 10−3  1.9 × 10−5  Comparison analytic model  IMF  TC  2CTC  CA  2CCA  IS  $$\mu _\Phi =1.56$$  3.2 × 10−4  0.094  1.7 × 10−3  0.69  1.0 × 10−3  1.1 × 10−7  $$\mu _\Phi =1.56$$ No wind  1.1 × 10−3  7.7 × 10−5  1.0 × 10−3  6.0 × 10−6  4.4 × 10−4  0.054  $$\mu _\Phi =2.17$$  0.19  0.92  0.27  0.56  0.25  0.013  $$\mu _\Phi =23.1$$  0.058  4.2 × 10−3  0.014  1.1 × 10−3  9.9 × 10−3  0.35  $$\mu _\Phi =\infty$$  0.50  1.3 × 10−3  0.011  1.4 × 10−4  9.2 × 10−3  0.35  $$\mu _\Phi =1.56$$ Driven  3.0 × 10−3  0.52  0.13  0.91  0.080  2.1 × 10−3  $$\mu _\Phi =2.17$$ Driven  3.5 × 10−3  0.72  0.086  0.94  0.050  1.6 × 10−4  $$\mu _\Phi =23.1$$ Driven  0.033  0.19  0.032  0.63  0.033  2.9 × 10−5  $$\mu _\Phi =\infty$$ Driven  4.4 × 10−3  0.43  0.036  0.60  4.4 × 10−3  1.9 × 10−5  View Large We find that our runs with weak initial fields ($$\mu _\Phi =23.1$$ and ∞) and no driving have peak masses that exceed the observed mode of the IMF (m* = 0.2  M⊙ – Chabrier 2005), and that are continuing to rise (Fig. 5), at the termination of the models. As discussed in Section 3.1, these models also approach a state of rapid, near global free-fall collapse, contrary to observational constraints. Furthermore, the shape of the mass distributions is closest to that of the IS case. The ever-increasing peak of the PMFs in these simulations appears to be a manifestation of the overheating problem identified by Krumholz et al. (2011), whereby values of εff close to unity produce such high accretion luminosities that fragmentation is strongly suppressed, leading to a mass function that shifts to ever-higher masses as collapse continues but no new stars can be created. Furthermore, these simulations have low K-S statistic p-values when compared against all analytic models except for the IS one. In contrast, the moderate and strong initial field cases all achieve stable or slowly increasing median masses in the range 0.05 M⊙ < m* < 0.15 M⊙, independent of the presence of external turbulent forcing. All of the models with continuously driven turbulence also show stable peak masses independent of the strength of the initial magnetic field. At the termination of these models, the median masses lie between the analytically predicted peak masses for the single CCA and single CTC models of 0.05 and 0.095  M⊙, respectively. Because of the median masses for the CA and TC models are close to those of the numerical simulations, these simulations also have high K-S p-values when compared to the CA and TC theoretical models. The shape of the mass functions are not amenable to reliably differentiating among the moderately magnetized models or the weakly magnetized models with continuous turbulent forcing to within the statistical limitations of the number of sources in these simulations. Fig. 7 shows the protostellar mass distribution at the termination of the strongly magnetized $$\mu _\Phi =1.56$$ without turbulent driving or protostellar outflow ejection. The effect of outflows on the protostellar mass distribution can be seen clearly by comparing this figure to the top left panel of Fig. 6, which shows the undriven $$\mu _\Phi =1.56$$ model with outflows. Clearly omission of outflows leads to a significant increase in the typical stellar mass and very low K-S test p-values. This difference is partly due to direct removal of mass by the outflows, but this effect alone cannot explain the difference in mass distribution, partly because it is too large (a factor of ≳10 in mass), and partly because there is also a dramatic drop in the total number of stars when we turn-off outflows, something that cannot easily be explained by mass ejection. Instead, the shift in the mass function here appears to result from the interaction between outflows and radiative heating. Outflows indirectly reduce the luminosity of individual stars by lowering their accretion rates (Krumholz et al. 2012; Hansen et al. 2012; Myers et al. 2014), and create paths of reduced optical depth, allowing more efficient escape of radiation from the gas immediately around the protostars (Krumholz, McKee & Klein 2005; Cunningham et al. 2011). Due to the omission of these two effects in the no-outflow model, radiative heating is much more effective at suppressing fragmentation, leading to the production of fewer, much more massive stars. We showed in Section 3.1 that the strong initial magnetic support in the no-outflow model was sufficient to limit the star formation efficiency to values consistent with observation, but clearly the shape of the mass distribution of this model is not consistent with observation. Figure 7. View largeDownload slide Histogram of the protostellar masses at the end of the $$\mu _\Phi =1.56$$ model without protostellar feedback. Figure 7. View largeDownload slide Histogram of the protostellar masses at the end of the $$\mu _\Phi =1.56$$ model without protostellar feedback. It is also interesting to compare our results to previous simulations making use of a similar set of physical processes. Our $$\mu _\Phi =1.56$$ model with continuous driving is comparable to the recent models by Li et al. (2018); the simulations have similar initial (βIC and $$\mathcal {M}_{{\rm A},\rm IC}$$) and post-driving (β>med, $$\mathcal {M}_{{\rm A},\rm >med}$$ and ρ>med) conditions. However, while in our case this model produces an IMF peak at slightly less than 0.1  M⊙, the Li et al. (2018) simulations yield at peak at 0.2  M⊙, comparable to the observed IMF peak. Given the similarities of the dimensionless conditions of the two cluster models, we attribute this difference to the outflow ejection parameters discussed in Section 2.3. The models presented here eject three times the momentum flux in the early phases of pre-stellar collapse and therefore have more disruptive effects on their parent cores. While Cunningham et al. (2011) point out that observations do not precisely constrain the outflow model ejection parameters, the fact that the model in Li et al. (2018) achieves better agreement with the well-constrained median IMF suggests that the outflow ejection parameters used in that model are probably more realistic. This conclusion, however, does not affect the relative differences between the various runs that we have identified; indeed, a shift to higher masses as a result of somewhat weaker outflows would worsen the agreement between the observed IMF and those produced in the weakly magnetized, non-driven models. The other obvious prior work to which we can compare are the simulations of Hansen et al. (2012). Indeed, our non-magnetized, decaying turbulence model is identical to theirs apart from the random phases used during the initial turbulent driving phase. Both models ran for comparable times as well. However, the model presented in this work exhibits rapidly increasing accretion rates and commensurately increasing source masses at late time, while the simulations of Hansen et al. do not. While this might at first seem quite surprising, the result can be understood by examining Fig. 1. This figure shows that the stellar mass at first increases relatively slowly, but then rapidly increases as the turbulence decays and free-fall collapse begins. The behaviour in the Hansen et al. simulation is consistent with what we see just before the onset of rapid collapse in our simulation. Since the exact time at which runaway collapse begins is almost certainly sensitive to the random driving, the most likely explanation for the difference is that the Hansen et al. simulation was simply not run quite long enough to reach the runaway collapse phase where the collapse begins to alter the IMF. 3.4 Stellar multiplicity The statistics of multiple systems provides another observational constraint against which our models can be tested. To compute the stellar system multiplicity, we use the algorithm of Bate (2009). The algorithm works recursively: one identifies the most strongly bound pair of stars in the simulation, and replaces it with a single object at the centre of mass position of the pair with the total mass and momentum of its constituent protostars, then repeats until there are no more bound pairs that can be combined. In cases where combining the most bound pair of objects would create a system with a multiplicity greater than four, one skips it and combines the next most bound combination instead. This is motivated by the fact that high-multiple bound systems are dynamically unstable and such systems are likely to break apart if the model were continued further in time. At the end of this procedure, we obtain S single star systems, B binary systems, T triple systems, and Q quadruple systems. For the set of bound systems, the multiplicity fraction is defined as (Bate 2009; Krumholz et al. 2012; Myers et al. 2014)   $$\textrm{MF} = \frac{B+T+Q}{S+B+T+Q}$$ (13) Fig. 8 shows the multiplicity fraction as a function of the primary star mass in each of our simulations. Given the limited number of systems in at least some of our simulations, computing the multiplicity fraction requires some care. We do so in two ways. First, we form running averages, meaning that, for each primary star (i.e. single star or the most massive member of a multiple), we evaluate equation (13) for the set of systems consisting of that primary and the next most massive and next least massive primaries. This quantity is shown as the blue curve in Fig. 8. Second, we can compute multiplicity in bins. To do so, we start with the least massive primary, and then place the primaries in bins of mass, with the bin width taken to be either sufficient to contain four primaries, or 0.2 dex, whichever is greater. For each bin, we regard the systems in it as samples drawn from a binomial distribution of multiples or singles, and from those samples, we compute a central estimate and a 68 per cent confidence interval on the true multiple system fraction in that bin (see Krumholz et al. 2012 for full details on how this computation is carried out). We plot the binned 68 per cent confidence intervals as the shaded regions in Fig. 8. Figure 8. View largeDownload slide Multiple stellar system fraction at the end of each of the models with turbulent driving (right) and with decaying turbulence (left), computed as a function of primary star mass. The four rows show the $$\mu _\Phi =1.56, 2.17, 23.1$$, and ∞, from top to bottom. For the simulations, the blue curve indicates evolution of a running average over three systems as a function of mass, while shaded regions indicate the mean multiplicity and central 68 per cent probability range on the multiplicity in bins – see the main text for details on how the running average and binned values are computed. Black crosses denote the mass range and mean multiplicity found in observational surveys. The two highest mass surveys are statistical lower bounds. The horizontal marker on these observations indicate the lower bound. The observational data, from low to high mass, are compiled from Basri & Reiners (2006) and Allen (2007) (combined into the left most point), Fischer & Marcy (1992), Raghavan et al. (2010), Preibisch et al. (1999), Mason et al. (2009), and these are the same data used in Bate (2012) and Krumholz et al. (2012). Figure 8. View largeDownload slide Multiple stellar system fraction at the end of each of the models with turbulent driving (right) and with decaying turbulence (left), computed as a function of primary star mass. The four rows show the $$\mu _\Phi =1.56, 2.17, 23.1$$, and ∞, from top to bottom. For the simulations, the blue curve indicates evolution of a running average over three systems as a function of mass, while shaded regions indicate the mean multiplicity and central 68 per cent probability range on the multiplicity in bins – see the main text for details on how the running average and binned values are computed. Black crosses denote the mass range and mean multiplicity found in observational surveys. The two highest mass surveys are statistical lower bounds. The horizontal marker on these observations indicate the lower bound. The observational data, from low to high mass, are compiled from Basri & Reiners (2006) and Allen (2007) (combined into the left most point), Fischer & Marcy (1992), Raghavan et al. (2010), Preibisch et al. (1999), Mason et al. (2009), and these are the same data used in Bate (2012) and Krumholz et al. (2012). Comparing either the running averages or the binned results to the observed multiplicity fractions that we also plot in Fig. 8, we find that all of the models generally agree with observations within the range of statistical uncertainties. It is noteworthy that the weakly magnetized models without turbulent forcing are consistent with observed multiplicity statistics even though their PMFs (Section 3.3) are strongly inconsistent with observations. The driven turbulence $$\mu _\Phi =2.17$$ and 23.1 models somewhat overpredict the multiplicity fraction for systems with primary mass Mp ∼ 1  M⊙, but this may be the result of poor statistics in this mass range in those models, and we note that the observational surveys are based on field stars, while there are strong hints that multiplicity fraction is higher for still-embedded systems (Tobin et al. 2016), which is the more relevant comparison for our simulations.5 We also note that, while the bins at low mass are generally consistent with the observations individually, taken together it is clear that we systematically underproduce multiple systems at low mass relative to the observational data. This is a likely a resolution effect. Our models do not resolve gas-particle gravity forces below the scale where sink source terms couple to the gas –4Δx. Low-mass multiple systems must necessarily be close and consequently more difficult to resolve to remain bound. This limitation is also noted in the multiplicity fractions extracted from the models of Bate (2009), Krumholz et al. (2012), and Myers et al. (2014). 3.5 Turbulent decay and protostellar feedback rates Fig. 9 shows the temporal evolution of the volume-weighted rms gas velocity for each model. As intended, the models that include steady turbulent forcing remain within 10 per cent of a constant rms velocity with αvir ∼ 1. For the runs without driving, at early time, t ≲ tff, turbulent decay causes the rms velocity to drop. However, after 1–2 free-fall times (depending on the level of magnetic support), the system begins to undergo a global collapse, and this causes the velocity dispersion to rise again. Figure 9. View largeDownload slide Volume-weighted root-mean-squared gas velocity, normalized by the gas sound speed, as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). Note the difference in scales between the left- and right-hand panels. Figure 9. View largeDownload slide Volume-weighted root-mean-squared gas velocity, normalized by the gas sound speed, as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). Note the difference in scales between the left- and right-hand panels. Fig. 10 shows the temporal evolution of the mass-weighted rms velocity dispersion. This averaging is more analogous to mass-biased observational line width–size measurements. However, it is also heavily biased towards dense, locally collapsing regions, which somewhat mutes the effect of turbulent decay relative to infall. Consequently, the clear drop in velocity dispersion seen in Fig. 9 for the non-driven runs is absent here. We infer that mass-weighted line-width diagnostics are surprisingly insensitive to interruptions in large-scale driving source of durations up to ≲ 1tcross. Figure 10. View largeDownload slide Same as Fig. 9, but now showing mass-weighted rather than volume-weighted velocities. As with Fig. 9, note the difference in scales between the two panels. Figure 10. View largeDownload slide Same as Fig. 9, but now showing mass-weighted rather than volume-weighted velocities. As with Fig. 9, note the difference in scales between the two panels. It is also interesting to compare the $$\mu _\Phi =1.56$$ cases with and without outflows. In principle, the outflows have more than enough mechanical power to alter the turbulence. We illustrate this in Fig. 11, which shows the outflow mechanical luminosity versus time, $$L_{\rm outflow} = (1/2) \dot{M}_{\rm wind} v_{\rm wind}^2$$. This is typically ∼1–10 L⊙ at late times in our simulations, For comparison, the turbulent decay rate is $${\sim } \bar{\rho } v_{\rm rms}^2 L^3/t_{\rm cross} \sim 10^{-3}\ {\rm L}_{{\odot }}$$, roughly three orders of magnitude smaller. However, when we examine the volume-weighted rms velocity, the $$\mu _\Phi =1.56$$ runs with and without outflows are nearly identical, and in the mass-weighted plot they are only slightly different. These results indicate that outflow mechanical power is efficiently dissipated by radiative shocks and does not quickly couple to the volume-filling turbulence of the surrounding molecular cloud, consistent with earlier works exploring the coupling of continuously driven jets through a stratified environments (Banerjee, Klessen & Fendt 2007). What effects the outflows do have are limited to the dense gas surrounding protostellar cores, which is why the mass-weighted average shows an effect from outflows but the volume-weighted one does not. However, we emphasize that even in a mass-weighted sense the effect is modest. In our simulations, outflows do not appear to be efficient drivers of turbulence. Figure 11. View largeDownload slide Total mechanical luminosity of the protostellar outflows as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). The data have been boxcar-smoothed with an interval of 0.05 tff. Figure 11. View largeDownload slide Total mechanical luminosity of the protostellar outflows as a function of time for the models with turbulent forcing (right) and with decaying turbulence (left). The data have been boxcar-smoothed with an interval of 0.05 tff. 3.6 Local magnetic support The turbulent driving phase introduces stretching and amplification of the magnetic field, relative to the initial magnetic field strength. The Alfvèn Mach number captures the strength of the turbulence relative to the initial magnetic field. As discussed in Section 2, the magnetized models in this work span the range $$1 < {\cal M}_{\rm A}< 15$$. However, the weakly magnetized models endure more significant stretching of the magnetic field and more small-scale magnetic field amplification in the turbulent driving phase than the more strongly magnetized models. The tangled magnetic field is further amplified due to collapse after self-gravity is switched on. Therefore, it is only the mean, large-scale field that is weak in our ‘weak-field’ model. On the scale of the pre-stellar dense cores, the magnetic field is significantly more tangled and stretched than in stronger field models. Similar field amplification due to turbulence and the early gravitational collapse phase into starless dark clumps was noted in Li et al. (2015, section 8 in particular). In this section, we examine the magnetic field structure in the vicinity of protostars, with the goal of understanding how the combination of turbulent and gravitational amplification maps from the initial, large-scale magnetic field to the final, small-scale one in the dense gas that is bound to and actively feeding accretion on to protostellar sources. We seek to characterize the mean mass-to-flux ratio, averaged over the bound gas around each protostar that is poised for imminent accretion. We take this gas to be that which satisfies the condition   $$\frac{(\boldsymbol {v}-\boldsymbol {v_*})\cdot (\boldsymbol {x}-\boldsymbol {x_*})}{|\boldsymbol {x}-\boldsymbol {x_*}| } + c_{\rm s} < \sqrt{\frac{Gm_*}{|\boldsymbol {x}-\boldsymbol {x_*}|}},$$ (14)where $$\boldsymbol {x}$$, $$\boldsymbol {v}$$, and cs are the position, velocity, and sound speed of the gas in the computational domain and $$\boldsymbol {x_*}$$, $$\boldsymbol {v_*}$$, and m* are the position, velocity, and mass of each protostar. This criterion selects gas with insufficient thermal or kinetic energy to escape eventual accretion to the protostar, but distinguishes between infalling gas (for which $$\boldsymbol {v}-\boldsymbol {v_*}$$ is anti-aligned with $$\boldsymbol {x} - \boldsymbol {x_*}$$, and thus the dot product is negative), which is more bound, and outflowing gas (vectors aligned), which is less so. We assign gas that is bound to more than one star under this criterion only to the star to which it is most strongly bound. For each star i, we have now defined a mass mgas, i of ‘bound gas.’ We next compute the volume-weighted mean magnetic field direction in this gas, $$\hat{\boldsymbol B}_{\rm i}$$, and we define the associated magnetic flux Φi as the net flux through the intersection of the plane defined by $$\hat{\boldsymbol B}_{\rm i}$$ and the volume of assigned gas. Finally, we characterize the local mass-to-flux ratio of the gas assigned to each star as μΦ, i = mgas, i/MΦ, i where MΦ, i is the magnetic critical mass computed from equation (1). It should be noted that the significance of $$\mu _\Phi$$ changes if a protostar is present in the cloud: The protostellar gravity also acts on the gas, so that accretion can occur normal to the field even for $$\mu _\Phi <1$$. Fig. 12 shows the temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The average is weighted by the mass of each protostar. Due to the computational burden associated with post-processing this result from our models, the plots in Fig. 12 are sampled only every 0.2tff, >med in time. In the figure, the curve for each model begins after the formation of the first sink particle. Figure 12. View largeDownload slide Temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The averaging is weighted by protostellar mass. The plot is presented in units of the magnetically critical mass, for the models with turbulent forcing (right) and with decaying turbulence (left). Figure 12. View largeDownload slide Temporal evolution of the mean mass-to-flux ratio averaged over the bound gas around each protostar. The averaging is weighted by protostellar mass. The plot is presented in units of the magnetically critical mass, for the models with turbulent forcing (right) and with decaying turbulence (left). In interpreting Fig. 12, it is important to bear in mind that this diagnostic characterizes the net impact of multiple competing processes. Under our sink particle prescription (Lee et al. 2014), infalling gas is released from its associated magnetic flux the instant that it deposited on to the sink particle. This prescription is intended to mimic resistive effects near the surface of the protostar at length scales below that which can be resolved by our models. Consequently, accretion causes the mass-to-flux ratio associated with locally bound dense gas clumps to decrease. The accreted mass is removed from the gas but the flux remains. However, the magnetic flux may subsequently escape, because it is no longer anchored to the gravitationally bound material – this is a form of magnetic interchange instability (Zhao et al. 2011; Krasnopolsky et al. 2012; Cunningham et al. 2012; Li et al. 2014). If the magnetic flux tubes escape the gravitationally bound region, this drives the mass-to-flux ratio back up. Similarly, turbulent magnetic reconnection acts to diffuse magnetic flux and increases the local mass-to-flux ratio in regions that undergo reconnection. In this case, the gas mass remains but the magnetic flux is displaced outward (Santos-Lima, de Gouveia Dal Pino & Lazarian 2012). Outflow ejection and turbulent forcing also impart mechanical energy to the system, and the resulting flows can further amplify, entangle, and ultimately reconnect magnetic flux tubes (Sur et al. 2012; Li et al. 2015). In this way, accretion, outflow, field tangling, and reconnection all influence the mass-to-flux ratio. Since accretion itself is episodic, so too is the evolution of the mean mass-to-flux ratio of the bound gas in Fig. 12. We anticipate the degree of tangling and episodic reconnection to be most pronounced in the weak-field models as weaker fields are more readily tangled and amplified by the turbulent flow. This too is borne out in Fig. 12. In the undriven models shown in the left-hand panel of Fig. 12, the mass-to-flux ratios in the bound gas systematically decreases with time for each of the models. This indicates that buoyant flux tubes rising away from the protostars as a consequence of accretion act enhance the degree of magnetic support to the bound gas, and this enhancement overwhelms the effect of turbulent reconnection. The $$\mu _\Phi =1.56$$ model without outflows indicates systematically stronger magnetic support in bound gas than the case with outflows. This indicates that outflows contribute to magnetic reconnection and consequent outward transport of magnetic flux, and explains why these models achieve comparable star formation efficiency (εff, IC = 12 per cent in case without outflows and εff, IC = 23 per cent in case with outflows). The enhanced magnetic support in case without outflows offsets, the loss of mechanical feedback support when the outflows are turned off. The $$\mu _\Phi =2.17$$ model is also strongly magnetically supported, with the mass-to-flux ratio in the bound gas approaching $$\mu _\Phi \sim 1$$ at late time, consistent with its accretion rate εff, IC = 14 per cent, slightly higher than the more strongly magnetized $$\mu _\Phi =1.56$$. The model that began with a field characterized globally as $$\mu _\Phi =23.1$$ attains a characteristic mass-to-flux ratio $$\mu _\Phi \approx 5$$ via field amplification and gravitational compression by the time the first protostellar particle forms. While the level of magnetic support in the bound gas increases with time, the degree of magnetic support achieved by the end of the model is insufficient to significantly reduce the net accretion rate from that of the purely hydrodynamical case. The evolution of the mass-to-flux ratio of the bound gas in the models with driven turbulence models is shown in the right-hand panel of Fig. 12. In this case, the model that began with a global $$\mu _\Phi =2.17$$ retains a roughly constant mass-to-flux ratio in the bound gas of $$\mu _\Phi \sim 2$$ throughout its evolution. The model that began with a global $$\mu _\Phi =1.56$$ attains stronger magnetic stabilization in the dense gas with $$\mu _\Phi < 1$$ by the time the first protostars have formed. Subsequent evolution indicates a slowly decreasing level of magnetic support in the bound gas, stabilizing to $$\mu _\Phi \gtrsim 1$$. The model that began with a global $$\mu _\Phi =23.1$$ enhances the degree of magnetic support in the bound gas with time, reaching to $$\mu _\Phi \sim 4$$ by the end of the model. Taken in aggregate, our results indicate that small-scale dissipation effects near the surface of a protostar provide a feedback loop for stabilizing the star formation efficiency. Bound gas cores that are undersupported initially undergo more rapid collapse than cores that are better supported, independent of whether the initial support is predominately turbulent or magnetic. The more rapidly collapsing regions reconnect and expel magnetic flux to the core from the newly formed protostar at a faster rate. The end result of this process over a time-scale >1tff is that cores with continuous mechanical support that remain in near virial equilibrium approach mass-to-flux ratios of $$\mu _\Phi < 5$$ independent of the large-scale mean field strength. Models with decaying or otherwise sub-virial levels of turbulent support will lead to greater local enhancement of the level of magnetic support in the bound gas, approaching $$\mu _\Phi < 2$$ at late time. 4 MAGNETIC AND RADIATIVE EFFECTS IN SHAPING THE IMF In Section 3.3, we demonstrated that our models have protostellar mass distributions that are stable or evolving slowly towards higher mass only slowly compared to cloud-collapse time-scales. The models with weak initial fields or a lack of driving all have median masses that are too large in comparison to the observed IMF, while the strong field and driven cases have smaller median masses that are at least roughly compatible with these models reproducing the observed IMF. It is therefore interesting to investigate in more depth what mechanisms shape this mass distributions. To explore this question, we turn to the techniques described in Krumholz et al. (2016), who analysed the simulations of Myers et al. (2014); these simulations use the same physics we include here, but focus on a much denser region appropriate to sites of high-mass star formation. Comparing the two sets of simulations therefore enables us to understand how this question depends on the star-forming environment. 4.1 Analysis methodology We begin with a brief summary of the Krumholz et al. (2016) analysis, and refer readers interested in the details to that paper for a full description. The central idea of this analysis is to determine what physical process is responsible for halting fragmentation in the vicinity of each protostar, thereby leaving the mass available for accretion and the build-up of more massive stars, rather than the production of additional low-mass stars. To this end, we examine spherical regions of gas within a $$4000\ {\rm {\rm au}}$$ radius around each protostellar sink particle at time intervals of 0.05tff, >med. We divide the volume of interest into 128 spherical shells of radius r concentric with the position of the sink particle. For each shell, we compute the total gas enclosed mgas, excluding the sink particle mass, the mean density ρ of the enclosed material, and its mass-weighted mean sound speed cs.6 From the resulting profiles in ρ and cs, we compute an effective temperature   $$T_{\rm eff} = \frac{\mu m_{\rm H}}{k_{\rm B}} c_{\rm s}^2 ,$$ (15)and, more importantly from the standpoint of determining whether gas will fragment, the Bonnor–Ebert mass (Ebert 1955; Bonnor 1956),   $$m_{\rm BE} = 1.86 \sqrt{\frac{c_{\rm s}^3}{G^3\rho }}.$$ (16)The Bonnor–Ebert mass characterizes the degree of thermal support. Absent other forces, objects less massive than mBE are stable against collapse and objects more massive than mBE are unstable against collapse. To elucidate the importance of heating by protostars in stabilizing the core, we repeat the computation with the sound speed cs = 0.19 cm s−1, the sound speed of molecular gas at 10 K used as the background isothermal temperature in the simulations. We shall denote this quantity as mBE, 10. Comparison of mBE with mBE, 10 enables us to determine the importance of radiative heating in limiting fragmentation. To characterize the importance of magnetic forces, we compute the magnetic flux Φ threading each spherical shell, defined so that $$\Phi = \int |\boldsymbol {B}\cdot \hat{n}|\, d^2x$$ is computed along a circular cross-section that is concentric with the spherical shell and has normal $$\hat{n}$$. The absolute value is taken in the integration on the assumption that oppositely directed fluxes do not reconnect and cancel. We consider 12 possible orientations for $$\hat{n}$$, uniformly distributed on the unit sphere, with one value aligned with the star's angular momentum vector, and use the largest value of Φ, we find for all the $$\hat{n}$$ values considered. From this, we can compute an effective mean magnetic field strength inside a radius r,   $$B_{\rm eff}(<\!r) = \frac{\Phi }{\pi r^2},$$ (17)a magnetic critical mass $$m_\Phi$$ (equation 1), and the magnetically supported mass (Mouschovias & Spitzer 1976; McKee & Ostriker 2007)   $$m_{\rm B}=\frac{m^3_\Phi }{m^2_{\rm gas}}.$$ (18)We prefer mB to $$m_\Phi$$ for this analysis as the former is an intensive quantity that does not depend on the size of the volume considered. It has been well established that magnetic fields influence the spectrum of stars formed and reduce the star formation efficiency (Price & Bate 2008). Here, we elucidate the relative contribution or radiative feedback to magnetic support by comparing the enclosed mass mgas with the thermally and magnetically supported masses mBE and mB. This provides direct insight into the role of thermal support (augmented by radiative heating) and magnetic fields in shaping the IMF. For sufficiently small shells around each protostar, we always find mgas ≪ mBE and mgas ≪ mB, i.e. the enclosed mass is too small to fragment and form another star given the level of thermal and magnetic supports. This gas will almost certainly be accreted on to the star that it surrounds, augmenting this star's mass. At sufficiently large radii, mgas ≫ mBE and mgas ≫ mB, i.e. the mass enclosed is large enough that is able to fragment and produce additional stars, rather than being accreted on to the star it surrounds. Thus, the shell at which mgas ≈ max (mBE, mB) marks a rough dividing line between gas that will end up in the existing star and gas that will end up elsewhere, and m* + mgas provides a rough estimate of the mass to which the star is likely to grow. We define the thermal critical mass mBE, crit as the smallest (and almost always sole) mass where mgas = mBE, and analogously for mB, crit and mBE, 10, crit. Thus these critical masses provide a rough estimate of the future mass supply available for each star. Krumholz et al. (2016) show that, in the simulations of Myers et al. (2014), stars with mass m* ≪ 0.2  M⊙ almost invariably have stabilized gaseous envelopes such that mBE, crit ≫ m*, explaining the relative rarity of brown dwarfs as compared to stars: brown dwarf mass fragments can form, but when they do, they stabilize enough gas around themselves that they usually continuing accreting well out of the brown dwarf mass regime. Moreover, Krumholz et al. (2016) find that their stabilized regions have mBE, crit ≫ mB, crit ≫ mBE, 10, crit, indicating that, while magnetic support is more important than thermal support in the absence of radiative heating, once one considers the effects of stellar radiation, it is the dominant factor in limiting fragmentation and thus establishing the location of the IMF peak. 4.2 Profiles of mean density, temperature, and magnetic field Fig. 13 shows the average density profile around central stars with three different masses; the averages plotted are taken over all protostars at all times. The density profiles always show a central dip at 130 au < r < 260 au, which is a result of our sink particle accretion zones; we therefore ignore this region in our analysis. At larger radii, we find that the moderately and weakly magnetized models ($$\mu _\Phi \ge$$ 2) without driving have density profiles ρ ∝ r−3/2. This profile is consistent with that expected for free-fall collapse on to a point mass. However, for Bondi-type flow, the density at a fixed distance from the central object should increase with mass. The fact that it is not inconsistent with pure Bondi flow, but is consistent with the TC model of McKee & Tan (2003) and with the model of Murray et al. (2015) for the for collapsing, self-gravitating, TCs. In contrast, the strongly magnetized ($$\mu _\Phi = 1.56$$) undriven cases, and all driven cases for a central star mass of m* = 0.15  M⊙, show slightly shallower density profiles. This is consistent with a build-up of magnetic flux resulting in a slower than free-fall collapse. Figure 13. View largeDownload slide Mean density profiles ρ(r) averaged in spherical shells around protostars. The panels show the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. The bottom row depicts results for the model with continued turbulence driving after the onset of self-gravity. Different colours show means for protostars of mass 0.025, 0.05, and 0.15  M⊙. For each mass bin, the central line shows the mean, while the shaded band shows the 1σ dispersion for protostars in that mass bin. The dashed black line shows the density scaling ρ ∝ r−3/2, expected for free-fall collapse. Figure 13. View largeDownload slide Mean density profiles ρ(r) averaged in spherical shells around protostars. The panels show the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. The bottom row depicts results for the model with continued turbulence driving after the onset of self-gravity. Different colours show means for protostars of mass 0.025, 0.05, and 0.15  M⊙. For each mass bin, the central line shows the mean, while the shaded band shows the 1σ dispersion for protostars in that mass bin. The dashed black line shows the density scaling ρ ∝ r−3/2, expected for free-fall collapse. Fig. 14 shows averaged profiles of Teff. Except for the $$\mu _\Phi = 1.56$$ undriven model, these profiles indicate increasing Teff with central star mass. As Krumholz et al. (2016) point out, to maintain a constant density profile, infall velocities and accretion rate must increase as the progenitor mass increases. Since accretion luminosity is the dominant source of radiative heating in low-mass sources, the temperature around them rises with their mass. The weakly magnetized models show the most rapid increase in heating with mass. Krumholz et al. (2016) noted a similar trend in models of denser high-mass star-forming regions, and also found that the temperature profiles were always close to Teff ∝ r−0.3. Here we find substantially shallower profiles, likely because the lower optical depth environment we are simulating is less effective at trapping the radiation (Krumholz & McKee 2008). Figure 14. View largeDownload slide Same as Fig. 13, but showing the mean effective temperature Teff(r) (equation 15). The black dashed line shows a Teff ∝ r−0.3. The initial temperature of the models is T = 10 K. Figure 14. View largeDownload slide Same as Fig. 13, but showing the mean effective temperature Teff(r) (equation 15). The black dashed line shows a Teff ∝ r−0.3. The initial temperature of the models is T = 10 K. In Fig. 15, we show the corresponding average profiles of Beff. Consistent with our findings in Section 3.6, despite the large differences in large-scale magnetic flux, there is very little difference between the small-scale magnetic field strengths around protostars in the different simulations. Nor does the effective field strength increase with star mass, as one might expect from a naive picture in which accretion drags magnetic flux towards a star. Figure 15. View largeDownload slide Same as Fig. 13, but showing the mean effective magnetic field Beff(<r) (equation 17). The black dashed horizontal line indicates the initial magnetic field strength before driving to the initial turbulent state. Figure 15. View largeDownload slide Same as Fig. 13, but showing the mean effective magnetic field Beff(<r) (equation 17). The black dashed horizontal line indicates the initial magnetic field strength before driving to the initial turbulent state. The non-dependence of the small-scale magnetic fields on either the large-scale flux or the amount of gas that has already been accreted strongly suggests that the magnetic field profiles we find are the result of a self-regulation mechanism operating on small scales. We can think of this self-regulation as due to two pairs of competing processes, although given our limited resolution and the complex nature of the flows in our simulations, we cannot easily identify which mechanisms dominate. First, there is a competition between dynamo amplification and reconnection. In a turbulent dynamo, magnetic fields amplify until these two effects balance (Sur et al. 2012; Li et al. 2015). Second, there is a competition between advection and magnetic interchange instabilities. As gas accretes, magnetic flux tubes are advected into the cores around the protostars, which tends to drive up Beff. The countervailing process is escape of flux through magnetic interchange instabilities. In reality, the interchange instabilities are triggered by resistive effects close to the surface of the protostar, which decouple flux tubes from the mass that anchors them and allows them to drift outward (Zhao et al. 2011; Krasnopolsky et al. 2012; Cunningham et al. 2012; Li et al. 2014); we model this by not accreting any magnetic flux into the star particles. If the flux liberated from the protostar diffused out to a radius R, then it would contribute a mean field   $$B_{\rm eff,*}=2G^{1/2}\left(\frac{L}{Z}\right)\frac{m_*}{\mu _\Phi R^2},$$ (19)where we have assumed that the protostar has formed from a cylinder of gas with a height Z less than the box height L. For m* = 0.15  M⊙ and $$\mu _\Phi =1.56$$, the results in Fig. 15 place a lower limit R > 1000 au on the distance that the field has diffused outward. 4.3 Critical mass In Fig. 16, we show the critical mass of gas supported by thermal pressure, mBE, thermal pressure fixed to the initial temperature, 10 K mBE, 10, and magnetic fields, mB, as a function of central stellar mass m* as described in Section 4.2. When mcrit > m*, the mass of the stabilized envelope around a protostar is much greater than the mass of the star itself, and thus is likely to greatly increase its mass. Figure 16. View largeDownload slide Critical mass of gas supported by thermal pressure mBE, thermal pressure fixed to the initial temperature 10 K mBE, 10 and magnetic fields mB as a function of central stellar mass m*. Shaded bands indicate 1σ dispersion over the protostars. The dashed black line indicates where the mass supported against collapse by pressure or magnetic forces is sufficient to double the present stellar mass, m* = mcrit. The panels shows the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. Figure 16. View largeDownload slide Critical mass of gas supported by thermal pressure mBE, thermal pressure fixed to the initial temperature 10 K mBE, 10 and magnetic fields mB as a function of central stellar mass m*. Shaded bands indicate 1σ dispersion over the protostars. The dashed black line indicates where the mass supported against collapse by pressure or magnetic forces is sufficient to double the present stellar mass, m* = mcrit. The panels shows the results for the strong (left), moderate (centre), and weak (right) magnetic field models. The top row depicts results for the models without turbulence driving after the onset of self-gravity at t = 0. First consider thermal support. Krumholz et al. (2016) found that mBE was typically a factor of ∼10 larger than mBE, 10, indicating that radiative heating increased the amount of mass that could be thermally supported by a factor of ∼10. In our lower mass cluster models, the difference between mBE and mBE, 10, is less pronounced, closer to 0.5 dex than 1 dex. However, the difference is smaller not because mBE is smaller, but rather because mBE, 10 is larger. That is to say, in our lower density simulation, the characteristic density is lower, and thus the Bonnor–Ebert mass in unheated gas is larger. In comparison, the mass that can be supported after heating, mBE, is almost unchanged and independent of environment – the increase in the temperature of the gas has been compensated by an increase in its density. Next consider magnetic support. The results in Fig. 16 are consistent with the findings of Krumholz et al. (2016) that the magnetic fields are relatively unimportant in supporting the cores compared to thermal pressure. Again, this points to the importance of local processes in regulating the magnetic field structure in the vicinity of protostars. The magnetic flux in the $$\mu _\Phi =1.56$$ case is sufficient to support a mass $$m_\Phi =119\ {\rm M}_{{\odot }}$$ against collapse, but the field threading the local region around the protostars is strong enough to support only 0.01–0.1  M⊙. This result explains why the mode of the IMFs (Section 3.3) in our models are relatively insensitive to the large-scale magnetic field. The critical mass supported against fragmentation is determined primarily by thermal effects enhanced by radiation. Magnetic fields thus have a significant effect on the rate of star formation, but not on the IMF. 5 CONCLUSIONS We present a numerical study of low-mass star cluster formation, and how this process is affected by the strength of the large-scale magnetic field, the presence or absence of turbulent energy input, and the effects of protostellar outflows. Our simulations include MHD, radiative transfer and protostellar radiation feedback, and protostellar outflows, a combination of physical processes that has received little attention in the literature to date. Our primary finding is that our models are able to reproduce the observed (low) efficiency of star formation and the observed location of the IMF peak (to within a factor 2) only in models that include outflows and that do not undergo global collapse, either because they are supported by external energy input that drives the turbulence, or because they are close enough to the line of magnetic criticality for magnetic fields to inhibit collapse. In contrast, multiplicity statistics are insensitive to all variations in the initial conditions. Krumholz et al. (2011, 2012) had previously conjectured that there might be a link between the low efficiency of star formation and the location of the IMF peak, but this was based on a single set of non-magnetized simulations. Our much more extensive suite of radiation-MHD simulations puts this conclusion on a much firmer footing, and suggests a general principle: the IMF peak location is directly linked to the inefficiency of star formation. The primary mechanism controlling the IMF peak appears to be accretion-powered radiation from forming protostars, which produces enhanced thermal pressure on small scales in the vicinity of forming stars, as noted earlier by Bate (2009) as well. In comparison, the small-scale magnetic field structure around forming protostars is both insensitive to the large-scale magnetic field strength, and dynamically sub-dominant compared to radiation when it comes to regulating gas fragmentation. These results strongly favour a picture in which the small-scale magnetic field structure around protostars is regulated by local competitions between accretion of magnetized material and dynamo amplification, which tend to enhanced the field strength, and magnetic interchange instability and turbulent reconnection, which tend to suppress it. The competition between these processes tends to force the local mass-to-flux ratio in gas bound to protostars to ∼2 (normalized to the critical value), regardless of the large-scale field. Our models lack resolution to capture the effects of thermal support against gravitational contraction deep within very low-mass cores and we lack non-ideal MHD effects that are most strongly coupled to turbulent density profiles on small scales. However, these effects should influence our modelled protostellar mass distributions only on the low-mass tail m* ≲ 0.05 M⊙ and our conclusions are not sensitive to these effects. We acknowledge that higher resolution models with non-ideal effects and multiple realizations to capture better statistics at protostellar mass extrema could improve the comparisons with observation made here by capturing a wider range of source masses with higher fidelity. Acknowledgements Support for this research was provided by NASA through NASA ATP grants NNX13AB84G (RIK and CFM), NNX15AT06G (MRK), the US Department of Energy at the Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 (AJC and RIK), the NSF through grant AST-1211729 (RIK and CFM), and the Australian Research Council's Discovery Projects funding scheme (project DP160100695, MRK). This research was also supported by grants of high-performance computing resources from the National Center of Supercomputing Application through grant TGMCA00N020, under the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation grant number OCI-1053575, the computing resources provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and resources and services from the National Computational Infrastructure, which is supported by the Australian Government. LLN-JRNL-737724. Footnotes 1 In contrast to these results, Murray (2011) and Lee, Miville-Deschênes & Murray (2016) argue εff ∼ 1 in the most luminous star clusters, based on their high ratios of ionizing luminosity to molecular gas mass. This interpretation has been disputed by Feldmann & Gnedin (2011) and Krumholz (2014), who suggest that these ratios are high not because εff is large, but rather because very massive clusters have strong feedback that destroys their parent clouds efficiently, thereby rendering invalid Murray's and Lee et al.'s assumption that one can identify and measure the gas mass of these clusters’ progenitors. However, this dispute is immaterial for the low-mass, low-density clusters with which we are concerned here, for which there is a strong observational consensus that εff is low. 2 We note that the time required to saturate the magnetic energy depends on the initial field strength. For models initialized with magnetic fields much weaker than those considered in this paper, significantly more than two dynamical times are required to fully saturate the magnetic energy of the system due to small-scale dynamos (Federrath et al. 2011; Tricco, Price & Federrath 2016). However, the models considered in this work are too strongly magnetized for small-scale eddies to generate significant field amplification, and the magnetic power spectra have been shown to saturate within a few crossing times in these regimes (Li et al. 2012). 3 The choice of the temperature at which one should evaluate the mean opacity for a grey method such as ours is in fact very subtle, even in calculations that separately track the dust, gas, and radiation temperatures. This is because tabulated opacity tables, such as the ones from Semenov et al. (2003) that we use, generally assume that all three temperatures are equal, and different temperatures matter for different physical effects – for example, condensation of ice mantles on grains depends on the gas temperature, evaporation of grains depends on the dust temperature, and the radiation spectrum depends on the radiation temperature. Absent tabulated opacities that consider out-of-equilibrium dust, gas, and radiation fields, there is no single choice of temperature at which the opacity can be computed that properly captures all these divergent effects. However, greater accuracy could be obtained by distinguishing the absorption opacity, which is primarily affected by the radiation temperature, and the emission opacity, which is primarily determined by the dust temperature. 4 We do not consider any of McKee & Offner's (2010) ‘tapered accretion’ analytic models, since the median accretion rate in our simulations shows no signs of tapering (Fig. 5). 5 We unfortunately cannot place the data from Tobin et al. (2016) on Fig. 8, because the masses of the individual stars are largely unconstrained. 6 As in Krumholz et al. (2016), we focus on cumulative rather than differential quantities because they are less subject to noise; see appendix A of Krumholz et al. (2016) for a demonstration that this choice does not alter the qualitative results. REFERENCES Allen P. R., 2007, ApJ , 668, 492 CrossRef Search ADS   Banerjee R., Klessen R. S., Fendt C., 2007, ApJ , 668, 1028 CrossRef Search ADS   Basri G., Reiners A., 2006, AJ , 132, 663 CrossRef Search ADS   Bate M. R., 2009, MNRAS , 392, 590 CrossRef Search ADS   Bate M. R., 2012, MNRAS , 419, 3115 CrossRef Search ADS   Bonnor W. B., 1956, MNRAS , 116, 351 CrossRef Search ADS   Chabrier G., 2005, in Corbelli E., Palle F., eds, Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later . Springer, Dordrecht, p. 41 Chakrabarti S., McKee C. F., 2005, ApJ , 631, 792 CrossRef Search ADS   Commerçon B., Hennebelle P., Henning T., 2011, ApJ , 742, L9 CrossRef Search ADS   Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ , 725, 466 CrossRef Search ADS   Cunningham A. J., Klein R. I., Krumholz M. R., McKee C. F., 2011, ApJ , 740, 107 CrossRef Search ADS   Cunningham A. J., McKee C. F., Klein R. I., Krumholz M. R., Teyssier R., 2012, ApJ , 744, 185 CrossRef Search ADS   Ebert R., 1955, Z. Astrophys. , 37, 217 Evans N. J., II, Heiderman A., Vutisalchavakul N., 2014, ApJ , 782, 114 CrossRef Search ADS   Falgarone E., Troland T. H., Crutcher R. M., Paubert G., 2008, A&A , 487, 247 CrossRef Search ADS   Federrath C., 2015, MNRAS , 450, 4035 CrossRef Search ADS   Federrath C., Klessen R. S., 2012, ApJ , 761, 156 CrossRef Search ADS   Federrath C., Klessen R. S., Schmidt W., 2009, ApJ , 692, 364 CrossRef Search ADS   Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010, A&A , 512, A81 CrossRef Search ADS   Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ , 731, 62 CrossRef Search ADS   Feldmann R., Gnedin N. Y., 2011, ApJ , 727, L12 CrossRef Search ADS   Fischer D. A., Marcy G. W., 1992, ApJ , 396, 178 CrossRef Search ADS   Gammie C. F., Ostriker E. C., 1996, ApJ , 466, 814 CrossRef Search ADS   García-Burillo S., Usero A., Alonso-Herrero A., Graciá-Carpio J., Pereira-Santaella M., Colina L., Planesas P., Arribas S., 2012, A&A , 539, A8 CrossRef Search ADS   Goldbaum N. J., Krumholz M. R., Matzner C. D., McKee C. F., 2011, ApJ , 738, 101 CrossRef Search ADS   Goldsmith P. F., 2001, MNRAS , 557, 736 Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, ApJ , 747, 22 CrossRef Search ADS   Hennebelle P., Commerçon B., Joos M. et al.  , 2011, A&A , 528, A72 CrossRef Search ADS   Heyer M., Gutermuth R., Urquhart J. S., Csengeri T., Wienen M., Leurini S., Menten K., Wyrowski F., 2016, A&A , 588, A29 CrossRef Search ADS   Klein R. I., Fisher R. T., McKee C. F., Truelove J. K., 1999, Numer. Astrophys. , 240, 131 CrossRef Search ADS   Klessen R. S., Hennebelle P., 2010, A&A , 520, A17 CrossRef Search ADS   Krasnopolsky R., Li Z.-Y., Shang H., Zhao B., 2012, ApJ , 757, 77 CrossRef Search ADS   Krumholz M. R., 2014, Phys. Rep. , 539, 49 CrossRef Search ADS   Krumholz M. R., 2015, preprint (arXiv:1511.03457) Krumholz M. R., McKee C. F., 2008, Nature , 451, 1082 CrossRef Search ADS PubMed  Krumholz M. R., Tan J. C., 2007, ApJ , 654, 304 CrossRef Search ADS   Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ , 611, 399 CrossRef Search ADS   Krumholz M. R., McKee C. F., Klein R. I., 2005, ApJ , 618, L33 CrossRef Search ADS   Krumholz M. R., Matzner C. D., McKee C. F., 2006, ApJ , 653, 361 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., Bolstad J., 2007, ApJ , 667, 626 CrossRef Search ADS   Krumholz M. R., Cunningham A. J., Klein R. I., McKee C. F., 2010, ApJ , 713, 1120 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ , 740, 74 CrossRef Search ADS   Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ , 745, 69; Erratum 2013, ApJ, 779, 89 CrossRef Search ADS   Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ , 754, 71 CrossRef Search ADS   Krumholz M. R., Myers A. T., Klein R. I., McKee C. F., 2016, MNRAS , 460, 3272 CrossRef Search ADS   Lada C. J., Lombardi M., Alves J. F., 2010, ApJ , 724, 687 CrossRef Search ADS   Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ , 745, 190 CrossRef Search ADS   Lada C. J., Roman-Zuniga C., Lombardi M., Forbrich J., Alves J. F., 2013, ApJ , 77, 133 CrossRef Search ADS   Larson R. B., 1981, MNRAS , 194, 809 CrossRef Search ADS   Lee A. T., Cunningham A. J., McKee C. F., Klein R. I., 2014, ApJ , 783, 50 CrossRef Search ADS   Lee E. J., Miville-Deschênes M.-A., Murray N. W., 2016, ApJ , 833, 229 CrossRef Search ADS   Lee Y.-N., Hennebelle P., 2016a, A&A , 591, A30 CrossRef Search ADS   Lee Y.-N., Hennebelle P., 2016b, A&A , 591, A31 CrossRef Search ADS   Li P. S., Martin D. F., Klein R. I., McKee C. F., 2012, ApJ , 745, 139 CrossRef Search ADS   Li P. S., Myers A., McKee C. F., 2012, ApJ , 760, 33 CrossRef Search ADS   Li P. S., McKee C. F., Klein R. I., 2015, MNRAS , 452, 2500 CrossRef Search ADS   Li P. S., McKee C. F., Klein R. I., 2018, MNRAS , 473, 4220 CrossRef Search ADS   Li Z.-Y., Nakamura F., 2006, ApJ , 640, L187 CrossRef Search ADS   Li Z.-Y., Banerjee R., Pudritz R. E., Jørgensen J. K., Shang H., Krasnopolsky R., Maury A., 2014, in Beuthe H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI . Univ. Arizona Press, Tucson, AZ, p. 173 Londrillo P., del Zanna L., 2004, J. Comput. Phys. , 195, 17 CrossRef Search ADS   Mac Low M.-M., 1999, ApJ , 524, 169 CrossRef Search ADS   Mason B. D., Hartkopf W. I., Gies D. R., Henry T. J., Helsel J. W., 2009, AJ , 137, 3358 CrossRef Search ADS   Matzner C. D., 2002, ApJ , 556, 302 CrossRef Search ADS   Matzner C. D., 2007, ApJ , 659, 1394 CrossRef Search ADS   McKee C. F., Offner S. S. R., 2010, ApJ , 716, 167 CrossRef Search ADS   McKee C. F., Ostriker E. C., 2007, ARA&A , 45, 565 CrossRef Search ADS   McKee C. F., Tan J. C., 2003, ApJ , 585, 850 CrossRef Search ADS   McKee C. F., Li P. S., Klein R. I., 2010, ApJ , 720, 1612 CrossRef Search ADS   Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS , 198, 7 CrossRef Search ADS   Miyoshi T., Kusano K., 2005, J. Comput. Phys. , 208, 315 CrossRef Search ADS   Mouschovias T. C., Spitzer L., Jr, 1976, ApJ , 210, 326 CrossRef Search ADS   Murray N. W., 2011, ApJ , 729, 133 CrossRef Search ADS   Murray D. W., Chang P., Murray N. W., Pittman J., 2015, MNRAS , 465, 1316 CrossRef Search ADS   Myers A. T., Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ , 735, 49 CrossRef Search ADS   Myers A. T., McKee C. F., Cunningham A. J., Klein R. I., Krumholz M. R., 2013, ApJ , 766, 97 CrossRef Search ADS   Myers A. T., Klein R. I., Krumholz M. R., McKee C. F., 2014, MNRAS , 439, 3420 CrossRef Search ADS   Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ , 703, 131 CrossRef Search ADS   Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI . Univ. Arizona Press, Tucson, AZ, p. 53 Ostriker E., Gammie C. F., Stone J. M., 1999, ApJ , 513, 259 CrossRef Search ADS   Padoan P., Nordlund Å., Kritsuk A. G., Norman M. L., Li P. S., 2007, ApJ , 661, 972 CrossRef Search ADS   Padoan P., Haugbølle T., Nordlund, Å., 2012, ApJ , 759, L27 CrossRef Search ADS   Padoan P., Pan L., Haugbølle T., Nordlund Å., 2016, ApJ , 822, 11 CrossRef Search ADS   Padoan P., Haugbølle T., Nordlund Å., Frimann S., 2017, ApJ , 840, 48 CrossRef Search ADS   Pan L., Padoan P., Haugbølle T., Nordlund Å., 2016, ApJ , 825, 30 CrossRef Search ADS   Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 2011, ApJ , 729, 72 CrossRef Search ADS   Peters T., Klaassen P. D., Mac Low M.-M., Schrön M., Federrath C., Smith M. D., Klessen R. S., 2014, ApJ , 788, 14 CrossRef Search ADS   Preibisch T., Balega Y., Hofmann K.-H., Weigelt G., Zinnecker H., 1999, New Astron. , 4, 531 CrossRef Search ADS   Price D. J., Bate M. R., 2008, MNRAS , 385, 1820 CrossRef Search ADS   Price D. J., Bate M. R., 2009, MNRAS , 398, 33 CrossRef Search ADS   Raghavan D. et al.  , 2010, ApJS , 190, 1 CrossRef Search ADS   Robertson B., Goldreich P., 2012, ApJ , 750, L31 CrossRef Search ADS   Rosen A. L., Krumholz M. R., McKee C. F., Klein R. I., 2016, MNRAS , 463, 2553 CrossRef Search ADS   Rosen A. L., Krumholz M. R., Oishi J. S., Lee A. T., Klein R. I., 2017, J. Comp. Phys. , 330, 924 CrossRef Search ADS   Salim D. M., Federrath C., Kewley L., 2015, ApJ , 806, L36 CrossRef Search ADS   Salpeter E. E., 1955, ApJ , 121, 161 CrossRef Search ADS   Santos-Lima R., de Gouveia Dal Pino E. M., Lazarian A., 2012, ApJ , 747, 21 CrossRef Search ADS   Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A , 410, 611 CrossRef Search ADS   Stone J. M., Ostriker E. C., Gammie C. F., 1998, ApJ , 508, L99 CrossRef Search ADS   Sur S., Federrath C., Schleicher D. R. G., Banerjee R., Klessen R. S., 2012, MNRAS , 423, 3148 CrossRef Search ADS   Tobin J. J. et al.  , 2016, ApJ , 818, 73 CrossRef Search ADS   Tricco T. S., Price D. J., Federrath C., 2016, MNRAS , 461, 1260 CrossRef Search ADS   Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., II, Howell L. H., Greenough J. A., 1997, ApJ , 489, L179 CrossRef Search ADS   Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., II, Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ , 495, 821 CrossRef Search ADS   Usero A. et al.  , 2015, AJ , 150, 115 CrossRef Search ADS   Wang P., Li Z.-Y., Abel T., Nakamura F., 2010, ApJ , 709, 27 CrossRef Search ADS   Wurster J., Price D. J., Bate M. R., 2017, MNRAS , 466, 1788 CrossRef Search ADS   Zhao B., Li Z.-Y., Nakamura F., Krasnopolsky R., Shang H., 2011, ApJ , 742, 10 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations