STARFORGE: Toward a comprehensive numerical model of star cluster formation and feedback

Requirements for a complete dynamical model of star formation and feedback

While some progress has been made with simpler models that consider only e.g. turbulence and gravity (Padoan, Nordlund & Jones 1997;Padoan & Nordlund 2002;Krumholz, McKee & Klein 2005;Hennebelle & Chabrier 2008;Padoan & Nordlund 2011;Hopkins 2012;Hennebelle & Chabrier 2013), other physics are likely to be important. In particular, feedback is important for understanding the end-point of star formation (the disruption of GMCs), and its implications for other questions such as the IMF and stellar multiplicity have only begun to be explored. Many analytic and semi-analytic calculations of the effects of feedback in GMCs have been performed (for reviews, see McKee & Ostriker 2007;Krumholz, McKee & Bland-Hawthorn 2019), yielding useful dimensional arguments and analytic insights. But GMCs are complex, turbulent, inherently threedimensional entities that evolve on their internal crossing time-scale (Larson 1981;Mac Low & Klessen 2004). Thus, even under the gross simplification of treating GMCs as isolated entities (i.e. neglecting galactic environment), (semi-)analytic predictions inevitably hinge on many highly uncertain assumptions. With so much parameter freedom it is difficult to say whether a given model is correct for the right reasons, limiting physical insight and ultimately predictive power. Direct numerical simulations of star formation are a necessary tool to resolve these uncertainties.

In the past two decades, great progress has been made incorporating stellar feedback into direct numerical simulations of star-forming GMCs (for reviews see Dale 2015;Krumholz et al. 2019). But these studies have shown that further progress on the key questions of star formation requires next-generation simulations that do all of the following:

(i) Resolve individual star formation: Many simulations of star cluster formation do not attempt to resolve the formation and ongoing accretion of individual stars across the entire stellar mass range, instead relying on a sub-grid SF prescription that either enforces a certain underlying IMF directly or is fine-tuned to recover the observed one (Colín, Vázquez-Semadeni & Gómez 2013;Howard, Pudritz & Harris 2017;Kim et al. 2017;Sormani et al. 2017;Su et al. 2018;Grudić et al. 2018a;He, Ricotti & Geen 2019;Lahén et al. 2019;Wall et al. 2020). But there is an infinite number of ways to do this, each with different implicit assumptions about how star formation works, and the choice of prescription can have a major effect upon simulation results (Grudić & Hopkins 2019), limiting predictive power. Simulations should ideally attempt to resolve the formation and accretion of individual stars (or sink particles), and to recover a realistic IMF self-consistently from physical (not numerical) processes. This is obviously necessary anyway if one wishes to study the physical origins of the IMF and stellar multiplicity.

(ii) Follow stellar dynamics: SF simulations that do not integrate stellar orbits explicitly generally discretize the stellar mass formed into a collisionless fluid represented by gravitationally softened particles (e.g. Grudić et al. 2018a;Lahén et al. 2019;Li et al. 2019), which can produce qualitatively correct star cluster density profiles (Grudić et al. 2018b;Lahén et al. 2020), but have the severe limitation that the collisionless description (and phase-space density conservation) breaks down on mass scales M ࣠ 100 M , so if cluster formation is a hierarchical assembly from smaller masses (e.g. Bonnell, Bate & Vine 2003), then individual stellar dynamics is always important at some stage in the process. A simulation must also treat dynamics on the scale of binary separations to accurately predict stellar multiplicity, let alone phenomena such as common disc accretion (e.g. Lee et al. 2019;Muñoz, Miranda & Lai 2019;Duffell et al. 2020).

(iii) Follow MHD, chemistry, and cooling: Obviously, following the dynamics of GMCs, star formation, and accretion requires gasdynamical simulations, and stars cannot form if gas cannot radiatively cool. Moreover, the ISM is magnetized, and this fact can easily have important implications for star formation. The magnetic field can act as an additional source of support, potentially stabilizing otherwiseunstable cores (Chandrasekhar 1951;Mouschovias & Spitzer 1976), affecting the IMF (Price & Bate 2007;Guszejnov et al. 2020b), the rate of star formation (Federrath & Klessen 2012;Federrath 2015), and altering the pressure balance, morphology, growth of instabilities, and transport of energy in feedback bubbles (Krumholz, Stone & Gardiner 2007;Offner & Liu 2018;Krumholz & Federrath 2019).

(iv) Scale up to massive GMCs: Current star formation simulations that do both (i) and (ii) have focused upon lower mass systems, simulating gas masses of 100-1000 M (Federrath 2015;Cunningham et al. 2018;Jones & Bate 2018;Lee & Hennebelle 2018;Li, Klein & McKee 2018;Wurster, Bate & Price 2019;Colman & Teyssier 2020), producing ∼10-100 M in stars. Low-mass clusters are important to model, as they can be readily compared to well-studied sites of star formation in the Solar neighbourhood (e.g. Evans, Heiderman & Vutisalchavakul 2014), but the overwhelming majority of SF in our Galaxy occurs in massive complexes with gas mass 10 5 M (McKee & Williams 1997;Murray & Rahman 2010). Simulated lowmass clusters are also less likely to host massive ( 10 M ) stars, and hence cannot be used to study massive SF. 1(v) Account for all major feedback channels: 3D MHD simulations of SF have not generally incorporated all known dynamically important feedback mechanisms (jets, winds, full-spectrum radiation, and SNe) acting in concert. A comprehensive treatment of feedback is needed because different feedback channels are effective on different scales, and can interact non-linearly. For example, direct radiation pressure from a massive star is ineffective if it couples deep within the star's potential well (Krumholz 2018), and radiation pressure in general may be subdominant to protostellar outflows for regulating the growth of individual massive stars (Rosen & Krumholz 2020). But by regulating accretion or punching optically thin holes, outflows could help photons to couple their momentum farther away from the star, eventually allowing them to disrupt the host GMC (Fall, Krumholz & Matzner 2010;Murray, Quataert & Thompson 2010;Hopkins, Quataert & Murray 2012;Raskutti, Ostriker & Skinner 2016;Kim, Kim & Ostriker 2018;Grudić et al. 2018a;Hopkins & Grudić 2019;Hopkins et al. 2020a). Meanwhile, jets can be a powerful feedback mechanism that can regulate star formation on the 1 pc scales of individual cores and dense clumps (Matzner & McKee 2000;Nakamura & Li 2007;Wang et al. 2010;Cunningham et al. 2011;Federrath 2015;Offner & Chaban 2017;Cunningham et al. 2018), but may have only weak effects upon the gas kinematics at larger (i.e. 10 pc) scales within the GMC (Murray, Goyal & Chang 2018). Many other synergies between feedback mechanisms can also be theorized.

Enter STARFORGE

In this work, we introduce the STARFORGE (STAR FORmation in Gaseous Environments) project, 2 a new initiative to perform next-generation 3D star cluster formation simulations in massive GMCs. The STARFORGE numerical framework that we have implemented in the GIZMO code (Hopkins 2015, hereafter H15) simultaneously follows the formation, accretion, and dynamics of individual stars in massive GMCs, with optional physics modules capable of accounting for all of the most widely discussed stellar feedback mechanisms (jets, radiative heating and momentum, stellar winds, and SNe), satisfying the requirements (i)-(v) laid out above. In Guszejnov et al. (2020b, hereafter Paper 0), we used numerical simulations (using an early version of the methods presented here) to show that the simple recipe of isothermal MHD turbulence and gravity fails to yield a realistic IMF and star formation history in Milky Way conditions, motivating the need for additional physics implemented in STARFORGE. In the present paper (Paper 1), we present and test the numerical methods of STARFORGE, permitting simulations that combine all of the important SF physics discussed here into a realistic simulation of GMC evolution and star cluster formation. In Guszejnov et al. (2021, hereafter Paper 2), we use the algorithms presented here to explore the effects of protostellar jets upon SF across an unprecedented parameter space of GMC properties and jet model parameters.

This paper is structured as follows. In Section 2, we present the 'core' algorithms that any 3D star cluster formation simulation must have in some form: MHD, gravity, sink particle methods, and an integration scheme that couples gas and stars stably and achieves acceptable accuracy in stellar dynamics. In Section 3, we describe the treatment of cooling, chemistry, and thermodynamics (treating the opacity limit and protostellar heating either with selfconsistent radiative transfer, or simple inexpensive approximations). In Section 4, we describe and test algorithms for the numerical coupling of stellar feedback in the form of jets, winds, SNe, and radiation. In Section 5, we explore various potential applications of these methods beyond isolated GMC simulations, and also enumerate the many caveats and uncertain assumptions inherent in simulating SF and feedback in this manner, motivating future work. In Section 6, we summarize our findings and outline the programme of the STARFORGE project.

C O R E A L G O R I T H M S F O R S TA R F O R M AT I O N

The STARFORGE framework is implemented in the GIZMO multimethod, multiphysics N-body and MHD simulation code (Hopkins 2015). 3 GIZMO was selected for the project for several reasons. It implements second-order, Galilean-invariant, Lagrangian meshless finite-volume MHD methods (Hopkins & Raives 2016, hereafter HR16), which have several useful advantages for SF problems (discussed in Section 2.1). It includes a gravity solver (Section 2.2) that is spatially adaptive (solved consistently with the MHD discretization) with near-ideal scaling up to 10 6 cores (Hopkins et al. 2018b). In addition to solving the MHD equations, GIZMO's meshless discretization and reconstruction schemes provide a flexible framework for solving additional, non-core physics such as diffusion, conduction, and non-ideal MHD terms (Hopkins 2017), radiative transfer (Hopkins & Grudić 2019;Hopkins et al. 2020a), and stellar feedback (Hopkins et al. 2018a). All equations are integrated according to a flexible, hierarchical powers-of-two time-stepping scheme (Section 2.3) that makes it possible to follow processes over a wide range of time-scales, from the ∼ 10 Myr lifetime of a GMC to a 1yr binary orbit.

Conceptually, our approach follows previous Lagrangian 3D star formation simulations (Klessen & Burkert 2000;Bate, Bonnell & Bromm 2003): we discretize the mass of the GMC and the surrounding medium into discrete elements of mass m, and integrate their evolution in time according to the MHD equations.

3 http://www.tapir.caltech.edu/ ∼ phopkins/Site/GIZMO.html Eventually the self-gravitating MHD equations can no longer be followed self-consistently at the centres of runaway core collapse, so we replace these centres with sink particles (Bate, Bonnell & Price 1995) nominally representing individual protostars. These sink particles interact with the gas via gravity, accretion, and optionally feedback, with feedback rates determined by a sub-grid model of (proto-)stellar evolution based upon that used in Offner et al. (2009). We target a MHD resolution scale on the order of a few 10 au, comparable to state-of-the-art low-mass star cluster formation simulations. We defer physics on smaller scales (e.g. disc formation, accretion, jet launching, and protostellar evolution) to a sub-grid approach, acknowledging the various caveats that this entails (Section 5.2).

We provide a glossary of the various numerical resolution-related quantities in Table 1.

Magnetohydrodynamics

The default MHD solver used by STARFORGE simulations is the Meshless Finite Mass (MFM) method presented in HR16, 4 which we will briefly summarize. This method discretizes the fluid into a finite number of gas cells of mass M i , each representing a domain of volume V i = M i /ρ i as determined by the kernel 5 volume partition described in H15. This partition defines the effective face areas A gg between each interacting pair of gas cells g and g , 6 between which the conservative MHD equations are evolved in standard finitevolume fashion:

where (V U) g gives the usual conserved quantities (mass, momentum, energy, etc.) integrated over the volumetric domain of the cell, and F gg is the tensor of their fluxes. The fluxes are obtained by solving the appropriate (HLLD) Riemann problem using the fluid states reconstructed at the interface according to a slope-limited, second-order least-squares gradient estimator, evaluated in the frame moving with the interface to ensure Galilean invariance. In MFM, the interfaces are defined and move such that the mass flux vanishes identically, so the method follows the motion of constant-mass, quasi-Lagrangian fluid elements. Cells exchange conserved quantities ensuring machine-precision conservation in this operation. Magnetic field divergence errors are controlled by augmenting equation (1) with the usual Powell et al. (1999) and Dedner et al. (2002) source terms and using the Hopkins (2016) constrained gradient method for obtaining the consistent fluid reconstruction operator that minimizes the numerically unstable terms.

