STARFORGE: the effects of protostellar outflows on the IMF

ABSTRACT

The initial mass function (IMF) of stars is a key quantity affecting almost every field of astrophysics, yet it remains unclear what physical mechanisms determine it. We present the first runs of the STAR FORmation in Gaseous Environments project, using a new numerical framework to follow the formation of individual stars in giant molecular clouds (GMCs) using the gizmo code. Our suite includes runs with increasingly complex physics, starting with isothermal ideal magnetohydrodynamics (MHD) and then adding non-isothermal thermodynamics and protostellar outflows. We show that without protostellar outflows the resulting stellar masses are an order of magnitude too high, similar to the result in the base isothermal MHD run. Outflows disrupt the accretion flow around the protostar, allowing gas to fragment and additional stars to form, thereby lowering the mean stellar mass to a value similar to that observed. The effect of jets upon global cloud evolution is most pronounced for lower mass GMCs and dense clumps, so while jets can disrupt low-mass clouds, they are unable to regulate star formation in massive GMCs, as they would turn an order unity fraction of the mass into stars before unbinding the cloud. Jets are also unable to stop the runaway accretion of massive stars, which could ultimately lead to the formation of stars with masses |${\gt}500\, \mathrm{M}_{\rm \odot }$|⁠. Although we find that the mass scale set by jets is insensitive to most cloud parameters (i.e. surface density, virial parameter), it is strongly dependent on the momentum loading of the jets (which is poorly constrained by observations) as well as the temperature of the parent cloud, which predicts slightly larger IMF variations than observed. We conclude that protostellar jets play a vital role in setting the mass scale of stars, but additional physics are necessary to reproduce the observed IMF.

MHD, turbulence, stars: formation, stars: jets, stars: luminosity function, mass function

1 INTRODUCTION

Star formation involves a large set of interconnected complex physical processes, including gravity, turbulence, magnetic fields, chemistry, and radiation (Girichidis et al. 2020). While each of these processes is necessary for a full picture of star formation, it is important to understand what role each of them plays and how they interact with each other.

Due to the complexity of the physics involved, star formation models often consider only a subset of the relevant physical processes to make the problem analytically (and even numerically) tractable. The simplest such model considers only the equations of isothermal hydrodynamics coupled to gravity, which models the dense, |${\sim}10\,\rm K$| interstellar medium (ISM) found in molecular clouds in our Galaxy (e.g. Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012). Recent numerical works have shown that the mass spectrum of collapsed fragments in such systems does not converge with numerical resolution (see e.g. Martel, Evans & Shapiro 2006; Kratter et al. 2010; Guszejnov, Krumholz & Hopkins 2016; Federrath, Krumholz & Hopkins 2017; Guszejnov et al. 2018b; Lee & Hennebelle 2018a), so additional physics must play a role.

Observations suggest that molecular clouds have significant support from magnetic fields (Crutcher 2012). In theoretical and numerical works, the addition of magnetic fields to isothermal star formation models have been shown to impose a resolution independent scale on the stellar mass spectrum (see e.g. Padoan et al. 2007; Padoan & Nordlund 2011; Haugbølle, Padoan & Nordlund 2018). While some of these studies claimed to reproduce the observed initial mass function (IMF), our recent study (Guszejnov et al. 2020) showed that, for clouds similar to giant molecular clouds (GMCs) in the Milky Way (MW), the mean stellar masses predicted by these magnetized, gravoturbulent models are an order of magnitude higher than observed (i.e. the mean stellar mass is |${\sim}4\, \mathrm{M}_{\rm \odot }$| in the simulations, while |${\sim}0.4\, \mathrm{M}_{\rm \odot }$| is observed). This study also found that stellar masses in isothermal magnetohydrodynamics (MHD) also increase with time and are sensitive to initial conditions (ICs; see analysis in section 4.2 of Guszejnov et al. 2020), leading to order of magnitude variations in the predicted characteristic scale of the IMF. Observations, however, have found the IMF to be near-universal within the MW, with variations in the IMF peak mass within a factor of <3 (see reviews of Bastian, Covey & Meyer 2010 and Offner et al. 2014, as well as analysis of Dib 2014).

Of course the ISM is not isothermal, one of the key assumptions of the above models is the gas can cool more rapidly than other relevant time-scales, making it effectively isothermal for this problem. This behaviour, however, is only a crude approximation of the real thermochemistry and radiative cooling, detailed calculations (e.g. Glover & Clark 2012) have shown significant temperature differences between low-density regions (⁠|${\sim}10^2\, \mathrm{cm}^{-3}$|⁠, |$T\sim 30\, \mathrm{K}$|⁠) and high-density regions where collapse occurs (⁠|${\sim}10^5\, \mathrm{cm}^{-3}$|⁠, |$T\sim 10\, \mathrm{K}$|⁠). Even at high densities, the isothermal assumption inevitably breaks down completely at high densities, when the cloud becomes opaque to its own cooling radiation, leading to an increase in temperature and thus a suppression of fragmentation (for the original idea, see Low & Lynden-Bell 1976; Rees 1976; for modern interpretations, see Colman & Teyssier 2020).

Another key feature of the simple models above is that they neglect feedback from the forming protostar and the stars that previously formed. These processes can dramatically affect the star formation process, as accreting protostars heat their surroundings (Offner et al. 2009; Krumholz 2011; Bate 2012; Myers et al. 2013; Guszejnov & Hopkins 2016; Guszejnov et al. 2016). Previously formed massive stars can also heat a large portion of their progenitor cloud and shut down star formation altogether (Grudić et al. 2018; Kim, Kim & Ostriker 2018; Li et al. 2019). The mass-loss of accreting protostars is dominated by high velocity bipolar outflows that can significantly affect their environment (see reviews of Frank et al. 2014; Bally 2016). These outflows are thought to be driven by highly collimated bipolar jets that entrain the ambient gas (Rosen & Krumholz 2020). These jets in turn are launched by MHD interactions between the protostar and the accretion disc (Shu et al. 1988; Pelletier & Pudritz 1992), with radiation pressure also contributing to their driving (Kuiper et al. 2010; Vaidya et al. 2011). These jets not only reduce the accretion rates of stars but also disrupt local accretion flows and drive turbulence on small scales (Matzner 2007a; Nakamura & Li 2007; Wang et al. 2010; Cunningham et al. 2011; Federrath et al. 2014a; Offner & Arce 2014; Offner & Chaban 2017; Murray, Goyal & Chang 2018).

Past work has shown that protostellar jets significantly reduce the global star formation rate (SFR) in a cloud (Hansen et al. 2012; Federrath et al. 2014a). Protostellar jets have been shown to play a role in setting the mass scale of stars, preventing ‘overaccretion’ from stars heating up their surroundings, thus preventing the gas from fragmenting and forming new stars (Krumholz, Klein & McKee 2012; Cunningham et al. 2018; Li, Klein & McKee 2018).

Simulations that take into account the above processes are necessary to understand the effects of each physical process, but so far these have generally been limited to simple physics or a very narrow range of cloud ICs. In this paper, we introduce the first results from the STAR FORmation in Gaseous Environments (STARFORGE) project.¹ These MHD simulations achieve a dynamic range in mass resolution that is an order of magnitude higher than any previous star cluster simulation, allowing us to simulate the detailed evolution of GMCs while following the formation of individual low-mass stars (see companion methods paper of Grudić et al. 2020b, henceforth referred to as Paper I). In this study, we perform and analyse a set of simulations with different ICs and levels of physics to identify the effects of non-isothermality and protostellar jets on the IMF.

We present our results in Section 3 with a focus on how the characteristic masses of sink particles (stars) change with the inclusion of additional physics and variations in the ICs (e.g. cloud temperature, surface density, level of turbulence). In Section 4, we introduce a simple toy model to explain the effects of protostellar outflows. The implications of these results as well as the potential role of further physics are discussed in Section 5. We summarize our conclusions in Section 6 and leave the details on how exactly our results vary with ICs to Appendix  A

2 METHODS

2.1 Physics

A full description and presentation of our methods including a variety of tests and algorithm details are given in a companion methods paper (Paper I), therefore we only briefly summarize them here.

2.1.1 Core physics

Similar to our previous studies of isothermal collapse with and without magnetic fields (Guszejnov et al. 2018b and Guszejnov et al. 2020, to the latter of which we will henceforth refer to as Paper 0), we simulate star-forming clouds with the gizmo code² (Hopkins 2015), using the Lagrangian meshless finite-mass method for MHD (Hopkins & Raives 2016), assuming ideal MHD (with the constrained gradient scheme of Hopkins 2016 to ensure that |$\nabla \cdot \mathbf {B}=0$|⁠).

