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 ∼4M in the simulations, while ∼0.4M 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 2010and Offner et al. 2014,aswellasanalysisofDib2014).

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 (∼10foot_1 cm -3 , T ∼ 30 K) and highdensity regions where collapse occurs (∼10 5 cm -3 , T ∼ 10 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;R e e s1976; for modern interpretations, see Colman&T eyssier2020).

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;Bate2012;Myersetal.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će ta l .2018; Kim, Kim & Ostriker 2018;Lietal.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;W a n ge ta l .2010; Cunningham et al. 2011;Federrath et al. 2014a;Offner&Arce2014;Offner&Chaban2017;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. 1 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će ta l .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

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.

Similar to our previous studies of isothermal collapse with and without magnetic fields (Guszejnov et al. 2018b andGuszejnov 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 2 (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 ∇•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 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 ∼10 au.

Once sinks form, they follow the protostellar evolution model from Offner et al. (2009) ,w h i c hi sa l s ou s e di nt h eORION 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)andPaper I.

We compare simulations with two different thermodynamics modules. Our 'isothermal' simulations enforce an isothermal equation of state (EOS) with c s = 0.2kms -1 (effective gas temperature T ∼ 10 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 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 selfconsistently 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 dust = 10 K and a temperature floor of T floor = 10 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).

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

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 K f 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 K f 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)a n d many other works, which puts f w f K in the middle of the observed range. It is useful to introduce the momentum loading parameter

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( s e ePaper 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,w h e r e m is the mass resolution element of our simulation, for which our fiducial value is m = 10 -foot_3 M , 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.

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 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. We generate our ICs using MAKECLOUD. 3 Unless otherwise specified our runs utilize 'Sphere' ICs, meaning that we initialize a spherical cloud (T = 10 K, radius R cloud , and mass M 0 ) 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 -2 (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 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 -2 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 ff ≡ 3π 32Gρ 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 )thatgives the desired 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).

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

where c s is the gas sound speed and v 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),

where R cloud and M 0 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 massto-flux ratio), which for a uniform magnetic field can be expressed as

where E grav and E mag are the gravitational and magnetic energy (assuming a uniform initial field), respectively, while the normalization constant is c 1 ≈ 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:

where ρ 0 is the density of the gas in the cloud. The initial turbulence also has a characteristic length-scale: the sonic length, L sonic ,o n which the turbulent dispersion becomes supersonic. The correspondingmassscaleisthesonic mass:

where we used the supersonic linewidth-size relation (σ 2 (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 ρ/ρ 0 ∼ 1 + 1 3 M 2 (Padoan, Nordlund & Jones 1997), which scales as

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 0 , size R 0 , turbulent virial parameter α turb , normalized magnetic flux μ,a n d 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 × 10 5 M that lie along a mass-size relation similar to observed GMCs in the MW (e.g. Larson 1981,s p e c i fically assuming ≡ M 0 /π R 2 cloud = 63M pc -2 ). These clouds are marginally bound (α turb = 2) and start out at T = 10 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 a2× 10 4 M 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 0 ) 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.

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.

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 3 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 2 , similar to previous results with driven turbulent boxes (e.g. Federrath & Klessen 2012;Murrayetal.2015Murrayetal. , 2018)).

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 max ∝ √ t power law, which steepens to M max ∝ t 3 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)canlead to the disruption of the parent cloud and subsequently preventing the formation of new stars. In more massive clouds (>10 4 M , similar to MW GMCs), protostellar jets show no sign of arresting star formation before the SFE exceeds ∼10 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 Table 2. ICs of clouds used in our runs, with M 0 , R cloud , α turb , μ,andT 0 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 0 , see Section 2.1.2). We also report the initial 3D sonic Mach number M, thermal virial parameter α th , total virial parameter α,Alfvén Mach number M 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).

Derived parameters Highest resolution run 1 and2) 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 ∼ 0.1M )toblue(M ∼ 10 M ). This simulation resolves a dynamic range from ∼50 pc down to ∼30 au.

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.

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 50 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). Fig. 5 shows the evolution of the mean and median sink masses along with that of M 50 as a function of SFE. For runs without 1 and2). 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 ff = 3π 32Gρ 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 3 for Sphere ICs and SFE ∝ t 2 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 3 once massive stars have formed.

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 ( ≈ 60 M pc -2 ) similar to MW GMCs (M2e2-M2e5,seeT able2). 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. jets, we find that the mean sink mass and M 50 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 50 due to the runaway accretion of massive sinks, similar to the M 50 ∝ SFE 1/3 scaling found in the isothermal case in Paper 0. Switching to Box ICs has little effect on the evolution of M 50 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 50 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, andA2), while M 50 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.

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-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 50 , 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 50 increases with time, diverging from observations at higher SFE values. M sink /M 0 ) for the runs shown in Fig. 3.W e 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).