Because the volume partition associated with each cell can have complicated shapes (see Hopkins 2015 for discussion), it is useful to define an effective cell size x g ≡ V 1/3 g ≡ ( m g /ρ g ) 1/3 (the equivalent cell side length for a cubic cell of the same volume and 4 MFM is our method of choice, but all STARFORGE methods are compatible with any quasi-Lagrangian MHD method implemented in GIZMO, including MFV and SPMHD, enabling easy comparisons. 5 Following HR16, we adopt the M4 cubic spline as the default kernel partition function W gg = W (|x gx g |, H g ), with kernel radius of compact support H g defined recursively by H g = 2 x g where x g is the kernel-weighted mean cell separation: 6 Throughout this work, we adopt index notation for gas cells and sinks where i and j denote any element regardless of type, g and g denote gas cells, and s and s denote sink particles.

Table 1. Glossary of numerical resolution-related quantities in STARFORGE simulations. M -3 is the mass resolution m in units of the fiducial 10 -3 M resolution, n 3 = X H ρ/ M p ≈ 0.7ρ/ M p is the local number density of H in units of 10 3 cm -3 , and c s, 0.2 is the minimum gas isothermal soundspeed in units of 0.2 km s -1 .

Symbol

and volume-equivalent spherical radius

where the latter expressions are given in terms of the typical STARFORGE mass resolution of 10 -3 M and the number density of H atoms n H,g ≈ 0.7ρ g /m p . We emphasize, as discussed in HR16, that MFM has little in common with smoothed-particle MHD (SPMHD) -MFM is formally a member of the class of Arbitrary Lagrangian-Eulerian (ALE) finite-volume Godunov methods, much more closely related to Voronoi moving-mesh methods (e.g. Springel 2010;Duffell & MacFadyen 2011), and in fact reduces to a Voronoimesh method in the limit of sharply peaked kernel functions with exact volume quadrature.

Meshless, Lagrangian, Galilean-invariant MHD methods have several advantages for simulating SF in GMCs with feedback. In Lagrangian methods, Galilean invariance implies that the time-step required does not depend upon the bulk flow velocity v as t ∝ x/v (as is required for stable advection in Eulerian fixed-grid methods), so larger time-steps are possible in the highly supersonic flows of GMCs, and the presence of very high ( 100 km s -1 ) bulk velocities in accretion flows or winds does not incur such a high cost, an issue often encountered by Eulerian simulations combining high velocities with high-spatial refinement levels.

Galilean and rotational invariance also ensure that structures formed in the simulation (e.g. dense cores and clumps) to evolve internally in a manner independent of their mean bulk velocity and orientation with respect to the coordinate axes (to machine precision). Numerical errors are velocity-independent, and a significant source numerical diffusivity in supersonic flows in Eulerian methods (the grid advection operation, Robertson et al. 2010;Pontzen et al. 2021) is absent. These advantages can enable more rapid convergence of phenomena involving highly supersonic flows, large density con-trasts, angular momentum conservation, and coupling to self-gravity, all of which are highly relevant in SF. We refer the reader to H15, HR16, and Hopkins (2016) for demonstrations of the performance of the MFM MHD method in a wide variety of standard test problems.

Gravity

We compute the gravitational field g = -∇ , tidal tensor T = -∇∇ (where = ∇ -2 4πGρ is the gravitational potential) and the gravitational jerk j at the location of every gas cell and sink particle in the simulation using a modified version of the massively parallel, approximate tree-force algorithm introduced in Springel (2005, hereafter S05). This algorithm recursively subdivides the simulation domain into an oct-tree structure, and uses the monopole approximation of the field contribution of the contents of a tree node, unless an the opening criterion is satisfied, in which case the opening criteria are re-evaluated recursively for all sub-nodes and forces evaluated accordingly. We use the acceleration-based opening criterion introduced in S05 (which requires that the quadrupole error term a Q ∼ GML 2 /r 4 of a node is less than a specified fraction of the total field g), but also always require the original Barnes & Hut (1986) opening criterion: a tree node is always opened if it subtends an angle θ ≡ L r > , where L is the side length of the node, r is the distance between the node centre of mass and the target point for field evaluation, and = 0.5 is the maximum opening angle. This ensures that a dense sub-system of a hierarchically structured system that dominates its own field (e.g. a dense clump or a binary) still has some control over the accuracy of the force contribution from surrounding material, which is still needed to predict its centreof-mass motion (Grudić et al. 2020). T and j are computed in the same pass through the gravity tree as g, summing the respective monopole contributions of tree nodes and particles according to the same opening criterion (Vogelsberger et al. 2008;Grudić & Hopkins 2020). Gravitational forces are updated for gas cells only as frequently as required per the Grudić (2020) adaptive criterion, using j to construct a predictor of g between updates. This generally decreases overall cost of calling the gravity solver by a factor of 2 or better.

Softening

We use a softened form of the gravitational force law for sink particles or gas cells that fall within each other's respective softening radii S i and S j , i.e. r ij < max S i , S j , ensuring that all interactions are antisymmetric, conserving total linear and angular momentum. Softening for gravitational interactions between gas cells is fully adaptive, i.e. we set S i = H i so that the gravitational force resolution is always scaled to be consistent with the cell volume partitioning assumed when solving the MHD equations. This prevents various unphysical effects that are seen if the hydro and gravitational resolution scales are mismatched (Bate & Burkert 1997). At the second-order consistency of our MHD method, H15 showed that this can be done by using the same spherically symmetric compact spline softening scheme as smoothed particle hydrodynamics (SPH), with the same additional terms and symmetrization to ensure conservation of momentum and energy (Price & Monaghan 2007). The full form of the pairwise force law between gas cells is given in H15 (equations H8-H10).

For interactions between sink particles, it would be ideal to use the full, unsoftened 1/r 2 force law to be able to follow stellar dynamics on all scales down to stellar radii. Unfortunately, this is not presently possible for our code, because binaries with arbitrarilyclose separations (down to surface contact) can form and harden dynamically (Heggie 1975), potentially requiring very short timesteps. In theory, this workload could be accomplished by our hierarchical individual time-stepping scheme (Section 2.3), but in practice the massively parallel architecture of GIZMO is such that global overheads tend to eventually bottleneck time-steps involving a small subset of the total particle number (e.g. two sink particles in a very short-period binary). Therefore, we adopt a finite, fixed softening radius S for sink particles in a given simulation, allowing us to follow collisional dynamics accurately on spatial scales S , while limiting the effective hardness of binaries having separation ࣠ S .

Lastly, softened interactions between fixed-softening sinks and adaptively-softened gas cells must be handled specially, because the respective softenings can differ by orders of magnitude -e.g. if a star with S ∼ 20 au is moving through a diffuse part of the GMC where n H ∼ 10 cm -3 , and hence a gas cell would have size ∼ 0.1 pc (equation 2). In such a case using the same Price & Monaghan (2007) symmetrization as gas-gas interactions (averaging the forces) would result in unphysical noise, because the interaction between the gas and star on the scale of S is totally unresolved, but the back-reaction on the star depends sensitively upon its position with respect to the cell centre. Instead, we take the maximum of S and the gas kernel radius H as the softening radius, in both directions of the pairwise interaction (thus conserving momentum). A natural choice for the fixed S is to match it to the finest possible Jeans-resolved spatial MHD resolution, which we will show to be ∼ 20 AU in Section 2.4 for our fiducial mass resolution of 10 -3 M .

In all pairwise interactions, the tidal tensor T and jerk j (for sinks) are summed using spatial derivatives of the same softened force kernel that is used for g, with the same symmetrization scheme used for that particular interaction.

Time-step criteria

Gas cells obey all of the standard local, Galilean-invariant, MHDspecific time-step criteria given in HR16, except that we neglect the gravitational component of the acceleration in the Power et al. (2003) acceleration criterion. Instead, both gas cells and sink particles obey the tidal time-step constraint (Grudić & Hopkins 2020):

where T denotes the Frobenius norm of the tidal tensor and η is a tolerance parameter controlling the overall accuracy of integration. In Grudić & Hopkins (2020), we showed that T encodes a reliable estimate of the local gravitational dynamical time t dyn ∼ -1 = r 3 /GM, respecting the equivalence principle (invariance to the addition of a uniform external field g ) and interpolating between appropriate limits for a continuous mass distribution and in the vicinity of a point mass (e.g. sinks) more accurately and robustly than the usual acceleration-based criterion.

Sink particles also obey their own special time-step criterion for ensuring orbital integration accuracy (von Hoerner 1960):

where

and

where j runs over all other sink particles, = h /2.8 is the Plummerequivalent sink softening radius, and r ss , v ss , M s , m ss are the separation, relative velocity, and respective masses of sink particles s and s . Note that t 2-body is simply the harmonic mean of a kinematic orbital crossing time-scale ∼r/v and an orbital dynamical timescale t dyn = -1 ∼ r 3 /GM tot , but replacing r with a softened version r 2 + 2 . We treat this as a single time-step criterion using the harmonic mean of the two time-scales because the smooth interpolation in the regime t c, min ∼ t dyn, min yields slightly better integration accuracy for eccentric binary orbits. Unlike the tidal criterion (equation 4), t 2-body is symmetric between pairs of sinks, ensuring that binaries are updated synchronously when this is the dominant time-step constraint, which can give better conservation of orbital parameters. The global min operations can be evaluated efficiently in the pass through the gravity tree, combining stellar masses that exist within the same tree node if it is not opened (consistent with the force approximation).

Sink particles also observe various time-step constraints derived from local gas conditions, to ensure the stability of local gas interactions occurring within the hydrodynamic stencil, such as accretion and feedback injection. First, it cannot time-step more than 4× the smallest time-step of a gas neighbour:

where g runs over all overlapping, potentially interacting gas neighbours, i.e. r sg < max (H s , H g ). This is analogous to the constraint imposed for neighbouring gas cells in GIZMO, following Saitoh & Makino (2009). A sink particle's time-step is also constrained to anticipate the infall and/or orbital motion of surrounding gas, via a gas free-fall time criterion:

where η is the parameter controlling overall integration accuracy and x s is the effective gas cell size in the vicinity of the sink. Sinks also obey a local Courant-Friedrichs-Lewy (CFL)-type time-step constraint:

where v s is the velocity of the sink, c s, s , v A, s , and v gas,s are the local gas sound speed, Alfvén speed, and gas velocity (reconstructed using a simple kernel-weighted interpolation), c is the (possibly reduced) speed of light (only included if radiative transfer is enabled), and v fb is an estimate of the maximum velocity of gas emerging from the sink due to feedback:

where v SN is the SN ejecta velocity given by equation (47) (or 0 if the sink is not currently going SN), v wind is the stellar wind velocity (equation 45), and v shell is the greater of the velocity of an energy-conserving (Weaver et al. 1977) or momentum-conserving (Steigman, Strittmatter & Williams 1975) shell as its radius reaches the resolution scale x:

where Ṁwind is the sink's wind mass-loss rate (equation 32), v wind is the wind velocity (equation 45), L wind = 1 2 Ṁwind v 2 wind is the mechanical luminosity of the wind, L is the bolometric luminosity of the sink, ρ is the local gas density. We have found that including something like the v fb term in t CFL, can be important to prevent the sink particle from 'overshooting' the amount of feedback it injects, i.e. injecting feedback over a time-step longer than the time required for the local gas cells to react to it, leading to an unstable solution. 7Reciprocally to the v fb term in equation ( 10), gas cells also obey a time-step constraint to anticipate the arrival of feedback from a star:

where s runs over all sink particles and the min can be evaluated efficiently in the gravity tree pass. If a hyperbolic RT solver (e.g. M1) is enabled, we also enforce a local radiation CFL condition:

Likewise, we enforce appropriate local time-step criteria in the relevant methods papers for the various optional physics (e.g. nonideal MHD and dust) described above.

Time integration

Given a choice of individual time-step as in Section 2.3.1, we require a time-integration scheme that achieves acceptable truncation error.

The error budget of a multiphysics SF simulation is dominated by errors in MHD, radiative transfer, stellar evolution/feedback, and gravity, all of which are necessarily approximate and/or have large modelling uncertainties. Some errors may not even converge away in the limit of infinite resolution: e.g. moments-based radiative transfer methods will not converge to the exact radiative transfer solution in general. Hence, for gas, high-order integration schemes are unlikely beat down the leading-order error terms in the global GMC and star cluster evolution. For all gas cells, and as a robust fall-back option for stars in special circumstances, we use the standard second-order Kick-Drift-Kick (KDK) integrator (S05):

where a i is the total gravitational+MHD+radiative acceleration of cell/particle i, which is re-evaluated after every initial half-step kick.

Some additional control on orbital integration accuracy for stars is needed to preserve the properties of binaries once formed. Any numerical integration scheme will incur a certain fractional energy error E/E per orbit (as well as an angular momentum error J/J and phase error φ/φ, but here we use E/E as an overall proxy for integration error, as is standard). A true symplectic integrator such as the leapfrog with constant time-steps will preserve orbital energy and angular momentum on average, but true symplecticity is lost once the adaptive KDK version is adopted and errors accumulate over time, causing the semimajor axis to change with each orbit (S05). If a fairly typical ∼ 100 yr, e = 0.9 binary is to survive the ∼ 10 Myr lifetime of its host GMC in a simulation, then we require | E/E| 1 and hence an energy error per orbit of 10 -5 . In Fig. 1, we show that this would require 2000 time-steps per orbit with the KDK integrator (and would demand a minimum time-step at periastron of ࣠ 1 d), which we have found to demand an excessivelylarge overhead. Instead, we adopt a modified version of the fourthorder Hermite integrator (Makino & Aarseth 1992) for stars. At the beginning of the time-step, we evaluate the initial accelerations a s,0 and jerk j s,0 of all sinks in a special sinks-only gravity tree pass. We then perform the initial prediction step:

re-evaluate a s and j ,s using the new positions and velocities, and then perform the correction step:

where the order here matters because the update to x s requires the updated version of v s . This is subtly different from the original implementation in Makino & Aarseth (1992), in that we perform two force/jerk evaluations per time-step, one at the beginning of the timestep and one after the prediction step, whereas the original Hermite scheme only re-evaluates the force/jerk after the prediction step. We discovered serendipitously that this small modification gives a scheme that converges at the same order, but can give order-ofmagnitude smaller energy errors at fixed time-step size in binary integration (Fig. 1). In a typical direct N-body application the entire cost of the simulation is force/jerk evaluation and there is not much parallelization overhead, so this advantage would be nullified by simply taking 2× smaller time-steps at equal cost. In GIZMO, the force/jerk comes relatively cheaply, but there can be significant global overheads involved in taking smaller time-steps, so our modified Hermite scheme is more suitable. For our standard choice of η = 0.01, this method achieves a relative energy error of <10 -6 per orbit for an e = 0.9 binary (and this decreases steeply for smaller e).

In a given time-step, a sink particle is first provisionally timestepped according to the KDK scheme, co-evolving it alongside the gas update cycle so that the gas-star coupling seen by the gas is unaltered by the Hermite integration (but saving the initial state of the time-step). At the end of the time-step, the sink particle is eligible to accrete gas cells. If it does, low-order integration errors are introduced and j ceases to be well-defined, so we simply keep the more-robust KDK result for that time-step. If it does not accrete, it is eligible to take a Hermite time-step, updating via equation ( 16) using the saved values x s,0 , v s,0 , a s,0 , and j s,0 , and performing the subsequent force evaluation and correction step (equation 17). Given Relative energy error accumulated per orbit integrating e = 0.9 binary motion with the second-order Kick-Drift-Kick (KDK) (S05), fourthorder Hermite (Makino & Aarseth 1992) and our modified Hermite integrator (equations 16-17), as a function of the time-step tolerance parameter η, which controls the number of steps taken per orbit per our adaptive time-stepping scheme (equations 4 and 5 and powers-of-two block scheduling, Section 2.3). Conserving binary properties in a ∼ 10 Myr GMC simulation ( 10 4 -10 7 binary orbits) is only practical with a higher order scheme. Both Hermite schemes happen to converge at fifth order in this problem when using our time-step criteria, and our modified version performs better at fixed η, for an extra force and jerk evaluation per time-step.

the order of the MHD reconstruction, and the inability to define the jerk given e.g. shocks, there would be no gain from using this integrator for gas.

Existing tests

The standard MHD and gravity algorithms in GIZMO have been extensively tested and applied to hundreds of different problems in the literature, so we will not repeat these. We do note these tests have demonstrated that our default implementation can simultaneously accurately evolve phenomena including gas in regular or warped Keplerian discs, strong interacting shocks, current sheets and flux tubes, supersonic and subsonic turbulence, fluid mixing instabilities (Kelvin-Helmholz, Rayleigh Taylor, etc.), multifluid dustgas dynamics, collisional+collisionless gravitational dynamics, and reproduces the correct linear growth rates of the magnetorotational instability (MRI) and non-ideal Hall MRI and anisotropic MHD instabilities (magnetothermal, heat-flux-bouyancy) (Hopkins 2015;Hopkins & Raives 2016;Zhu & Li 2016;Lupi, Volonteri & Silk 2017;Deng et al. 2019a, b;Moseley et al. 2019;Rennehan et al. 2019;Hu & Chiang 2020;Panuelos, Wadsley & Kevlahan 2020). Tests of idealized problems involving self-gravitating MHD including the Evrard (1988) problem (spherical collapse of a self-gravitating polytrope), the MHD Zel'dovich (1970) pancake (self-gravitating collapse of an initially linear density perturbation along one dimension in a 3D Hubble flow) demonstrate that the MFM/MFV methods in GIZMO (as well as related moving-mesh methods) converge much more rapidly than popular AMR or SPH methods applied to the same problem (Hopkins 2015;Hopkins & Raives 2016;Hubber, Rosotti & Booth 2018).

Several studies have argued that the most notable advantages of MFM compared to SPH or AMR methods may come in astrophysical discs, which are crucial for the physics of stellar accretion but are often marginally resolved in our simulations (meaning that more-rapid convergence at fixed resolution is especially valuable). For example, (1) MFM accurately conserves angular momentum and prevents both unphysical disc 'spreading' and/or clumping/fragmentation via artificial viscous instabilities in SPH or catastrophic angular momentum loss from spurious coordinatealignment torques which are inescapable in AMR (Hopkins 2015;Few et al. 2016;Zhu & Li 2016;Lupi et al. 2017;Panuelos et al. 2020;Deng, Ogilvie & Mayer 2021). ( 2) Few et al. (2016) found MFM more rapidly converges to correct linear growth rates for spiral arms and other disc instabilities, compared to AMR or SPH, while Deng et al. (2021) found a similar result for physical parametric instabilities of warped discs. (3) Deng, Mayer & Meru (2017) showed MFM was the only method surveyed which exhibited convergence to exact solutions for gravitoturbulent fragmentation in cooling discs. (4) MFM, at fixed resolution, has been shown to more accurately capture boundary-layer mixing in discs (especially those formed via collisions), avoiding artificial suppression of subsonic turbulence and mixing common in SPH (Zhu & Li 2016;Deng et al. 2019b). ( 5) Hubber et al. (2018) demonstrated that MFM simulations of 'gap opening' via massive planets or binaries in discs converged more rapidly and maintained the gaps more accurately than equivalent SPH or AMR simulations (which tended to produce artificially-high torques and therefore stellar accretion rates in this regime).

Isothermal collapse tests

Next, we consider a variant of the 'isothermal test case' from Boss & Bodenheimer (1979): a uniform-density, un-magnetized, spherical solar-mass core with initial radius 5 × 10 16 cm, in uniform rotation with = 7.2 × 10 -13 rad s -1 and a 10 per cent m = 2 azimuthal density perturbation with an isothermal equation of state P = c 2 s ρ, c s = 0.166 km s -1 (Burkert & Bodenheimer 1993;Bate & Burkert 1997;S05). We use the MFM hydrodynamics solver with the default STARFORGE gravity and time-stepping set-up (Sections 2.2-2.3), and initialize the cells in a glass configuration with the density perturbation imposed by rescaling cell masses. In Fig. 2, we plot the maximum density in the simulation as a function of time while varying the average cell mass m from 10 -3 -10 -7 M , and compare with SPH results from Bate & Burkert (1997) andS05. The solution appears to converge to an answer in fair agreement with the highest resolution SPH results in S05. Moreover, our solution with N gas = 10 4 cell resolution is closer to its respective converged limit than SPH simulations with 3.34 × 10 4 and 8 × 10 4 particles, respectively. However, at low enough resolution numerical effects become apparent, as evidenced by the ∼ 10 per cent delay of the collapse of our lowest resolution run with only 10 3 gas cells. It is important to assess the effects of resolution upon SF simulations, as this will inform our sink particle prescription. A common convergence parameter for self-gravitating isothermal hydrodynamics simulations is the number of Jeans lengths per cell (Bate & Burkert 1997;Truelove et al. 1997;Hubber, Goodwin & Whitworth 2006):

The consequences of underresolving λ J (i.e. allowing x λ J ) vary from method to method, and have been the subject of extensive study. Truelove et al. (1997, hereafter T97) found that Eulerian grid simulations that do not enforce f J < 1 4 are subject to artificial fragmentation (AF), wherein fragments of unphysical origin can form even in a smooth, symmetric collapse. A similar effect is seen in SPH simulations if care is not taken to match the gravitational resolution to the hydrodynamic resolution (Section 2.2.1), e.g. adopting a constant gas softening length that is much smaller than the particle spacing (Bate & Burkert 1997). Clearly, AF is undesirable, so a variety of approaches have been developed to prevent it, e.g. by fine-tuning the sink particle formation, accretion, and merger criteria in conjunction with the refinement scheme (e.g. Krumholz, McKee & Klein 2004;Haugbølle, Padoan & Nordlund 2018). AF does not occur in SPH simulations that maintain consistency between gravitational and hydrodynamic resolution (Bate & Burkert 1997;Whitworth, Boffin & Francis 1998;Hubber et al. 2006), and more recently it has been confirmed that this is true for MFM as well in the linear Jeans problem (Hubber et al. 2018;Yamamoto et al. 2021 in preparation). With these methods, fragments that should physically collapse but are insufficiently Jeansresolved either do not collapse, or simply collapse more slowly. We plot a surface density map of the filament formed in the T97 self-gravitating isothermal hydrodynamics test problem, at the time that the maximum density is ρ max = 10 -9.5 g cm -3 (cf. T97 Fig. 4). From top left to bottom right, we vary the mass resolution m from 8 × 10 -5 -1.6 × 10 -7 M , which varies the maximum number of Jeans wavelengths per cell

max from 1.1 to 0.14. Failure to resolve the Jeans length simply coarse-grains the structure of the filament -there is no evidence of artificial fragmentation when the Jeans length is poorly resolved.

Here, we also check for AF in the exact test problem simulated in T97, a variant of the Boss & Bodenheimer (1979) problem, using an initial Gaussian density profile. With a 10 per cent m = 2 initial density perturbation, T97 found that the converged solution is the formation of a single collapsing filament, but if the Jeans resolution criterion x > λ J /4 was violated then they would obtain an unphysical solution containing two filaments instead. In Fig. 3, we plot the structure formed in the simulation at the time that the maximum density exceeds 10 -9.5 g cm -3 , at a variety of mass resolutions such that the T97 criterion is strongly violated at our lowest resolution (8 × 10 -3 M , x ≈ 1.1λ J ), and is well satisfied at our highest (1.56 × 10 -7 M , x ≈ 0.14λ J ). No additional filament or fragment forms even when the T97 criterion is strongly violatedthe effect of poor resolution appears consistent with a simple spatial coarse graining of the structure of the filament. T97 also found that the version of the problem with no initial density perturbation resulted in the formation of numerical fragments, unless f J < 1/4 was enforced. We have verified that this is not the case for MFM: axisymmetry is preserved accurately throughout the collapse, even when the Jeans resolution criterion is strongly violated.

Our findings for MFM appear consistent with previous results in the linear Jeans problem (Hubber et al. 2018;Yamamoto et al. in preparation): unstable scales that are well-resolved (f J 1) collapse as they should, and scales that should be stable are stable. Marginally resolved (f J ∼ 1) unstable wavelengths are either artificially stabilized, or collapse more slowly than is physical (e.g. the lowest resolution run in Fig. 2), and these effects converge away with sufficient resolution.

Resolution criteria

What density-and length-scales should then be considered 'resolved' in isothermal self-gravitating flows? This depends on what threshold value of f J is considered acceptable for the question at hand, which is generally problem-dependent with no one straightforward answer. But assuming we do adopt a certain maximum f J, max to demarcate the boundary of 'trusting' results in a certain problem, the maximum Jeans-resolved density is,

and the minimum Jeans-resolved cell length is

We caution that direct comparisons of the 'resolved' density-or length-scale between SF simulations should ideally be made at fixed f J, max (i.e. correcting by appropriate f J, max factors), and that even then there can be many confounding factors when comparing across different methods.

This discussion of Jeans resolution neglects magnetic fields, which can supplement thermal pressure as a source of support against gravi-tational collapse. For the purposes of gravitational stability analyses, it effectively adds the Alfvén speed v A = |B|/ √ μ 0 ρ to the thermal sound speed c s in quadrature, i.e. c s → c 2 s + v 2 A = 1 + 2 β c s , modulo some geometry-specific O (1) factors in the β term (Chandrasekhar 1951;Mouschovias & Spitzer 1976). Assuming that the convergence parameter for isothermal, self-gravitating MHD is instead the magnetic Jeans number obtained by substituting the above into equation ( 18), as has been argued in various works (Federrath et al. 2010;Myers et al. 2013), our assessment of the resolving power of the simulations (equations 19 and 20) is overly conservative. However, the densest gas in isothermal MHD core collapse attracts towards β ∼ 1 (Mocz et al. 2017;Wurster et al. 2019;Guszejnov et al. 2020b), so the corrections to our analysis from magnetic fields are expected to be modest for the present purposes. Even if not, this would merely make our effective f J threshold more conservative, so e.g. our sink algorithm would not follow gas as far into the marginally resolved regime, and our simulations are better-resolved than as quoted in equations ( 19) and (20).

Sink particles

We use sink particles to model the accretion, dynamics, and feedback of individual stars and protostars (e.g. Bate et al. 1995;Krumholz et al. 2004;Federrath et al. 2010;Hubber, Walch & Whitworth 2013;Bleuler & Teyssier 2014). A sink particle represents a designated region in the domain of the simulation in which physical processes are considered unresolved, and are relegated to sub-grid prescriptions. The general strategy is to put a sink in the centre of a collapsing core once the collapse process can no longer be followed self-consistently by the MHD scheme, and to allow this sink to accrete subsequent infalling material according to certain physically motivated rules.

Our sink implementation formally distinguishes between resolved accretion, i.e. the actual mass transfer from the gas in the simulation domain to the sink particle, and unresolved accretion: the transfer of mass from the sink's internal gas reservoir (comprising unresolved gas in the envelope or the protostellar disc) on to the protostar itself (and potentially into the protostellar outflow). Other works equate the two types of accretion, often assuming that gas removed from the simulation domain arrives at the protostar immediately (e.g. Krumholz et al. 2004), or using a detailed sub-grid model to decide how rapidly resolved accretion should occur (Hubber et al. 2013). For us it is important to model accretion on to the protostar distinctly from resolved accretion into the sink region, because we discretize resolved accretion into quanta -the mass resolution m -but would like a continuous estimate of the protostellar accretion rate for modelling the protostellar evolution and accretion luminosity. 8 Our algorithm is most similar to that of Bate et al. (1995), with some additional rules for formation 8 We have experimented with our own implementation of the algorithm of Hubber et al. (2013), which uses an estimate of Ṁ that interpolates between disc-like and Bondi-like regimes based on local gas kinematics, and uses that estimator to directly determine how much gas should be removed. However, we have found that in some problems the estimator of Ṁ used to set the rate of resolved accretion can underestimate the actual accretion rate of the surrounding flow, so mass piles up within the softening radius of the sink particle and the actual accretion rate ends up being set by the need to remove gas cells on too small a time-step (circumventing the normal criteria), defeating the purpose of trying to estimate and enforce the proper accretion rate as determined by physical processes. One potential issue is that the α-disc parameter used in the disc-like regime must be known a priori, otherwise the accretion rate will not match the boundary flow. This will generally vary with Figure 4. Diagram of the flow of mass due to accretion and feedback, as managed by our sink particle algorithm. We follow gas cells of mass m until they satisfy all sink particle accretion criteria (Section 2.5.2) and they are transferred to the sink's accretion reservoir representing the envelope or disc gas mass present on scales < R sink . Mass is accreted from the reservoir towards the protostar according to the smoothed accretion prescription (equation 32), and if protostellar jets are enabled a fraction of this mass f w is diverted to the jet reservoir. The rest arrives at the star, and mass is transferred from the star to the wind reservoir according to the wind mass loss rate (which is set to an extremely large value (with appropriate velocity) if the star goes SN, equations 48 and 47). The jet and wind reservoirs return gas to the simulation domain via their respective feedback channels (waiting until a sufficient mass is available to inject or spawn their respective mass quanta). and accretion, and some additional modelling of protostellar accretion and feedback. We sketch the flow of mass dictated by our algorithm in Fig. 4, and describe the algorithm in detail in this section.

Formation

A gas cell is eligible to turn into a sink particle if and only if it satisfies the following criteria:

(i) Density threshold: The gas cell is denser than a density threshold ρ th , which we take to be the maximum density of marginal Jeans resolution, ρ J (equation 19), assuming f J = 1/2.

(ii) Density maximum/no overlapping sink: The gas cell is the densest of all neighbouring gas cells or sink particles with overlapping kernel radii (with r gi < max (H g , H i )). For the purposes of this criterion we take sinks to have infinite density, i.e. overlapping with a pre-existing sink always prevents sink formation.

(iii) Increasing density: The gas cell's density is increasing: ∇ • v < 0, according to the same least-squares matrix gradient estimator of ∇v used for reconstructing fluid quantities for the MHD solver.

(iv) Virial/Jeans criterion: The gas cell is gravitationally unstable/self-gravitating at the resolution scale, as determined by a local Jeans analysis including contributions from thermal pressure, magnetic fields, and velocity dispersion (Federrath et al. 2010;Hopkins, Narayanan & Murray 2013a). We evaluate a local virial parameter for the gas cell:

where x = ( m/ρ) 1/3 is the local cell length, and • denotes the Frobenius norm. We permit sink formation only if α g < 2. It is easy to turbulent and numerical viscosity, magnetic torques, gravitational torques, etc., and cannot generally be fit by a single choice of α.

verify that this reduces to the usual requirement that the cell is Jeansunstable at the resolution limit where kinetic energy is negligible, and recovers the Hopkins et al. (2013a) kinematic virial criterion when the ∇v term dominates (e.g. preventing sink formation in Toomre-stable flows stabilized by shear near a star).

(v) Tidal criterion: The tidal tensor T at the position of the gas cell is fully compressive (possesses three negative eigenvalues). Note that the linearization of the gravitational field about a point x i is g (x)g (x i ) ≈ T • (xx i ), so if T has three negative eigenvalues then the gravitational force seen in the frame comoving with a ballistic trajectory starting at x i will always be directed back towards the origin. This is similar in motivation to the potential minimum requirement in Federrath et al. (2010) and the Hill sphere criterion in Hubber et al. (2013), intended to pick out actual centres of collapse from the shape of the gravitational landscape. Unlike a potential minimum criterion, the tidal criterion respects the equivalence principle, i.e. it is invariant to the transformation g → g + g for a spatially uniform g , which should not physically affect the internal dynamics of the simulation in any way, but would displace the location of a potential minimum. However, it is less strict than a potential minimum criterion, e.g. it is satisfied at every point in a uniform sphere (in which T is constant and negative-definite), whereas the potential minimum criterion is satisfied at one point, or none if the external field g exceeds the internal field.

(vi) Can collapse before accretion: The local gas free-fall time t ff = 3π 32Gρg is shorter than both the time-scale for approaching a sink particle and the orbital time-scales around that sink particle, as estimated by evaluating equations ( 6) and ( 7) for gas cells.

When a gas cell is converted to a sink particle, it is removed from the simulation domain, and the volume it occupied is reassigned to surrounding cells when they re-compute their volume partitions.

Accretion criteria

Gas cells are accreted by a sink particle if they satisfy the following criteria:

(i) Sink radius: A gas cell is only eligible for accretion if its centre approaches within sink radius R sink . We take R sink to be the greater of the sink particle softening radius S or volume-equivalent radius of a gas cell at the density of marginal Jeans resolution ρ J (equation 19, assuming f J = 1/2):

where c s denotes the isothermal sound speed at sink formation (i.e. it is set at formation, and kept constant thereafter).

(ii) Boundedness criterion: The gas cell satisfies

where u g is the specific internal energy of the gas, v A, g is its Alfvén speed, and (r gs ) is the softened gravitational potential of the sink, evaluated at the separation between the gas and sink r gs . This checks that the sink-gas system is gravitationally bound, and could not escape to infinity in isolation.

(iii) Angular momentum criterion: The gas cell possesses less angular momentum than a circular Keplerian orbit around the sink at r gs (Bate et al. 1995):

(24)

In the limit of ballistic flow, this ensures that the orbit of the gas cell lies within R sink (so we do not capture e.g. a gas cell that only makes a single close passage but then interacts outside R sink and escapes).

(iv) Size/density criterion: The volume of the gas cell is less than the volume within R sink :

This has the effect of ensuring that only gas having spatial resolution on the scale of R sink can be accreted, which may be necessary for the other criteria to be reliable predictors of the gas's dynamics (and whether it is legitimately being accreted). In any true resolved accretion flow, gas will pile up around the sink until this criterion is eventually satisfied. It is analogous to maintaining the maximum refinement level in the vicinity of a sink in an AMR simulation (Krumholz et al. 2004).

It is possible for a gas cell to satisfy all of these criteria for more than one sink. In such instances, the gas is accreted by the sink s with which it has shortest mutual dynamical time t dyn = -1 = r 3 gs G(mg +ms ) . The quantization of resolved accretion into parcels of mass m has certain important limitations. Clearly, a sufficiently slow accretion flow with Ṁ m/t, where t is some time-scale of interest, cannot be captured. In the limit of a ballistic, Bondi-like flow, we can take t = t dyn (< R) = R 3 /GM at some radius of interest R. Then, assuming the physical accretion flow has a certain Ṁ , the most optimistic radius down to which the flow can be resolved is

where we insert typical values for m, Ṁ , and M . Hence, the accretion flow becomes less well-resolved for smaller accretion rates and greater stellar masses. This may impose some numerical bias towards higher accretion rates in the accretion histories of sink particles, and underestimate more extended periods Bondi-Hoyle accretion from low-density gas. However, the effect does converge to the correct solution with sufficient mass resolution.

Accretion from reservoir on to protostar

To model the continuous accretion of the protostar for the purposes of modelling protostellar evolution and feedback, we use a simple prescription:

where Ṁ ,s is rate at which mass arrives at the protostar, f w is the fraction of gas mass transferred into the protostellar outflows instead (Section 4.2), and t acc is the accretion time-scale. Both f w and t acc are variable, prescription-dependent quantities (we discuss f w further in Section 4.2 and in Paper 2), but by default we take t acc to be the mean time interval between the arrival of gas cells of mass m, assuming the accretion rate is c 3 s /G:

which is dimensionally the same as the free-fall time at the maximum Jeans-resolved density, t J ∼ (Gρ J ) -1/2 . This feeds the protostar with an exponentially smoothed version of the discrete resolved accretion rate, with a 1/e-folding time equal to t acc . For the smallest plausible continuous accretion rate in the initial core collapse, Ṁ ∼ c 3 s /G (Shu 1977), our choice of t acc is simply the mean time interval between the accretion of mass quanta m, which guarantees that it limits unphysical discreteness noise without 'oversmoothing' accretion.

Note that the prescription in equation ( 32) is not meant to model the physical accretion rate at the protostellar surface in detail, and is merely a numerical scheme to obtain a continuous version of the resolved accretion rate with a smoothing time-scale adapted to the mass resolution. If the accretion flow is a direct radial infall (e.g. Bondi accretion) then the relevant physical accretion timescale is the free-fall time (generally shorter than t acc ). In the regime where the gas hits an angular momentum barrier before reaching 9 The raw accreted angular momentum of a sink particle is typically of order √ GM s R sink , which depends on the numerical parameter R sink , and is typically orders of magnitude greater than the angular momentum of a star (Hubber et al. 2013). To determine the actual stellar angular momentum evolution one must model unresolved AM transfer processes. 10 One possible solution for Lagrangian MHD codes (not explored here) would be to introduce a numerical resistivity η sink that interpolates between ∼0 when r R sink and η sink ∼ √ GM R sink when r ∼ R sink , which would diffuse flux away from the star as mass is carried into the sink, modelling the physical non-ideal processes that occur near protostars.

the protostar, accretion will generally take many orbits, and might be better described by e.g. a Shakura & Sunyaev (1973) α-disc type model (in which the dimensionless parameter α encodes the net effect of gravitational torques, magnetic fields, outflows, and viscosity upon angular momentum transport). In principle, our continuous accretion rate estimator could be fed into a physical model to obtain a more realistic estimate of the rate at which mass arrives at the protostar. However, protostellar accretion on sub-10 au scales is subject to a wide variety of poorly understood complex microphysics (e.g. making the specific choice of α an open problem), so we do not attempt to model such processes here.

Stellar evolution

In simulations with feedback, it is necessary to model the evolution of the protostar or star in the sink particle to inform the emergent luminosity, spectral energy distribution (SED), mass-loss rate, and wind/outflow velocity. We model the star or protostar according to a one-zone sub-grid model whose sole input is the present protostellar mass and the accretion rate (equation 32), originally following Nakano et al. (2000) and based upon the particular implementation of Offner et al. (2009). The model evolves the protostellar radius R explicitly, and is calibrated to recover the results of detailed numerical simulations of individual protostellar evolution. This model has been used in many subsequent works by different groups with different codes (e.g. Myers et al. 2014;Federrath, Krumholz & Hopkins 2017;Murray et al. 2018), so we describe it only briefly and refer the reader to Offner et al. (2009) for full details. The evolution is separated into distinct phases:

(i) Pre-collapse: If M < 0.01 M then the protostar is presumed to be a ≈ 4 au first Larson core that has yet to undergo the second collapse phase (Larson 1969;Masunaga, Miyama & Inutsuka 1998).

(ii) No burning: Once M > 0.01 M the core undergoes the second collapse to protostellar density, but deuterium has yet to ignite.

(iii) Deuterium burning at fixed core temperature: D burning has started, fixing the core of the protostar at ≈ 1.5 × 10 6 K.

(iv) Core burning at variable core temperature: The core temperature has begun to rise and D is convected to the core on short time-scales (burning it roughly as rapidly as it arrives at the protostar).

(v) Shell deuterium burning: If D is still arriving rapidly enough, the protostar swells and forms an outer convective zone where the D ignites.

(vi) Main sequence: The star has reached a central core temperature sufficient to ignite H.

At each time-step, the state of the protostar is updated based upon the present mass, accretion rate, and evolutionary phase, and dictates the evolution of the stellar radius R and the emergent luminosity L (which includes terms from accretion, Kelvin-Helmholtz contraction, D burning, and H burning, as given in Offner et al. 2009). We use the Tout et al. (1996) fits for the mass-dependent zero-age main sequence luminosity L MS and radius R MS , and neglect the effects of stellar evolution beyond the main sequence (apart from modelling a Wolf -Rayet phase for winds, Section 4.3, and an eventual SN for >8 M stars, Section 4.4). For the purposes of modelling SNe, stars >8 M have a mass-dependent stellar lifetime: We plot the main-sequence luminosity L MS , radius R MS , and resulting effective temperature T eff from Tout et al. (1996), the stellar lifetime t per equation ( 34), the flux of H-ionizing > 13.6eV photons Q HI (assuming a black-body spectrum of temperature T eff , Section 4.5), the wind mass-loss rates Ṁwind and Ṁwind,WR for main-sequence and Wolf-Rayet stars, and the wind velocity v wind (with wind quantities assuming solar metallicity, see Section 4.3).

which interpolates between ∼ 10 Gyr for solar-type stars and ∼ 40 Myr for 8 M stars, and asymptotes to ∼ 3.7 Myr for the most massive ( 100 M ) stars. In Fig. 5, we plot L MS , R MS , t , and various other useful derived quantities for stellar feedback (Section 4) as a function of the zero-age main-sequence mass M ZAMS .

We do not model the end of life of stars less massive than 8 M (i.e. planetary nebulae), but this could be important for calculations that run for much longer than a GMC lifetime (e.g. Section 5.1.2). We also presently neglect the formation of relic compact objects, but this would be a trivial modification to the inputs of the SN/wind module (simply reserving a certain relic mass), given a more-detailed stellar evolution prescription.

Singular isothermal collapse test

We first validate the formation and resolved accretion criteria of sink algorithm in the Shu (1977) singular isothermal sphere problem, the collapse of a core with an initial density profile ρ = Ac 2 s 4πr 2 , where A parametrizes the family of solutions and collapse occurs for A > 2. This problem possesses a single central singularity (to be represented by the sink), and admits a semi-analytic, spherically-symmetric solution for all fluid quantities [from the numerical solution of Shu (1977) equations 11 and12]. This unambiguous reference solution allows it to quickly expose numerical quirks and bugs, whereas testing the sink particle algorithm on e.g. a full turbulent GMC collapse problem is both more expensive and less conclusive because the 'correct' solution (or whether it exists for a given physics setup) is unknown a priori. An insufficiently-strict sink formation prescription (or an overly strict accretion prescription) can result in the formation of multiple spurious sinks when there should be a single singularity. Errors in momentum conservation or gravity can cause the sink to drift from the centre of collapse, causing subsequent gas to arrive off-centre and form spurious discs or sinks. Passing this test does not prove that a sink algorithm is valid for all problems, but failing this test is a strong indicator that the algorithm is flawed.

For reliable test results, the initial conditions should represent the analytic initial density field with equal-mass elements as we use in our simulations, but this is non-trivial. For MFM, we initialize 125 000 equal-mass gas cells on a uniform radial grid (producing the desired r -2 density profile) with random initial angular positions, and relax the resulting Poisson sampling noise in the IC to a glass by reversing the sign of gravity and allowing cells to slide around on their respective initial radial shells, with an artificial drag force to damp out the motions towards equilibrium. We then rescale to survey various values of A. Exactly one sink forms in each test, and we plot its mass accretion rate in Fig. 6 for A values ranging from 3 to 1000. Agreement with the semi-analytic solution is excellent across the entire range of A surveyed. We also verify Galilean invariance by running a version of the A = 4 setup with a velocity boost of 100c s : even at this extreme bulk Mach number, the solution is preserved owing to the machine-precision Galilean invariance of the hydro, gravity, and sink particle algorithms. The error bars in Fig. 6 plot the ±σ variations of our continuous accretion rate estimator (equation 32), showing that its average error is at most a factor of ∼2 for the lowest A values and accretion rates, and generally much less for higher accretion rates.

Effect of sink prescriptions

Because we wish to use the properties of sink particles to predict the IMF that emerges from a given set of physics, it is important to ensure that the results of STARFORGE simulations are insensitive to the specific parameter choices made in our sink algorithm for e.g. the density threshold and sink radius, and ideally have some robustness to the specific choice of sink formation and accretion Shu (1977) singular isothermal collapse problem, in units of c 3 s /G, as a function of the instability parameter A, such that the initial density profile is ρ = Ac 2 s 4πGr 2 and A = 2 is the threshold of stability. To the analytic solution (dashed) we compare results simulated with our default hydro, gravity, and sink particle algorithms simulated with the ICs at rest (squares) and with the initial gas moving at Mach 100 to verify Galilean invariance (star). Error bars indicate the ±σ quantiles of the accretion rate estimator used to feed the sub-grid protostar (equation 32) -variance is driven by the discretization of resolved accretion into chunks of mass m. Agreement with the analytic solution is excellent, and in all instances exactly one sink particle is formed (replacing the central singularity). criteria as well. For this we re-run the M2e4 GMC set-up in Paper 0, a 2 × 10 4 M initially spherical GMC of radius 10 pc at 10 -3 M resolution with numerous variations from the prescription described in this section, listed in full in Appendix A. By including only minimal physics (isothermal MHD and gravity), the incremental effects of sink numerics are expected to be more pronounced than in a more complex setup with realistic thermodynamics and feedback. In this sense this test could be considered a worst-case assessment of the sensitivity of STARFORGE results to sink prescriptions.

The reference numerical parameters in this setup are m = 10 -3 M , R sink = S = 18 au, and ρ J = 2.6 × 10 -14 g cm -3 , and some tests vary these quantities (Appendix A). We plot the results of this sink parameter study in Fig. 7: the SFE, number of sinks N * , mass-weighted median sink mass M 50 , median sink mass M med , and maximum sink mass M max . Our SFE and IMF results are remarkably robust to wide variations in the sink particle prescription and parameters, including ρ th over a factor of 10 6 , and R sink and S over a factor of 200. The SFE is particularly robust, with very good agreement across all tests (the one outlier is consistent with a simple delay in accretion). The only setups that produced markedly different results in the IMF were ones with obvious flaws, such as ignoring the density maximum criterion (making the choice of which gas cell to turn into a sink generally non-unique), making R sink much smaller than S (making the accretion criteria unreasonably difficult to satisfy, because gravity is unresolved at the sink radius), and increasing R sink and S by a factor of >100 (a gross mismatch with the simulation resolution). Moreover our results do not hinge on any one particular ingredient or assumption -neglecting each formation criterion in our prescription in turn had small effects (apart from the density maximum criterion). Stripping down our accretion criteria to simpler versions also made a negligible difference. Hence, in practice there is a fair amount of redundancy between the different elements of our prescription.

We conclude from this experiment that the results of STARFORGE simulations are unlikely to have any strong dependence upon the details of our sink implementation, at least within the space of Bate et al. (1995)-like algorithms we have explored. Hence, to our knowledge, our sink implementation lacks parameter freedom for 'fine tuning' to ensure a particular desired result -rather, simulation results are mainly sensitive to physical processes as desired. We generally recommend that ρ th and R sink be matched to the nominal density and spatial resolution limits of the simulations (e.g. equations 19 and 20), but e.g. the exact numerical prefactors we have adopted for these quantities are not important, within reasonable limits. Our results hint that simpler prescriptions can perform just as well as ours, but we err on the side of redundancy because the cost of evaluating the various sink formation and accretion criteria is small, and no criterion appears to be unreasonably strict (or else the sink algorithm would allow the simulation to crash due to a runaway gas pile-up).

T H E R M O DY NA M I C S

Although often idealized as such, GMCs are not truly isothermal, and many potentially important effects in star formation require an explicit treatment of the thermal structure of the ISM, such as the dynamics of fragmentation (Bate et al. 2003;Larson 2005) and the evolution of bubbles driven by wind, radiative, and SN feedback. We evolve the internal energy of the gas according to the MHD equations explicitly (HR16), accounting for all gravitational and MHD work terms with heating and cooling. We explicitly evolve the species Z, He, C, N,O, Ne, Mg, Si, S, Ca, and Fe, and by default assume initial solar abundances (Z, He, C, N,O, Ne, Mg, Si, S, Ca, Fe) = (0.02, 0.28, 3.26 × 10 -3 , 1.32 × 10 -3 , 8.65 × 10 -3 , 2.22 × 10 -3 , 9.31 × 10 -4 , 1.08 × 10 -3 , 6.44 × 10 -4 , 1.01 × 10 -4 , 1.73 × 10 -3 ), and re-scale abundances appropriately to the desired initial metallicity (but this can be freely varied).

We use a gas equation of state (EOS) with a variable adiabatic index γ , to account for variations in the equilibrium mixture of paraand ortho-hydrogen and the collisional dissociation of H 2 above ∼ 2000 K (Vaidya et al. 2015). However, we do not roll the heat of ionization into the EOS as in Vaidya et al. (2015), because this is handled separately by our cooling/chemistry solver. We fit to the values of γ given in Vaidya et al. (2015) (neglecting the feature corresponding to ionization) as a function of internal energy:

where

is a sigmoid function, δ k = (-0.38,0.22,-0.068,-0.42,0.65), a k = (5. 95,6.18, 10.26,7.71,98.87), b k = (9. 25, 9.89,10.24,11.13, 14.28), and u is the specific internal energy in cm 2 s -2 . We operator-split the adiabatic MHD evolution with a standard implicit cooling algorithm, which solves for equilibrium internal energy, temperature, net cooling/heating rate, mean molecular weight, and ionization state of the gas (treating the adiabatic heating rate from the MHD solver as an additional heating term). Our treatment of cooling and heating terms largely follows the FIRE-2 simulations (described fully in Hopkins et al. 2018b Appendix B), in accounting for free-free, photoionization/recombination, Compton, photoelectric, metal-line, molecular, fine-structure, dust collisional, Figure 7. Effect of variations in sink particle formation and accretion prescriptions and properties upon the results of a simulation of a 2 × 10 4 M GMC of radius 10 pc, including just gravity and isothermal MHD (i.e. re-simulating M2e4 from Guszejnov et al. (2020b) at the same 10 -3 M resolution). We plot the SFE (top left), number of sinks N (top centre), and the mass-weighted median ( M 50 ), number-weighted median ( M med ), mean ( M mean ), and maximum ( M max ) mass statistics of the stellar mass function as a function of time from the beginning of SF. We highlight the most 'extreme' variations: neglecting the density maximum and virial sink formation criteria in turn, reducing R sink by a factor of 1/4 (to ∼ 5 au) without reducing the softening, reducing the minimum density for sink formation ρ th by a factor of 10 -3 , and increasing R sink and the stellar softening S by a factor of 100 (to 1800 au). A variety of other overlapping sink parameter/prescription variations are plotted in grey, listed in full in Appendix A. Our predictions are fairly insensitive to most of these variations, provided they are within reasonable physical limits. and cosmic ray heating and cooling processes. This cooling module has had various evolutionary updates since Hopkins et al. (2018b) that are not important for our results here [e.g. updating to the Faucher-Giguère (2020) UV background, which is similar to the previously-used Faucher-Giguère et al. ( 2009) UVB at z = 0], but will be detailed in full in an upcoming paper (Hopkins et al. 2021, in preparation). Note that our treatment of hot ( > 10 6 K) gas considers the dominant radiative cooling mechanisms (i.e. freefree emission and metal lines) as described in Hopkins et al. (2018b), but neglects heat conduction by thermal electrons by default, which may moderate expansion of wind and SN bubbles.

Simulations without explicit RHD

If an explicit RHD solver is not enabled, we either make the minimal assumption that T dust is constant, or use a simple, inexpensive RT approximation similar to Guszejnov, Krumholz & Hopkins (2016) and Federrath et al. (2017). This approximation uses the LEBRON radiative transfer algorithm (Hopkins et al. 2018b) to estimate the IR radiation energy density from local sources at the position of a gas cell in the gravity tree pass, summing over contributions from all stars:

where L s are the respective bolometric luminosities of the sink particles, r gs are the distances from the gas cell to the sinks. We then solve for T dust assuming local equilibrium between absorption and emission according to a β = 1 opacity law [i.e. κ(ν) ∝ ν, e.g. Draine 2006]. T dust is then the solution to the quintic equation

where the index k runs over three radiation field components with respective energy densities and effective temperatures: the above component from local sources, which is assumed to have T rad, IR = T dust, IR , the CMB with e CMB = 0.262 (1 + z) 4 eV cm -3 and T rad, CMB = 2.73K(1 + z), and the dust-reprocessed component of the interstellar radiation field (ISRF) with e ISRF = 0.31eV cm -3 and T rad,ISRF = max 20K, T rad,CMB , with fiducial values appropriate for solar neighborhood conditions, and adjustable depending upon the simulated environment. Note that this differs slightly from Guszejnov et al. (2016) and Federrath et al. (2017), who adopted the optically thick grey-opacity radiative transfer solution in the static diffusion limit T ∝ e 1/4 r , as was found in Offner et al. (2009).

Tests

To test the code's ability to capture the transition from isothermal to adiabatic behaviours as the ISM gets optically thick to cooling radiation (important e.g. for the opacity limit for fragmentation), we simulate collapse of a 1 M , uniform-density Jeans-unstable core with both methods described here: explicit RHD with the M1 solver and the simpler TreeCol-based prescription (equation 40). We initialize the cloud with a mass resolution of 10 -5 M (sufficient to marginally resolve the first Larson core, Bate et al. 2003), arranging the cells in a uniform-density glass configuration with density 5.3 × 10 -18 g cm -3 , and assume solar metallicity with a dust-to-gas ratio of 0.01, as has been simulated by many other RHD studies (Larson 1969;Masunaga et al. 1998;Vaytet & Haugbølle 2017). In Fig. 8, we plot the thermal evolution of the centre of the core as a function of its density, showing good agreement with previous calculations: in all instances, the transition from isothermal (T ∼ const.) to adiabatic (T ∝ ρ 2/5 ) evolution occurs at ∼ 10 -13 g cm -3 . Small Local injection adds momentum, energy, and either mass or photons to pre-existing interacting gas cells (coloured domains with circles representing the mesh-generating points) around the sink particle (red star), in proportion to the solid angle subtended by each cell according to the 'effective' faces constructed around each source from the neighbouring mesh-generating points (thick black lines). We use a weighting scheme that ensures statistical isotropy and exact momentum and energy conservation, and in the case of photons accurately accounts for unresolved extinction between the star and cell centre (Section 4.5).

residual differences between the four calculations at the level seen here are expected, due to varied assumptions about radiation initial and boundary conditions, dust properties and opacities, molecular equation of state, etc.

Local injection

When coupling feedback in the form of mass, momentum, and energy from stars, a natural approach is analogous to the manner in which the MHD equations are solved: simply determine the fluxes of those respective quantities at the faces of surrounding gas cells, or more generally distribute those quantities in some weighted fashion in the local hydrodynamic stencil -we refer to this technique as 'local injection' (illustrated for a meshless/unstructured mesh code in Fig. 9). This is what is done in virtually all grid codes and is the technique typically used in GIZMO simulations for local radiative feedback, SNe, and stellar winds (Hopkins et al. 2018a(Hopkins et al. , b, 2020a)). We adopt this algorithm for feedback coupling (described in full in Hopkins et al. 2018a) where appropriate: for stellar winds when the free-expansion radius is unresolved (Section 4.3) and 2). We create new gas cells of mass M w (blue circles) at the sink radius (or at a fraction of the local inter-cell spacing, whichever is smaller), at antipodal positions and velocities to ensure conservation of centre of mass and momentum. The angle θ from the sink angular momentum vector J sink can be controlled to allow arbitrarilycollimated injection. M w can be set smaller than m (the 'normal' mass resolution of ambient gas) to improve resolution in diffuse feedback-driven cavities. Note that this depicts the launching of a pair of gas cells (as for jets and winds), but for SN we launch 24 at once in an angular grid (Section 4.4), designed to ensure exact conservation and isotropy.

