DOES GOD PLAY DICE WITH STAR CLUSTERS?

,


INTRODUCTION
Astrophysicists aim to interpret light from galaxies and star clusters in terms of their stellar populations and the physical processes that affect them.More often than not, detailed knowledge of the individual stars is not available, so inferences about the system rely on a model of their stellar demographics.The cornerstone of such models is the stellar initial mass function (IMF): the distribution of stellar birth masses (Salpeter 1955;Hopkins 2018).The standard approach to modeling a galaxy or star cluster is to sample stellar masses from the IMF and evolve the stars forward in time according to a stellar evolution model (e.g.Leitherer et al. 1999;Bruzual & Charlot 2003).
A common model assumption is that stellar masses are sampled independently and identically from an invariant IMF ("stochastic sampling"), i.e., the probability of a star forming with mass m 1 does not depend on the mass of another star m 2 , nor the macroscopic properties of the system the star is forming in.Or, in other words, 10 little star clusters are statistically identical to one big one.This is the simplest approach and is motivated by the observed lack of strong, systematic variations in the IMF in normal star-forming environments (Bastian et al. 2010;Offner et al. 2014).
Stochastic sampling should not necessarily be interpreted physically: the initial gas flows in molecular clouds that produce individual stars cannot be truly random and independent of each other.The formation of †NASA Hubble Fellow Corresponding author: mgrudic@carnegiescience.edu one star can affect another via gravity and feedback in the form of protostellar outflows, winds, radiation, supernovae, and cosmic rays.So stochastic sampling cannot be the ground-truth star formation process, and if star clusters deviate from it significantly then so will any resulting inferences from the integrated light of galaxies (Weidner & Kroupa 2006;Kroupa et al. 2013;Kroupa & Jerabkova 2021).
Of course, all models are wrong at some level, but some are useful.Stochastic sampling from an invariant IMF has enjoyed wide application due to its simplicity and the lack of an alternative model with compelling physical motivation.It has been used to characterize star clusters and their demographics from photometry (Fumagalli et al. 2011;da Silva et al. 2012), where it gives a marginally-better fit to the data than IMF-averaged photometric models without stochasticity (Krumholz et al. 2015) and can enable inferences about the highly-incomplete low-mass cluster population, if accurate (Tang et al. 2023).Stochastic sampling has also seen increasing use in high-resolution galaxy and molecular cloud simulations (Fujii & Portegies Zwart 2015;Hu et al. 2017;Su et al. 2018;Lahén et al. 2020;Smith 2021;Liow et al. 2022).On molecular cloud (≲ 100pc) scales, the granularity of feedback and gravity from individual stars is quite important; the specific timing, sampling, and distribution of massive star formation has a large effect on subsequent star formation (Grudić & Hopkins 2019;Dinnbier & Walch 2020;Lewis et al. 2023).However, lacking a detailed physical picture of star cluster formation, it has been difficult to assess the fidelity of stochastic sampling fully.
Recently, the STARFORGE framework1 (Grudić et al. 2021) has enabled radiation magnetohydrodynamics (RMHD) star formation simulations with significantly enhanced realism and self-consistency.By accounting for key feedback processes and resolving individual star formation, the simulations can follow the entire sequence of star formation from initial cloud collapse, through star formation and accretion, feedback, to the eventual cloud disruption (Grudić et al. 2022).This allows us to predict the final masses of stars and hence the IMF in a consistent framework (Guszejnov et al. 2022).
In this paper we present the results of a statistical ensemble of 100 STARFORGE simulations and measure their combined and individual stellar mass statistics.The statistical power of this suite allows us to compare stochastic sampling models to a physical model of star cluster formation with all relevant processes for the first time.In particular we aim to test whether statistical sampling, independent of cloud mass, provides a good representation of star formation process realized in nature.We will focus on the statistics of massive (≳ 8M ⊙ ) stars, which dominate the light and stellar feedback from young star clusters, and hence are most relevant for photometric modeling and numerical simulations.

