Effects of the environment and feedback physics on the initial mass function of stars in the STARFORGE simulations

ABSTRACT

One of the key mysteries of star formation is the origin of the stellar initial mass function (IMF). The IMF is observed to be nearly universal in the Milky Way and its satellites, and significant variations are only inferred in extreme environments, such as the cores of massive elliptical galaxies and the Central Molecular Zone. In this work, we present simulations from the STARFORGE project that are the first cloud-scale radiation-magnetohydrodynamic simulations that follow individual stars and include all relevant physical processes. The simulations include detailed gas thermodynamics, as well as stellar feedback in the form of protostellar jets, stellar radiation, winds, and supernovae. In this work, we focus on how stellar radiation, winds, and supernovae impact star-forming clouds. Radiative feedback plays a major role in quenching star formation and disrupting the cloud; however, the IMF peak is predominantly set by protostellar jet physics. We find that the effect of stellar winds is minor, and supernovae ‘occur too late’ to affect the IMF or quench star formation. We also investigate the effects of initial conditions on the IMF. We find that the IMF is insensitive to the initial turbulence, cloud mass, and cloud surface density, even though these parameters significantly shape the star formation history of the cloud, including the final star formation efficiency. Meanwhile, the characteristic stellar mass depends weakly on metallicity and the interstellar radiation field, which essentially set the average gas temperature. Finally, while turbulent driving and the level of magnetization strongly influence the star formation history, they only influence the high-mass slope of the IMF.

MHD, turbulence, stars: formation, stars: luminosity function, mass function, stars: winds, outflows

1 INTRODUCTION

Although star formation is one of the most fundamental processes of astrophysics, there is no widely accepted theory of star formation, despite decades of intensive work from both observers and theorists (McKee & Ostriker 2007; Krumholz 2014). The primary reason for this is the large set of interconnected, complex physical processes, including gravity, turbulence, magnetic fields, chemistry, and radiation (Girichidis et al. 2020). Furthermore, these processes interact in a non-linear way that also create interactions between vastly different scales (e.g. feedback from massive stars affecting their progenitor cloud). Thus, in order to understand star formation, it is vital to investigate the role each physical process plays and how it modifies the outcome. A key question is how these processes affect the at-formation stellar mass spectrum, i.e. the initial mass function (IMF).

Since the large set of complex physical processes involved prevents a direct treatment, star formation models (both analytical and numerical) can only include a limited subset of physical processes. Analytic models and early simulations modelled the dense, star-forming clouds of the Milky Way (MW) as isothermal, turbulent objects collapsing under self-gravity (e.g. Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012). Recent numerical works have shown that this set of physics is inadequate to produce a converged mass spectrum for collapsed fragments (i.e. stars; see Martel, Evans & Shapiro 2006; Kratter et al. 2010; Guszejnov, Krumholz & Hopkins 2016; Federrath, Krumholz & Hopkins 2017; Guszejnov et al. 2018b; Lee & Hennebelle 2018), so additional physics must play a role.

Star-forming clouds are observed to have significant support from magnetic fields (Crutcher 2012). Both theoretical and numerical works have found that the addition of magnetic fields 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). Several of these studies claimed to reproduce the observed IMF. However, a larger parameter study of high-resolution simulations (Guszejnov et al. 2020b, henceforth referred to as Paper I) showed that the mass scale imposed by magnetic fields is both too large and too sensitive to initial conditions (ICs), in a way that would violate the observed near-universality of the IMF (Bastian, Covey & Meyer 2010).

While clouds are close to isothermal for a wide range of densities, there can be significant deviations even at low densities around |$\sim 10^2\, \mathrm{cm}^{-3}$| (see Glover & Clark 2012). In high-density regions (⁠|$\sim 10^5\, \mathrm{cm}^{-3}$|⁠), the isothermal assumption breaks down as the cloud becomes opaque to its own cooling radiation, leading to increased temperatures and suppressed fragmentation (Low & Lynden-Bell 1976; Rees 1976; Lee & Hennebelle 2018; Colman & Teyssier 2020). In our recent work (Guszejnov et al. 2021, henceforth referred to as Paper II), we showed that non-isothermal effects by themselves are insufficient to lower the characteristic stellar mass to the observed value, and these effects are most significant in providing a minimum mass scale for star formation (see Low & Lynden-Bell 1976).

Recent works show that the energy and momentum injected by newly formed stars to their environment (i.e. stellar feedback) can dramatically affect the star formation process (Offner et al. 2009; Krumholz 2011; Bate 2012; Myers et al. 2013; Guszejnov & Hopkins 2016; Guszejnov et al. 2016; Rosen et al. 2020). The first of these processes to act during star formation is the launching of high-velocity bipolar outflows by accreting protostars (see reviews of Frank et al. 2014; Bally 2016) that are likely driven by collimated bipolar jets launched along the rotational axis of the accreting protostar (Rosen & Krumholz 2020), which are produced by the magnetic interaction between the protostar and its accretion disc (Shu et al. 1988; Pelletier & Pudritz 1992) and radiation pressure in the case of massive protostars (Kuiper et al. 2010; Vaidya et al. 2011; Rosen et al. 2016). Previous work has shown that these jets both reduce accretion rates and drive turbulence on small scales (Matzner 2007; Nakamura & Li 2007; Wang et al. 2010; Cunningham et al. 2011; Offner & Arce 2014; Federrath et al. 2014a; Offner & Chaban 2017; Murray, Goyal & Chang 2018; Rohde et al. 2021; Appel et al. 2022). Recently in Paper II, we showed that jets can dramatically lower the characteristic stellar mass as they disrupt the accretion flow around stars, allowing the nearby gas to fragment (similar to the results of Cunningham et al. 2018; Li, Klein & McKee 2018; Mathew & Federrath 2021). Overall, jets reduced stellar masses by an order of magnitude compared to magnetohydrodynamic (MHD) simulations without stellar feedback, thereby bringing the simulated stellar mass spectrum in agreement with the observed IMF. However, Paper II also found that jets by themselves are insufficient to regulate star formation and that massive stars undergo runaway accretion without additional feedback processes.

Radiation from accreting protostars is thought to be a crucial ingredient of the star formation process for both low- and high-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; Rosen et al. 2016, 2020; Cunningham et al. 2018) as they can heat their surroundings, preventing fragmentation. Once stars reach the main sequence, they start emitting ionizing radiation as well as isotropic line-driven stellar winds, both of which can dramatically affect their surroundings, potentially halting stellar accretion (Krumholz, Klein & McKee 2012; Rosen et al. 2016, 2020; Cunningham et al. 2018; Li et al. 2018; Rosen 2022). Massive stars in particular provide feedback powerful enough to affect the IMF in their entire natal cloud (Gavagnin et al. 2017), as well as completely quench star formation (see Krumholz, McKee & Bland-Hawthorn 2019 for review and fig. 1 of Grudić et al. 2020 for a literature compilation of theoretical predictions), helping to limit the star formation efficiency (SFE) of clouds to a few per cent (Grudić et al. 2018, 2019, 2021b; Kim, Kim & Ostriker 2018; Li et al. 2019).

Figure 1.

Row 1: surface density maps for M2e4_C_M_J_RT with M₀/Δm = 2 × 10⁷ initial gas cells (see Table 1) at different times for the Sphere IC. 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{20\, pc}$| to |$\sim \!\mathrm{30\, au}$| and is run until stellar feedback quenches star formation and disrupts the cloud (see the right-hand column). Row 2: 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). These kinematic maps can highlight feedback processes that would be invisible in surface density maps (i.e. protostellar jets).

The lifetime of massive stars is several Myr, comparable to the lifetime of star-forming clouds. The resulting supernovae (SNe) dominate the momentum input by stellar feedback in the interstellar medium (ISM) and are critical to regulate star formation on galactic scales (Somerville & Davé 2015; Naab & Ostriker 2017; Vogelsberger et al. 2020) and could be a vital ingredient to the cloud-scale star formation process as well. Their effect, however, is reduced by the fact that they act ‘late’ in the star formation process, so it is unclear how much they either affect the IMF or regulate star formation. Simulations of star cluster formation have found that they have negligible impact upon SFE and bound cluster masses compared to early feedback (i.e. radiation), even in massive Giant Molecular Clouds (GMCs) that survive long enough to host an SN before disruption (Grudić et al. 2021b). Nevertheless, SNe must at least play an indirect role in cloud-scale star formation as they are thought to be one of the main drivers of galactic turbulence and thus set the properties of GMCs (e.g. Hopkins, Quataert & Murray 2011, 2012; Ostriker & Shetty 2011; Faucher-Giguère, Quataert & Hopkins 2013; Padoan et al. 2017; Seifried et al. 2018; Gurvich et al. 2020; Guszejnov et al. 2020a).

Simulations that take into account the above processes are necessary to understand the effects of each physical process, but so far such studies have generally been limited to simple physics or to a very narrow range of cloud ICs. In this paper, we introduce a suite of results from the STAR FORmation in Gaseous Environments (STARFORGE) project,¹ which include all of the above physical processes. These radiation-magnetohydrodynamic (RMHD) 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 molecular clouds while following the formation of individual stars with stellar masses as low as ∼0.1 M_⊙ (see methods paper of Grudić et al. 2021a, henceforth referred to as Methods Paper). In this study, we perform and analyse a set of simulations with different ICs and levels of physics to identify the impact of outflows, stellar radiation, winds, and SNe on the IMF and the star formation history of clouds.

First, we provide a brief overview of the STARFORGE code base and the various parameters and metrics we use in Section 2 (for a more detailed discussion, see the Methods Paper). We present our results in Section 3 with a focus on how the star formation history and the characteristic masses of sink particles (stars) change with the inclusion of additional physics. In Section 4, we explore variations in the ICs (e.g. cloud surface density and level of turbulence) and physical parameters (e.g. turbulent driving). The implications of these results as well as the potential role of further, not-yet included physics are discussed in Section 5. We summarize our conclusions in Section 6.

Note that for brevity we are only showing figures essential for the main points of this paper. For additional figures, we refer the reader to the online supplementary materials of this paper, which can also be found in a GitHub repository.²

2 NUMERICAL METHODS

2.1 The STARFORGE simulations

For this work, we utilize simulations from the STARFORGE project, which are run with the gizmo simulation code.³ A full description and presentation of the STARFORGE methods including a variety of tests and algorithm details are given in the Methods Paper; therefore, we only briefly summarize the key points here. Readers familiar with the STARFORGE simulations should skip ahead to Section 3.

2.1.1 Physics

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.

Sink particles represent individual stars. Once they form, they follow the protostellar evolution model from Offner et al. (2009).