Gravity is solved with an improved version of the Barnes–Hut tree method from Springel (2005) with high-order integration of sink particle trajectories to accurately follow multiple sink systems (see Paper I). Force softening is fully adaptive for gas cells (Price & Monaghan 2007). Sink particles (representing stars) have a fixed Plummer-equivalent softening radius of |$7.56\, \rm au$|⁠. We adopt the sink formation and accretion algorithm from Bate, Bonnell & Price (1995), while accurately accounting for thermal, magnetic, kinetic, and gravitational energies and angular momentum, again described in Paper I. As such we are able to follow the formation and evolution of binaries and multiples with separations larger than |${\sim}10\, \rm au$|⁠.

Once sinks form, they follow the protostellar evolution model from Offner et al. (2009), which is also used in the orion code. In this model, the protostar is treated as a collapsing polytrope: the collapse is divided into distinct phases during which the qualitative behaviour changes. These phases are ‘pre-collapse’, ‘no burning’, ‘core deuterium burning at fixed temperature’, ‘core deuterium burning at variable temperature’, ‘shell deuterium burning’ and ‘main sequence’. This module dynamically evolves stellar properties (e.g. radius, accretion, and internal luminosities) throughout the simulation. For details, see appendix B of Offner et al. (2009) and Paper I.

2.1.2 Thermodynamics

We compare simulations with two different thermodynamics modules. Our ‘isothermal’ simulations enforce an isothermal equation of state (EOS) with |$c_\mathrm{ s}=0.2\, \rm km\,s^{ -1}$| (effective gas temperature |$T\sim 10\, \mathrm{K}$|⁠). Our ‘non-isothermal’ or ‘cooling’ simulation runs utilize the radiative cooling and thermochemistry module presented in Hopkins et al. (2018) that contains detailed metallicity-dependent cooling and heating physics from |$T=10\,\mathrm{ to}\,10^{10}\,$| K, including recombination, thermal bremsstrahlung, metal lines (following Wiersma, Schaye & Smith 2009), molecular lines, fine structure (following Ferland et al. 2013), and dust collisional processes. The cooling module self-consistently solves for the internal energy and ionization state of the gas (see appendix B of Hopkins et al. 2018). The gas adiabatic index is calculated from a fit to density based on the results of Vaidya et al. (2015). Note that a constant dust temperature of |$T_\mathrm{dust}=10\, \mathrm{K}$| and a temperature floor of |$T_\mathrm{floor}=10\, \mathrm{K}$| are assumed here. As detailed in Paper I, this module does not explicitly evolve radiation hydrodynamics (RHD), but it does attempt to approximately capture the transition between optically thick and optically thin cooling regimes. It does so following Rafikov (2007) and modelling each gas cell as a plane-parallel atmosphere with optical depth to escape integrated using the TreeCol algorithm (Clark, Glover & Klessen 2012).

2.1.3 Protostellar jets

Protostars eject a significant portion of the accreting material in bipolar jets. To represent this process, we adopt the following jet model: each accreting protostar launches an f_w fraction of the accreting mass in bipolar jets along its rotational axis with a velocity of

$$\begin{eqnarray*} v_\mathrm{jet} = f_K \sqrt{G M_\ast /R_*}, \end{eqnarray*}$$

(1)

which is just f_K times the Keplerian velocity at the surface of the star, where the R_* stellar radius is evolved using the protostellar evolution model of Offner et al. (2009). Observations estimate the f_w mass loading parameter to be in the range of 0.1–0.4 (see review by Frank et al. 2014), while simulations found values of 0.1–0.6 (e.g. Seifried et al. 2012). The f_K velocity scaling parameter is not observed directly, however, f_Kf_w can be derived from the observed momentum injection rate by assuming a constant protostellar radius (see section 2.4 of Cunningham et al. 2011), which yields the constraint f_Kf_w ∼ 0.05–0.4. In our runs, we adopt f_w = 0.3 and f_K = 0.3, similar to the values used by Cunningham et al. (2011) and many other works, which puts f_wf_K in the middle of the observed range. It is useful to introduce the Γ momentum loading parameter  

$$\begin{eqnarray*} \Gamma = \frac{2}{3} \frac{f_w f_K}{(1-f_w)}, \end{eqnarray*}$$

(2)

which describes how the momentum output of the jets per unit accreted mass scales with these parameters (see Section 4 for a derivation).

The numerical implementation of jets is described in Paper I, briefly we spawn new gas cells around the sink particle and launch them along the sink particle’s angular momentum axis using the same angular distribution model as Cunningham et al. (2011), which corresponds to a vanishingly small opening angle. We find that the exact value of the opening angle has little effect on the results, provided that it is <1 (see Paper I for details). These gas cells are spawned in pairs (to conserve momentum and centre of mass exactly) and in mass quanta of Δm_jet = 0.1Δm, where Δm is the mass resolution element of our simulation, for which our fiducial value is |$\Delta m=10^{-3}\, \mathrm{M}_{\rm \odot }$|⁠, sufficient to predict the shape of the IMF in the stellar (≳0.1 M_⊙) mass range (see Paper I for resolution study).

In Paper I, we find that for Δm_jet/Δm ≤ 1 the sink mass spectrum is insensitive to our choice of Δm_jet, so we adopt Δm_jet/Δm = 0.1 in our simulations. We will show that the effects from the jet module are primarily determined by the Γ momentum loading parameter, see Section 3.2.3 for a details.

2.2 Initial conditions and parameters of clouds

2.2.1 Initial conditions

The main aim of the STARFORGE project is to identify the roles different physical processes play in star formation from the protostellar to the GMC scale (au to |$100\, \mathrm{pc}$|⁠). This investigation requires simulations of GMC scale clouds with individual star formation and progressively more complicated physics: starting with magnetized, isothermal gas (see Paper 0), then enabling gas thermodynamics without stellar feedback and finally adding protostellar outflows, see Table 1 for the different ‘rungs’ of this ‘physics ladder’. To explore the dependence of our results on ICs and simulation parameters, we also carry out a detailed parameter study. Note that the STARFORGE numerical framework can incorporate many other important feedback physics (e.g. radiative heating, winds, and supernovae: for methods see Paper I), which will be explored in future papers.

Table 1.

Labels used by throughout this paper to identify simulations with different physics. See Section 2.1 and Paper I for details on the individual physics modules.

Physics label	MHD	Thermodynamics	Protostellar jets
I_M	Ideal (M)	Isothermal (I)	Not included
C_M	Ideal (M)	ApproxRad (C)	Not included
C_M_J	Ideal (M)	ApproxRad (C)	Included (J)

We generate our ICs using makecloud.³ Unless otherwise specified our runs utilize ‘Sphere’ ICs, meaning that we initialize a spherical cloud (⁠|$T=10\, \mathrm{K}$|⁠, radius R_cloud, and mass M₀) with uniform density, surrounded by diffuse gas with a density contrast of 1000. The cloud is placed at the centre of a periodic 10R_cloud box. The initial velocity field is a Gaussian random field with power spectrum E_k ∝ k⁻² (Ostriker, Stone & Gammie 2001), scaled to the value prescribed by α_turb. The initial clouds have a uniform B_z magnetic field whose strength is set by the parameter μ. There is no external driving in these simulations.

We also run simulations using ‘Box’ ICs, similar to the driven boxes used in Federrath et al. (2014a) and Cunningham et al. (2018). These are initialized as constant density, zero velocity periodic cubic box with |$T=10\, \mathrm{K}$|⁠. This periodic box is then ‘stirred’ using the driving algorithm from Federrath et al. (2010) and Bauer & Springel (2012). This involves a spectrum of E_k ∝ k⁻² of driving modes in Fourier space at wavenumbers 1/2 - 1 the box size, with an appropriate decay time for driving mode correlations (t_decay ∼ t_cross ∼ L_box/σ_3D). This stirring is initially performed without gravity for five global freefall times |$(t_{\mathrm{ff}}\equiv \sqrt{\frac{3 \pi }{32 G \rho _0}})$|⁠, to achieve saturated MHD turbulence. The normalization of the driving spectrum is set so that in equilibrium the gas in the box has a turbulent velocity dispersion (σ_3D) that gives the desired |$\mathcal {M}$| and α_turb. We use purely solenoidal driving, which remains active throughout the simulation after gravity is switched on. We take the box side length L_box to give a box of equal volume to the associated Sphere cloud model, and thus define α_turb using the volume-equivalent R_cloud in equation (4). An important difference between the Sphere and Box runs is that in the case of driven boxes, the magnetic field is enhanced by a turbulent dynamo (Federrath et al. 2014b) and saturates at about α_B ∼ 0.1 (see Paper 0), so for Box runs the ‘pre-stirring’ magnetic field strength (defined by μ) does not directly specify the actual initial magnetic field strength when gravity is turned on (however, the ‘pre-stirring’ flux in the box will still affect the large-scale geometry of the magnetic field).

2.2.2 Parameters surveyed

To describe our ICs, we introduce several parameters, such as the 3D sonic Mach number  

$$\begin{eqnarray*} \mathcal {M}^2\equiv \left\langle ||\mathbf {v}_\mathrm{turb}||^2/c_{\rm s}^2\right\rangle , \end{eqnarray*}$$

(3)