for photon injection (Section 4.5). This method ensures machineprecision conservation of mass, momentum, and energy, and ensures that feedback is injected isotropically (or according to the desired angular weighting) even when the local spatial arrangement of cells is anisotropic (unlike e.g. a simple kernel weighting).

Cell spawning

Local injection with Lagrangian methods can run into a major challenge: the resolution is concentrated where gas dens, but feedback structures such as jet cavities, SN remnants, and wind bubbles can be very diffuse. Moreover, if feedback is driving mass away then this problem grows worse over time. And if the inter-cell spacing becomes sufficiently large, it may cease to be a good approximation to instantaneously dump mass, momentum, and energy, because the gas has finite traveltime. Moreover, if we restrict injection to the nearest neighbour cells we cannot inject outflows more collimated than the solid angle subtended by a neighbouring cell (which can be large). So where local injection is not feasible or appropriate, we instead create new gas cells, in a procedure we refer to as 'cell spawning' (Fig. 10). A similar technique has been used previously in SPH simulations of stellar winds (Dale & Bonnell 2008) and protostellar outflows (Rohde et al. 2019), and recently in GIZMO to simulate AGN jets (Torrey et al. 2020). We adapt it here to simulate protostellar jets, stellar winds (when resolvable), and SN ejecta.

Cell spawning can be viewed as the inverse of the gas cell deletion operation that occurs during sink accretion: a new cell is created at a certain position with a certain mass, velocity, and internal energy, and the volume partition in its vicinity is re-computed to accommodate it in the next density iteration. We take the distance between the centreof-mass mesh-generating point of the spanned cell and the sink to be