‘Non-isothermal’ or ‘cooling’ STARFORGE runs utilize the radiative cooling and thermochemistry module from Hopkins et al. (2022) that contains detailed metallicity-dependent cooling and heating physics from |$T=10\!-\!10^{10}\,$|K, including recombination, thermal bremsstrahlung, metal lines (following Wiersma, Schaye & Smith 2009), molecular lines, fine structure (following Glover & Abel 2008), and dust collisional processes. The cooling module self-consistently solves for the internal energy and ionization state of the gas (see Hopkins et al. 2022 and appendix B of Hopkins et al. 2018b). STARFORGE simulations use two different treatments of radiation transport: In this work, we present radiation-hydrodynamical (RHD) simulations that co-evolve the gas, dust, and radiation temperature self-consistently (unlike in Paper II), including the stellar luminosity in various bands accounting for photon transport, absorption, and emission using dust opacity. In addition to local sources (i.e. stars), we include an external heating source that represents the interstellar radiation field (ISRF).

As shown in Paper II, protostellar jets represent a crucial feedback mechanism. We model their effects by having sink particles launch a fixed fraction (f_w = 0.3) of the accreted material along their rotational axis with a fixed fraction (f_K = 0.3) of the Keplerian velocity at the protostellar radius.

In addition to their radiative feedback, massive main-sequence stars inject a significant amount of momentum and energy into their surroundings through stellar winds. We calculate the mass-loss rates and wind velocities based on Smith (2014) and Lamers, Snow & Lindholm (1995), respectively. Winds are implemented either through local mass, momentum, and energy injection or direct gas cell spawning, while SNe are spawned at the end of the lifetime of all >8 M_⊙ stars.

2.1.2 Cloud parameters

To describe our ICs, we introduce several parameters, using the same definitions as in Papers I and II. First, we introduce the 3D sonic Mach number

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

(1)

where c_s is the gas sound speed and |$\boldsymbol {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 ||\boldsymbol {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}$$

(2)

where R_cloud and M₀ are the cloud (spherical-equivalent) radius and total mass. 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}$$

(3)

where the normalization constant c₁ ≈ 0.4 (Mouschovias & Spitzer 1976).

For a detailed definition of these quantities and the others listed in Table 1, see section 2 in Paper I.

Table 1.

Simulations used in this paper described with STARFORGE label conventions. Top: physics modules included, see Section 2.1.1 and Methods Paper for details on the individual physics modules. Bottom: ICs of clouds used in our runs, with M₀, R_cloud, α_turb, and μ being the initial cloud mass, size, virial parameter, mass-to-magnetic flux ratio, and temperature, respectively (note that the initial gas–dust temperature is set by ISRF). We also report the initial 3D turbulent velocity dispersion σ, thermal virial parameter α_th assuming |$T=10\, \mathrm{K}$|⁠, 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 I 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.

Physics label	Thermodynamics	MHD	Protostellar jets	Stellar radiation	Stellar winds and SNe
‘Physics ladder’ of star formation
I_M	Isothermal (I)	Ideal (M)	Not included
C_M	Non-isothermal, RHD (RHD)	Ideal (M)	Not included
C_M_J	Non-isothermal, RHD (RHD)	Ideal (M)	Included (J)	Not included
C_M_J_R	Non-isothermal, RHD (RHD)	Ideal (M)	Included (J)	Included (R)	Not included
C_M_J_R_W	Non-isothermal, RHD (RHD)	Ideal (M)	Included (J)	Included (R)	Included (W)
Physics variation tests
ISRFx10	Includes all like C_M_J_R_W, but the background ISRF is 10 times the solar circle value
ISRFx100	Includes all like C_M_J_R_W, but the background ISRF is 100 times the solar circle value
Z01	Includes all like C_M_J_R_W, but metallicity is 10 per cent of the solar value
Z001	Includes all like C_M_J_R_W, but metallicity is 1 per cent of the solar value

	Input parameters						Derived parameters								Resolution
Cloud label	M0 (M⊙)	Rcloud (pc)	Lbox (pc)	αturb	μ	σ (km s−1)	α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	1.0	0.02	2.04	10	7.8	0.02	6 × 10−2	7 × 10−3	0.1	2 × 105	36
M2e3	2 × 103	3		2	4.2	1.9	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	3.2	0.008	2.03	10	0.78	0.02	3 × 10−3	7 × 10−5	0.1	2 × 107	36
Parameter variation tests
M2e4_R3	2 × 104	3		2	4.2	5.8	0.008	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	1.9	0.02	2.04	10	2.3	0.02	1 × 10−2	6 × 10−4	0.1	2 × 107	36
M2e4_a1	2 × 104	10		1	4.2	2.3	0.008	1.03	10	0.78	0.02	3 × 10−3	4 × 10−5	0.1	2 × 107	36
M2e4_a4	2 × 104	10		4	4.2	4.5	0.008	4.03	10	0.78	0.02	3 × 10−3	1 × 10−4	0.1	2 × 107	36
M2e4_mu1.3	2 × 104	10		2	1.3	3.2	0.008	2.21	3.1	0.078	0.2	3 × 10−3	7 × 10−5	0.4	2 × 107	36
M2e4_mu0.4	2 × 104	10		2	0.42	3.2	0.008	4.01	3.1	0.0078	2	3 × 10−3	7 × 10−5	4	2 × 107	36

2.1.3 Initial conditions

We generate our ICs using MakeCloud,⁴ identical to Paper II. Unless otherwise specified, our runs utilize ‘Sphere’ ICs, meaning that we initialize a spherical cloud (radius R_cloud and mass M₀) with uniform density, surrounded by diffuse gas with a density contrast of 1000 (see Fig. 1). The cloud is placed at the centre of a periodic 10 R_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. Although the gas is initialized at |$T=10\, \mathrm{K}$|⁠, but the gas–dust mixture quickly reaches equilibrium with the ISRF, for which we assume solar neighbourhood conditions. Also, the gas is initially fully atomic and reaches an equilibrium molecular fraction by the time star formation begins (note that in Box runs the equilibrium is reached during the ‘stirring’ phase).

We also run simulations using ‘Box’ ICs, similar to the driven boxes used in e.g. Federrath et al. (2014a) and Cunningham et al. (2018), see Fig. 2. These are initialized as a constant density, zero velocity periodic cubic box with the same temperature prescription as the ‘Sphere’ ICs. 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 times the box size, with an appropriate decay time for driving mode correlations (t_decay ∼ t_cross). This stirring is initially performed without gravity for five global freefall times (t_ff, see equation 5), 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 that gives the desired Mach number |$\mathcal {M}$| and virial parameter α_turb. We use purely solenoidal driving, which remains active throughout the simulation after gravity is switched on, unless specified otherwise. We take the box side length L_box to give a box of equal volume to the associated Sphere cloud model. 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 a relative magnetic energy level of α_B ∼ 0.1 (see Paper I), so for Box runs, the ‘pre-stirring’ magnetic field strength (defined by the normalized flux μ) 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).

Figure 2.

Same as Fig. 1 but for Box IC.

Table 1 shows the target parameters for the runs we present in this paper. The input parameters are the cloud mass M₀, radius R₀, turbulent virial parameter α_turb, and normalized magnetic mass-to-flux ratio μ. Similar to Paper II, we set up our clouds to 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 either |$T=10\, \mathrm{K}$| or in equilibrium with the ISRF (RHD runs only). For the initial magnetization, we assume α_B = −2E_mag/E_grav = 0.02, which translates to μ = 0.4. Note that due to the much higher computational cost of RHD runs, only clouds up to |$2\times 10^4\, \mathrm{M}_{\rm \odot }$| are simulated with explicitly evolved radiation. Also, in runs where the star formation is not quenched (i.e. those without stellar radiation), we restrict ourselves to times when the SFE (= M_⋆/M₀) is below 10 per cent since most MW GMCs achieve an 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 an SFE of <1 per cent; see Federrath & Klessen 2013). Note that the STARFORGE simulations have an effective mass resolution of |$\Delta m=10^{-3}\, \mathrm{M}_{\rm \odot }$|⁠, making the mass function incomplete for |$M\lesssim M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$|⁠, which are thus omitted from our analysis. See Appendix  A for a detailed explanation for our choice of M_complete.

2.1.4 Star formation metrics

To describe the different aspects of the star formation process in our simulations, we use a series of variables, including the star formation efficiency SFE

$$\begin{eqnarray} \mathrm{SFE}\equiv M_\mathrm{sink}/M_0 = M_{\star }/M_0, \end{eqnarray}$$

(4)

where M₀ is the initial mass of the cloud, while M_sink and M_⋆ are the total mass in sink particles and stars, respectively. Note that this metric is normalized by the total cloud mass, so the SFE of individual self-gravitating subregions can be higher. A characteristic time-scale of the problem is the initial freefall time t_ff of the cloud:

$$\begin{eqnarray} t_\mathrm{ff}\equiv \sqrt{\frac{3 \pi }{32 G \rho _0}}, \end{eqnarray}$$

(5)

where ρ₀ is the initial cloud density. Combining these two leads to the SFE per freefall time ϵ_ff, a common metric in the literature (Krumholz & McKee 2005):

$$\begin{eqnarray} \epsilon _\mathrm{ff}(t) \equiv \dot{M}_\mathrm{sink}(t) t_\mathrm{ff}/ M_{\rm gas}, \end{eqnarray}$$

(6)

where M_gas is the available initial gas mass. To better describe the time dependence of these values, we also introduce |$\tilde{t}$|⁠, the number of initial freefall times elapsed since the cloud started star formation:

$$\begin{eqnarray} \tilde{t}(t) \equiv \left(t-t_\mathrm{SF\, starts}\right)/t_\mathrm{ff}, \end{eqnarray}$$

(7)

where |$t_\mathrm{SF\, starts}$| is the time the first star forms in the simulations.

A key prediction of any star formation model is the initial mass spectrum of stars (i.e. the IMF), which we assume to be identical to the mass spectrum of sink particles at the end of the simulation (for caveats, see Section 5.3). To illustrate the evolution of the mass spectrum and to make comparisons easier, it is useful to derive a characteristic mass scale of the stellar population. Note that all such scales are only calculated for the stars above our |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$| completeness limit. In this work and its appendices, we use several different summary statistics of the mass function.

• M_med: The number-weighted median stellar mass. Half of the stars will be more massive than this value. In strongly peaked distributions, it provides a good estimate for the peak of the IMF (i.e. the peak of dN/dlog M). Note that by including only stars above our M_complete completeness limit, the median stellar mass is not sensitive to mass resolution (see Appendix  A).

• M_mean: The number-weighted mean stellar mass. Similar to M_med, it provides an estimate for the IMF peak.