where c_s is the gas sound speed and |$\mathbf {v}_\mathrm{turb}$| is the turbulent velocity field, while 〈...〉 denotes mass-weighted averaging. It is also useful to introduce the turbulent virial parameter α_turb, which measures the relative importance of turbulence to gravity, following the convention in the literature (e.g. Bertoldi & McKee 1992; Federrath & Klessen 2012),

$$\begin{eqnarray*} \alpha _{\mathrm{turb}}\equiv \frac{5 ||\mathbf {v}_\mathrm{turb}||^2 R_\mathrm{cloud}}{3 G M_\mathrm{0}} = \frac{5 \mathcal {M}^2 c_{\rm s}^2 R_\mathrm{cloud}}{3 G M_\mathrm{0}}, \end{eqnarray*}$$

(4)

where R_cloud and M₀ are the cloud (spherical-equivalent) radius and total mass, respectively. The relative importance of the magnetic field is commonly described by the normalized magnetic flux (or mass-to-flux ratio), which for a uniform magnetic field can be expressed as

$$\begin{eqnarray*} \mu = c_1\sqrt{\frac{-E_{\rm grav}}{E_{\rm mag}}}, \end{eqnarray*}$$

(5)

where E_grav and E_mag are the gravitational and magnetic energy (assuming a uniform initial field), respectively, while the normalization constant is c₁ ≈ 0.4. With this normalization μ = 1 corresponds to the critical point in the stability of a homogeneous sphere in a uniform magnetic field (Mouschovias & Spitzer 1976).

Clouds have several characteristics mass scales defined by ICs. Such a scale is the Jeans mass, representing the scale below which thermal pressure can prevent the gravitational collapse of a fluid element:

$$\begin{eqnarray*} M_{\rm Jeans}\equiv \frac{4 \pi c_{\rm s}^3}{3\sqrt{G^3\rho _0}}, \end{eqnarray*}$$

(6)

where ρ₀ is the density of the gas in the cloud. The initial turbulence also has a characteristic length-scale: the sonic length, L_sonic, on which the turbulent dispersion becomes supersonic. The corresponding mass scale is the sonic mass:

$$\begin{eqnarray*} M_{\rm sonic}\equiv \frac{c_{\rm s}^2 L_{\rm sonic}}{G}=\frac{c_{\rm s}^2 R_\mathrm{cloud}}{\mathcal {M}^2 G}, \end{eqnarray*}$$

(7)

where we used the supersonic linewidth–size relation (σ²(L) ∝ L). Another mass scale of an isothermal turbulent flow is the turbulent Bonnor–Ebert mass, the maximum gas mass that can support itself against its own self-gravity plus external pressure in post-shock compressed gas with |$\rho /\rho _0 \sim 1+\frac{1}{3}\mathcal {M}^2$| (Padoan, Nordlund & Jones 1997), which scales as

$$\begin{eqnarray*} M_{\rm BE}^{\rm turb}\sim 2 M_{\rm Jeans}\left(1+\frac{1}{3}\mathcal {M}^2\right)^{-1/2}=\frac{8 \pi c_{\rm s}^3}{3\sqrt{G^3\rho _0\left(1+\frac{1}{3}\mathcal {M}^2\right)}}. \end{eqnarray*}$$

(8)

Note that many other parameters are also used in the literature that can be expressed in terms of the ones introduced in this subsection (see section 2 in Paper 0 for how they relate to each other).

Table 2 shows the target parameters for the runs we present in this paper. The input parameters are the cloud mass M₀, size R₀, turbulent virial parameter α_turb, normalized magnetic flux μ, and initial temperature. Since our primary goal is to study the IMF in similar environments to the MW, we set up our fiducial runs as clouds between 2000 and |$2\times 10^5\, \mathrm{M}_{\rm \odot }$| that lie along a mass–size relation similar to observed GMCs in the MW (e.g. Larson 1981, specifically assuming |$\Sigma \equiv M_\mathrm{0}/ \pi R_\mathrm{cloud}^2 = 63 \mathrm{M}_{\rm \odot }\, \mathrm{pc}^{-2}$|⁠). These clouds are marginally bound (α_turb = 2) and start out at |$T=10\, \mathrm{K}$|⁠, the temperature of the cold ISM. For the initial magnetization, we assumed −E_mag/E_grav = 0.01, which translates to μ = 0.4 (note that this choice has little effect on the results, see Section 3.2.4). For the treatment of protostellar jets, we use our fiducial parameters of f_w = 0.3 and f_K = 0.3 (see Section 2.1.3). Since observed clouds can deviate from the observed linewidth–size relation (Heyer et al. 2009), we also simulate several clouds with different surface densities, turbulence, and magnetic support. Note that for these studies, we use clouds with a |$2\times 10^4\, \mathrm{M}_{\rm \odot }$| initial mass (M2e4), due to the high computational cost of larger runs. Also, since most MW GMCs achieve a star formation efficiency (SFE = M_⋆/M₀) of 1 per cent–10 per cent over their lifetime (see Krumholz 2014 for a discussion, and note that some clouds have <1 per cent, see Federrath & Klessen 2013), we restrict our analysis to times when SFE is below 10 per cent.

Table 2.

ICs of clouds used in our runs, with M₀, R_cloud, α_turb, μ, and T₀ being the initial cloud mass, size, virial parameter, mass to magnetic flux ratio, and temperature, respectively (note that in all our runs T_floor = T₀, see Section 2.1.2). We also report the initial 3D sonic Mach number |$\mathcal {M}$|⁠, thermal virial parameter α_th, total virial parameter α, Alfvén Mach number |$\mathcal {M}_{\rm A}$|⁠, plasma β, magnetic virial parameter α_B, as well as the relative Jeans, sonic, and magnetic mass scales (see section 2 in Paper 0 for definitions). Note that the parameters in this table apply to both Box and Sphere runs as they are set up to have identical initial global parameters, with L_box being the box size for Box runs and R_cloud being the cloud radius for Sphere runs. Note that Box runs have slightly different initial parameters (e.g. Mach number, virial parameter) due to the non-exact scaling of the driving, so the values shown here are the target values. Many of the above clouds have been run at different mass resolutions as part of the resolution study in Paper I, in the table we note for each the highest resolution that was run (Δm mass resolution and Δx_J minimum resolved Jeans length; see section 2 in Paper I for details).

	Input parameters							Derived parameters								Highest resolution run
Cloud label	M0 (M⊙)	Rcloud (pc)	Lbox (pc)	αturb	μ	T0 (K)	|$\mathcal {M}$|	αth	α	|$\mathcal {M}_{\rm A}$|	β	αB	|$\frac{M_{\rm Jeans}}{M_0}$|	|$\frac{M_{\rm sonic}}{M_0}$|	|$\frac{M_{\Phi }}{M_0}$|	M0/Δm	ΔxJ (au)
MW cloud analogues
M2e2	2 × 102	1		2	4.2	10	5	0.02	2.04	10	7.8	0.02	6 × 10−2	7 × 10−3	0.1	2 × 105	36
M2e3	2 × 103	3	4.8	2	4.2	10	9.3	0.02	2.04	10	2.3	0.02	1 × 10−2	6 × 10−4	0.1	2 × 107	3.6
M2e4	2 × 104	10	16	2	4.2	10	16	0.008	2.03	10	0.78	0.02	3 × 10−3	7 × 10−5	0.1	2 × 107	36
M2e5	2 × 105	30		2	4.2	10	29	0.002	2.02	10	0.23	0.02	5 × 10−4	7 × 10−6	0.1	2 × 108	36
Parameter variation tests
M2e4_a4	2 × 104	10		4	4.2	10	22.6	0.008	4.03	10	0.78	0.02	3 × 10−3	4 × 10−5	0.1	2 × 107	36
M2e4_a1	2 × 104	10		1	4.2	10	11.3	0.008	1.01	10	0.78	0.02	3 × 10−3	1 × 10−4	0.1	2 × 107	36
M2e4_a05	2 × 104	10		0.5	4.2	10	8	0.008	0.51	10	0.78	0.02	3 × 10−3	5 × 10−5	0.1	2 × 107	36
M2e4_R3	2 × 104	3		2	4.2	10	29	0.002	2.02	10	0.23	0.02	5 × 10−4	7 × 10−6	0.1	2 × 107	36
M2e4_R30	2 × 104	30		2	4.2	10	9.3	0.02	2.04	10	2.3	0.02	1 × 10−2	6 × 10−4	0.1	2 × 107	36
M2e4_mu13	2 × 104	10		2	13	10	16	0.008	2.01	31	7.8	0.002	3 × 10−3	7 × 10−5	0.04	2 × 107	36
M2e4_mu1.3	2 × 104	10		2	1.3	10	16	0.008	2.21	3.1	0.078	0.2	3 × 10−3	7 × 10−5	0.4	2 × 107	36
M2e4_T30	2 × 104	10		2	4.2	30	9.3	0.024	2.04	10	2.3	0.02	1 × 10−2	6 × 10−4	0.1	2 × 107	36
M2e4_T60	2 × 104	10		2	4.2	60	6.6	0.048	2.07	10	4.6	0.02	5 × 10−2	2 × 10−3	0.1	2 × 107	36