where x s is the average inter-cell spacing in the vicinity of the sink. We prescribe the initial radial direction and velocity according to the desired angular pattern of the feedback mechanism being realized (Sections 4.2-4.4). We assign purely radial initial velocities, but non-radial velocities could potentially be used to model angular momentum transport (Federrath et al. 2010). Spawned cells are assigned an initial temperature of 10 4 K, and a very small, random initial magnetic field scaled such that the initial plasma β is 10 6 . Spawning is allowed to occur when the sink's internal respective feedback reservoir (Fig. 4) contains enough mass to produce at least N spawn × M w cells, where the number of cells spawned at a time N spawn and spawned cell mass resolution M w are specified for the respective feedback channel. Note that cells do not necessarily have to be spawned with a mass resolution equal to the nominal average 'ambient' gas cell mass m : rather M w can be chosen to be smaller to achieve better time resolution of feedback and to improve spatial resolution within diffuse feedback cavities. For jets and winds, we have generally found the choice M w = 0.1 m to be a good compromise between computational cost and resolution. For SNe, we simply take M w = m. With these choices, we note the spatial resolution in diffuse feedback bubbles will generally be fairly coarse (∼ 1 pc) for a typical m ≈ 10 -3 M , possibly making it challenging to resolve channels and leakage of hot gas (e.g. Rogers & Pittard 2013).