• M₅₀: The mass-weighted median stellar mass. Half of the stellar mass is in stars more massive than this value. This measures ‘where the mass is’ in the IMF. This metric is insensitive to low-mass objects and probes the high-mass tail of the IMF. For example, in Paper II, we found this metric to be strongly affected by the runaway accretion of massive stars that happens if stellar feedback is insufficient.

• γ_{1, 10}: The effective slope of the IMF between |$1\hbox{ and }10\, \mathrm{M}_{\rm \odot }$| that is derived from the mean stellar mass within the same range. Assuming the IMF in that range to be a pure power law of the form |$\mathrm{d}N/\mathrm{d}M \propto M^{-\gamma _{1,10}}$|⁠, there is a one-to-one mapping between γ_{1, 10} and the mean mass within that range. This metric probes the most well-constrained part of the IMF: Observations find this slope to be near-universal in the Local Group with a value of −2.35 ± 0.25 (Salpeter 1955; Kroupa 2002; Offner et al. 2014).

• M_max: The mass of the most massive star. Since stellar feedback is strongly non-linear with mass, M_max is expected to have a major impact on the evolution of the cloud. Note, however, that M_max is partially stochastic in our simulations, both due to truly random elements of the simulation and the chaotic nature of the problem (Geen et al. 2018; Grudić et al. 2021a).

• L/M: The light-to-mass ratio that is the ratio of the total bolometric luminosity from stars divided by the total stellar mass (in units of L_⊙ M_⊙⁻¹). Since stellar luminosity is |$L_{*}\mathrel {\propto \atop \sim } M_*^{7/2}$|⁠, this quantity is set by the most massive stars in the simulation and has the same stochasticity issue as M_max.

Finally, we compare the stellar mass distributions predicted by our simulations with the observed IMF. In simulations where star formation is quenched and the cloud is disrupted, we use the post-disruption sink mass distribution as the IMF.⁵ In runs that neglect the feedback processes required to quench star formation, we take the sink mass distribution at |$\mathrm{SFE}\sim 10{{\ \rm per\ cent}}$| (similar to the values runs with quenched star formation achieve) to be the IMF. We compare these simulated IMFs with the fitting function of Kroupa (2002) that is widely used in the literature,⁶ henceforth referred to as K02. For the mass scales above, we calculate the 95 per cent confidence interval for the individual values at different total stellar masses by varying the IMF with the uncertainty values reported in K02 and calculating the various mass scales for each realization. Similarly, for the IMF, we take the 95 per cent confidence intervals in each individual mass bin of the distribution.

3 ROLE OF DIFFERENT PHYSICS IN SETTING THE IMF AND STAR FORMATION HISTORY

We carried out a suite of simulations using combinations of the ICs and physics from Table 1. In this section, we concentrate on a set of runs with increasingly complex physics, forming a set of runs we denote as the ‘physics ladder’ of star formation. The ‘ladder’ starts from a base model of isothermal MHD and gravity (I_M, explored in Paper I), then adds non-isothermal thermodynamics and protostellar jets (C_M and C_M_J, respectively). Note that these have been explored in Paper II with a different thermodynamics treatment that does not directly evolve the radiation fields (ApproxRad, see Methods Paper). In our analysis, we found no qualitative difference in the star formation history or the IMF between those and the new RHD results. The next step is including stellar radiation (C_M_J_R) and, finally, stellar winds and SNe (C_M_J_R_W), which we recently explored in a single example simulation in Grudić et al. (2022). We also include a run with all physical processes except protostellar jets to showcase their importance (C_M_R_W). Due to the computational costs of these simulations, the M2e4 clouds were the largest ones the full ‘physics ladder’ was run on, and we use these runs in the subsequent analysis (see Table 1 for details on the IC). In the following subsections, we investigate the effects of each ‘rung’ of the physics ladder on the star formation history of the cloud (Section 3.1) and the stellar mass spectrum (Section 3.2).

3.1 Star formation history

Fig. 3 shows the star formation histories of simulations with identical Sphere ICs (M2e4) for different rungs of the ‘physics ladder’ suite. In all runs, the SFE evolves as |$\mathrm{SFE}\propto \tilde{t}^3$|⁠, similar to the findings for runs without feedback or with protostellar jets only (Papers I and II, respectively). Note that these results are sensitive to the ICs, see Section 4.8 for details. In runs without feedback (I_M,C_M), star formation continues unregulated. With the addition of protostellar jets (C_M_J), star formation is partially suppressed at later times, but continues without cloud disruption. The addition of stellar radiation (C_M_J_R) allows massive stars to greatly influence the cloud evolution, first by blowing away nearby gas (and thus stopping accretion) and by eventually quenching star formation in the cloud. Winds do not qualitatively alter this picture. Meanwhile, SNe go off at the end of the full physics run (C_M_J_R_W), but it occurs after the cloud is disrupted, too late to meaningfully alter the star formation history of the cloud.

Figure 3.

The evolution of the SFE (left-hand panel), SFE per freefall-time ϵ_ff (centre panel), and the number of sink particles N_sink (right-hand panel) as a function of time for a set of runs with increasingly complex physics in M2e4 clouds, corresponding to labels I_M, C_M, C_M_J, C_M_J_R, C_M_J_R_W, as well as C_M_R_W (see Table 1). The times of the formation of the first O star and the first SN are marked for the full physics (C_M_J_R_W) case. Note that t = 0 denotes the start of the simulation, while |$\tilde{t}=0$| is set to the start of star formation (see equation 7). Without stellar feedback, the cloud undergoes runaway star formation. The inclusion of jets suppresses star formation at later times, preventing new stars from forming but still allowing massive stars to accrete. The addition of stellar radiation is required to completely quench star formation. Note that in runs without radiative feedback, star formation does not quench before reaching non-physical values (>10 per cent), where the runs are terminated.

3.1.1 Cloud disruption

One key question of star formation is the regulation of star formation within clouds to achieve an SFE of a few per cent (Krumholz 2014). Figs 3 and 4 show that the simulations from the ‘physics ladder’ suite attain three qualitatively different endings. Runs without any stellar feedback (I_M, C_M) continue to form stars at an accelerated rate as the parent cloud undergoes global gravitational collapse. With the inclusion of protostellar jets (C_M_J), global collapse is slowed at high SFE values due to the excess momentum provided by jets to the gas, but star formation continues at a reduced rate. With the inclusion of radiative feedback, massive stars are able to disrupt the cloud and blow away the remaining gas. Note that this blown-out gas is still able to form stars, but at a much reduced rate (see late times in Fig. 3). This late-stage star formation happens late enough in the simulation that SNe from previously formed massive stars can affect it. SNe completely shut down any remaining star formation and blow away all remaining gas.

Figure 4.

Surface density maps at 5 Myr and at the end of simulations for runs with different levels of stellar feedback: no feedback (C_M, left-hand panels), protostellar jets only (C_M_J, middle panels), and with jets, radiation, winds, and SNe enabled (C_M_J_R_W, right-hand panels). The inclusion of additional feedback physics dramatically changes both the sink mass distribution and cloud morphology. Note that the end of the simulations is at different times, as runs that do not experience cloud disruption (C_M and C_M_J) are stopped at 10 per cent SFE.

Fig. 5 shows the evolution of the cloud energetics within the ‘physics ladder’ suite. Since we use Sphere ICs (see Section 2.1.3), turbulence initially decays in the simulations, until gravitational collapse provides enough kinetic energy to saturate roughly to the virialized values (α_turb ∼ 1). During this time, the initial magnetic field is amplified by the turbulent dynamo, increasing the relative magnetic energy to gravitational from |$\sim 1{{\ \rm per\ cent}}$| to |$\sim 10{{\ \rm per\ cent}}$| (see Paper I and references therein for a discussion). The inclusion of non-isothermal physics does not significantly affect the overall energy evolution of the cloud. Note that the virial parameter of the gas in the late stages of this run does go above the α = 2 boundedness limit, but the cloud is still bound due to the sink particles (stars). Stellar feedback in the form of protostellar jets dramatically alters the cloud’s energetics, as jets entrain nearby gas, creating outflows. The overall effect stellar feedback has on the cloud energetics is a dramatic increase in kinetic and thermal energy, raising the virial parameter to well above the boundedness limit. However, this does not mean that the cloud is fully disrupted and SF is quenched, as a large amount of heated gas remains (see Fig. 4) that can still be accreted by the central cluster, fuelling star formation and the growth of existing massive stars. The addition of radiative heating produces initially similar trends as previous runs, but the formation of the first massive, main-sequence O star dramatically affects the evolution of the cloud. Massive stars emit an enormous amount of radiation, outproducing the luminosity of all other stars. Since massive stars form in the dense, central regions of our simulated clouds, a significant portion of their radiation cannot escape, leading to a marked increase in thermal and kinetic energy. This surge in thermal and kinetic energy completely unbinds the cloud relatively quickly, leading to α > 10 values. Without jets, stellar masses are significantly higher, leading to more massive stars and a significantly earlier disruption of the cloud. The addition of stellar winds makes it easier for massive stars to unbind the cloud, but they have no qualitative effect on the simulations in our suite. Finally, SNe occur only after the cloud has been unbound by radiative feedback. It should be noted that for more massive clouds, which have longer dynamical times, SNe could occur early enough to play a role in cloud disruption.

Figure 5.

Evolution of cloud energetics as a function of time in the ‘physics ladder’ suite. Vertical lines mark the times the cloud reaches the boundedness limit (α = 2) and when the cloud is considered completely disrupted (α = 10). On the energy evolution trends, we mark the appearance of the first main-sequence O star (purple circle), the point where a significant portion of the cloud has been turned into stars (⁠|$\mathrm{SFE}=5{{\ \rm per\ cent}}, 10{{\ \rm per\ cent}}$|⁠, marked with black X and +, respectively) and all SNe (red stars). Overall radiative feedback plays a key role in the disruption of the cloud, with the formation of the first O star marking the turning point in the cloud’s energy evolution regardless of the presence of jet feedback. Once these massive stars formed, the kinetic and thermal energy dramatically increases, which leads to the unbinding of the cloud (apparent in the decreasing gravitational binding energy).

Note that the disruption time of the cloud is highly sensitive to the time and mass of the first O-type star, which itself is subject to the initial turbulent realization as well as stochastic effects (see Section 4.1). Thus, the effects of different physics or parameters on the cloud disruption time can only be studied in a statistical sense for which we currently have too few simulations.

3.2 Sink particle IMF