SIMULATIONS
We use the GIZMO code (Hopkins 2015) to run an ensemble of 100 STARFORGE RMHD star cluster formation simulations.The simulations use the "full" STAR-FORGE physics setup introduced in Grudić et al. (2022), with numerical methods for sink particles, stellar evolution, and feedback coupling described in full in Grudić et al. (2021).The simulations solve the equations of ideal magnetohydrodynamics (MHD) using the GIZMO Meshless Finite Mass MHD solver (Hopkins & Raives 2016) with constrained-gradient divergence minimization (Hopkins 2016) and a fixed Lagrangian mass resolution of 10 −3 M ⊙ .
Gas cells turn into sink particles if they satisfy all of a list of criteria intended to identify centers of runaway gas collapse.They then accrete other gas cells that have sufficient density, proximity to the sink, low angular momentum, and gravitational boundedness (e.g.Bate et al. 1995).There is not stochasticity in these criteria.Each sink particle hosts a star that accretes the sink's gas reservoir smoothly, evolving in radius and deuterium abundance according to the Offner et al. (2009) protostellar evolution model.Stars on the main sequence may still accrete and move along it, and their bolometric luminosity and radius are given by the Tout et al. (1996) fits.Stars' emergent spectra are approximated as a black-body with and their metallicity-and luminosity-dependent massloss rate and wind velocity are given by the fitting formulae in Grudić et al. (2021Grudić et al. ( , 2022)).Massive stars can progress to a Wolf-Rayet phase with strong winds and an eventual supernova, but this does not occur during the star-forming phase of any of the clouds simulated here, which are disrupted primarily by protostellar outflows and ionizing radiation.
We solve the time-dependent, frequency-integrated radiative transfer (RT) equation with the GIZMO mesh-less frequency-integrated M1 solver (Hopkins & Grudić 2019;Hopkins et al. 2020) with a reduced speed of light of 30 km s −1 , in 5 frequency bands ranging from Lyman continuum to far IR.The different bands couple to matter through photoionization, photodissociation, photoelectric heating, and dust absorption, accounting for photon momentum.We account for radiative cooling and heating from all major molecular, atomic, nebular, and continuum processes using the cooling module shared with the FIRE-3 simulations (Hopkins et al. 2023), except that we perform dust radiative transfer self-consistently with the IR band of the RT solver, evolving the dust, gas, and radiation temperatures independently.

Initial conditions
All 100 simulations start with a uniform-density molecular cloud of mass M = 2000M ⊙ and radius R = 3pc, in a 30pc periodic box filled with an ambient medium 10 3 times less dense.The initial magnetic field is uniform and normalized to give a mass-to-flux ratio of 4.2 times the critical value inside the cloud.The simulations differ only in their initial random velocity fields, which are different statistical realizations of a Gaussian random field with a ∝ k −2 power spectrum, all normalized so that the initial virial parameter is α turb = 5σ 2 R/(3GM ) = 2, with a natural mix of compressive and solenoidal modes.We use a fixed quasi-Lagrangian mass resolution of 10 −3 M ⊙ for normal gas cells and 10 −4 M ⊙ for protostellar jet and stellar wind cells, hence each cloud initially consists of 2 × 10 6 gas cells.