Care must be taken when handling MHD interactions between cells of greatly differing masses and sizes in MFM, particularly for new cells which change the local volume partition and can perturb ∇ • B (interacting with out ∇ • B cleaning algorithm). Spawned cells with M w < m/2 default to a lower-order but morerobust reconstruction in the Riemann problem. We also limit the magnitude of the oriented effective face area A gg between cells (all interacting cells, wind or not) to the lesser of their geometric areas, A max = min πh 2 g , πh 2 g . A spawned cell g is merged into a normal cell g if they are hydrodynamically interacting neighbours, they are moving towards each other, and |v g -v g | < min c s,g , c s,g . We have found that this is a good indicator that the spawned cell has merged with the surrounding ISM, and hence its additional resolution is no longer required.

Lastly, to ensure physical and convergent results, it is important for any feedback algorithm to ensure conservation of momentum and centre of mass. We achieve this by always spawning cells in multiples of 2, such that each cell has an antipodal counterpart in the opposite direction, giving machine-precision conservation.

Physics prescription

Collimated, bipolar protostellar outflows (jets) have special importance as a feedback mechanism. Many works have demonstrated their potential importance both in setting the IMF and SFE, because they are the only channel can be both prompt and omnipresent, immedi-ately regulating the growth of individual stars without requiring e.g. massive stars with dynamically relevant wind or radiative fluxes to be present (Matzner & McKee 1999a;Krumholz et al. 2019). We find that they are likely to be an important ingredient for the IMF in Paper 2, consistent with previous studies (Krumholz, Klein & McKee 2012;Myers et al. 2014;Federrath et al. 2014;Li et al. 2018;Cunningham et al. 2018). The details of the MHD jet launching mechanism (e.g. whether it is better-described by the canonical 'X-wind' or 'D-wind' models, Pudritz & Norman 1983;Shu et al. 1994) are the subject of active research, and depend upon a variety of complex microphysics operating at sub-au scales that are not practical to resolve in our simulations (although the simulations may provide important context for subsequent 'zoom-in' studies of individual stars). Thus, we model jets according to a simple phenomenological prescription following Cunningham et al. (2011), parametrizing the jets' properties in three parameters: the fraction f w of mass accreted by the envelope-discstar system that is diverted to the jet (see Fig. 4), the fraction f K of the Keplerian velocity at the protostellar radius R at which jets are launched, such that