In this subsection, we analyse the effects the different physical processes in the ‘physics ladder’ have on the stellar mass spectrum. In the simulation, sink particles represent stars (or star systems with unresolved stellar companions when the at-formation separation is below the |$\mathrm{d}x\sim 30\, \mathrm{au}$| Jeans resolution of the simulation). Note that an analysis of a subset of these simulations has already been presented in Papers I and II, so here we concentrate on the new results from the simulations that include radiation, winds, and SN feedback processes.

3.2.1 Evolution of stellar mass scales

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 limit. While in most cases, these objects represent a vanishingly small fraction of the total sink mass (see Paper I for an example and Guszejnov et al. 2018b for a counterexample), their large number skews the mean and median sink masses. We mitigate these effects by taking the mean and median only for stars more massive than the completeness limit of the simulation. The mass-weighted median mass of sinks M₅₀ does not suffer from this effect (see Krumholz et al. 2012 and Paper I), but this choice makes the characteristic mass scale overly sensitive to the most massive sinks, leading to significant variations due to low number statistics (see Section 4.1 for details). Thus, to give a more holistic picture of the evolution of stellar masses, we analyse the evolution of all four mass scales (M_mean, M_med, M₅₀, and M_max, see Section 2.1.4 for definitions).

In Fig. 6, we plot the evolution of these mass scales as a function of SFE (all simulations are run to |$\mathrm{SFE}\gt 10{{\ \rm per\ cent}}$| unless SF is quenched earlier). As in Papers I and II, without feedback, all mass scales of the sink particles are significantly higher than observed, although the switch to explicit RHD without feedback does allow high-density regions to cool below 10 K, leading to somewhat lower mass stars. Note that this was not the case in Paper II, where radiation was not explicitly evolved and a 10 K floor was enforced (see Table 1). Introducing protostellar jets brings low and intermediate stellar masses in line with observed values, but SF is only suppressed, not quenched (Paper II), while massive stars undergo runaway accretion, which is apparent in both M₅₀ and M_max. With the introduction of radiative feedback, the overall mass scales are not affected, but massive stars no longer undergo runaway accretion and SF is quenched, leading to a final |$\mathrm{SFE}\sim 7{{\ \rm per\ cent}}$|⁠.

Figure 6.

The evolution of the number-weighted mean (M_mean = ∑M_sink/N_sink, top left-hand panel), number-weighted median (defined such that N_sink(M > M_med) = N_sink/2, bottom left-hand panel), mass-weighted median (M₅₀, the mass scale above which half of the total sink mass resides, top right-hand panel), and the maximum (M_max, bottom right-hand panel) sink mass as a function of SFE for the ‘physics ladder’ suite shown in Fig. 3. We also show with a shaded region the 95 per cent confidence interval for these values using the fitting function of K02. These are obtained by constructing a large set of IMFs whose parameters are sampled around the fiducial values with the uncertainties described by K02, then sampling each of these distributions up to the current total stellar mass in the simulation. For these sampled populations, we calculate the 95 per cent confidence interval of the various stellar mass metrics. All simulations are run to |$\mathrm{SFE}\gt 10{{\ \rm per\ cent}}$| unless SF is quenched earlier. 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.

We find that all four mass scales for all rungs of the ‘physics ladder’ exhibit an evolution with SFE; thus, the final SFE value where SF is quenched plays an important role in setting the final mass scales of the stellar mass spectrum. For the clouds in this suite (M2e4, see Table 1), stellar winds and SNe did not significantly alter the final SFE value (⁠|$\approx 7{{\ \rm per\ cent}}$| in both cases). Overall, for the runs with both jets and radiative feedback, we find that

$$\begin{eqnarray} M_{\rm mean}\propto \mathrm{SFE}^{0.5\pm 0.2} \\ M_{\rm 50}\propto \mathrm{SFE}^{0.8\pm 0.2} \\ M_{\rm med}\propto \mathrm{SFE}^{0.3\pm 0.2} \\ M_{\rm max}\propto \mathrm{SFE}^{0.8\pm 0.2}, \end{eqnarray}$$

(8)

where the uncertainty is the mean-squared fitting error.

3.2.2 Initial mass function

While the various characteristic mass scales provide some information on the sink (stellar) mass distribution, a holistic view of the sink mass spectrum (IMF) is necessary to understand the effects of each physical process. Fig. 7 shows the mass distribution of sinks at |$\mathrm{SFE}=7{{\ \rm per\ cent}}$|⁠, corresponding to the final SFE of the clouds with radiative feedback. The base isothermal MHD + gravity model produces an extremely top-heavy IMF with stellar masses a factor of 20 higher than observed. The introduction of detailed thermodynamics allows the gas to cool below the 10 K isothermal temperature limit in dense regions, leading to lower stellar masses (see Fig. 13 later for details). With the addition of protostellar jets, stellar masses are dramatically reduced (significantly more than the mass-loss from outflows) and the sink mass spectrum takes on a similar shape as the observed IMF. Due to the finite mass resolution of the simulations (⁠|$\mathrm{d}m = 10^{-3}\, \mathrm{M}_{\rm \odot }$|⁠), our IMF is incomplete in the brown dwarf regime (M < 0.08 M_⊙). Nevertheless, the IMF peak is essentially identical between all runs with protostellar jets at the same SFE value. Jets are essential in setting the IMF peak. Stellar radiation, winds, and SNe do not directly affect this mass scale; instead, they are mechanisms to quench SF, thus setting the final value of the SFE (see Section 3.1.1).

Figure 7.

Right-hand panel: distribution of sink particle masses measured in each simulation at 5 per cent SFE (= ∑M_sink/M₀) for the runs shown in Fig. 6. We also show the K02 fitting function for the IMF with a shaded region illustrating the uncertainties. Left-hand panel: the effective slope γ_{1, 10} of the sink mass distribution between |$1\hbox{ and }10\, \mathrm{M}_{\rm \odot }$| for the same runs as a function of SFE along with the 95 per cent observational confidence interval from K02. This high-mass slope is fairly stable for runs with jet feedback and produces values close to −2. Additional feedback (i.e. radiation, winds) steepen the slope as they are more dominant for massive stars.

The high-mass end of the IMF shows apparent deviations between models. However, this region is very sensitive to small number statistics (see Section 4.1), so we use the effective slope between |$1\hbox{ and }10\, \mathrm{M}_{\rm \odot }$| as a proxy to compare the high-mass end of the IMF (see Section 2.1.4). Fig. 7 shows that all runs with protostellar jets produce an effective slope close to −2. The addition of radiative and wind feedback suppresses the accretion of massive stars and prevents runaway accretion and thus the eventual flattening of the high-mass IMF; they, however, do not change the IMF slope value from −2. Note that the shallow effective slope for the no-jet runs (I_M, C_M) is due to the peak of the distribution being within |$1\!-\!10\, \mathrm{M}_{\rm \odot }$|⁠; at higher masses, they also produce a power-law tail of −2 slope (see Paper I and Guszejnov, Hopkins & Grudić 2018a for details).

4 SENSITIVITY TO INITIAL CONDITIONS AND PARAMETERS

In this section, we analyse how the results of our full physics runs (C_M_J_R_W, see Table 1) depend on changes in the ICs. We test for variations in the following initial parameters: the initial cloud surface density (Σ), virial parameter (α_turb), magnetization (μ), metallicity (Z), as well as the ISRF and turbulent driving; see Table 2 for specifics. Our aim is to formulate a general expression for how the IMF is affected by variations in different initial parameters and how these variations influence the star formation history of the cloud.

Table 2.

List of parameter variations investigated in Section 4 and the relevant IC/physics labels from Table 1.

Parameter	Default value	Tested variations	Labels in Table 1
Initial turbulence	αturb = 2 (Marginal boundedness)	x0.5, x2	M2e4_a1, M2e4_a4
Surface density	|$\Sigma =63\, \mathrm{M}_{\rm \odot }\,\mathrm{pc}^{-2}$| (MW average)	x10, x0.1	M2e4_R3, M2e4_R30
Cloud mass	|$M_0=2\times 10^4\, \mathrm{M}_{\rm \odot }$|	x0.01, x0.1	M2e2, M2e3
Mass-to-flux ratio	μ = 4.2 (1 per cent relative magnetic energy)	x0.3, x0.1	M2e4_mu1.3, M2e4_mu0.4
Interstellar radiation (ISRF)	Solar circle values (Habing 1968; Draine 1978)	x10, x100	ISRFx10, ISRFx100
Metallicity	Z = Z⊙	x0.1, x0.01	Z01, Z001

4.1 Difference between realizations

Before analysing runs with different initial parameters, we need to first examine what kinds of variations are possible for the same initial parameters. We ran three full physics (C_M_J_R_W) versions of an M2e4 cloud that had the same global parameters but used different initial turbulent realizations, i.e. the runs had different initial velocity fields even though the global turbulent parameters (velocity dispersion and power spectrum) were kept identical.

The grey lines in Fig. 8 show that the qualitative star formation history is similar between the runs; however, the evolution of the cloud’s virial state can show large variations between runs, as it is mainly set by the feedback of massive O stars, making the evolution of α highly sensitive to the formation times and masses of the most massive stars. We note that even though the runs had identical initial, global parameters, the final SFE values varied mildly between 6 per cent and 8 per cent. The cloud disruption time (i.e. time when it reaches α = 10) varies dramatically between runs, between 1.2 and 1.8 freefall times.

Figure 8.

The evolution of the SFE (left-hand panel) and the final sink mass distribution (IMF, right-hand panel) for runs with identical physics (C_M_J_R_W). We denote our fiducial M2e4 run (α_turb = 2, |$\Sigma =63\, \mathrm{M}_{\rm \odot }\,\mathrm{pc}^{-2}$|⁠, |$M=2\times 10^4\, \mathrm{M}_{\rm \odot }$|⁠) with a solid black line. Grey lines show runs with identical parameters but different turbulent realizations. We also show the results of runs with 2× higher and lower virial parameter α_turb and 10× higher and lower surface density Σ and for lower mass clouds (M2e2, M2e3). With the exception of the low surface density and low-mass cases, the star formation history is well described by a rising power law that flattens at different final values, ranging between 1 per cent and 15 per cent. Meanwhile, the stellar mass distribution (IMF) appears to be nearly invariant to variations in these ICs and agrees with the MW IMF within observed uncertainties (shaded area, K02) for stellar masses above M_complete.

Fig. 9 shows that although stellar masses start out similar in the simulations, variations of up to a factor of 2 can develop in M_med at fixed time. Similarly, there are small variations in the effective high-mass slope. The variations in M_max cause significant differences in cloud evolution (i.e. cloud disruption time), as cloud disruption is highly sensitive to the formation history of the most massive stars.

Figure 9.