3 RESULTS

We carried out a suite of simulations using the ICs from Table 2 and the different physics combinations from Table 1.

All simulations develop filaments, clumps, and cores, and begin global collapse (see Fig. 1 for the case with protostellar jets). In the runs with protostellar jets, once star formation begins jets disrupt the flow around newly formed stars (see Fig. 2), reducing their accretion rates and allowing new stars to form. In the following subsections, we investigate different aspects of star formation with different physics enabled.

Figure 1.

Surface density maps for M2e5_C_M_J with M₀/Δm = 2 × 10⁸ initial gas cells (see Tables 1 and 2) at different times. The colour scale is logarithmic and the circles represent sink particles (stars) that form in high-density regions where fragmentation can no longer be resolved, their size increasing with mass as well as their colour changing from red (⁠|$M\sim 0.1\, \mathrm{M}_{\rm \odot }$|⁠) to blue (⁠|$M\sim 10\, \mathrm{M}_{\rm \odot }$|⁠). This simulation resolves a dynamic range from |${\sim}\!\mathrm{50\, pc}$| down to |${\sim}\!\mathrm{30\, au}$|⁠.

Figure 2.

Zoomed-in surface density maps for a medium-sized cloud with protostellar jets enabled (M2e4_C_M_J). Symbols and colour maps are similar to Fig. 1. The final image (top, right) shows the kinetic energy weighted surface density (weight = m(1 + [v/v₀]²), where |$v_0=1\, \mathrm{km\,s^{ -1}}$|⁠), as well as the local velocity field (white arrows), whose length scales with velocity, to highlight the jet (which has high velocity but low density, making it challenging to see in surface density maps).

3.1 Star formation history

Fig. 3 shows the star formation history of several clouds with identical ICs (M2e4) but with different physics modules and turbulent driving (see Sphere versus Box ICs in Section 2.2.1). For the Sphere runs, we find that the SFE in all cases follows a similar broken power law, which starts linearly (note that this ‘early time’ slope is potentially sensitive to the definition of the time zero-point) and transitions to SFE(t) ∝ t³ at later times, similar to the findings of Paper 0 for the isothermal case and other simulations without turbulent driving from the literature (e.g. Myers et al. 2014). Note that while protostellar jets do reduce the SFR, their net effect is only a shift in the curve, delaying the onset of the cubic regime from roughly 10 per cent of the freefall time to about 20 per cent. The results for the Box runs are qualitatively similar, but their SFRs are slower: they scale as SFE ∝ t², similar to previous results with driven turbulent boxes (e.g. Federrath & Klessen 2012; Murray et al. 2015, 2018).

Figure 3.

Left: Evolution of the SFE (SFE(t) = ∑M_sink(t)/M₀) as function of time for a subset of runs with M2e4 runs but different physics (see Tables 1 and 2). A run with Box ICs is also included for comparison. Note that t = 0 is set to the start of star formation and that the right-hand panel (α_turb) uses a linear x-axis for time t, as opposed to a log axis, with negative values t < 0 representing time before the first sink forms. The SFE rises as a broken power law of time and reaches about 10 per cent in 1–2 freefall times (⁠|$t_{\mathrm{ff}}=\sqrt{\frac{3 \pi }{32 G \rho _0}}$|⁠, which is about 3.5 Myr in these runs). Note that the overall shape of the SFE is unchanged by enabling feedback physics, however, jets shift the curve, effectively delaying the star formation process by a factor of 2 in t/t_ff. Meanwhile, the slope of the curve is sensitive to the type of IC used, we find SFE ∝ t³ for Sphere ICs and SFE ∝ t² for Box ICs. Middle: Number of sink particles in the simulations as a function of time. Non-isothermal thermodynamics suppresses the formation of low-mass sink particles while the additional turbulence on small scales enhances it. Switching to Box ICs leads to a shallower exponent, similar to the SFE case above. Right: Maximum sink mass in the same simulations as a function of time. In all cases, the maximum mass asymptotes to ∝ t³ once massive stars have formed.

Fig. 3 also shows the number of sink particles, N_sink, over time. For most runs, N_sink follows a similar trend to the SFE, which produces a roughly time-invariant mean sink mass (see Fig. 5). Note that even though switching to driven turbulence (Box IC) reduces the SFR, the mean sink mass remains roughly similar (SFE ∝ N_sink). This implies that the sink mass distribution (IMF) in the simulation is determined by local physics (e.g. jets) instead of large-scale boundary conditions (i.e. turbulent driving spectrum).

We find that the maximum sink mass increases over time, starting as a |$M_\mathrm{max}\propto \sqrt{t}$| power law, which steepens to M_max ∝ t³ once massive sinks (stars) form, as they undergo runaway accretion. This plays out qualitatively similarly in all runs here, regardless of physics or turbulent driving. The main effect of protostellar jets is that they reduce the maximum sink mass by about an order of magnitude at fixed total sink mass in the simulation.

Fig. 4 shows that the inclusion of protostellar jets (C_M_J) can lead to the disruption of the parent cloud and subsequently preventing the formation of new stars. In more massive clouds (>10⁴ M_⊙, similar to MW GMCs), protostellar jets show no sign of arresting star formation before the SFE exceeds |${\sim}10{{\ \rm per\ cent}}$|⁠. Note that SFE is challenging to measure observationally, but observed clouds in the range of sizes and masses we have simulated are generally believed to have a typical SFE of only a few  per cent (Lee, Miville-Deschênes & Murray 2016; Vutisalchavakul, Evans & Heyer 2016; Grudić et al. 2019b; Kruijssen et al. 2019; Chevance et al. 2020).

Figure 4.

The evolution of the SFE (left), number of sink particles N_sink (middle), and turbulent virial parameter α_turb (α_turb = 2 being equivalent marginal gravitational boundedness) as function of time for a subset of runs with protostellar jets enabled (C_M_J) that have clouds with increasing initial masses and a constant surface density (⁠|$\Sigma \approx 60\, \mathrm{M}_{\rm \odot }\,\mathrm{pc}^{-2}$|⁠) similar to MW GMCs (M2e2–M2e5, see Table 2). Note that here t = 0 is set to the start of star formation. We find that after reaching a sufficiently high SFE, protostellar jets are able to unbind low-mass clouds the time of which is marked with a vertical dashed line. After this jets are able to quench star formation, almost completely stopping the formation of new sink particles.

3.2 Sink mass distribution (IMF)

Sink particles represent stars (or systems with separations below the resolution limit) in our simulations, so we use their mass spectrum as an analogue of the IMF. Since it is possible for the sink mass spectrum (IMF) not to converge numerically at the lowest masses, while still converging on shape at higher masses or providing characteristic mass scales, we investigate the effects of different physics on both the various characteristic mass scales and the shape of the sink mass spectrum.

3.2.1 Characteristic mass of stars

A common issue in numerical simulations is that the low-mass end of the sink mass spectrum is sensitive to numerical resolution and simulations often have a large number of very low mass objects near their resolution. While in most cases these objects represent a vanishingly small fraction of the total sink mass (see Paper 0 for an example and Guszejnov et al. 2018b for a counterexample), their large number skews the mean and median sink masses. Adopting the mass-weighted median mass of sinks M₅₀ as the characteristic mass scale mitigates this effect (see Krumholz et al. 2012 and Paper 0), but this choice makes the mass scale overly sensitive to the most massive sinks that can undergo runaway accretion (see Fig. 5).

Figure 5.

The evolution of the number-weighted mean (M_mean = ∑M_sink/N_sink, left), number-weighted median (defined such that N_sink(M > M_med) = N_sink/2, centre), and mass-weighted median (M₅₀, the mass scale above which half the total sink mass resides, right) sink mass as a function of SFE for the runs shown in Fig. 3. We also show with a shaded region the 95 per cent confidence interval for these values if one sampled the Kroupa (2002) IMF at the current SFE value in the cloud. Protostellar jets reduce sink masses and bring all three mass scales closer to those of the observed IMF, however, M₅₀ increases with time, diverging from observations at higher SFE values.

Fig. 5 shows the evolution of the mean and median sink masses along with that of M₅₀ as a function of SFE. For runs without jets, we find that the mean sink mass and M₅₀ both increase with time due to the runaway accretion of the massive sinks. Note the introduction of non-isothermal physics has little effect on the three mass scales and without jets they are all significantly larger than those observed in the MW. The introduction of jets allows low-mass stars to form again, such that the mean mass is roughly time invariant while all three mass scales are near their observed values. But as star formation progresses (SFE > 1 per cent), we find that all simulations show an increasing trend in M₅₀ due to the runaway accretion of massive sinks, similar to the M₅₀ ∝ SFE^1/3 scaling found in the isothermal case in Paper 0. Switching to Box ICs has little effect on the evolution of M₅₀ or the mean sink mass, except for a delay in the runaway accretion of massive stars. For the median mass, however, turbulent driving appears to suppress the formation of very low mass stars.