and the collimation angle θ 0 , such that the angular distribution of injected wind momentum is given by (Matzner & McKee 1999b):

where θ is the angle with respect to the angular momentum axis of the sink J s of the sink. This concentrates injection in a narrow cone of angular size ≈θ 0 about the angular momentum axis of the star (see Fig. 11). Our default choices are f w = f K = 0.3 and θ 0 = 0.01, following Cunningham et al. (2011) and other authors. These choices put the momentum loading f w f K roughly in the middle of the observed range (see Cunningham et al. 2011, section 2.4;Federrath et al. 2014, section 3.5;and references therein). The values of f w and f K do matter: in Paper 2, we find that the product f w f K affects the IMF peak. The specific value of θ 0 is likely to be unimportant because in any realistic turbulent accretion scenario J s will generally tend to precess during accretion over an angular region much larger than θ 0 (Rosen & Krumholz 2020), and even without precession jet cavities will expand in the perpendicular direction, opening up an ever-increasing solid angle (Arce & Sargent 2006;Offner et al. 2011). In our tests, we will show that our results are insensitive to variations in θ 0 of at least a factor of 10, which is consistent with prior hydrodynamic outflow simulations carried out by Offner & Arce (2014).

Numerical methods

We always couple jets via cell spawning (Section 4.1.2), waiting until sufficient mass is available in the jet reservoir to spawn two cells of mass M w . The angular direction with respect to the sink angular momentum J s is sampled randomly for the first cell from equation ( 43), and the second cell is pointed in the opposite direction, conserving momentum and centre of mass. Our prescription ignores both the angular momentum and magnetic flux content of the jet material. Although outflows can be the dominant mechanism of angular momentum transport within the disc (i.e. on scales smaller than the sink radius), the material in the disc already had to have very low angular momentum to get to the base of the jet, so it is unlikely that this angular momentum has important effects on larger We plot some examples of the effects of the jet module in simulations in Fig. 11. We generally observe realistic-looking structures in the simulations, with broad bipolar cavities penetrated by a narrow jet, surrounded by a disc (if angular momentum support is important) or a pseudo-disc (of infalling material funnelled by the bipolar cavities). v jet tends to increase as the star accretes (because R varies only weakly in equation 42), so the jet tends to catch up with itself, piling up and cooling in a plume-like region reminiscent of Herbig-Haro objects (e.g. Bally 2016).

Lacking a test problem for our jet module that admits an analytic or universally agreed-upon solution, we resort to heuristic methods to validate it: checking for numerical convergence, checking for robustness to uncertain or arbitrary numerical parameters, and finally comparing with another published solution from a different code.

Resolution tests

Here we consider the effects of numerical resolution upon a simulation of the M2e3 GMC model introduced in Paper 0 with gravity, MHD, the cooling module without explicit RT or protostellar radiation (Section 3), and the jet module enabled. The initial condition is a spherical, uniform-density GMC with mass 2 × 10 3 M , radius 3 pc, and virial parameter 2 (for an initial turbulent Mach number of ∼9). The GMC is surrounded by a diffuse medium with 1/1000 the density filling a 30 pc box, and the initial magnetic field is uniform throughout the box with a strength of 2.3 μG. The normal mass resolution m varies from 0.1 to 10 -4 M , and the jet mass resolution M w varies from 0.01 to 10 -5 M . The sink formation density threshold, sink radius, and sink softening scale are varied according to the mass resolution (section 2.5), scaling ρ th ∝ m 2 over a factor of 10 6 between 3 × 10 -17 and 3 × 10 -11 g cm -3 , and scaling R sink = S ∝ m from 1800 to 1.8 au. Hence, there is no purely-numerical, resolution-related quantity that is held fixed in our resolution study, so the possibility of inferring 'false' convergence is ruled out. 12 In the top left-hand panel of Fig. 12, we plot the evolution of the SFE for the 10 simulations in our resolution study. Star formation tends to start sooner at lower resolution, opposite to what is seen in the collapse of Jeans-mass clumps (Fig. 2), possibly owing to increased numerical dissipation of turbulence at low resolution. The star formation history ceases to appear sensitive to resolution below ≈0.01 M . This corresponds to the resolution criterion we derived in Paper 0, that the sonic mass M sonic ≈ M 0 M -4 ≈ 0.3M be resolved in at least ≈30 gas cells.

We examine the resolution dependence of the stellar mass spectrum at fixed total stellar mass in panels 2-6 of Fig. 12: the 12 Scaling all purely numerical, dimensional sink-related quantities with resolution is important for resolution studies in SF simulations. Otherwise, it is possible that the constant value of e.g. the density threshold or sink radius imprints a characteristic Jeans mass or length, leading one to falsely infer convergence. In near-isothermal problems with Lagrangian codes, this entails scaling ρ th ∝ m 2 and R sink ∝ m. For AMR codes, one must scale ρ th ∝ x -1/2 min and R sink ∝ x min , where x min is the spatial resolution at the finest refinement level (and assume a Truelove et al. (1997) Jeans refinement scheme). AMR resolution studies may also need to scale the based-level grid resolution to resolve turbulent fragmentation (Haugbølle et al. 2018), i.e. fixing the number of refinement levels.

Figure 12.

Resolution study of a simulation of a 2000 M GMC with initial radius 3 pc, with MHD, gravity, cooling physics, and protostellar jets (Section 4.2). We vary the mass resolution m from 10 -4 -0.1 M (darker is finer), scaling sink accretion and softening radii ∝ m from 2-2000 au (equation 22) and scaling the sink density threshold ∝ m -2 (equation 19), so that there is no characteristic numerical scale that could cause false convergence. We plot the resulting star formation history and the mass-weighted median, mean, and median stellar masses M 50 , M mean , and M med , and black dashed lines correspond to the value for a Kroupa (2002) IMF. SFE and IMF statistics cease to scale systematically with resolution past a certain threshold. number of stars N , the mass-weighted median stellar mass M 50 , and the median, mean, and maximum stellar masses. All IMF statistics, whether mass-or number-weighted, eventually cease to change systematically with resolution. As in the isothermal case, the resolution threshold for M 50 and M max to stabilize is ≈0.01 M . Number-weighted statistics require somewhat higher resolution, with ≈10 -3 M being the marginal value for accurately predicting the mean stellar mass and 4 × 10 -4 M for the median. We also plot the effects of resolution upon statistics taken over different stellar mass cuts in Appendix B, finding that the resolution required depends upon the mass cut (consistent with a simple resolution-dependent low-mass incompleteness effect, with incompletness starting below <100 m).

By design, there is no purely numerical dimensional quantity that could imprint a characteristic stellar mass here, so it is likely that the predicted IMF is shaped largely by the physical processes modelled in the simulation, i.e. there may exist a well-defined, physical IMF that emerges from the combined physics of cooling, MHD, gravity, stellar dynamics, and protostellar outflows, and this IMF resembles the observed one. We explore this IMF prediction across a wide parameter space in Paper 2. Assuming that other feedback processes do not demand further resolution requirements, this experiment gives some idea of the mass resolution needed to predict e.g. the mean stellar mass in STARFORGE simulations. Note that some incompleteness in the IMF may persist to higher resolution -to obtain a complete IMF, we may require a resolution of ∼10 -5 M to fully resolve the collapse of clumps at the opacity limit, forming the smallest brown dwarfs (Bate et al. 2003). But brown dwarfs contain only a small fraction of the total number and mass in stars and are not expected to exert significant feedback. Hence, many major questions involving cluster formation, feedback, and the physics underlying the typical mass of stars can be addressed at much lower resolution. We adopt 10 -3 M as our standard resolution for these purposes, but note that the coincidence of our results here with the Paper 0 sonic mass resolution criterion m 0.03M sonic ∼ 0.03M 0 M -4 ≈ 0.01M (T /10 K) 2 /100M pc -2 -1 suggests that this is the more-general convergence criterion. This resolution criterion can be more demanding for e.g. high-surface density GMCs found in the Galactic centre (Oka et al. 2001;Longmore et al. 2012) but not necessarily because such clouds can also be warmer.

We take this opportunity to comment on the computational cost and scaling of these simulations with feedback. Our fiducial resolution (10 -3 M , 2 × 10 6 cells) run cost roughly 10 000 corehours run on the Frontera supercomputer at the Texas Advanced Computing Center equipped with 2.7 GHz Intel Xeon 'Cascade Lake' processors (56 cores per node), and the simulations in the following section at the same resolution all had comparable cost. The largest simulation shown in Fig. 12 (mass resolution 10 -4 M , 2 × 10 7 gas cells) required roughly 100 000 core-hours, running for 17 wall-clock days on 4 Frontera nodes. The largest STARFORGE simulation with jet feedback run so far ( M GMC = 2 × 10 5 M with 2 × 10 8 cells, Paper 2) required roughly 4.8M core-hours in 70 wall-clock days on the Stampede-2 machine at the Texas Advanced Computing Center with 2.1 Ghz Intel Xeon 'Skylake' processors (24 cores per node). Note that these are numbers for simulations without explicit RHD; we anticipate that our full RHD STARFORGE simulations currently in progress will be a factor of ∼5-10 more expensive than their non-RHD Figure 13. Effects of various numerical and physical variations upon the star formation history, kinematics, and IMF evolution of the same 2000 M GMC simulated in Fig. 12, at 10 -3 M resolution. Quantities are as described in Fig. 12, plus the evolution of the 3D mass-weighted RMS Mach number M (bottom left). Reference: Baseline settings, with cooling, MHD, and protostellar jets enabled. No jets: Jet module disabled. No jets; f acc = 0.5: Jet module disabled, and sinks delete half the mass that they accrete. Protostellar heating: Including the approximate protostellar heating prescription described in Section 4.1.2, with physical and artificially large (× 10) luminosities. Isotropic jets: spawning jet cells in random (versus highly collimated, equation 43) directions. θ 0 = X: Changing the jet collimation angle θ 0 (equation 43) from the standard value of 0.01. No MHD: Setting the magnetic field to 0.t acc × 10: Scaling the accretion smoothing time-scale (33) × 10. No jets if <X: Jets are disabled for stars with mass < X. M w = m: Setting the jet cell mass resolution to the nominal mass resolution m, as opposed to the standard m/10. Ang. mom return: Returning accreted angular momentum to surrounding gas as in Hubber et al. (2013). counterparts, mainly due to their more stringent time-step constraints.

Effect of physics variations and numerical details

In Fig. 13, we explore the effects of 13 other variations on this set-up (both numerical and physical) upon the same SFH and IMF statistics, as well as the GMC kinematics. Neglecting jet feedback altogether results in much higher terminal SFE, even if we delete half the accreted mass from the simulation, a simple prescription used by previous works to model jet feedback, deleting either on-the-fly or in post-processing (Padoan, Haugbølle & Nordlund 2012;Federrath & Klessen 2012;Haugbølle et al. 2018). Models that do not explicitly treat feedback also seriously underestimate the level of turbulence in the GMC, and predict a much more top-heavy IMF (e.g. greater values of M 50 ). Haugbølle et al. (2018) found that deleting half the accreted mass gave a good fit to the observed IMF, but simulated a somewhat different GMC set-up, and it can easily be reasoned on dimensional grounds that the accretion efficiency one needs to emulate the effect of jets could generally be problem-dependent. Overall, runs with jet feedback behave dramatically different to the baseline run with feedback: both the total stellar mass formed and average individual stellar masses are roughly an order of magnitude smaller, because the momentum content of the jets disrupts both local protostellar accretion flows and eventually the cloud itself (see Paper 2).