The evolution of the median stellar mass M_med as a function of SFE (left-hand panel) and the evolution of the effective IMF slope γ_{1, 10} (right-hand panel) for the same runs as in Fig. 8 using the same notation. Although M_med varies significantly at fixed SFE, the different runs have different final SFE values, leading to similar final M_med values for all runs, well within the observational uncertainties. The effective IMF slopes in most runs are shallower than the canonical Salpeter (1955) value, but are also within observational uncertainties without a significant trend in any of the varied parameters.

Table 3 shows the summary IMF statistics from Section 2.1.4 in the three simulations that have identical global parameters but different initial turbulent realizations. We compare these with the values we obtain by sampling the K02 IMF fitting function between M_complete and |$150\, \mathrm{M}_{\rm \odot }$|⁠, while varying the IMF parameters within the uncertainties reported in K02. For all IMF statistics, we find that the runs with our fiducial parameters fall within the 95 per cent confidence intervals we get from sampling K02. In other words, the IMF produced by runs with our fiducial parameters is within observational uncertainties with the K02 IMF, although the resulting IMF slope is consistently shallower than the canonical value for all realizations. This means that the canonical −2.3 slope cannot be reproduced just by varying the initial turbulent realization.

Table 3.

Summary statistics of runs with identical parameters to the fiducial run but with three different turbulent realizations as well as the values expected from K02. Note that the statistics are all calculated between |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$| and |$150\, \mathrm{M}_{\rm \odot }$|⁠. For K02, the values and their errors are obtained by sampling the IMF at a fixed total stellar mass of |$1000\, \mathrm{M}_{\rm \odot }$| while varying its parameters, then taking the median value and the 95 per cent confidence intervals, respectively. For the simulations, we simply take the mean and standard variation of the values in the three runs. For all statistics, the simulated values fall within the confidence intervals we get from random sampling K02.

IMF statistics	Sampling K02 IMF	Simulations
Mmed (M⊙)	|$0.27_{-0.05}^{+0.10}$|	0.32 ± 0.02
Mmean (M⊙)	|$0.6_{-0.2}^{+1.0}$|	1.1 ± 0.1
M50 (M⊙)	|$1.4_{-0.8}^{+15.5}$|	4.3 ± 0.5
Mmax (M⊙)	|$52_{-40}^{+84}$|	47 ± 7
γ1, 10	|$-2.3_{-0.59}^{+0.54}$|	−1.93 ± 0.04
L/M (L⊙ M⊙−1)	|$900_{-800}^{+3900}$|	900 ± 100

4.2 Initial level of turbulence

In a turbulent medium, shocks can create self-gravitating overdensities that ultimately form stars. In a globally collapsing medium (i.e. our Sphere ICs, see Fig. 1), gravitational compression also triggers star formation. Since the cloud starts without a global infall motion, the initial level of turbulence (set by α_turb) determines how long the turbulent star formation channel dominates global collapse, as shown in Fig. 8. This is due to higher α_turb both enhancing the turbulent SF channel and slowing down the global collapse of the cloud. Previous work has found that in a turbulent medium without global collapse (i.e. with periodic boundary conditions), |$\mathrm{SFE}\propto \tilde{t}^2$| (Federrath & Klessen 2012; Murray et al. 2017, 2018 and Paper II), while our previous results showed |$\mathrm{SFE}\propto \tilde{t}^3$| for global collapse-dominated simulations (i.e. in an isolated cloud without external turbulent driving; see Paper II). So the net effect of a higher α_turb is delaying the |$\mathrm{SFE}\propto \tilde{t}^3$| regime, effectively lowering ϵ_ff. This delays star formation, and leads to an ultimately lower final SFE value, 10 per cent, 8 per cent, and 4 per cent for α_turb values of 1, 2, and 4, respectively.

Figs 8 and 9 show that although varying the initial α_turb turbulent virial parameter significantly affects the star formation evolution of the cloud, the final stellar mass scales and the IMF are insensitive to the initial level of turbulence in the cloud.

4.3 Cloud surface density

Surface density is considered to be a key parameter in determining the star formation history of a cloud (e.g. Fall, Krumholz & Matzner 2010; Grudić et al. 2018; Li et al. 2019). Fig. 8 shows that to be the case for our simulations as well. The average and high surface density runs produce similar star formation histories (within the uncertainties of Section 4.1); however, the low surface density run is dramatically different as it only has a single burst of star formation, because feedback from the newly formed stars easily disrupts the low-density cloud. Σ has a strong effect on the final SFE values giving 1 per cent, 8 per cent, and 14 per cent for Σ values of 6.3, 63, and |$630\, \mathrm{M}_{\rm \odot }\,\mathrm{pc}^{-2}$|⁠, respectively.

Despite the vastly different star formation history between the low surface density run and the others, the final stellar mass scales and spectra are essentially identical (Fig. 9). This is highly desirable in simulations as observed clouds have order-of-magnitude variations in surface density (Heyer & Dame 2015; Miville-Deschênes, Murray & Lee 2017), while the observed IMF is near universal (Offner et al. 2014).

4.4 Cloud mass

Observed molecular clouds have a large variety of masses from a few thousand to a million M_⊙ (Rice et al. 2016), but, due to computational costs, we are only able to probe the |$\le 2\times 10^4\, \mathrm{M}_{\rm \odot }$| range. Fig. 8 shows that for these relatively low-mass clouds, the star formation history is insensitive to the initial mass for the majority of their lifetime. We find that M₀ has a negligible effect on the final SFE values giving 7 per cent–8 per cent values for initial masses of 200, 2000, and |$2\times 10^4\, \mathrm{M}_{\rm \odot }$|⁠, respectively. It should be noted that SF is much more stochastic in the lowest mass cloud due to sampling effects.

Having such similar star formation histories, it is not surprising that the sink mass scales and spectra are also similar (Figs 8 and 9). The lowest cloud mass run (M2e2) does deviate from the higher mass ones, but the difference appears consistent with missing massive stars due to sampling effects. Overall, within the probed mass range, the initial cloud mass has no significant effect on any part of the star formation process. Note that this might not be true for more massive clouds, as their longer freefall times (assuming fixed surface density) could allow SN feedback to affect the cloud evolution (see discussion in Section 5).

4.5 Cloud magnetization

Magnetic fields provide support against collapse (Mouschovias & Spitzer 1976; Shu, Adams & Lizano 1987), and can affect the dynamics of turbulence and feedback in GMCs (Mac Low 1999; Krumholz, Stone & Gardiner 2007; Offner & Liu 2018). Fig. 10 shows that in our simulations, increasing the strength of the initial magnetic field (corresponding to a decrease in the mass-to-flux ratio μ) slows down global collapse and significantly reduces the star formation rate of the cloud. Due to the slower star formation rate, massive stars in the highly magnetized run have more time to unbind the cloud, resulting in the lower final SFE values of 8 per cent, 6 per cent, and 4 per cent for the 4.2, 1.3, and 0.4 mass-to-flux ratio runs, respectively.

Figure 10.

The evolution of the SFE (left-hand panel) and the final sink mass distribution (IMF, right-hand panel) for runs with identical physics (C_M_J_R_W). We denote our fiducial M2e4 run (μ = 4.2, Z/Z_⊙ = 1, e_ISRF = e_{ISRF, Solar}) with a solid black line. Grey lines show runs with identical parameters but different turbulent realizations. We also show the results of runs with 3× and 10× higher μ magnetic fluxes α_turb and 10× higher and lower surface density. The star formation history is well described by a rising power law that flattens at different final values, ranging between 3 per cent and 10 per cent. Unlike in Fig. 8, the stellar mass distribution (IMF) appears somewhat sensitive to these variations, specifically the number of high-mass stars (i.e. high-mass slope of the IMF). However, even these variations are within the observed uncertainties of the MW IMF (shaded area, K02) for stellar masses above M_complete.

Higher magnetization leads to a slower star formation and has a significant effect on the IMF. Figs 10 and 11 show that high magnetic fields significantly suppresses the formation of massive stars, steepening the high-mass slope to the canonical −2.35 value of Salpeter (1955) above |$10\, \mathrm{M}_{\rm \odot }$| in the μ = 0.4 case. This most magnetized run is also the only simulation in this work where jets, radiation, and winds are insufficient to unbind the cloud. Star formation is only quenched once the first SN explodes and disrupts the cloud.

Figure 11.

The evolution of the median stellar mass M_med as a function of SFE (left-hand panel) and the evolution of the effective IMF slope γ_{1, 10} (right-hand panel) for the same runs as in Fig. 10 using the same notation. While M_med is insensitive to changes in the initial magnetic field, it is higher at all times if the gas has a lower metallicity or a higher background radiation. The effective IMF slopes show mild variations, but only the |$z=0.01\, Z_\odot$| run produces slopes shallower than allowed by observational uncertainties.

The relatively minor changes in the IMF peak can be explained by Fig. 12, which shows how increasing the initial cloud mass-to-flux ratio only increases the magnetization in low-density (⁠|$\lt 10^3\, \mathrm{cm}^{-3}$|⁠) gas. Despite the different ICs, all runs saturate to the same B ∝ ρ^1/2 line, similar to the results of Mocz et al. (2017), Wurster, Bate & Price (2019), and Paper I. This can be understood as the effect of a turbulent dynamo enhancing the magnetic fields and driving the systems towards a common B–ρ relation at high densities (Federrath et al. 2014b). This relation roughly corresponds to v_A(ρ) ∼ c_{s, 0}, where v_A(ρ) is the local Alfvén velocity at density ρ, while c_{s, 0} is the isothermal sound speed, which is the relation we would expect from equipartition. A possible explanation is that the normalization of the B–ρ relation is enforced by a local dynamo effect (similar to the global α_B saturating in driven boxes; see Federrath et al. 2011) that is driven by the local gravitational collapse. In numerical experiments, β ∼ 1 is generally achieved for trans- or modestly super-sonic turbulence (Stone, Ostriker & Gammie 1998), which was indeed found on all scales in individual collapsed cores by Mocz et al. (2017). The fact that this B–ρ relation at high densities is independent from the initial magnetization explains why there are no differences in the mass scale of low-mass stars (which dominate M_med), as they accrete most of their material from their natal core that has n > 10⁵ cm⁻³. Massive stars in these simulations accrete material from scales much larger than their natal cores (Grudić et al. 2022), so changes in magnetic support at lower densities can affect their accretion flow (Lee et al. 2014).

Figure 12.

Magnetic field strength as a function of gas density one freefall time after the start of the simulation for runs with different initial normalized mass-to-flux ratios μ. Solid lines show the mass-weighted median values, while shaded regions denote the 1σ (68 per cent) intervals. To achieve satisfactory statistics, we stacked the distribution from five snapshots around the target simulation time. We also show B ∝ ρ^1/2 scaling law (in the isothermal regime, this corresponds to the Alfvén velocity v_A being an order unity times the thermal sound speed). Note that the simulations have a density threshold ρ_SF for sink particle formation that we mark with a shaded region, see Methods Paper for definition.