RESULTS
All simulations are run for 10 Myr (5 cloud freefall times) and follow the usual qualitative sequence predicted by molecular cloud simulations with star formation and feedback, e.g. as documented in detail in Grudić et al. (2022).Turbulent motions seed dense substructures that go on to collapse into stars, which then proceed to accrete and disperse the cloud via feedback, eventually halting star formation.The clouds form a total of 19,009 stars totaling 1.26 × 10 4 M ⊙ , for an average final cluster mass of 126M ⊙ , i.e. a final star formation efficiency of 6%, broadly consistent with cloud-scale SFEs inferred in the Solar neighborhood (Evans et al. 2009(Evans et al. , 2014;;Megeath et al. 2016).Roughly 80% of the stellar systems appear to be expanding freely at the end of their respective simulations, and the others have a steady halfmass radius and are likely bound and destined to dissolve over a longer timescale.For the purposes of this paper we adopt a broad definition of "cluster" and do not distinguish between expanding associations and bound clusters.
The spectrum of zero-age main sequence masses of stars formed in all 100 simulations is plotted in Figure 1, with the Kroupa (2001) and Chabrier (2005) Solar neighborhood field star IMFs for comparison.This IMF is similar to the one that we have presented for clouds with these parameters and physics in previous work (Grudić et al. 2022;Guszejnov et al. 2022), with a power-law slope of ∼ −2 in the 1 − 10M ⊙ range and a turnover or break around 0.2M ⊙ .The median stellar mass is 0.22M ⊙ and the log-dispersion is σ logM = 0.52 dex, similar to estimates for the Solar neighborhood field population (Chabrier 2005).The high-mass slope α = −2 For comparison we plot estimated ±σ probable ranges of the Kroupa (2001) and Chabrier (2005) Milky Way field IMFs, using parameter uncertainties quoted in Kroupa (2001) and assuming log mc/M ⊙ = log 0.2 ± 0.1dex and σ = 0.55 ± 0.10dex for the Chabrier (2005) form, controlling all IMFs for the number of > 0.1M ⊙ stars.The dashed line plots the best-fit ∝ M −2 powerlaw, and the shaded region plots the ±3σ confidence region of a modified Schechter function fit, sampled from the posterior distribution of the model parameters given the > 1M ⊙ sample.The two lower histograms show the IMF obtained from one 2 × 10 4 M ⊙ cloud run with the same physics model, and a run with the same initial conditions but no feedback.
is shallower than the normally-assumed -2.35, but well within the variation measured between individual clusters and within the uncertain range quoted for the Milky Way field population (Kroupa 2001).

The high-mass IMF truncation
Stacking 100 realizations unveils a new feature of the IMF: a high-mass truncation.Figure 1 shows that there is an unambiguous high-mass break around the sample maximum stellar mass of 28M ⊙ , for the 100 low-mass clusters 2 .To assess this statistically, note that if the ∝ M −2 power-law portion extrapolated to M up = 150M ⊙ , we would expect ∼75 stars more massive than 28M ⊙ in our sample.The probability of obtaining zero stars above the cutoff is 2 × 10 −33 assuming Poisson statistics.Even a steeper ∝ M −2.35 power-law up to 150M ⊙ can be discarded with p = 5 × 10 −12 by this argument.The maximum stellar mass found in these simulations is confidently inconsistent with the size-of-sample effect from sampling a pure power-law extending to the largestknown stellar masses (e.g.Crowther et al. 2016), a stan-dard model assumption.The sample maximum is also significantly smaller than even the smallest cluster mass in the sample (∼ 80M ⊙ ), so it is not directly explained by the limited mass budget.
There is no physical reason to expect a perfectly-sharp truncation, so we model the high-mass (> 1M ⊙ ) IMF as a power-law with a smooth truncation factor.A modified Schechter function with varying sharpness provides a better fit than the standard Schechter (1976) form with p = 1, e.g., its additional parameter is strongly preferred by the Schwarz (1978) information criterion (∆BIC = −17.5).Assuming flat priors with α ∈ [0, −4], log M * ∈ [0, 4], and p ∈ [0, 10] the posterior distribution gives marginalized parameters α = −2.00± 0.02, M * = 24 ± 4M ⊙ , and p = 7 ± 2, and we plot the 3σ confidence interval of the parameter distribution in Figure 1.This sample IMF is also statistically distinct from those found in various previous simulations using the same setup with different initial conditions (Grudić et al. 2022;Guszejnov et al. 2022).Most of those exhibited the same high-mass slope of -2 but produced more-massive stars in clusters with fewer members.For example, in Figure 1 we also plot the IMF from a cloud run with an identical physics setup, surface density, virial parameter, and mass-to-flux ratio, but 10 times greater mass (2 × 10 4 M ⊙ ).This single cloud produces about 10× fewer stars than the 100 smaller clouds.Despite this, the massive cloud produces 3 stars more massive than the 28M ⊙ maximum of the smaller clouds, with a maximum of 44M ⊙ .This strongly suggests that the upper cutoff predicted by the STARFORGE models scales in some way on the properties of the host GMC, in agreement with other recent simulations (Smith et al. 2023).We do emphasize that our model space can still account for the existence of ≳ 100M ⊙ stars -e.g. the 10× higher surface-density cloud in Guszejnov et al. (2022) formed a 107M ⊙ star.
Figure 1 also shows the IMF from a 2000M ⊙ cloud with the same bulk properties as the others, but evolved without any stellar feedback physics.Despite forming only 69 stars, this cloud produced 8 stars more massive than the most massive star in the feedback simulations, with a maximum stellar mass of 224M ⊙ .Clearly then the upper cutoff is highly sensitive to feedback.And because feedback generally operates with variable effectiveness depending on the environment (Fall et al. 2010), the upper cutoff is naturally sensitive to cloud properties.
Due to this variation in the upper truncation in the IMF, no sampling model employing the usual constant ∼ 100 − 150M ⊙ upper cutoff -stochastic or not -adequately describes the statistics of the simulated star clusters.Instead, a model in which the truncation depends on the environment is required.-Expected number of massive (> 8M ⊙ ) stars in a cluster as a function of the number of stars in a cluster, measured over time in all simulations.We count both the number of stars that will eventually exceed 8M ⊙ and the number of stars that currently exceed 8M ⊙ at a given number of cluster members.Horizontal bars show bin widths, vertical bars show counting uncertainties, downward arrows indicate bins with no stars.For comparison we plot the expected number of massive stars when drawing from the total sample IMF (Fig. 1) at random (dashed).
IMF.We count massive stars in two different ways: (1) we include a protostar in the count after its formation if it will eventually be massive, and (2) we include a star in the count only after it has reached 8M ⊙ .
The results, binned over the full set of simulation snapshots (over all stages of cluster formation), are summarized in Figure 2. The average number of eventuallymassive stars tends to be ∼ 3× greater than that expected from random sampling -i.e., massive stars in the simulations tend to start forming earlier than expected if stars of different masses formed in random order.The opposite is true for the counts of presently-massive stars: there tend to be fewer than expected from random sampling as the cluster forms, while it has ≲ 100 stars.Stars of different masses do not finish forming (i.e., reach their maximum mass) at random times.A massive star is never among the first 20 fully-formed stars in any of these simulations.This is likely due to the finite formation time required for the massive stars to accrete their mass, which can be as long as a few Myr (Grudić et al. 2022), comparable to the the total lifetime of the simulated clouds.

Final stellar demographics
Next we examine the statistics of the most important properties for the fully-formed cluster's protometric and feedback properties: the maximum stellar mass M max and the cluster luminosity L bol .In Figure 3 we plot the relation between the number of member stars and L bol and M max respectively for the simulated clusters.Note that the two plots are qualitatively identical because L bol is strongly correlated with M max due to the steep scaling of stellar luminosity with stellar mass.
The L bol and M max relations exhibit ∼ 0.4dex and ∼ 0.15dex of scatter, respectively; the number of stars in the cluster does not uniquely predict either quantity.Remarkably, the degree of scatter agrees well with that predicted from stochastic sampling from the sample IMF, or from a modified Kroupa (2002) IMF with slope α = −2 and hard cutoff M up = 28M ⊙ .However, a standard IMF with α = −2.3 and M up = 150M ⊙ exhibits significantly more scatter in both directions due to the different shape and cutoff.
It is worth re-emphasizing why this result is not necessarily expected: a priori, the data plotted in Figure 3 did not have to exhibit any significant scatter whatsoever, let alone resemble random sampling.In principle star formation could have proceeded as a completely welldetermined, organized process in which the macroscopic system properties tightly constrain the outcome (e.g.Weidner & Kroupa 2006).The highly-nonlinear, dissipative physics at play might have damped out any particulars of the initial seed and regressed toward an emergent attractor behavior, perhaps a detailed self-regulation by feedback.But this is not the outcome realized by the physics of star formation, as captured by the simulations.

DISCUSSION
4.1.The variable high-mass truncation of IMF What does it mean that the most-massive known stars in the Local Group tend to be found in the most-massive young star clusters?Is it a size-of-sample effect from sampling a universal power-law (Larson 1982), or do certain clusters form under special conditions that promote the formation of very massive stars?Distinct high-mass truncations are evident in the raw present-day mass functions of 30 Dor and Carina OB2 (Weidner & Kroupa 2004;Wright et al. 2010), but with the severe caveat that the present-day mass function of these > 4 Myr old systems will also be truncated from the death of the earliest-forming massive stars.This sensitivity to the star formation history introduces some nontrivial modeldependence to the measurement in evolved clusters and associations.Nevertheless, various observational works restricted to < 4Myr old clusters have found evidence against the pure size-of-sample explanation (Hsu et al. 2012;Stephens et al. 2017;Yan et al. 2023).
The situation in the simulations is clearer.In §3.1 we showed that our 100 low-mass cluster formation simulations, taken together with previous STARFORGE simulations of other conditions, suggest an environmentallyvarying high-mass truncation in the IMF that is generally distinct from the usually-assumed universal 100−150M ⊙ cutoff.We do not yet have sufficient coverage of parameter space to map out the detailed scaling, but the 10× more-massive cloud IMF shown in Figure 1 indicates that at least total cloud mass is a significant factor.
One factor affecting the maximum stellar mass in a given cloud could simply be the time that a star has to draw gas from the cloud.These simulated clouds tend to disperse due to feedback after a time proportional to the free-fall timescale t ff , and at fixed mean surface density t ff ∝ M 1/4 , so more-massive clouds tend to live longer, and more-extended accretion is possible.At first glance this explanation seems tempting for the cases shown in Fig. 1, because the ∼ 10× more-massive cloud produces a star ∼ 10 1/4 = 1.8 times more massive.
But the cloud lifetime alone does not explain all variation, and predicts the wrong trend in other cases: Guszejnov et al. ( 2022) found their densest, shortest-lived cloud produced more-massive stars than the fiducial case.The other facet of the problem is feedback: simulations of these same clouds without feedback were easily able to produce 100+M ⊙ stars (Guszejnov et al. 2020(Guszejnov et al. , 2021)).The relation between the cutoff and extensive cloud properties may therefore be understood in terms of the way massive stars tend to form in these simulations: from accretion of extended gas flows over timescales as long as several Myr (Grudić et al. 2022), which are hardly separated in length-or time-scale from the parent cloud size and crossing time.The accretion of massive stars and the resulting influence of feedback on cloud dynamics are therefore coupled.Denser clouds are generally more resistant to feedback, and may therefore produce more-massive stars if all else is equal.

Analytic estimates
To build some analytic intuition for how feedback and cloud properties might influence the upper cutoff of the IMF, let us consider how the stellar mass required to disrupt a cloud might scale under some simplifying assumptions.We estimate an upper bound on M up by determining the mass that a star must have to singlehandedly disrupt a cloud within some fraction of a freefall time N ff ≲ 1, i.e. soon enough to prevent significant additional star formation.
If the cloud is disrupted by ionizing radiation, we may adopt the Hosokawa & Inutsuka (2006) D-type HII region expansion solution and ask what H ionizing photon production rate Q H 0 is required for the bubble to sweep out the cloud radius R within N ff freefall times, t ff = 3π/ (32Gρ).Assuming R is much larger than the Strömgren radius (and hence the self-similar ∝ t 4/7 expansion law applies), we obtain where M cl is the mass of the cloud, X H the mass fraction of H, µ i ∼ 0.61 the mean molecular weight of ionized gas, Σ = M cl /πR 2 the mean cloud surface density, T i ∼ 10 4 K the temperature of ionized gas, and α B (T i ) ≈ 2.6 × 10 −13 T i /10 4 K −0.85 cm 3 s −1 is the case B recombination coefficient (Osterbrock & Ferland 2006, Table 2.1).More conveniently, where Hence the maximum stellar mass scales within increasing surface density and cloud mass.For the 100 simulations with M 4 = 0.2 and Σ 2 = 0.7, the median maximum 3 In identifying Q H 0 with a certain stellar mass we have assumed that all of the ionizing flux originates in the most massive star.In practice, the most massive star accounts for > 50% of the total ionizing flux in 95 of our 100 simulations, so this is not an especially bad assumption.In the stellar mass range ∼ 10 − 30M ⊙ Q H 0 scales so steeply with stellar mass that the most massive star tends to dominate.
stellar mass was ∼ 18M ⊙ , implying N ff ∼ 0.4, fairly consistent with the cloud disruption timescales documented for these simulations (Grudić et al. 2022;Guszejnov et al. 2022).Also notable is the dual role of T i , setting both the gas pressure driving bubble expansion and the efficiency of recombination.Dramatically lower ionizing fluxes are required to disrupt a cloud if T i is greater, e.g. in very low-metallicity clouds lacking metal line coolants.This agrees qualitatively with the 0.01Z ⊙ simulations presented in Guszejnov et al. (2022), which did generally have significantly lower maximum stellar masses (and star formation efficiency) than their Z ⊙ counterparts. 4hotoionization feedback is not universally effective across all of cloud parameter space, so other mechanisms could also affect the maximum stellar mass.The picture may change qualitatively when scaling to sufficiently massive (e.g.≳ 10 6 M ⊙ ) clouds, with freefall times exceeding 10Myr, both because the greater escape speed relative to the ionized gas sound speed will make a given amount of feedback less efficient, and because massive stars can die over this longer timescale.
Consider now the Steigman et al. (1975) solution for an expanding bubble driven by a radial momentum source Ṗ : R = 3 Ṗ /2πρ 1/4 t 1/2 .This could model direct (single-scattered) radiation pressure on dust grains, or a radiatively-efficient stellar wind bubble.The required Ṗ to disrupt a cloud or clump is For photon momentum from a massive star we have assuming a photon escape fraction f esc , so Eq. 5 gives , (7) so again the maximum stellar mass scales positively with cloud mass and surface density, and the numerical value may plausibly account for the upper IMF cutoff we find for M 4 = 0.2, Σ 2 = 0.7 at the order-of-magnitude level, assuming that (1 − f esc )N 2 ff ≈ 1.For plausible values of f esc and metallicity, the momentum present in main-sequence OB stellar winds is generally less than that of photons, so the corresponding expression for stellar wind momentum predicts a much higher stellar mass.Stellar winds could be more relevant if they instead maintained a radiatively-inefficient interior (Weaver et al. 1977), but this appears difficult to realize in turbulent GMCs due to venting and strong mixing at the bubble interface (Lancaster et al. 2021).
Lastly, as stars approach the Eddington limit (≳ 30M ⊙ ) multiple-scattered infrared radiation pressure should also become important (Larson & Starrfield 1971;Wolfire & Cassinelli 1987;Kuiper et al. 2010), likely tapering off the scaling law in a manner that is sensitive to dust properties.But how the bound on the maximum stellar mass from IR radiation pressure connects to the large-scale cloud properties is less clear, because it operates mainly on much smaller scales: within the dense protostellar envelope, just beyond the dust destruction front at a radius of a few AU.

Broader implications
If different cloud properties imprint different upper truncations in the IMF, then there are important implications for the stellar populations of galaxies.GMCs span a broad range of masses and densities within even a single galaxy (Heyer & Dame 2015;Freeman et al. 2017;Sun et al. 2018), so the GMC-level IMFs would vary in their upper truncation and the resulting galaxy-wide IMF would have a different shape.In particular, if M up correlates positively with the host cloud mass then the galactic IMF will generally be steeper than the individual cluster IMFs.This is one way in which our models may be able to reproduce the observed high-mass slope in the field population, which is probably steeper than our measured α = −2.This is also a key prediction of the integrated galactic IMF (IGIMF) theory developed by Weidner & Kroupa (2006) and many subsequent works (see Kroupa et al. 2013 for review), and it has been argued that it could explain differences between FUV and Hα brightness profiles (Pflamm-Altenburg & Kroupa 2008).
However, our eventual star cluster demographics are more random than e.g. the Kroupa et al. (2013) "optimal sampling" model.Crucially, the simulations have a scatter of 0.13 dex in the relation between cluster mass and the most-massive star, whereas optimal sampling has none.This variation compounds to the ∼ 0.4 dex of scatter in luminosity shown in Figure 3. Nevertheless, it is the systematic variation in M up -and not so much the specifics of the sampling process -that is required for the composite IMF to have a different shape from the cluster-level IMF, and our simulations give physical motivation to this hypothesis.Whether the corrections to normal population synthesis are large or small in practice remains to be seen, as it will depend on how specifically M up varies across the vast parameter space of starforming cloud properties.
Because the fully-formed clusters' photometric properties are similar to stochastic sampling from a truncated IMF in the simulations, our result is similar to the one tested observationally by Andrews et al. (2013).They compared photometric models stochastically sampled from a Kroupa (2001) IMF assuming M up = 30M ⊙ and 120M ⊙ respectively, finding that M up = 120M ⊙ is better able to account for the photometry of ∼ 10 3 M ⊙ young star clusters in NGC 4214 5 .However, we emphasize that 1. our clusters are an order of magnitude less massive, and 2. the simulations suggest that the M up should vary with the properties of the cluster or its progenitor cloud, and α can also vary from the exact Salpeter (1955) value.Thus it is not surprising that a model with a hard cutoff of 30M ⊙ with α = −2.35exactly across the entire galaxy is in tension with the data.Given the widely-varying properties of giant molecular clouds, we expect there should be conditions in most galaxies capable of producing stars more massive than 30M ⊙ .
In principle, star cluster photometry may constrain the variation of α and/or M up in greater detail, e.g. by incorporating a model for their scaling and intrinsic scatter in a full Bayesian forward-model (e.g.Krumholz et al. 2019).However in practice α and M up tend to be degenerate because both parameters mainly affect the ratio of ionizing (and hence Hα) luminosity to bolometric luminosity.Photometric models allowing even α alone to vary freely tend to be degenerate (Ashworth et al. 2017), so simulations may be useful for establishing physicallymotivated priors to break degeneracies.

Implications for simulations and feedback modeling
Stochastic sampling from a fixed standard IMF has become increasingly popular in numerical simulations of galaxies and star cluster formation, with the aim of modeling feedback, stellar dynamics, and chemical enrichment more accurately on small scales.But in §3.1 we found that the conditions we simulated do not allow the formation of ≳ 30M ⊙ stars, requiring an environmentally-dependent upper IMF truncation to account for the range of stellar masses observed.Ours is likely an extreme case for very low-mass clusters, but nevertheless the variation of M up could be important for feedback modeling in general (Weidner & Kroupa 2005;Yan et al. 2021;Steyrleithner & Hensler 2023).Fortunately the specific choice of M up will not have a large effect on the rate of core-collapse supernovae from a given stellar population, because the number of progenitors is dominated by∼ 8M ⊙ stars, well below the cutoff seen even in these low-mass clusters.Therefore the overall energy budget of supernova feedback will not be highly sensitive if .But a lower value of M up would produce significantly less radiative feedback from photoionization and earlytime radiation pressure.Compared to a Kroupa (2001) IMF with α = −2.3 and M up = 150M ⊙ , a cluster with our sample IMF with α = −2 and M up ∼ 28M ⊙ would have a zero-age bolometric luminosity 5× less, and produce 10× fewer H-ionizing photons.There would also be significantly fewer stars producing long-lived, powerful Wolf-Rayet winds (Meynet & Maeder 2005).Since these early feedback mechanisms are key for regulating star formation on cloud and star cluster scales, cloud and star cluster properties and the pre-processing of the ISM for supernova explosions may be affected.Variations in the relative abundances of SN progenitors of different masses would also result in measurably-different chemical abundance patterns from supernova ejecta (e.g.McWilliam et al. 2013), and different mass-metallicity relations (Köppen et al. 2007;Yan et al. 2021).
Even if the correct IMF is known, we also found in §3.2 that the way in which stars of different masses are sampled as the cluster formed is not well-modeled by random sampling.The discrepancy between the number of currently-massive stars and random sampling is partly due to the non-negligible time required for massive stars to form.Hence it is possible that stochastic models could be augmented with a model accounting for the finite time required to assemble a star.Without accounting for the time delay, massive stars would be predicted to form too early in a cluster's lifetime, and their feedback would influence subsequent star formation disproportionately.If the scatter in Fig. 3 is less than permitted by observations, as argued by Yan et al. (2023), then some element is missing from the simulations.Intuitively, if star formation proceeds in roughly a turbulent crossing time (e.g.Elmegreen 2000), as it does in these simulations, then there is insufficient time for different parts of the cloud to communicate, and some of the randomness in the initial conditions will survive.Slowing down the star formation process over multiple crossing times would allow the system to become well-mixed and more self-regulated.In STARFORGE, we see something like this in our most highly-magnetized (mass-to-flux ratio of 0.4) cloud, which has a much smoother, more-quiescent star formation history (Guszejnov et al. 2022).
More generally, a potentially-important uncertainty is the initial and boundary conditions.We have simulated an ensemble of isolated, initially-uniform, spherical clouds with Gaussian random initial turbulence, and real GMCs are unlikely to be well-described by any of these properties.Some studies have found that such spherical models may exhibit realistic properties and star formation dynamics once clumpy substructure has had time to form (Rey-Raposo et al. 2015;Priestley et al. 2023).But any dense cloud must inevitably have some continuity with the converging galactic gas flow that formed it, so isolated cloud evolution is unlikely to capture the full dynamics of star formation.The problem of the IMF and its sampling process should be revisited with alternative initial/boundary conditions setups (e.g., Lane et al. 2022), and ideally a larger-scale setup embedded within a galaxy.
And lastly, we re-emphasize that we have intentionally simulated a very particular case of low-mass star formation to stress-test the usual assumptions of stellar population synthesis, and that this is not representative of the massive complexes that form most stars.We have shown that discrepancies with stochastic sampling exist in detail, and this should not change qualitatively with the scale of the star-forming region.But it remains to be seen whether the differences from stochastic sampling are large in general.For instance, if the sampling always regresses to the mean once ∼ 100 stars have formed then the correction would be negligible in any massive starforming complex.But 100 just happens to be the typical cluster size here, so it could also be that clusters regress to the mean on Fig. 2 only once most stars have formed, in which case the difference from stochastic sampling will be important at all scales.It will be possible to clarify this issue with an upcoming statistics suite of larger clouds, currently in progress.

CONCLUSION
In this work we use numerical simulations to test two common assumptions when modeling stellar populations: that stellar masses are randomly, independently sampled, and that the underling IMF is invariant.Within our models, both assumptions break down to some extent: • There is a high-mass truncation in the IMF M up that can generally be significantly lower than largest known stellar masses (Fig. 1), at least in low-mass clusters.This truncation varies with cloud properties such as mass, so different environmental conditions are required to account for the full observed range of stellar masses, and in general it is not valid to assume that all clouds can form arbitrarily-massive stars.
• Stars of different masses do not form in a random order.Massive stars in the simulated clouds begin forming earlier and finish forming later than average, so assuming random sampling during star formation overestimates the number of massive stars that exist at a given time.
Although star formation is not random and uncorrelated in our deterministic simulations, stochastic sampling does appear to provide an adequate description of the end result of star formation, with the major caveat that the high-mass slope and truncation for the system in question must be known, or generally allowed to vary as a model parameter.Stars do not form through a roll of the dice, but the laws of physics may conspire to make it appear so.
Fig. 1.-Histogram of the 19,009 ZAMS stellar masses resulting from the stacked 100 star cluster formation simulations.For comparison we plot estimated ±σ probable ranges of theKroupa (2001) andChabrier (2005) Milky Way field IMFs, using parameter uncertainties quoted inKroupa (2001) and assuming log mc/M ⊙ = log 0.2 ± 0.1dex and σ = 0.55 ± 0.10dex for theChabrier (2005) form, controlling all IMFs for the number of > 0.1M ⊙ stars.The dashed line plots the best-fit ∝ M −2 powerlaw, and the shaded region plots the ±3σ confidence region of a modified Schechter function fit, sampled from the posterior distribution of the model parameters given the > 1M ⊙ sample.The two lower histograms show the IMF obtained from one 2 × 10 4 M ⊙ cloud run with the same physics model, and a run with the same initial conditions but no feedback.
Fig.2.-Expected number of massive (> 8M ⊙ ) stars in a cluster as a function of the number of stars in a cluster, measured over time in all simulations.We count both the number of stars that will eventually exceed 8M ⊙ and the number of stars that currently exceed 8M ⊙ at a given number of cluster members.Horizontal bars show bin widths, vertical bars show counting uncertainties, downward arrows indicate bins with no stars.For comparison we plot the expected number of massive stars when drawing from the total sample IMF (Fig.1) at random (dashed).
Fig. 3.-Relation between number of stars and maximum stellar mass (left) and bolometric luminosity (right) at the end of star formation in the 100 simulations; each datapoint represents one simulated star cluster.Error bars plot the binned median and ±σ range of the simulation results (horizontal bars show bin widths).For comparison we plot the ±σ quantiles of clusters randomly sampled from sample IMF, a Kroupa (2002) with standard slope α = −2.3 and cutoff Mup = 150M ⊙ , and a modified Kroupa IMF with slope α = −2 and hard cutoff Mup = 28M ⊙ .
and T i,4 = T /10 4 K. Approximating the Diaz-Miller et al. (1998) scaling of ionizing flux with stellar mass 3 over the range 10 − 30M ⊙ as a power-law log Q H 0 = 44.7 + 7.3 log (M ⋆ /10M ⊙ ), we find 4.3.CaveatsStar formation simulations require many physical approximations and can only treat a segment of the physical dynamic range self-consistently, so their results should be interpreted carefully.A detailed discussion of known and potential caveats and uncertainties of the STARFORGE numerical model is given inGrudić et al. (2021) §5.2.