All results are fairly insensitive to variations in the collimation angle θ 0 from 0.01 to 1, the jet mass resolution M w from 0.1 to 1 m, and the accretion smoothing time-scale t acc from 1 to 10 × mc 3 s /G (equation 33). Results are also insensitive to whether we allow jets from stars <0.1 M , whether we include protostellar heating, and whether we return accreted angular momentum to the surrounding gas as in Hubber et al. (2013) (which would influence the stellar angular momenta and jet directions in turn). Results only differ significantly for variations that are ruled out observationally and/or unphysical: neglecting magnetic fields (resulting in a much more bottom-heavy IMF, Guszejnov et al. 2018), assuming jets are emitted isotropically (reducing the typical stellar masses and increasing the overall SFE), artificially increasing the coupled protostellar heating luminosity by a factor of 10 (which made the IMF noticeably more top-heavy, as in Krumholz et al. 2012), and disabling jets for all stars <1 M (which also made the IMF more top-heavy).

Comparison with AMR simulations

To conclude our testing of the jet module, we test its results against a code that differs from ours in all regards except for the physical equations solved and the physical assumptions underlying our feedback models. Our objective here is to verify that the results of simulations with jets are robust to such details.

We have reproduced the set-up described in Cunningham et al. (2018), who simulated low-mass cluster formation in a periodic box of side length 0.65 pc containing 185 M with the ORION2 AMR constrained-transport MHD code with radiative transfer and protostellar jets (Cunningham et al. 2011;Li et al. 2012). We re-run their driven turbulence simulation with an initial mass-toflux ratio μ = 2.17, initially driving turbulence at M ∼ 6.6 for ∼ 1 Myr and then switching on gravity. We did not use exactly the same turbulent driving pattern, so our results should only be compared statistically, hence we ran an ensemble of simulations from Figure 14. Comparison between the results of the driven μ = 2.17 simulation in Cunningham et al. (2018) with MHD and protostellar heating and jets, and three different STARFORGE replications at standard 10 -3 M resolution, using our simple approximate RT treatment of dust heating and our protostellar jet module (modified so that v jet = min √ GM /R , 60 km s -1 , as in Cunningham et al. 2018). Top: Star formation efficiency versus time. Bottom: Stellar mass function at the final SFE of ∼ 8 per cent, comparing the stacked mass function from the three STARFORGE runs with Cunningham et al. (2018), showing the respective median stellar masses as vertical lines. three different initial turbulence realizations. We adopted the same modified prescription for the jet speed as Cunningham et al. (2018), v jet = min √ GM /R , 60 km s -1 , and the two codes' respective protostellar evolution modules setting R both follow Offner et al. (2009). Cunningham et al. (2018) account for protostellar radiation with a flux-limited diffusion solver, while we use the inexpensive tree-based approximation described in Section 3. Our simulations adopt a mass resolution of 10 -3 M and cost roughly 200 core-hours each when run on a single Frontera CLX node.

We plot the resulting star formation history and IMF at 8 per cent SFE in Fig. 14. The specific shape of the SF history does appear to be a function of the initial turbulent driving, but we had two seeds that matched that of Cunningham et al. (2018) for an appreciable fraction of the SF history, both modulo a small time difference due to the different initial turbulent states. Therefore predictions regarding the regulation of SF due to protostellar feedback appear similar for the two codes.

The final IMFs are also in fair agreement: the median stellar mass predicted by both codes is ∼0.15 M (shown as vertical lines). The only statistically significant difference is between the predicted numbers of 0.01-0.03 M brown dwarfs, with STARFORGE runs finding only ∼1/4 as many. This may be due any number of details in the cooling and RT modules that differ (with more higher temperatures or more radiative heating suppressing brown dwarfs, Bate 2009b; Offner et al. 2009), or a mere resolution effect (as we expect some numerical IMF incompleteness at masses ࣠ 30 m).

In summary, the respective implementations of STARFORGE and ORION2 find very similar results for the regulation of SF and the stellar mass range of the IMF in low-mass star cluster formation, despite these two codes' detailed numerical differing in every regard. We scale a set-up similar to this to GMCs as much as >1000× more massive in Paper 2, exploring the broader implications of protostellar feedback in massive GMCs.

Tests