4.6 Cloud metallicity

The thermodynamic behaviour of the gas strongly depends on its elemental and dust abundances (metallicity), which is expected to strongly affect the mass scale of stars (Larson 2005; Sharda & Krumholz 2022). Also, metal line cooling becomes weaker at low Z, which increases the temperature of H ii regions (Osterbrock 1989), making feedback more mechanically efficient and reducing the SFE (He, Ricotti & Geen 2019). Fig. 10 shows that in our simulations with different metallicity values, the star formation histories are qualitatively similar, i.e. an initial |$\propto \tilde{t}^2$| phase followed by |$\propto \tilde{t}^3$|⁠, but the transition is delayed at low metallicity. Abundances also affect the final SFE, giving 8 per cent, 5 per cent, and 3 per cent for metallicity (Z/Z_⊙) values of 1, 0.1, and 0.01, respectively.

Despite the similar (but delayed) star formation rates, the stellar mass scales are significantly different. Fig. 11 shows that lower metallicities lead to increased stellar masses and the final stellar mass spectrum (IMF) is consistently more top heavy (i.e. has a shallower slope) for lower metallicity values. This is due to the less efficient cooling of the gas with the absence of metals.

Fig. 13 shows that lowering the cloud metallicity increases the gas temperature at densities above |$\sim 10^3\, \mathrm{cm}^{-3}$|⁠. This can be understood as lowering metallicity suppresses molecular line cooling, the dominant cooling channel at densities |$\sim 10^3\, \mathrm{cm}^{-3}$|⁠, as well as reducing dust density, which in turn reduces dust cooling, the dominant cooling channel above |$\sim 10^5\, \mathrm{cm}^{-3}$|⁠. We find the resulting temperature to roughly follow

$$\begin{eqnarray} T \propto Z^{-1/4}, \end{eqnarray}$$

(9)

which is roughly consistent with what one would expect from blackbody radiation from dust in an optically thin medium that is in equilibrium with the ISRF (i.e. e_ISRF × const. = j_cool ∝ ZT⁴).

Figure 13.

Median temperature of the gas at different densities one freefall time after the start of the simulation for different initial Z metallicity values. Solid lines show the mass-weighted median values, while shaded regions denote the 1σ (68 per cent) intervals. To achieve satisfactory statistics, we stacked the distribution from five snapshots around the target simulation time. Note that the simulations have a ρ_SF density threshold for sink particle formation that we mark with a shaded region.

4.7 Interstellar radiation field

The ISRF is the only external heating source in the simulation. The ISRF sets the temperature of the gas in the absence of other sources (i.e. before stars form). Increasing the ISRF increases thermal support in the cloud but has little effect on the speed of global collapse (see Fig. 10). Star formation rates between the fiducial and the 10 times higher ISRF runs are similar, but the highest, 100 times larger ISRF run has significantly higher star formation rates. This, in turn, means mildly higher final SFE values with higher ISRF: 8 per cent, 10 per cent, and 11 per cent for solar circle, 10 times higher, and 100 times higher values, respectively.

Fig. 11 shows that increasing the ISRF leads to higher stellar masses, due to higher gas temperatures, which lead to an increase in most characteristic mass scales (e.g. M_Jeans and |$M_{\rm BE}^{\rm turb}$|⁠). This effectively shifts the IMF to higher masses, leading to a shallower effective high-mass slope γ_{1, 10}, even though the actual shape (e.g. high-mass slope) of the IMF is largely unchanged (see Fig. 10).

Fig. 14 shows that increasing the ISRF increases the gas temperature at densities above |$\sim 10^5\, \mathrm{cm}^{-3}$|⁠. This can be understood as increasing the ISRF directly increases the heating radiation with which dust grains interact, effectively raising the dust temperature. This, in turn, raises the gas temperature above the density where gas becomes strongly coupled to dust (⁠|$n\gt 10^5\, \mathrm{cm}^{-3}$|⁠). We find the resulting gas temperature roughly as follows:

$$\begin{eqnarray} T \propto e_\mathrm{ISRF}^{1/4}, \end{eqnarray}$$

(10)

similar to equation (9), and it is roughly consistent with what one would expect from blackbody radiation from dust in an optically thin medium that is in equilibrium with the ISRF.

Figure 14.

Median temperature of the gas at different densities one freefall time after the start of the simulation for runs with ISRFs of different strengths. Solid lines show the mass-weighted median values, while shaded regions denote the 1σ (68 per cent) intervals. To achieve satisfactory statistics, we stacked the distribution from five snapshots around the target simulation time. Note that the simulations have a density threshold ρ_SF for sink particle formation that we mark with a shaded region.

4.8 Role of turbulent driving

In our previous works (Papers I and II), we found that the star formation history of a cloud significantly differed between the globally collapsing Sphere ICs (similar to Bate 2009) and driven, periodic Box ICs (similar to Federrath et al. 2014a; Cunningham et al. 2018). Here, we investigate the effects of turbulent driving and the periodic boundary condition by comparing our default full physics (C_M_J_R_W) Sphere run with two Box runs, one with continuously driven turbulence and one where turbulence is allowed to decay after SF starts.

Fig. 15 shows that the Sphere and Box runs exhibit the same star formation scaling that was found in the literature, |$\mathrm{SFE}\propto \tilde{t}^3$| and |$\mathrm{SFE}\propto \tilde{t}^2$|⁠, respectively (see e.g. Paper II and Murray et al. 2018). The Box run without driving follows an intermediate behaviour, initially following an |$\mathrm{SFE}\propto \tilde{t}^2$| trend, then switching over to |$\mathrm{SFE}\propto \tilde{t}^3$| as turbulence decays and global collapse starts, similar to the high-α_turb Sphere runs in Fig. 8. Note that due to the periodic boundary conditions in the Box runs, neither gas nor radiation can escape the cloud. This means that Box runs do not experience cloud disruption; instead, radiation rises in them to unphysical levels as star formation progresses. Thus, there is no physically meaningful ‘final SFE’ value for Box runs.

Figure 15.

The evolution of the SFE and the star formation rate per freefall time (ϵ_ff) as a function of time for C_M_J_R_W runs using Sphere IC, Box IC, and Box IC with decaying turbulence. Similar to Fig. 8, we denote the results for different turbulent realizations for the Sphere run with grey lines. Note that due to its periodic geometry, Box runs do not experience cloud disruption, and the runs are terminated once the simulated volume is filled with unphysical levels of radiation. In our fiducial Sphere runs, SF progresses as |$\mathrm{SFE}\propto \tilde{t}^3$|⁠; however, in the driven Box run, it only rises as |$\mathrm{SFE}\propto \tilde{t}^2$|⁠. We attribute this to the external driving and weaker gravitational focusing, as the decaying Box run transitions between the two regimes as its turbulence decays. Note that in the driven Box case, the star formation rate ϵ_ff is roughly steady, while in the other cases, it varies orders of magnitude over the cloud lifetime.

Figs 15 and 16 show that there is a significant discrepancy between the stellar spectra of the Sphere and Box runs, with Sphere runs producing about a factor of 2 higher stellar masses. Thus, the effective slope γ_{1, 10} is steeper for the Box runs than observed but still within the limits of sampling uncertainty. Stellar masses in the Box run without turbulent driving start out similar to the driven case, but quickly switch over to the same track as the Sphere run. Note that observed molecular clouds experience both global collapse and external driving (Heyer & Dame 2015), so the behaviour of a realistic cloud would likely lie between the tracks of the Sphere and driven Box runs. So we conclude that an accurate modelling of external driving of turbulence in clouds is necessary for any simulation to reproduce the observed IMF (see recent work by Lane et al. 2022 for such a model).

Figure 16.

The evolution of the number-weighted median stellar mass (M_med, left-hand panel), the effective IMF slope (γ_{1, 10}, middle panel) as a function of SFE, and the final sink mass distribution (IMF, right-hand panel). The driven Box run exhibits lower stellar masses and a significantly steeper high-mass slope.

4.9 Scaling relations

From the previous results, we can formulate a general expression for the dependence of the IMF and the star formation history of the cloud as a function of initial global parameters. Specifically, we attempt to formulate scaling relations for the M_med median stellar mass (above |$0.1\, \mathrm{M}_{\rm \odot }$|⁠), the γ_{1, 10} effective IMF slope, and the terminal SFE of the cloud. This is done by carrying out least-squares fitting for the dependence of each individual global parameter while marginalized over the rest. We estimate the error in M_med and γ_{1, 10} by bootstrapping: We resample the sink mass distribution at fixed total sink mass and calculate the 66 per cent (1σ) confidence intervals. We obtain the scaling relations shown in Table 4. We can simplify this scaling relation by ignoring all low and insignificant exponents, i.e. all exponents γ with fitting error σ, where either |γ| < 2σ or |γ| < 0.1 + σ. We find that

$$\begin{eqnarray} M_{\rm med}\sim \mathrm{const.}, \end{eqnarray}$$

(11)

in other words, the median stellar mass is insensitive to all varied parameters. For γ_{1, 10}, we find that it varies as

$$\begin{eqnarray} e^{\gamma _{1,10}} \propto \mu ^{0.28\pm 0.16} Z^{-0.13\pm 0.10} \end{eqnarray}$$

(12)

while for the terminal SFE, we get

$$\begin{eqnarray} \mathrm{SFE}_\mathrm{final} \propto \alpha _{\mathrm{turb}}^{-0.72\pm 0.24}\, \Sigma ^{0.59\pm 0.17}\, \mu ^{0.19\pm 0.06}\, Z^{0.26\pm 0.04}. \end{eqnarray}$$

(13)

Table 4.

List of exponents obtained by least-squares fitting for the final SFE value, median stellar mass M_med, and the γ_{1, 10} effective high-mass slope IMF in Section 4.

Parameter	Final SFE exponent	Mmed exponent	|$e^{\gamma _{1,10}}$| exponent
Initial turbulence (αturb)	−0.72 ± 0.24	−0.12 ± 0.03	0.06 ± 0.27
Surface density (Σ)	0.59 ± 0.17	0.01 ± 0.01	−0.07 ± 0.10
Cloud mass (M0)	0.03 ± 0.02	0.0 ± 0.03	0.23 ± 0.19
Mass-to-flux ratio (μ)	0.19 ± 0.06	0.10 ± 0.01	0.28 ± 0.16
Interstellar radiation (ISRF)	0.06 ± 0.02	0.11 ± 0.01	−0.03 ± 0.06
Metallicity (Z)	0.26 ± 0.04	−0.08 ± 0.02	−0.13 ± 0.10