isothermal physics alone has little effect on the IMFfoot_5 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.

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 w f K /(1f 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).

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

which can also be expressed as M mean ∝ -0.65±0.05 c 2.5±0.5 s,min ρ -0.2±0.07 α 0.15±0.19 turb

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 ,andM 0 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 ∼ 60 M pc -2 ;seeLarson1981), we can simplify equation ( 9) as M mean ∝ -0.65±0.15 c 2.5±0.5 s,min α 0.15±0.19 turb M -0.12±0.07 0 .

(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.

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 <0.1 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 ∼10 pc in scale (bottom row of Fig. 9).

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

where Ṁacc is the mass accretion rate. Let us further assume that Ṁacc = const. and R * = const.foot_6 so that M * = Ṁacc t. This will simplify the above equation to

which we can integrate to get the total amount of momentum injected by jets over time t. Replacing t = M * / Ṁacc , we obtain

whereweha vealsousedthe 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 esc ∼ G(M * + M g )/R,whereR = (M g /(4π /3ρ)) 1/3 is the radius of the reservoir. This means Note that we exclude sink particles that do not reach 0.5M 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 <0.1 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.

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

Substituting in typical values for GMCs this becomes

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 50 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

turb , while M Jeans is independent, so we adopt M g = M Jeans in this model. Plugging it into equation ( 18), we get

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 ,andM turb BE ) 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 .

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 ṖJ > |F 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 Ṁacc follows a Bondi-like scaling, then (following similar logic as before) the core would be unbound when

.T h i s 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 momentumfeedback-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.

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 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.1kms -1 ) to orange (10 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.

(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 -2 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).

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 &Clark2012). 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;Rees1976), 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. 6 Including magnetic fields, however, in MW-like cloud conditions shifts the turnover mass of the IMF to much larger >20 M 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 (<0.1M ,seeFigs5 and 6).

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;W a n g et al. 2010;O f f n e r&A r c e2014; Offner & Chaban 2017;Murray et al. 2018). We do find that jets disrupt the local accretion flow, which greatly changes gas dynamics on <0.1pcscales,butthishas 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;Hansenetal.2012;Federrath et al. 2014a;Murrayetal.2018) 2) with the initial Jeans mass M Jeans (equation 6), sonic mass M sonic (equation 7), turbulent Bonnor-Ebert mass M turb BE (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 50 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. 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;L ie ta l .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 floor = T dust = 10 K. In future work, we will replace this with more realistic assumptions, but in Appendix B we show the IMF shape at 1M is quite sensitive to this value (this is essentially c s, 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 ∼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&Bate2016), 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 4 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 M ), by the time jets begin to unbind the parent cloud (causing a sharp rise to α turb 2) the integrated SFE has already reached 10 per cent values, much larger than observed in MW clouds that motivate our ICs (Krumholz 2014). Thus, for clouds with masses >1000 M 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ćetal. 2018Grudićetal. , 2019bGrudićetal. , 2020a) should be revisited.

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 M ), densities (10-1000 cm -3 ), virial parameters (0.1-10), and temperatures (10-30 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 >10 5 cm -3 , surface densities at ∼1000 times higher values, molecular temperatures ∼70-100 K; see Dame, Hartmann & Thaddeus 2001;Gao & Solomon 2004;Longmore et al. 2012). Meanwhile the IMF is observed to be nearuniversal, 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.

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 and6); (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;B a t e2012;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;Lietal.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;L ie ta l . 2019).

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 and 10 M ,m a s s i v es t a r s (>10 M ) 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 (∼10 4 M ) so for massive GMCs (>10 5 M ) 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.

. Evolution of M 50 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 and2) with variations in the initial turbulent virial parameter α turb , initial cloud mass M 0 , 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.

. Evolution of the mass-weighted median sink mass M 50 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 and2) 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).