At our fiducial resolution both the mean sink mass and M₅₀ are insensitive to numerical resolution (see Paper I). We also find that the mean sink mass exhibits a nearly time invariant trend between 1 per cent and 10 per cent SFE in most simulations (see Figs 5, A1, and A2), while M₅₀ increases with time in nearly all cases, so we adopt it as a proxy for the characteristic scale of the IMF for the remainder of the paper.

3.2.2 The IMF

While the various characteristic masses provide some information on the sink mass distribution, a holistic view of the IMF is necessary to understand the effects of each physical process. Fig. 6 shows the mass distribution of sink particles at 5 per cent SFE, which we will use as a proxy for the IMF. We find that the addition of non-isothermal physics alone has little effect on the IMF⁴ leaving the IMF top-heavy (see Paper 0). We find that the inclusion of protostellar jets dramatically changes the distribution, shifting the turnover to mass scales comparable to that observed in the MW. Switching to Box ICs does not qualitatively change the IMF apart from slightly suppressing the formation of very low mass objects. Driven turbulence also delays the runaway accretion of massive stars; that is why the IMF is not yet top-heavy at 5 per cent SFE for the Box run in Fig. 6.

Figure 6.

Distribution of sink particle masses measured in each simulation at 5 per cent SFE (SFE = ∑M_sink/M₀) for the runs shown in Fig. 3. We also show the Salpeter (1955), Kroupa (2002), and Chabrier (2005) fitting functions for the IMF. Again, jets greatly improve the agreement with the IMF for our chosen jet parameters (see Appendix  A for results with different values).

3.2.3 Role of jet momentum loading

Since jets have a dramatic effect on the IMF (see Fig. 6), we examine how our results depend on the f_w and f_K jet parameters (see Section 2.1.3). Fig. 7 shows the results of varying these parameters for an M2e4_C_M_J run. We find that the evolution of the cloud and the sink mass spectra depend primarily on Γ = f_wf_K/(1 − f_w), which determines the momentum loading of the jets (e.g. the results obtained for f_w = 0.1 and f_K = 1 are very similar to the results for our fiducial f_w = 0.3 and f_K = 0.3). Furthermore, we find that the number of sink particles appears to be insensitive to the values of the jet parameters, but there is a factor of 2–3 difference between jet and non-jet runs (see Fig. 3).

Figure 7.

SFE (top left), number of sink particles N_sink (top centre), mass-weighted median sink mass M₅₀ (top right), median sink mass M_med (bottom left), mean sink mass M_mean (bottom centre), and the maximum sink mass M_max (bottom right) in M2e4_C_M_J runs with different jet mass loading (⁠|$f_w=\dot{M}_\mathrm{jet}/\dot{M}_\mathrm{acc}$|⁠) and velocity scaling (⁠|$f_K=v_\mathrm{jet}/\sqrt{G M_\ast /R_\ast }$|⁠) parameters for protostellar jets (see Section 2.1.3). The shaded region represents the minimum and maximum values at a given time over three different initial turbulence realizations. Dashed lines indicate statistics of the Kroupa (2002) IMF. The IMF shape is sensitive to jet model parameters, particularly to Γ ∝ f_wf_K, which determines the momentum loading of jets.

3.2.4 Sensitivity to initial conditions

We investigate the sensitivity of the predicted M_mean in our C_M_J runs (as these produce the most realistic IMF) to ICs by systematically varying cloud parameters around our M2e4 reference cloud, as shown in Table 2. We also vary the momentum loading of protostellar jets. Using a least-squares fit for M_mean as a function of each varied parameter (at fixed 4 per cent SFE), we obtain

$$\begin{eqnarray*} M_{\rm mean}\propto \Gamma ^{-0.65\pm 0.15}\, c_{\rm s, min}^{2.5\pm 0.5}\, \Sigma ^{-0.3\pm 0.1}\, \alpha _{\mathrm{turb}}^{0.15\pm 0.19}\, M_0^{-0.12\pm 0.07}, \end{eqnarray*}$$

(9)

which can also be expressed as

$$\begin{eqnarray*} M_{\rm mean}\propto \Gamma ^{-0.65\pm 0.05}\, c_{\rm s, min}^{2.5\pm 0.5}\, \rho ^{-0.2\pm 0.07}\, \alpha _{\mathrm{turb}}^{0.15\pm 0.19}\, M_0^{-0.22\pm 0.10}, \end{eqnarray*}$$

(10)

where Γ is the momentum loading of jets (see equation 2 and Section 4), c_s,min is the adiabatic sound speed at the T_floor temperature floor, while ρ, α_turb, and M₀ are the initial density, virial parameter, and mass of the parent cloud, see Appendix  A for a detailed presentation of the results and the derivation of the exponents and their errors. Assuming a mass–size relation similar to that in the MW (corresponding to |$\Sigma \sim 60\, \mathrm{M}_{\rm \odot }\,\mathrm{pc}^{-2}$|⁠; see Larson 1981), we can simplify equation (9) as

$$\begin{eqnarray*} M_{\rm mean}\propto \Gamma ^{-0.65\pm 0.15}\, c_{\rm s, min}^{2.5\pm 0.5}\, \alpha _{\mathrm{turb}}^{0.15\pm 0.19}\, M_0^{-0.12\pm 0.07}. \end{eqnarray*}$$

(11)

Equations (9)–(11) imply that the number-weighted mean sink mass for clouds is only weakly dependent on most cloud properties and is primarily set by the jet momentum loading factor Γ and the c_s,min sound speed at the cloud temperature floor.

3.3 Effects of jets on the accretion flow

Fig. 8 shows that protostellar jets dramatically change the accretion history of sink particles. Their effects are more than just removing some fraction of the accreted gas (i.e. multiplying the accretion rates by a constant factor), as the ejected jets entrain local gas and thus disrupt the accretion flow. This dramatically reduces the mass flux towards the sink particles on |${\lt}0.1\, \mathrm{pc}$| scales, slowing their growth (but not preventing the runaway accretion of massive stars, see Fig. 10). The nature of jet feedback is also showcased by Fig. 9. Looking at the surface density map, we find that the large-scale (>0.1 pc) gas structure is almost identical between runs with and without jets (C_M and C_M_J), but the sink mass spectrum is dramatically different (see Fig. 6). This is due to the dramatic effect jets have on gas kinematics, disrupting accretion flows around stars and creating outflows that extend up to |${\sim}10\, \mathrm{pc}$| in scale (bottom row of Fig. 9).

Figure 8.

Top: Evolution of the mass of a sink particles that form within the first Myr of star formation in a run with and without jets (M2e4_C_M_J and M2e4_C_M). Note that we exclude sink particles that do not reach |$0.5\, \mathrm{M}_{\rm \odot }$| within the first Myr of their lifetime to avoid including ‘failed sinks’ that form around a massive sink particle that prevents them from growing. The addition of protostellar jets greatly reduces the growth rate of the sink, in some cases shutting down accretion. Bottom: The average radial velocity in various radial shells around the same sink particles as above, 100 kyr after their formation. The solid line shows the value averaged over sinks, while the shaded area shows the interquartile range of values. On pc scales and above, the two runs are essentially identical as the mass flux is set by the global collapse of the cloud. On |${\lt}0.1\, \mathrm{pc}$| scales, jets disrupt the local accretion flow, dramatically reducing the magnitude of the mean radial velocity, which in turn leads to a dramatically reduced mass flux.

Figure 9.