4.10 Variations in the light-to-mass ratio

We note that even though the IMF appears to vary mildly between the runs presented here, other summary statistics of star formation, like the light-to-mass ratio L/M, can vary significantly. Fig. 17 shows that although the set of runs shown in Sections 4.2 and 4.4 exhibit statistically indistinguishable IMFs, their final light-to-mass ratios vary by orders of magnitude. This is because L/M is highly sensitive to the most massive stars due to the steep scaling of stellar luminosity with stellar mass. Overall, |$L/M\mathrel {\propto }_{\sim } M_{\rm max}$|⁠, and their values fall within the 95 per cent confidence interval we obtain by sampling the K02 IMF. We also find that in our simulations, clouds with a higher final total stellar mass have a higher M_max and thus a higher L/M. We do not have sufficient statistics to rule out the existence of a high-mass cut-off for the IMF that depends on ICs (see e.g. Weidner & Kroupa 2006).

Figure 17.

Left-hand panel: the evolution of the light-to-mass ratio as a function of time for the same set of runs as in Fig. 8. The shaded region denotes the 1σ (68 per cent) intervals for L/M values obtained by random sampling the K02 IMF while varying its parameters within observed uncertainties. While these runs have statistically indistinguishable IMFs (see Fig. 8), their L/M values vary significantly between runs. Right-hand panel: the final L/M values and M_max values for all runs presented in Section 4. The shaded region shows the 95 per cent confidence intervals for sampling the canonical K02 IMF. L/M and M_max are highly correlated, so L/M is probing the highest masses of the IMF, making it sensitive to sampling effects.

5 DISCUSSION

A key goal of the STARFORGE project is to understand the roles different physical processes play in star formation by carrying out a suite of simulations with increasingly complex physics. A significant advantage of this suite, relative to comparing results from various groups in the literature, is that they use the same code base and ICs, allowing for a cleaner comparison.

5.1 Role of magnetic fields, thermodynamics, and stellar feedback

Our suite proceeded through the ‘physics ladder’ of star formation, starting from isothermal MHD + gravity, then adding non-isothermal gas thermodynamics, protostellar jets, stellar radiation, winds, and SNe.

Similar to other recent work in the literature (e.g. Haugbølle et al. 2018), we found that in STARFORGE runs, the magnetic fields impose a well-defined mass scale on the stellar mass spectrum (Paper I), preventing the runaway fragmentation found in non-magnetized, isothermal runs (Guszejnov et al. 2018b). This mass scale, however, is an order of magnitude higher than observed for MW-like conditions. It is also sensitive to ICs (e.g. surface density) in a way that could not be reconciled with the apparent universality of the IMF in the MW (Offner et al. 2014; Guszejnov, Hopkins & Ma 2017).

Due to the highly efficient cooling of molecular gas, most star formation models assume the gas to be isothermal (Girichidis et al. 2020). Detailed studies, however, showed that there is a significant scatter in the gas temperature, with a clear density dependence (see Glover & Clark 2012). At high densities, the isothermal assumption inevitably breaks down as the gas becomes opaque to its own cooling radiation, forming a hydrostatic Larson core (Larson 1969). This transition to an adiabatic behaviour can impose a mass scale on the stellar mass spectrum (Low & Lynden-Bell 1976; Rees 1976), but it was found to be significantly below the observed characteristic stellar mass scales. But recent works proposed that tidal screening around Larson cores can raise the relevant mass scales to be close to observed values (Lee & Hennebelle 2018; Colman & Teyssier 2020). Note that these works all pertain to non-magnetized clouds without feedback, so the only unique mass scale is imposed by the isothermal–adiabatic transition. In Paper II, we investigated these effects for magnetized clouds and found that transitioning to non-isothermal thermodynamics had little effect on the stellar mass spectrum, which was still predominantly set by magnetic effects, leading to stellar masses significantly above those observed. In this work, we carried out RMHD simulations where heating and cooling radiation is explicitly followed, unlike in Paper II. This allows the gas to self-shield, leading to temperatures below 10 K that is the fiducial temperature in isothermal approximations and was used as a floor temperature in Paper II. This leads to a reduction of stellar masses, but the overall mass scales (with only these physics) remain significantly above the observed values (see Fig. 7).

Previous work in the literature has shown that protostellar jets significantly affect star formation. Jets directly reduce stellar masses and slow down star formation in the cloud (Cunningham et al. 2011; Hansen et al. 2012; Federrath et al. 2014a) and can potentially drive turbulence on small scales (e.g. Nakamura & Li 2007; Wang et al. 2010; Offner & Arce 2014; Offner & Chaban 2017; Murray et al. 2018). In Paper II, we showed that protostellar jets play a vital role in setting stellar masses, reducing them by an order of magnitude. This reduction is significantly larger than what one would expect if the jets simply removed of the order of 1/3 of the accreted material, as jets not only remove accreted material, but also disrupt the accretion flows around protostars, allowing the nearby ISM to fragment and form new stars. Overall, in Paper II, we found that jets bring the stellar mass spectrum in line with observations in the MW, with the exception of the most massive stars. In this work, we reran the same simulations with explicitly evolved radiative heating and cooling (RHD) and found the results to be largely the same. In both works, protostellar jets significantly affect the virial state of the cloud and can suppress the formation of new stars. However, they are unable to expel the remaining gas and to prevent massive stars from accreting it, leading to runaway accretion in high-mass stars.

Radiative feedback from stars has long been theorized to play an important role in setting stellar mass scales (e.g. Krumholz 2011; Bate 2012; Myers et al. 2013; Guszejnov et al. 2016; Li et al. 2018, see Hennebelle et al. 2020 for a counterexample) and regulating star formation (e.g. Offner et al. 2009; Krumholz et al. 2012; Cunningham et al. 2018). Many of these previous works have produced mass spectra similar to the observed IMF, generally by including most of the relevant feedback processes (jets, radiation) in magnetized clouds (e.g. Cunningham et al. 2018). This work expands upon these results in several ways. First, the clouds in our simulations are an order of magnitude more massive than those in previous works with similar physics (e.g. Cunningham et al. 2018). Secondly, the ‘physics ladder’ suite allows us to disentangle the effects of the individual feedback processes. We see that radiative feedback plays a key role in disrupting the cloud and quenching star formation. In all our simulations, the formation of the first main-sequence O star marks a turning point in the global evolution of the cloud and the beginning of the disruption process. Note that most of the aforementioned previous works in the literature simulated much smaller clouds (⁠|$200\!-\!1000\, \mathrm{M}_{\rm \odot }$|⁠), so no massive stars formed in them. Stellar winds further enhance feedback from massive stars, but do not significantly alter it as they are only a significant channel of momentum feedback for massive O stars. Radiative feedback and winds both counteract the runaway accretion on to massive stars found in jet runs, both by shutting off accretion and by generally expelling gas from the cloud. Therefore, jets are responsible for setting the peak of the IMF, but radiative and wind feedback are responsible for preventing the high-mass end of the IMF from flattening. Note, however, that even with stellar feedback the resulting high-mass IMF slopes are consistent with the −2 value expected from scale-free fragmentation.

Finally, massive stars end their lives as SNe. These explosions are generally agreed to be critical for regulating star formation on galactic scales, and in particular to dominate the overall momentum input in the ISM by stellar feedback (Somerville & Davé 2015; Naab & Ostriker 2017; Vogelsberger et al. 2020), although non-linear interactions between different processes are also important (e.g. Hopkins et al. 2014, 2018a). However, they occur fairly late in the star formation process; simulations of cluster formation have found that they have negligible effects upon SFE and bound cluster masses compared to early feedback, even in massive GMCs that survive long enough to host SNe before disruption (Grudić et al. 2021b). In this work, we similarly find that SNe occur after the cloud has been completely disrupted by earlier feedback processes. In almost all of our simulations, the cloud is disrupted within 0.5–1.0 freefall times after the first O star forms, the only exception being the μ = 0.4 highly magnetized run, which does not completely disrupt until the first SNe go off. Note that for massive clouds (⁠|$10^5\!-\!10^6\, \mathrm{M}_{\rm \odot }$|⁠), we could plausibly expect SNe to trigger before the cloud is disrupted, even for less magnetized clouds. Even if SNe turn out to have no direct role in regulating SF in massive clouds, they are likely to have a major indirect role as they are thought to be one of the main drivers of galactic turbulence and thus set the properties of GMCs (e.g. Hopkins et al. 2011, 2012; Ostriker & Shetty 2011; Faucher-Giguère et al. 2013; Walch et al. 2015; Martizzi et al. 2016; Padoan et al. 2017; Seifried et al. 2018; Gurvich et al. 2020; Guszejnov et al. 2020a).

5.2 Environmental variations

In Section 4, we analysed the star formation history and stellar mass spectrum in simulations that include all levels of feedback physics (C_M_J_R_W) while varying various initial cloud parameters (see Papers I and II for similar studies for the lower rungs of the ‘physics ladder’).

We find that the final SFE values of the clouds significantly depend on initial cloud parameters, specifically the initial surface density, level of turbulence, magnetization, and metallicity (i.e. gas temperature; see equation 13), varying between 1 per cent and 12 per cent in the simulated clouds. These trends are in agreement with expectations from simple rule-of-thumb considerations, such as higher surface density making it harder for feedback to unbind the cloud, thus leading to a higher SFE, just as increased initial turbulence or magnetic support makes it easier to unbind the cloud, lowering the final SFE. Lowering the metallicity of the initial gas or lowering the ISRF lowers the final SFE, although the exact mechanism is unclear.

Regarding the sink mass spectrum (IMF), the initial cloud surface density and virial parameter have little effect on the final median stellar mass scale (see equation 11). This appears to contradict Papers I and II, where we found both parameters to significantly affect the stellar mass spectrum. The apparent contradiction is resolved by taking into account that those works lacked the relevant feedback physics to quench star formation and disrupt the cloud; thus, all comparisons were done at fixed SFE. For fixed SFE, the runs presented in this work show similar variations. However, when comparing the final, post-disruption IMFs, we find that the dependence of the final SFE on cloud properties effectively cancels these variations. Insensitivity to both the α_turb virial parameter and Σ surface density is required if we are to reproduce the near-universal IMF of the MW (Offner et al. 2014), as observed molecular clouds in similar regions exhibit an order-of-magnitude scatter in α_turb and a factor of 5 in Σ (Heyer et al. 2009; Kauffmann, Pillai & Goldsmith 2013; Heyer & Dame 2015). Following Table 4, we predict mild variations less than few tens of per cent in M_med, well within the observational uncertainties of the MW IMF. Decreasing metallicity or increasing the local ISRF raises the gas temperature, which in turn increases the relevant mass scales of star formation (e.g. Jeans mass, sonic mass). This leads to an increase of the median stellar mass; however, the shift is very small, roughly consistent with |$M_{\rm med}\mathrel {\propto }_{\sim } Z^{-1/10} e_\mathrm{ISRF}^{1/10}$|⁠. This is consistent with the weak trend predicted by Sharda & Krumholz (2022) for the characteristic stellar mass for |$Z\gt 0.01\, Z_\odot$|⁠. Previous simulations that included only a subset of the physics presented here (e.g. no MHD, jets, or winds) found similarly no significant IMF variations with metallicity (Bate 2019). Overall, M_med varies very weakly with initial gas parameters, consistent with the observed limited variations in the stellar IMF. The high-mass slope of the IMF is more sensitive to ICs, steepening with increasing initial magnetization and becoming more shallow for lower metallicity values.