We test this module on the problem of the self-similar expansion of a wind bubble propagating into a uniform medium with an Figure 15. Self-similar expansion of a stellar wind bubble with an adiabatic interior and radiative outer shell, as simulated with our stellar wind module (Section 4.3) for a star with Ṁwind = 10 -5 M and v wind = 3000 km s -1 in a 16 pc box containing 2 × 10 4 M , using different numerical methods (local injection (Section 4.1.1), cell spawning (Section 4.1.2), and a hybrid method that switches between them adaptively (Section 4.3). Top: Comparison of the evolution of the radius of the swept-up shell R shell (corresponding to the number of swept-up gas cells N shell ) to the known similarity solution, R shell = 0.763 L wind ρ 1/5 t 3/5 (Weaver et al. 1977). Bottom: Numerical velocity anisotropy of the bubble (= 0 in the similarity solution), generally 10 per cent except in very early phases when the bubble is not well-resolved (contains just a few cells). The switch between different coupling methods at 0.02 Myr is apparent when the 'Hybrid' curve deviates from the 'Local Injection' curve.

adiabatic interior and radiative exterior with negligible exterior pressure, which has the analytic solution R shell = 0.763 L wind ρ 1/5 t 3/5 (Weaver et al. 1977). We place a star with Ṁwind = 10 -5 M and v wind = 3000 km s -1 in a 16 pc box containing 2 × 10 4 M , with m = 0.01 M and M w = 10 -3 M , with initial temperature 10K. In Fig. 15, we plot the bubble expansion and find that agreement with the similarity solution is good whether we use pure local injection, pure cell spawning, or our hybrid method. 13 We also examine the velocity anisotropy σ v, max /σ v, min -1 where σ v, max and σ v, min are the maximum and minimum gas velocity dispersions along the principal axes of the gas momentum distribution, which is 0 in the exact solution. This is typically 10 per cent, except in the very early phase of the run with pure cell spawning (because the free expansion radius is not yet well resolved, so shot noise from individual injection steps is still apparent). For the hybrid method, the transition between methods occurs smoothly, with no clear spurious numerical artefacts. 13 We have also found negligible differences between the different wind methods in full star cluster formation simulations including winds as the only feedback mechanism.

Radiation

Following Hopkins et al. (2020a), STARFORGE simulations with the radiative transfer module enabled follow the emission, transport, and absorption of photons in five different bands in wavelength λ:

(i) Hydrogen ionizing (λ < 912 Å): Ionizing photons emitted by stars and responsible for the dynamics of H II regions, which are widely theorized to be the most important feedback effect from massive stars in typical Galactic conditions on global GMC scales (McKee, van Buren & Lazareff 1984;Krumholz & Matzner 2009;Dale, Ercolano & Bonnell 2012;Geen, Soler & Hennebelle 2017;Kim et al. 2018;Grudić et al. 2019;Olivier et al. 2021).

(ii) Far-UV/photoelectric (912 Å< λ < 1550 Å): Responsible for heating the ISM via the photoelectric effect on dust grains, and likely an important component of the thermal balance of the cold and warm neutral media in the outer parts of galaxies (Wolfire et al. 1995;Ostriker, McKee & Leroy 2010).

(iii) Near-UV (1550 Å < λ < 3600 Å): Contains most of the photon energy and momentum emitted by a young stellar population, and hence is the most important term for direct stellar radiation pressure (Fall et al. 2010;Murray et al. 2010;Raskutti et al. 2016;Kim et al. 2018;Hopkins & Grudić 2019).

(iv) Optical/near-IR (3600 Å < λ < 3 μm): Contains most of the light from old stellar populations, and carries a non-negligible fraction of photon momentum that can potentially couple on larger scales in a GMC due to reduced dust opacity compared to NUV.

(v) Mid/far-IR (mainly λ > 3 μm): Radiation absorbed and reradiated by dust, which is the primary cooling mechanism in the densest gas and can dominate the radiation pressure near massive protostars or in ULIRGs (Krumholz et al. 2009;Kuiper et al. 2011;Davis et al. 2014;Rosen et al. 2016;Tsang & Milosavljević 2018). Thi s is treated specially, as a component of the radiation field having a blackbody SED with a local effective temperature T rad , which is evolved self-consistently.

We can optionally further 'fine grain' these bands into narrower bins: for example, splitting ionizing and FUV photons into photoelectric, Lyman-Werner, H ionizing, He-ionizing, He-secondaryionizing, soft X-ray, etc., as described in Hopkins et al. (2020a), but this generally produces second-order effects. Note that, unlike most RHD SF simulations, we independently and explicitly evolve the dust temperature T dust , radiation temperature T rad (of the IR band),foot_6 and gas temperature T gas .

All sinks/stars are treated as potential sources for all bands above. In our default simulations, we calculate the emitted flux in each band by treating each sink/star as a blackbody with effective temperature T eff, ≈ 5780 K (L /L ) 1/4 (R /R ) -1/2 with L and R given by our stellar evolution module (Section 2.5.5), integrated over the relevant wavelengths. We ignore 'primary' gas emission at other wavelengths, as this is generally negligible in the problems of interest. Secondary gas/dust (re)-emission is treated as follows: recombination emission from absorbed ionizing photons are re-emitted into the optical/near-IR (OIR) band; the absorbed radiation energy in other bands (where the opacity is dust-dominated) is re-emitted by dust in the IR band with the evolved dust temperature T dust .

The limitation of this band-integrated treatment of radiation is that line emission and absorption are not followed explicitly. As such, it does not capture molecular line cooling explicitly (which we treat using approximate formulae, see Section 3). It is also not applicable to dust-free conditions where multiply-scattered Ly α photons are dynamically important (e.g. Smith, Bromm & Loeb 2017).

Absorption and scattering cross-sections/opacities and coupling to our thermo-chemistry (gas heating/cooling) routines largely follow Hopkins et al. (2020a). For ionizing bands, this employs standard photoionizing absorption cross sections which scale with the neutral H and neutral or partially ionized fractions for He, and absorbed ionizing photons directly couple to our detailed photoionization heating rates (see Hopkins et al. 2018b; note these are always calculated, our RHD simulations simply include local sources in addition to the UVB and ISRF). The opacities in FUV, NUV, OIR are given by the grey expressions: (κ FUV , κ NUV , κ OIR ) = (0.2 + 2000 f d , 1800 f d , 180 f d ) cm 2 g -1 where f d = f d /0.01 is the local dust-to-gas ratio relative to solar, defined as in Section 3. The FUV intensity directly enters the photo-electric heating rate in our thermochemistry calculation (Hopkins et al. 2018b;appendix B); absorbed radiation in FUV+NUV+OIR + IR bands contributes to determine the dust temperature T dust which interacts via dust-gas collisions (Section 3). For the IR band, we calculate the opacity as a function of ionized fraction, local dust-to-gas ratio, dust temperature, and radiation temperature, as κ IR = κ gas + κ 0 dust f d where κ gas ≈ 0.35 x e cm 2 g -1 from Thompson scattering with x e the free electron fraction, and κ 0 dust (T dust , T rad ) calculated from the tables of Semenov et al. (2003). Specifically, we take the 'standard' model with the 'porous 5-layered sphere' composition in Semenov et al. (2003) and for each dust temperature T dust (which gives a different dust composition) we explicitly calculate the Rosseland-mean opacity for each radiation temperature T rad assuming a blackbody-like spectral shape. We provide detailed fits in Appendix C. We assume a sublimation temperature of T sub dust = 1500 K. Finally, because our RHD methods account for both absorption and scattering we must define the albedo A. For simplicity we assume A ion = 0 (pure absorption) for ionizing bands; A FUV, NUV, OIR = 1/2, i.e. equal absorption and scattering opacities which is roughly appropriate for dust grains in FUV through OIR bands (see e.g. Weingartner & Draine 2001); and for IR we assume the Thompson portion of the opacity is pure-scattering while for the dust albedo we can interpolate reasonably accurately between the short-wavelength (A ∼ 1/2) and long-wavelength (Rayleigh scattering, A → 1) regimes by taking A IR ≈ ( Tr /2 + 1)/( Tr + 1) with Tr ≡ (T rad, IR /725 K) 2 .

Photon momentum (radiation pressure) is always transferred appropriately to the gas + dust when radiation is absorbed or scattered (from any band).

Photon transport

GIZMO employs modular RHD solvers, so in principle we can adopt and compare various methods for photon transport. But in our default explicit-RHD simulations we adopt the first-moment or M1 (Levermore 1984) method, which has the advantages of being computationally efficient (well-adapted to hierarchical time-stepping and multiphysics simulations), manifestly momentum and energy conserving in finite-volume form, able interpolate between opticallythick and optically thin limits, and well-tested in simulations of star cluster formation (Geen et al. 2015;Gavagnin et al. 2017;Geen et al. 2017;He et al. 2019). In particular, for questions involving radiation pressure forces on gas, shadowing in an inhomogeneous medium, and the transition between optically thin-thick regimes, it is (by construction) able to capture phenomena which cannot appear in the 0th-order flux-limited-diffusion (FLD) method (see references above and e.g. Davis et al. 2014;Rosdahl et al. 2015;Zhang & Davis 2017;Kannan et al. 2019). For each band i we explicitly evolve the first two moments of the intensity equation in the usual mixedframe approximation keeping all terms to O(v 2 /c 2 ), with all terms appropriately integrated over the relevant bands (Mihalas & Mihalas 1984;Lowrie, Morel & Hittinger 1999). This gives

where u is the local gas/dust velocity, ėabs ≡ c2 ψ a e r and ėem are the volumetric absorption and emission rates, g r ≡ f ru • (e r I + P r ) with e r and f r the radiation energy and flux densities, ψ a, s ≡ ρ κ a, s / c are the absorption+scattering coefficients, c and c the reduced (RSOL) and true speed of light, and P r ≡ e r D r the radiation pressure tensor (Eddington tensor D r ). 15 These are discretized and inter-cell fluxes are integrated in the same finite-volume (conservative) form as the MHD equations (Hopkins et al. 2020a).

Note that equations ( 49)-( 50) include terms up to formal O(v 2 /c 2 ) such as the u • (e r I + P r ) term differentiating g r and f r or the (ψ aψ s ) u • g r term representing the work done by radiation pressure which are often dropped in SF simulations where typical speeds v c. However, as many authors have pointed out, these terms actually dominate the behaviour in the infinite-optical-depth, tight-coupling or photon-trapped limit, where their actual order scales closer as O(v/v rad, diff ) (with v rad, diff the effective bulk speed of photon diffusion). Without this terms, the RHD equations in the trapped limit will give unphysical solutions (photons will not properly be advected and arbitrarily large radiation-pressure forces can arise). The terms in ( ėi absėi em ) u/c 2 , on the other hand, represent true relativistic beaming, and are negligible for our problems of interest.

In the gas/dust momentum equation we add the terms i (ψ i a + ψ i s ) g i r + ( ėi absėi em ) u/c 2 , the former representing the normal radiation pressure acceleration and the latter accounting for beaming. In the gas/dust total energy equation we add the terms i ( ėi absėi em ) + (ψ i sψ i a ) u • g i r , the former representing the energy of absorption + emission (handled with our thermochemistry as described above) and the latter work terms. The gas temperature is then evolved according to the thermodynamics modules in Section 3, coupled also to the dust temperature by way of dust-gas collisions (equation 36). The dust temperature is set assuming grain absorption-emission equilibrium, giving T 4 dust = ( Q abs / Q em ) e tot r c/(4 σ T ) where e tot r = i e i r , and Q abs, em are the appropriate absorption and emission efficiencies. 16 As noted above, the IR radiation field is treated as a blackbody shape with local effective temperature T rad, IR and total energy integrated over the cell domain E gamma, IR , which evolves as new photons are emitted or when radiation is exchanged between cells of different T rad . In emission, sinks emit with T em rad, IR = T eff and dust has T em rad, IR = T dust ; given a total emitted radiation energy E em in the cell timestep, the effective T rad, IR is then updated to guarantee both radiation energy and photon number conservation, giving T rad, IR (t

The same update scheme is used when cells exchange radiation energy.

Finally, we follow common practice and adopt a reduced speed of light (RSOL) with c < c to enable larger time-steps. In general, c should be larger than the bulk speeds of radiative diffusion or ionizing front expansion to capture the dynamics; Geen et al. (2015) argue this is satisfied for c 30 km s -1 in our problems of interest, and we verify this below.

Tests

First, we remind the reader of the test in Section 3, which demonstrates the accuracy of our IR RHD + thermochemistry models evolving T dust , T rad , and T gas appropriately in the collapse of a Jeansunstable core. To test other aspects of our RHD models, we next repeat the test setup of an O star in a box described in the previous sections (Sections 4.3-4.4) for radiation. First, we examine the ability of the photon injection, transport, and absorption schemes to capture radiation pressure in two regimes: where the photon mean free path is well-resolved ( x λ mfp ) and totally unresolved ( x > λ mfp ). We inject radiation in a single band with opacity scaled so that the cell opacity τ cell = x/λ mfp ranges from 0.015 to 1.5, and the global optical depth τ box = L box /λ mfp ranges from 1.9 to 190. Because the box is always optically-thick and thermal pressure is negligible, in all cases we expect the expanding shell solution to approach the momentum-conserving similarity solution R shell = L 6πρ 0 c 1/4 t 1/2 at late times with radial momentum approaching the total emitted photon momentum L t/c, and the solution should be spherically symmetric.

Fig. 18 shows that the dynamics of a radiation pressure-driven shell are captured accurately: all three solutions eventually approach the 16 For IR-IR band interactions, we assume Q abs ≈ Q em , though this could lead to small differences in T dust .

Figure 18. Radiation pressure-only test with an O star in a homogeneous box, accounting only for radiation pressure feedback with a range of opacities, varying the global and cell optical depths τ box and τ cell and injecting and transporting photons as described in Section 4.5. Top: position of the spherical shell swept up by the radiatively-driven bubble. All solutions approach the analytic similarity solution for momentum-driven bubbles ∝ t 1/2 , with the time lag determined by the absorption time-scale as expected. Middle: Radial gas momentum in units of the total emitted photon momentum Lt/c. The correct momentum is always coupled even when τ cell > 1, by accounting for unresolved local extinction in the injection phase (Hopkins & Grudić 2019). Bottom: Numerical velocity anisotropy, which is < 1 per cent in all but the earliest (i.e. worst-resolved, few-cell) phases of the expansion. similarity solution (with the least optically-thick run having a time delay λ mfp / c ≈ 0.3 Myr, owing to the finite light traveltime before absorption). The correct radial momentum is imparted whether or not λ mfp is well-resolved (due to the face-integrated injection method, detailed in Hopkins & Grudić 2019), and numerical anisotropy falls rapidly below ∼ 1 per cent once the bubble becomes well-resolved.

We repeat this experiment with our ionizing radiation band enabled, without radiation pressure, to test the ability of the code to follow the dynamics of expanding H II regions (and also to determine an appropriate value for c to accomplish this). We compare with the approximate Hosokawa & Inutsuka (2006) solution for the position of the ionization front:

where c i ≈ 11 km s -1 is the isothermal sound speed in ionized gas and the R St is the Strömgren radius:

where ṄLyC is the emission rate of H ionizing photons and α B is the case-B H II recombination coefficient (invoking the 'on-the-spot' approximation, Osterbrock 1989). We fix the ratio ṄLyC /α B so that R St is initially resolved in only two cell lengths, R St = 2( m/ρ 0 ) 1/3 , and survey c values between 5 and 300 km s -1 . We also run a version with c = 30 km s -1 with the gas fixed in place, to compare with the static Strömgren sphere solution and isolate artefacts of the RSOL approximation.

We plot the expansion of the ionization front for the different runs in Fig. 19. The 'No Hydro' frozen solution relaxes to the static Strömgren sphere solution after a time ∼ t St = R St / c as expected, and remains statistically (but not exactly) static thereafter. The solutions with c > 30 km s -1 relax to the Strömgren solution similarly, but then start to expand after a time ∼R St /c i and agree well with the equation (51) solution. But for c = 5 km s -1 , t St is longer than the physical sound crossing time of the bubble, so the bubble expansion is delayed artificially, confirming the finding of Geen et al. (2015) that c ∼ 30 km s -1 is roughly the marginal RSOL value for following the dynamics of H II regions accurately. Numerical anisotropy is again small (a few per cent ) once the bubble actually expands and becomes well resolved.

Local and global galaxy simulations

Stratified and/or shearing-box simulations have been used to simulate the evolution of a patch of the ISM within a galaxy at a resolution that is generally higher than what is attainable in global galaxy simulations (e.g. Hennebelle & Iffrig 2014;Walch et al. 2015;Martizzi et al. 2016;Kim & Ostriker 2017). These can be used to follow the formation and dispersal of GMCs self-consistently. All the algorithms presented here translate directly to this type of setup -the only differences are the initial conditions, boundary conditions, and additional inertial forces (all currently implemented in GIZMO). Resolved SF simulations could also be performed in a galactic context via a 'zoom-in' re-simulation of a GMC, taking a coarsely resolved GMC in a simulated galaxy (e.g. Guszejnov et al. 2020a) and up-sampling it to higher resolution (Rey-Raposo, Dobbs & Duarte-Cabral 2015).

The total stellar masses we form in the largest simulations in Paper 2 ( 10 4 M ) is within an order of magnitude of the total stellar mass in the faintest known dwarf galaxies (Simon 2019;Wheeler et al. 2019), so it may even be possible to perform a global cosmological galactic zoom-in simulation of an ultra-faint dwarf (UFD) with individually-resolved star formation. State-of-theart simulations like Wheeler et al. (2019) simulate UFD formation at a mass resolution of ∼30 M , so much higher resolution would be required, but this could potentially be achieved by using a T97-like refinement criterion to reach the required resolution in the dense gas, while the more diffuse gas not engaged in star formation could be kept at coarser resolution.

C O N C L U S I O N

We have presented and demonstrated the methods of the STAR-FORGE, combining the physics of MHD, gas self-gravity, stellar dynamics, thermodynamics, and all major dynamically important (proto-) stellar feedback mechanisms into a detailed numerical model of star formation. We have shown that the respective techniques for each mechanism give satisfactory results in test problems with known solutions. We also discovered a remarkable degree of robustness in the sink particle prescription (Section 2.5.8), and found good agreement when comparing with results in similar problems from a code that implements the same physics with completely different numerical methods (Section 4.2.3).

We found stable numerical results for the IMF down to a completeness limit of ∼ 0.1M at the modest (by SF simulation standards) mass resolution of ≈10 -3 M (Section 4.2.3, Fig. 12), and in Paper 2 we scale this setup up to GMCs as massive as 2 × 10 5 M , mapping out exploring the effects of protostellar jets upon the IMF at this scale for the first time. In subsequent works, we will present the full results of radiation MHD simulations of those same massive GMC models, with all feedback modules described in this work acting in concert.

We anticipate that STARFORGE will be a useful theoretical laboratory for disentangling the many physical mechanisms at work in GMCs. By starting with a realistic picture and switching different physics on and off in controlled experiments, it can help distil the essential elements of a working theory of star formation. It can also be used to calibrate sub-resolution prescriptions for effects such as stellar kinematics and protostellar feedback for use in lowerresolution star cluster and galaxy formation simulations, increasing the predictive power of such simulations in the densest gas and star clusters, where the details of prescriptions become important (e.g. Hopkins et al. 2013b;Grudić & Hopkins 2019;Li et al. 2020). It should also be useful as interpretive tool for observations, mapping out the effects of different physics upon the relations between observed gas tracer properties and star formation.

An important goal of this project is to reduce the dependence of SF simulations upon sub-grid prescriptions, which must make highly uncertain assumptions about how individual stars form. Our setup helps to accomplish this, but it only peels back one layer. To simulate feedback and the emergence of the IMF, we must make assumptions about how various sub-au physics (accretion, stellar winds, jet launching, protostellar and stellar evolution/death) proceed, and we list various ways in which incorrect assumptions about these processes could affect our results (Section 5.2). Hence, it is crucial to continue to advance our understanding of the processes governing the formation and internal evolution of individual stars and star systems.

upper envelope of the predicted IMF mass scale statistics because accretion is made much easier and IMF incompleteness is introduced by accreting independently-collapsing cores within the sink radius. However, these effects are small.

(iii) Increasing R sink and ×10 and ×100 without rescaling S . These had similar results to tests in which we scaled S as well.

(iv) Decreasing S from 18 to 1.8 au. This had negligible effects. (v) Changing the sink formation density threshold ρ th by a factor of 10 -3 and 10 3 . The ρ th × 10 -3 test also neglected the thermal term in the virial criterion (equation 21), which would otherwise have effectively imposed a density threshold of its own. This has the effect of rescaling the threshold number of Jeans wavelengths per cell at the sink formation threshold from 0.5 to 1.58 and 0.15, respectively [i.e. strongly violating versus satisfying the Truelove et al. (1997) criterion]. All results of the ×10 3 run were difficult to distinguish from the fiducial run, possibly because following collapse far beyond ρ J is unlikely to reveal any new fragments (Section 2.4). The ρ th × 10 -3 run formed a slightly larger number of sinks than the fiducial run, because the other sink formation criteria are not perfect predictors of runaway collapse when looking at gas of modest density, but effects were still quite small.

(vi) Rescaling R sink from 18 to 9 au and 4.5 au, respectively, while keeping the softening radius fixed at 18 au. The R sink = 9 au run had no clear systematic difference from the fiducial run. The 4.5 au run (highlighted in Fig. 7) was the largest outlier of our survey, forming stars with a slight delay with respect to the modal SFE, and producing a noticeably larger (×2) number of sinks, affecting the median and mean sink masses in turn. Mass-weighted statistics such as M 50 and M max were affected more modestly, but were also systematically lower than the modal solution. Note that this is not a particularly reasonable prescription: it forces sink particles to have a volume 1/64 the volume of a gas cell at the resolution limit, resulting in a gross mismatch between R sink , S , and the gas resolution scale at the time of sink formation, making it very difficult to satisfy all accretion criteria.

(vii) A minimal accretion prescription requiring only r < R sink . This had negligible effects.

(viii) A simple accretion prescription requiring boundedness (equation 23) and r < R sink . This had negligible effects.

(ix) Rescaling R sink by a factor of 1/4 while also rescaling S by the same factor and ρ th by a factor of 64 so that the cell spacing at sink formation matched R sink . This was much closer to the reference run than reducing R sink alone.

(x) Running our fiducial sink formation prescription in conjunction with the simpler sink accretion prescription given in Bate et al. (1995) (except neglecting the correction terms to the hydro force). This amounts to neglecting thermal and magnetic pressure in the boundedness calculation (equation 23), and not requiring that gas cells physically fit inside the sink volume. This had negligible effects because there is considerable redundancy between the different checks.

(xi) Ignoring the ∇ • v, virial, density maximum, tidal, and infall sink formation criteria in turn. These all had negligible effects except for neglecting the density maximum and virial criteria, which were at the upper envelope of number of sinks formed.

(xii) Including a version of the Hubber et al. (2013) angular momentum return prescription, exerting a net torque τ = Js tacc upon the surrounding gas to transfer angular momentum from the sink back to the gas (taking t acc to be 500 yr, likely much faster than the actual angular momentum transfer time-scale in a protoplanetary disc). This had negligible effects, but may be more pronounced in problems where protostellar discs are well-resolved and angular momentum support is important.

(xiii) Enforcing the additional criterion of 'collapse along all three axes', a stricter version of the ∇ • v criterion. We check that all three eigenvalues of the symmetric component of ∇v are negative (as opposed to merely their sum ∇ • v), similar to prescriptions used in Federrath et al. (2010) and Gong & Ostriker (2013). This had negligible effects. In Section 4.2, we perform a resolution study of a 2000 M GMC with mass resolution ranging from 0.1 to 10 -4 M with cooling, MHD, and protostellar jet physics enabled, and found that the predicted IMF statistics stabilized at sufficient resolution (Fig. 12). In Figs B1 and B2, we remake the relevant panels from Fig. 12 while cutting stars <0.1 M and <1 M , respectively, to determine the resolution requirements for statistics computed on different mass ranges of the IMF. Cutting at <0.1 M (Fig. B1), a mass resolution of ≈2 × 10 -3 M ap-