(Top row) Surface density maps for a very small (M2e2, left) and a large cloud (M2e5, right) with and without jets (C_M and C_M_J). Colours (encoding projected gas surface density) and symbols (representing sink particles) are similar to Fig. 1. While the gas and star distribution on larger scales is almost identical between the runs, the masses of sink particles are different as illustrated by the relative sizes and colours of the circles. (Bottom row) Same as above, but now shown with a colour map that encodes the 1D line-of-sight velocity dispersion (increasing from purple (⁠|$0.1\,\rm km\, s^{-1}$|⁠) to orange (⁠|$10\,\rm km\, s^{-1}$|⁠) and encodes surface density information in lightness (lighter is denser). While the protostellar jets are almost invisible in surface density maps (and hence dust and CO maps as well), due to their low density, they stand out in these kinematic maps.

4 A SIMPLE MODEL FOR THE CHARACTERISTIC MASS SCALE SET BY JETS

In this section, we present a simple, plausible (but not necessarily unique) model that may explain the scaling of the mean sink (stellar) mass in our simulations (see equation 10 and Appendix  B). The jet model in our simulation launches an f_w fraction of the accreted mass at f_K times the Keplerian velocity (see equation 1 and Section 2.1.3). The total momentum output by the jet per unit time is therefore

$$\begin{eqnarray*} \dot{P_J} = v_\mathrm{jet} \dot{M}_\mathrm{jet} = f_K f_w \dot{M}_\mathrm{acc} \sqrt{G M_\ast /R_*}, \end{eqnarray*}$$

(12)

where |$\dot{M}_\mathrm{acc}$| is the mass accretion rate. Let us further assume that |$\dot{M}_\mathrm{acc}=\mathrm{const.}$| and R_* = const.⁵ so that |$M_\ast =\dot{M}_\mathrm{acc} t$|⁠. This will simplify the above equation to

$$\begin{eqnarray*} \dot{P_J} = f_K f_w \left(1-f_w\right)^{1/2} G^{1/2}\dot{M}_\mathrm{acc}^{3/2} R_*^{-1/2}t^{1/2}, \end{eqnarray*}$$

(13)

which we can integrate to get the total amount of momentum injected by jets over time t. Replacing |$t=M_\ast /\dot{M}_\mathrm{acc}$|⁠, we obtain

$$\begin{eqnarray*} P_J(M_\ast) = \frac{2}{3}\frac{f_w f_K}{1-f_w} G^{1/2}M_\ast ^{3/2} R_*^{-1/2} = \Gamma G^{1/2}M_\ast ^{3/2} R_*^{-1/2}, \end{eqnarray*}$$

(14)

where we have also used the Γ momentum loading parameter from equation (2).

Let us assume that this protostar forms in a cloud of uniform density (ρ) that is much larger than the jet (i.e. GMC) and that there is a spherical gas reservoir of mass (M_g) around the protostar that would eventually be accreted on to it without feedback. Let us also assume that protostellar jets are the only feedback process and that all the momentum injected by jets is deposited uniformly in the mass reservoir. The reservoir will become unbound if enough momentum is injected for its gas to reach escape velocity |$v_\mathrm{esc}\sim \sqrt{G (M_\ast +M_g)/R}$|⁠, where R = (M_g/(4π/3ρ))^1/3 is the radius of the reservoir. This means

$$\begin{eqnarray*} \left(\frac{P_J(M_\ast)}{M_g}\right)^2 = \frac{G(M_g+M_\ast)}{R}, \end{eqnarray*}$$

(15)

so

$$\begin{eqnarray*} \Gamma ^2 G M_\ast ^{3} R_*^{-1} M_g^{-2} = G(M_g+M_\ast)M_g^{-1/3}\left(\frac{4\pi \rho }{3}\right)^{1/3}. \end{eqnarray*}$$

(16)

Assuming M_g ≫ M_* and a fixed ρ, we can solve for the SFE of the gas reservoir before it becomes unbound:

$$\begin{eqnarray*} \frac{M_\ast }{M_g} = \Gamma ^{-2/3} \left(\frac{4\pi }{3}\right)^{1/9} \rho ^{1/9} M_g^{-1/9} R_*^{1/3}. \end{eqnarray*}$$

(17)

Substituting in typical values for GMCs this becomes

$$\begin{eqnarray*} \frac{M_\ast}{M_g} &=& 0.02 \left(\frac{\Gamma }{0.085}\right)^{-2/3} \left(\frac{n}{100\, \mathrm{cm^{-3}}}\right)^{1/9} \left(\frac{M_g}{\mathrm{M}_{\rm \odot }}\right)^{-1/9}\nonumber\\ &&\times \,\left(\frac{R_*}{2 \mathrm{R_\odot }}\right)^{1/3}, \end{eqnarray*}$$

(18)

where we used the n number density instead of ρ for convenience, as well as our fiducial parameters of f_w = 0.3 and f_K = 0.3 to normalize Γ.

To get the mass scale of the IMF, we formulate an ansatz for the M_g gas reservoir mass. Possible candidates are the Jeans, sonic, and turbulent Bonnor–Ebert masses. Based on our scaling results from Section 3.2.4 and Appendix  A, we know that the characteristic mass scales of the IMF (M₅₀ and the M_mean) both show weak dependence with the cloud virial parameter α_turb, consistent with an exponent between 0 and 1/3. Of these mass scales |$M_{\rm sonic}\propto \mathcal {M}^{-2} \propto \alpha _{\mathrm{turb}}^{-1}$| and |$\mathrm{M_{BE}}\propto \mathcal {M}^{-1}\propto \alpha _{\mathrm{turb}}^{-1/2}$|⁠, while M_Jeans is independent, so we adopt M_g = M_Jeans in this model. Plugging it into equation (18), we get

$$\begin{eqnarray*} M_\ast &=& 0.12 \mathrm{M}_{\rm \odot }\left(\frac{\Gamma }{0.085}\right)^{-2/3} \left(\frac{c_{\rm s, min}}{200\, \mathrm{m\,s^{ -1}}}\right)^{24/9} \left(\frac{n}{100\, \mathrm{cm^{-3}}}\right)^{-1/3}\nonumber\\ &&\times \,\left(\frac{R_*}{2 \mathrm{R_\odot }}\right)^{1/3}. \end{eqnarray*}$$

(19)

We find that the parameters of this model all fall within the uncertainty thresholds we found by fitting in equation (10). Fig. 10 shows how that the mass scales commonly used in the literature (M_Jeans, M_sonic, and |$M_{\rm BE}^{\rm turb}$|⁠) are all correlated with the mean sink mass M_mean in our simulations with jets. Meanwhile, our toy model from equation (19) provides a surprisingly good fit to the results with only a few outliers. Of course, it is only a toy model and makes several strong assumptions (e.g. constant R_*, M_g ∼ M_Jeans). Essentially, in this model M_* is set by the characteristic reservoir mass (i.e. core mass) with a feedback efficiency factor that varies only weakly with gas properties and primarily depends on the jet momentum loading as M_*/M_g ∝ Γ^−2/3.

Figure 10.

Comparison of the mean sink mass M_mean measured in different simulations using C_M_J physics with different initial cloud conditions (see Table 2) with the initial Jeans mass M_Jeans (equation 6), sonic mass M_sonic (equation 7), turbulent Bonnor–Ebert mass |$M_{\rm BE}^{\rm turb}$| (equation 8) as well as the M_theory mass scale predicted by our toy model (equation 19) and the result M_fit obtained by arbitrary least-squares fit marginalized over these parameters (equation 9). A dashed line shows the best linear fit between these mass scales and the mean sink mass. Note that we chose 4 per cent as the reference SFE as some runs never reach higher values as jets disrupt the cloud. The errors are estimated by bootstrapping: we resample the sink mass distribution at fixed total sink mass and calculate the 95 per cent confidence interval of M₅₀ over these realizations, which we denote with errorbars. Note that simulation runs with variable momentum loading are only shown in the bottom row and are slightly offset to make the plot easier to parse.

We stress this particular model is not unique and should not be overinterpreted. For example, the time integral above implies jets accelerate gas slowly on time-scales long compared to core dynamical times. If this is not true, the criterion for unbinding gas becomes |$\dot{P}_{J} \gt |{\bf F}_{\rm grav}|$| where F_grav is the gravitational force. If we assume also (unlike our derivation above) that the protostar is sufficiently massive that its gravity is important in the envelope so |$\dot{M}_{\rm acc}$| follows a Bondi-like scaling, then (following similar logic as before) the core would be unbound when |${M}_{\ast } \sim c_{\mathrm{ s}}^{2}\, \Gamma ^{-2/3}\, G^{-1}\, R_{\ast }^{1/3}\, (M_{g}/\rho)^{2/9}$|⁠. This gives a comparably good fit to the scaling we empirically extract from the simulations, but without reference to the Jeans mass (in fact it depends quite weakly on whatever physics sets M_g). Instead, the c_s dependence in this model comes from the fact that higher c_s (all else equal) slows accretion and therefore reduces the instantaneous strength of feedback. What is robust is that in any momentum-feedback-regulated model, we expect M_* to scale inversely with Γ. We also note that the above toy model is not unique to protostellar jet and can easily be adapted to derive the characteristic stellar mass for other feedback mechanisms.

5 DISCUSSION

In isothermal MHD runs, Paper 0 found that magnetic fields impose a well-defined characteristic mass (related to the initial sonic mass) on the sink mass distribution that is insensitive to numerical resolution (unlike the non-magnetized isothermal hydrodynamics case; see Guszejnov et al. 2018b), similar to the results of Haugbølle et al. (2018). Above this mass scale, the sink mass distribution roughly follows a dN/dM ∝ M⁻² trend, similar to the observed IMF (Salpeter 1955), and likely arising as a general consequence of scale free physics on this dynamic range (Guszejnov, Hopkins & Grudić 2018a). Paper 0 found that this characteristic mass of stars is an order of magnitude higher than what is observed, and is sensitive to ICs in a way that violates the apparent near-universality of the IMF in the MW (Offner et al. 2014; Guszejnov, Hopkins & Ma 2017).

5.1 Role of non-isothermal thermodynamics

Isothermality is often assumed in star formation theories and simulations due to the highly efficient cooling of molecular gas (Girichidis et al. 2020), even though there is a significant scatter in the gas temperature with a clear density dependence (see Glover & Clark 2012). At high densities, the isothermality assumption must eventually break down, allowing for the formation of hydrostatic cores (Larson 1969) that are the progenitors of protostars. This transition from near-isothermal to adiabatic behaviour was originally proposed to be responsible for setting the peak of the IMF (see Low & Lynden-Bell 1976; Rees 1976), but the corresponding mass scale (∼0.008 M_⊙) was too low to explain observations. The idea has recently been revived by taking into account the tidal screening effect around the first Larson core (Colman & Teyssier 2020), which increases the relevant mass scale to be comparable to the observed IMF peak.

It is important to note that most of these simulations have been run on non-magnetized clouds, so the only unique mass scale in the sink mass spectrum arises from non-isothermal physics at high densities.⁶ Including magnetic fields, however, in MW-like cloud conditions shifts the turnover mass of the IMF to much larger |${\gt}20\, \mathrm{M}_{\rm \odot }$| scales (see Fig. 6 and Paper 0). Thus, for MW-like clouds, gas thermodynamics (i.e. the opacity limit) do not set the ‘mean’ characteristic or turnover mass scale of the IMF (which is of order ∼M_⊙), their effects are likely limited to the lowest mass scales of the IMF (⁠|${\lt}0.1\, \mathrm{M}_{\rm \odot }$|⁠, see Figs 5 and 6).

5.2 Role of protostellar jets

Previous work has shown that protostellar jets can expel a significant portion of accreting material, directly reducing stellar masses (e.g. Federrath et al. 2014a; Offner & Chaban 2017) and potentially driving small-scale turbulence (e.g. Nakamura & Li 2007; Wang et al. 2010; Offner & Arce 2014; Offner & Chaban 2017; Murray et al. 2018). We do find that jets disrupt the local accretion flow, which greatly changes gas dynamics on |${\lt}0.1\, \mathrm{pc}$| scales, but this has little effect on the global evolution of a massive GMC. Previous work has shown that protostellar outflows reduce the SFR of the parent cloud (Cunningham et al. 2011; Hansen et al. 2012; Federrath et al. 2014a; Murray et al. 2018), which we confirm.

Previous non-MHD simulations (e.g. Bate 2009; Krumholz et al. 2012) argued that radiation (specifically radiative heating by local protostars) and jets are the key ingredients to the IMF, where radiation heats the gas surrounding the star, preventing it from fragmenting and forming new stars, thus creating a mass reservoir that the protostar can almost fully accrete. This, however, can lead to an ‘overaccretion’ problem that is resolved by the addition of protostellar jets (Hansen et al. 2012; Krumholz et al. 2012). Later works also included MHD processes and produced IMFs similar to that observed (Cunningham et al. 2018; Li et al. 2018), but the combination of protostellar jets and MHD without radiation on cloud scales (>pc) was not investigated. Our simulations suggest that radiation may not be necessary to reproduce the observed IMF, as magnetic fields naturally provide support against fragmentation near newly formed stars. This is true regardless of the initial magnetization of the cloud as the turbulent dynamo drives the system towards a common B–ρ relation at high densities (see fig. 7 in Paper 0 and Appendix  A). However, several caveats are in order. (1) We focus primarily on statistics insensitive to the lowest mass stars (which may be most sensitive to radiation), as long as the IMF is shallower than Salpeter at low masses. We have not rigorously demonstrated that the low-mass IMF is numerically converged in our C_M_J simulations, even if M_mean is (see Paper I). (2) Our cooling/non-isothermal simulations include simple approximations to account for the transition between optically thin and thick cooling, rather than explicit radiation MHD; if these underestimate the cooling rates at high densities we might underestimate the need for radiative heating. (3) We enforce a constant dust and ‘floor’ temperature |$T_{\rm floor} = T_{\rm dust} = 10\,$| K. In future work, we will replace this with more realistic assumptions, but in Appendix  B we show the IMF shape at |${\lesssim}1\, \mathrm{ M}_{\odot }$| is quite sensitive to this value (this is essentially |$c_{s,\, \rm min}$|⁠, in our equation 10). So the IMF is sensitive to thermodynamics, and it remains to be seen whether more physical models for cooling and dust temperatures below |${\sim}100\,$| K can robustly reproduce the observed IMF without local radiative heating. (4) These simulations do not resolve protostellar discs (let alone disc fragmentation), whose stability may be critically impacted by radiative feedback. Also, our treatment neglects non-ideal MHD terms, so we see discs lose their angular momentum rapidly and simply accrete entirely on to the central sink owing to strong magnetic braking (Hennebelle & Fromang 2008; Wurster, Price & Bate 2016), artificially avoiding fragmentation (see Wurster, Bate & Price 2019 for a counterargument).

In addition to their effects upon sub-pc accretion flows and the IMF, we find that jets can have a significant global effect upon GMC kinematics and evolution in smaller clouds (Fig. 4). Specifically, protostellar jets alone appear sufficient to unbind initially bound clouds at least as massive as 2 × 10⁴ M_⊙, once a sufficiently high SFR and momentum injection rate are achieved. However, in Fig. 4 we see that for all but our least-massive clouds (⁠|$M=200\, \mathrm{M}_{\rm \odot }$|⁠), by the time jets begin to unbind the parent cloud (causing a sharp rise to α_turb ≫ 2) the integrated SFE has already reached |${\gg}10{{\ \rm per\ cent}}$| values, much larger than observed in MW clouds that motivate our ICs (Krumholz 2014). Thus, for clouds with masses |${\gt}1000\, \mathrm{M}_{\rm \odot }$| some other process (e.g. radiation from massive stars) must dominate cloud disruption. Even if some other feedback mechanism is ultimately responsible for GMC disruption (Section 5.4), the contribution of jets alone to the cloud kinematics can be significant. Therefore, it is likely that jet feedback has important non-linear interactions with other feedback mechanisms, potentially making it easier for stellar radiation to disrupt the cloud by increasing the initial turbulence or reducing the initial density at the time that massive stars break out from their envelopes. For this reason, previous simulations of feedback and cluster formation on GMC scales that neglected jets (including previous works by the present authors, e.g. Grudić et al. 2018, 2019b, 2020a) should be revisited.

5.3 Apparent sensitivity to initial conditions

We find that the mass scale set by protostellar jets exhibits significant sensitivity to variations in the ICs (see equation 10). Even if one could argue that parameters like the Γ momentum loading factor are set by atomic and nuclear physics in a way that they vary little between star-forming regions (despite differences in local metallicities), observed clouds, even in the Solar neighbourhood, have a wide range of masses (⁠|$10^4-10^6\, \mathrm{M}_{\rm \odot }$|⁠), densities (⁠|$10{-}1000\, \mathrm{cm}^{-3}$|⁠), virial parameters (0.1–10), and temperatures (⁠|$10{-}30\, \mathrm{K}$|⁠), see Kauffmann, Pillai & Goldsmith (2013), Heyer & Dame (2015), and Miville-Deschênes, Murray & Lee (2017). The properties of star-forming gas in more extreme environments (e.g. Galactic Centre, ULIRGS) can vary much more wildly (e.g. densities |${\gt}10^5\, \mathrm{cm}^{-3}$|⁠, surface densities at ∼1000 times higher values, molecular temperatures |${\sim}70{-}100\, \mathrm{K}$|⁠; see Dame, Hartmann & Thaddeus 2001; Gao & Solomon 2004; Longmore et al. 2012). Meanwhile the IMF is observed to be near-universal, with variations, even in extragalactic sources, within a factor of 3 or less in both the IMF peak and mass-to-light ratio. In equation (10), the strong dependence on temperature is perhaps most concerning here, as that alone would predict ∼4 variations among local clouds and factor ∼20–80 variations between the Solar neighbourhood and more extreme galactic environments.

5.4 Potential role of additional feedback physics

While we find that protostellar jets dramatically reduce the stellar mass scales to values similar to those observed, these models still have several shortcomings that only additional physics can address. The most significant issues with the current model are that (1) massive stars undergo runaway accretion, creating a top-heavy IMF (see Figs 5 and 6); (2) star formation continues potentially up to SFE of unity for massive GMCs; and (3) the stellar mass scale set by jets is sensitive to the temperature of the parent cloud, which may potentially violate the observed near-universality of the IMF.

As discussed in Section 5.2, one obvious step is the inclusion of radiative heating, which has been argued to be crucial in setting the mass scale of low-mass stars (Offner et al. 2009; Krumholz 2011; Krumholz, Klein & McKee 2011; Bate 2012; Myers et al. 2013; Guszejnov & Hopkins 2016; Guszejnov et al. 2016; Cunningham et al. 2018). Ionizing radiation of main-sequence stars as well as the stellar winds they emit could also potentially solve the runaway accretion of massive stars (Krumholz et al. 2012; Cunningham et al. 2018; Li et al. 2018). Furthermore, these feedback processes (along with supernovae) could allow massive stars to disrupt their natal cloud and quench star formation at the observed SFE levels (Grudić et al. 2019a; Krumholz, McKee & Bland-Hawthorn 2019; Li et al. 2019).

6 CONCLUSIONS

In this paper, we presented simulations from the STARFORGE project, which are high-resolution MHD simulations of the collapse of a GMC that also follow the evolution of individual stars. The runs include progressively more complex physics, starting from isothermal MHD, then adding cooling physics, then feedback in the form of protostellar jets. We found that the inclusion of jets dramatically alters the mass spectrum of sink particles (the simulation analogue of the observed stellar IMF). The resulting mass distribution is broadly similar to the observed IMF in both shape and scale, but additional physics is needed for a complete IMF theory.

We carried out a large suite of tests to determine the sensitivity of our results to variations in both ICs and input physics parameters. We found that the mean sink particle mass set by jets is insensitive to many parameters, but sensitive to the momentum loading of jets and the cold, dense gas and dust temperatures, and potentially the surface density. Based on observed variations in cloud properties these would lead to larger variations in the IMF than observed in the Solar neighbourhood and much larger variations in extreme environments (e.g. Galactic Centre, starburst galaxies).

While protostellar jets allowed our simulations to produce a realistic IMF at masses between |$0.1\,\mathrm{ and}\,10\, \mathrm{M}_{\rm \odot }$|⁠, massive stars (⁠|${\gt}10\, \mathrm{M}_{\rm \odot }$|⁠) undergo runaway accretion, leading to an increasingly top-heavy IMF with time, in increasing conflict with the observed IMF slope. Even though jets can ultimately quench star formation it requires >10 per cent of the cloud mass to be turned into stars for even low-mass GMCs (⁠|${\sim}10^4\, \mathrm{M}_{\rm \odot }$|⁠) so for massive GMCs (⁠|${\gt}10^5\, \mathrm{M}_{\rm \odot }$|⁠) star formation would likely continue until an order unity fraction of the gas turns into stars. Meanwhile, observed nearby clouds, whose properties motivate our ICs, achieve terminal SFE values of only a few per cent. We conclude that additional physics is required to stabilize the IMF and regulate star formation. Candidates for these processes will be explored in future work.

ACKNOWLEDGEMENTS

DG is supported by the Harlan J. Smith McDonald Observatory Postdoctoral Fellowship. MYG is supported by a CIERA Postdoctoral Fellowship. Support for PFH was provided by NSF Collaborative Research Grants 1715847 and 1911233, NSF CAREER grant 1455342, and NASA grants 80NSSC18K0562 and JPL 1589742. SSRO is supported by NSF Career Award AST-1650486 and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. C-AFG is supported by NSF through grant AST-1715216 and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. This work used computational resources provided by XSEDE allocation AST-190018, the Frontera allocation AST-20019, and additional resources provided by the University of Texas at Austin and the Texas Advanced Computing Center (TACC; http://www.tacc.utexas.edu).

DATA AVAILABILITY

The data supporting the plots within this article are available on reasonable request to the corresponding authors. A public version of the gizmo code is available at http://www.tapir.caltech.edu/∼phopkins/Site/GIZMO.html.

Footnotes

1

http://www.starforge.space

2

http://www.tapir.caltech.edu/∼phopkins/Site/GIZMO.html

3

https://github.com/mikegrudic/MakeCloud

4

Note that since Paper 0 improvements on the sink formation and accretion algorithms (see Paper I) have reduced the population of very low mass sinks in isothermal MHD runs compared to Paper 0, suggesting that a subpopulation of these was unphysical in origin (strengthening our conclusions about the necessity of additional physics to prevent an overly top-heavy IMF).

5

Note that our results in Fig. A2 show that the results are insensitive to whether we have an evolving or a constant R_*.

6

Note that the runs in Lee & Hennebelle (2019) did include magnetic fields and did not produce a top-heavy IMF. This is due to dense, highly turbulent ICs, which dramatically lowers the magnetic mass scale compared to what it would be in MW-like clouds (see section 4.3 in Paper 0), hence the opacity limit does dominate in this regime.

APPENDIX A: DEPENDENCE OF THE IMF ON INITIAL CONDITIONS

In this appendix, we present in detail the results of various test runs (see Table 2) with protostellar jets enabled. In Fig. A1, we find that both the mean and mass-weighted median sink masses are sensitive to the initial properties of the cloud (mass, virial parameter, surface density). The one exception is the initial level of magnetization, which appears to have negligible effects, similar to the isothermal case in Paper 0.

Figure A1.

Evolution of M₅₀ and the mean sink masses (left and centre columns, similar to Fig. 5) as well as the distribution of sink particle masses at 5 per cent SFE (right column) for M2e4_C_M_J (see Tables 1 and 2) with variations in the initial turbulent virial parameter α_turb, initial cloud mass M₀, cloud surface density Σ, and normalized magnetic mass-to-flux ratio μ. Note that we plot the IMF at a lower SFE for the surface density test, as the lowest surface density cloud becomes unbound before reaching 5 per cent SFE.

Fig. A2 shows the results of further tests where the parameters of the underlying physical models were varied. Fig. A2 shows that the mass spectrum is especially sensitive to the floor temperature of the simulation. We carried out an additional test where we varied the critical surface density Σ_crit where the cooling module transitions between optically thin and thick regimes. We found that varying Σ_crit (i.e. the opacity limit) by a factor of 10 in either direction has little effect on the sink mass spectrum. Transitioning to an isothermal EOS also has only minor effects that arise from the formation of very low mass sinks, which were previously suppressed by the EOS.

Figure A2.

Evolution of the mass-weighted median sink mass M₅₀ and the number-weighted mean sink mass M_mean (left and centre columns, similar to Fig. 5) as well as the distribution of sink particle masses at 5 per cent SFE (right column) for M2e4_C_M_J (see Tables 1 and 2) with variations in the floor temperature T_floor, the thermodynamics of the simulation (varying Σ_crit, the transition surface density between optically thin and thick cooling regimes), the parameters of the jet module (different f_K values as well as using a fixed R_* stellar radius, see Section 2.1.3), and the type of IC (Sphere versus Box, see Section 2.2.1).

As expected, changing the parameters of the jet module has significant effects, we find that the results are sensitive to the momentum loading of the jets, which is set by Γ, see Section 3.2.3 for details. Note that we also find that launching jets from a constant stellar radius, instead of the one set by the protostellar evolution model of Section 2.1.1, produces qualitatively similar results (see Fig. A2).

APPENDIX B: SCALING RELATIONS

In this appendix, we examine in detail how the mean sink mass depends on the ICs of the cloud (turbulent virial parameter α_turb, minimum sound speed c_s,min, normalized magnetic flux ratio μ, surface density Σ, and initial cloud mass M₀), by examining how the characteristic sink mass depends on each of them independently. We assume that the relation between the mean sink mass and the ICs is described by a multivariate power law. Using subsets of our runs from Table 2 where only one of these parameters is varied, we carry out least-squares fits to the individual exponents in turn, each at a fixed fiducial SFE value (4 per cent). To estimate the errors of the fitted exponents, we first estimate the errors in the mean sink mass using bootstrapping, which means resampling the sink mass distribution at fixed SFE and calculating the 95 per cent confidence interval of the mean mass over these new samples (see Fig. B1). We find the following fitting parameters and errors:

$$\begin{eqnarray*} M_{\rm mean}\propto \Gamma ^{-0.65\pm 0.15}\, c_{\rm s}^{2.5\pm 0.5}\, \Sigma ^{-0.3\pm 0.1}\, \alpha _{\mathrm{turb}}^{0.15\pm 0.19}\, M_0^{-0.12\pm 0.07}, \end{eqnarray*}$$

(B1)

which can be also expressed as

$$\begin{eqnarray*} M_{\rm mean}\propto \Gamma ^{-0.65\pm 0.05}\, c_{\rm s, min}^{2.5\pm 0.5}\, \rho ^{-0.2\pm 0.07}\, \alpha _{\mathrm{turb}}^{0.15\pm 0.19}\, M_0^{-0.22\pm 0.10}. \end{eqnarray*}$$

(B2)

See Fig. 10 for a visual representation of the goodness of the fit.

Figure B1.

Dependence of the mean sink mass at 4 per cent SFE on the initial turbulent virial parameter α_turb (top, left), minimum sound speed c_s,min (set by the floor temperature T_floor, top, middle), cloud mass M₀ (top, right), cloud surface density Σ (bottom, left), normalized magnetic mass-to-flux ratio μ (bottom, middle), and the Γ momentum loading of the jets (bottom, right). We chose 4 per cent as the reference SFE because some low mass and surface density runs are disrupted by jets before reaching 5 per cent SFE. Note that in the case of the jet momentum loading, we used runs with different initial turbulent realizations, which are shown slightly offset to make the plot easier to parse. The errors are estimated by bootstrapping: we resample the sink mass distribution at fixed total stellar mass and calculate the 95 per cent confidence interval of M₅₀ over these realizations, which we denote with errorbars. These scalings are discussed in the main text in Section 3.2.4.