Although our simulations only tested the effects of mild variations in initial parameters, we can extrapolate them to the more extreme star-forming regions, such as the Central Molecular Zone of the MW, starburst galaxies, or high-redshift galaxies. These regions have surface densities a factor of 100–1000 higher than in the MW (Solomon et al. 1997; Swinbank et al. 2011) and an ISRF that is a factor of 100 higher. While we plan to simulate star formation in such environments in the future, for now we can make a rough estimate of the IMF with Table 4 and find M_med to be within a factor of 2 of the MW value.

5.3 Caveats

While the simulations presented here are the current state of the art for simulating star-forming clouds, like other simulations in the literature, STARFORGE employs a large number of significant approximations and assumptions to make the simulations computationally tractable (see Methods Paper for detailed discussions). In particular, the runs used here have an |$\sim 30\, \mathrm{au}$| Jeans resolution, i.e. fragmentation on scales smaller than this is not resolved. This has a dramatic effect on the formation of protostellar discs and their fragmentation, causing the simulation to potentially miss closely formed binaries and overestimate stellar masses. However, we do not expect it to qualitatively affect the IMF above M_complete (see Appendix  A), except for potentially steepening the high-mass slope as massive stars are broken up through disc fragmentation.

Recent observations of dwarf galaxies (Hunter et al. 2021, 2022; Elmegreen, Martinez & Hunter 2022) showed that on galactic scales (∼400 pc), the ISM velocity dispersion correlates best with the local star formation rate with a 100 Myr delay. Meanwhile, in our simulations, clouds are destroyed in about two freefall times (corresponding to roughly 10 Myr) and the kinetic energy of the gas increases by an order of magnitude (see Fig. 5). Due to the relatively small size of the simulated volume and the simplified modelling of the surrounding ISM, the distribution and dissipation of kinetic energy in the ISM after cloud disruption is not captured accurately, and will be revisited in future work.

6 CONCLUSIONS

In this work, we presented simulations from the STARFORGE project, which are high-resolution RMHD simulations following the evolution of star-forming molecular clouds. The runs include progressively more complex physics, starting from isothermal MHD, then enabling explicitly solved heating and cooling radiation and adding stellar feedback in the form of protostellar jets, radiation, stellar winds, and SNe. Building on our past work, we investigate each rung of this ‘physics ladder’ of star formation to identify the role each process plays in star formation.

In previous works, we showed that isothermal MHD leads to a well-defined stellar mass spectrum (Paper I), and that the addition of protostellar jets is necessary to bring these scales in line with observations (Paper II). The runs presented in this paper reinforce those conclusions: Stellar mass scales are set by MHD turbulence that both creates the self-gravitating structures and prevents their runaway fragmentation (see the non-magnetized case in Guszejnov et al. 2018b). Protostellar jets dramatically reduce stellar mass scales by both directly removing accreted material and by disrupting the accretion flow around stars; however, they cannot prevent the most massive stars from undergoing runaway accretion. In these runs, radiation was explicitly evolved, allowing gas to cool below the isothermal temperature limit in dense regions. This leads to a significant reduction in stellar masses, which was not captured in Paper II. The addition of stellar radiation, winds, and SNe has little direct effect on the stellar mass spectrum, apart from preventing the runaway accretion of massive stars. They, however, play a dominant role in regulating star formation. In the presented runs, stellar radiation and protostellar jets are the dominant forms of feedback that quench star formation and disrupt the cloud, with the formation of the first main-sequence O star marking the turning point in the cloud’s evolution. While SNe do go off in these simulations, these exclusively happen at the end of the runs when the cloud has already been disrupted by radiative feedback. It should be noted that our simulations followed |$\le 2\times 10^4\, \mathrm{M}_{\rm \odot }$| clouds with lifetimes of |$\sim 7\, \mathrm{Myr}$|⁠, so it is possible that SNe play a significant role in more massive clouds whose lifetimes are longer, which we plan to explore in future work.

In addition to the ‘physics ladder’ suite, we present a suite of full physics simulations with varied initial parameters to determine how the star formation history and stellar mass spectrum (IMF) depend on ICs. The characteristic stellar masses are insensitive to the initial cloud mass, surface density, and level of turbulence. Note that this only applies to the final, post-disruption mass spectrum; comparisons at fixed times or SFE values show significant differences. Of the parameters probed in this study, the IMF peak is only affected by the cloud metallicity and the strength of the ISRF. Since both significantly alter the thermodynamics of the cloud, we conjecture that their effects can be attributed to a change in the mean temperature of star-forming gas. Meanwhile, the high-mass slope of the IMF becomes steeper with decreasing metallicity or increasing ISRF and magnetization. The scaling relations derived from our parameter study predict IMF variations that are within the observational uncertainties of the near-universal IMF observed in the MW.

SUPPORTING INFORMATION

STARFORGE_IMF_paper_extra_plots.zip

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

DG is supported by the Harlan J. Smith McDonald Observatory Postdoctoral Fellowship and the Cottrell Fellowships Award (#27982) from the Research Corporation for Science Advancement. Support for MYG was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51479 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. SSRO was supported by NSF CAREER Award AST-1748571, NASA grant 80NSSC20K0507, NSF grant 2107942, and by a Cottrell Scholar Award (#24400) from the Research Corporation for Science Advancement. CAFG was supported by NSF through grants AST-1715216, AST-2108230, and CAREER award AST-1652522; by NASA through grants 17-ATP17-0067 and 21-ATP21-0036; by STScI through grants HST-AR-16124.001-A and HST-GO-16730.016-A; by CXO through grant TM2-23005X; and by the Research Corporation for Science Advancement through a Cottrell Scholar Award. Support for PFH was provided by NSF Collaborative Research Grants 1715847 and 1911233, NSF CAREER grant 1455342, and NASA grants 80NSSC18K0562 and JPL 1589742. ALR acknowledges support from Harvard University through the ITC Post-doctoral Fellowship. This work used computational resources provided by XSEDE allocations AST-190018 and AST-140023, 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 paper are available on reasonable request to the corresponding authors. Additional figures can be found at our GitHub repository (https://github.com/guszejnovdavid/STARFORGE_IMF_paper_extra_plots). 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

https://github.com/guszejnovdavid/STARFORGE_IMF_paper_extra_plots

3

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

4

https://github.com/mikegrudic/MakeCloud

5

Note that due to the finite resolution of the simulation, sink particles may not represent individual stars, but we expect that to be the case above our |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$| completeness limit.

6

Note that there are several other fitting functions in the literature (e.g. Chabrier 2005; De Marchi, Paresce & Portegies Zwart 2010; Parravano, McKee & Hollenbach 2011; Dib, Schmeja & Hony 2017), but the differences between them are not significant for the purposes of this paper.

APPENDIX A: RESOLUTION EFFECTS ON THE SINK MASS SPECTRUM

The fiducial resolution level in our simulations is |$\mathrm{d}m=10^{-3}\, \mathrm{M}_{\rm \odot }$| (equivalent to |$\mathrm{d}x_\mathrm{Jeans}\sim 20\, \mathrm{au}$|⁠, see Methods Paper), a choice based on our previous work in Papers I and II, where this value was sufficient for the sink mass spectrum to be complete down to |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$|⁠. With the transition to the new RHD-based thermodynamics module and the inclusion of stellar radiation and winds, it is worth re-examining this choice. We do so by running a ‘full physics’ (C_M_J_R_W) simulation at different resolution levels. Since RHD simulations with M₀/dm ≫ 2 × 10⁷ are prohibitively expensive, we choose to do this resolution study on a smaller M2e3 cloud (see Table 1) within a resolution range of dm ∈ [10⁻², 10⁻⁴].

Fig. A1 shows that the star formation history of the cloud is sensitive to numerical resolution in the examined range. The final sink particle number and the evolution of the cloud virial state are virtually identical between our fiducial resolution and the 10 times higher value.

Figure A1.

The evolution of the star formation rate per freefall time (ϵ_ff, left-hand panel), number of sink particles (N_sink, middle panel), and the virial parameter (α, right-hand panel) as a function of time for an M2e3 cloud with all physics included (C_M_J_R_W, see Table 1).

Fig. A2 shows that the mean, median, and maximum sink masses are essentially identical between the fiducial |$\mathrm{d}m=10^{-3}\, \mathrm{M}_{\rm \odot }$| and the higher resolution run. Note that these metrics are calculated above the same |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$| as in the rest of the paper, regardless of the mass resolution of the simulation.

Figure A2.

The evolution of the number-weighted mean (M_mean, top left-hand panel), number-weighted median (M_med, top right-hand panel), mass-weighted median (M₅₀, bottom left-hand panel), and maximum (M_max, bottom right-hand panel) sink mass as a function of time for an M2e3 cloud with all physics included (C_M_J_R_W, see Table 1). Shaded regions show the 95 per cent confidence intervals. Note that all mass scales are calculated for sinks above our chosen limit of |$M_{\rm complete}=0.1\, \mathrm{M}_{\rm \odot }$|⁠.

Fig. A3 shows the mass distribution of stars at the end of the simulations. As expected, the spectrum extends to lower masses with higher resolution; however, the part above |$\sim 0.1\, \mathrm{M}_{\rm \odot }$| is identical between the fiducial and the high-resolution run. This is similar to the results we obtained in our previous works (Paper II), leading to a conservative, rule-of-thumb estimate of our completeness limit as M_complete ∼ 100dm. For our fiducial resolution this yields M_complete = 0.1 M_⊙, which we adopt as our completeness limit for this work.

Figure A3.

The sink mass spectrum at the end of the simulation for an M2e3 cloud with all physics included (C_M_J_R_W, see Table 1) at different mass resolutions.

Author notes

NASA Hubble Fellow.