Star Formation Regulation and Self-Pollution by Stellar Wind Feedback

Stellar winds contain enough energy to easily disrupt the parent cloud surrounding a nascent star cluster, and for this reason have been considered candidates for regulating star formation. However, direct observations suggest most wind power is lost, and Lancaster21a,b recently proposed that this is due to efficient mixing and cooling processes. Here, we simulate star formation with wind feedback in turbulent, self-gravitating clouds, extending our previous work. Our simulations cover clouds with initial surface density $10^2-10^4$ $M_{\odot} \, {\rm pc}^{-2}$, and show that star formation and residual gas dispersal is complete within 2 - 8 initial cloud free-fall times. The"Efficiently Cooled"model for stellar wind bubble evolution predicts enough energy is lost for the bubbles to become momentum-driven, we find this is satisfied in our simulations. We also find that wind energy losses from turbulent, radiative mixing layers dominate losses by"cloud leakage"over the timescales relevant for star formation. We show that the net star formation efficiency (SFE) in our simulations can be explained by theories that apply wind momentum to disperse cloud gas, allowing for highly inhomogeneous internal cloud structure. For very dense clouds, the SFE is similar to those observed in extreme star-forming environments. Finally, we find that, while self-pollution by wind material is insignificant in cloud conditions with moderate density (only $\lesssim 10^{-4}$ of the stellar mass originated in winds), our simulations with conditions more typical of a super star cluster have star particles that form with as much as 1\% of their mass in wind material.


INTRODUCTION
Self-gravitating molecular clouds are observed to be star-forming in a wide range of environments, from dwarf galaxies (Lopez et al. 2011;Cheng et al. 2021), to the Solar neighborhood (Schneider et al. 2020;Fujita et al. 2021), to galactic center starbursts (Johnson et al. 2015;Leroy et al. 2018;Levy et al. 2021;Emig et al. 2020). On the scale of molecular clouds, the lifetime star formation efficiency (SFE) depends on the processes that truncate star formation (SF) by destroying the cloud, some of which are external and some of which are internal (e.g. Chevance et al. 2020). For giant molecular clouds (GMCs) that create massive stars, internal processes associated with SF feedback dominate the reg-injection of kinetic energy and mass (so-called "mechanical" feedback) may also be important -sources are stellar winds at early stages followed by supernovae (SNe) (Agertz et al. 2013;Rogers & Pittard 2013;Walch & Naab 2015;Iffrig & Hennebelle 2015;Geen et al. 2016). For clouds that are relatively diffuse and have long evolutionary timescales, SNe are effective, and indeed SFRs in large-scale simulations with SN feedback are consistent with observations for these conditions even when dynamical effects of radiation are omitted (e.g. Kim & Ostriker 2017). For very high surface density clouds, where super star clusters (SSCs) form, the SF timescales are too short for SNe to affect evolution significantly, but winds can potentially play a key role.
In Lancaster et al. (2021a) we proposed a new theoretical model for the expansion of stellar wind bubbles under conditions present in the turbulent molecular clouds where stars are born. This model appeals to efficient loss of energy through turbulent, fractal, radiative mixing layers between the hot, shocked wind and the surrounding cloud. This results in bubble evolution that is momentum-driven (ṗ b = const., equivalent bubble radius R b ∝ t 1/2 , also as in Steigman et al. (1975)) rather than energy-driven (Ė b = const., R b ∝ t 3/5 ) as in the classical Weaver et al. (1977) solution. In Lancaster et al. (2021b), we presented a large suite of threedimensional (3D) hydrodynamical simulations to validate this theory. However, our previous work did not take into account the effects of gravity, either for dynamical evolution of the large-scale cloud, or for establishing a self-consistent star formation history (SFH) by inducing small-scale collapse.
Investigating the interplay between gravity, turbulence, and stellar wind feedback is the main goal of the present work. A key motivation is to confirm that energy losses are still large, consistent with the Lancaster et al. (2021a) solution. If not, it would be difficult to understand the observed high SFEs of SSC-forming clouds (Johnson et al. 2015;Leroy et al. 2018;Emig et al. 2020;Levy et al. 2021). If evolution were to follow the winddriven bubble solution from Weaver et al. (1977), the large momentum injection from an inefficiently-cooling bubble would easily destroy GMCs even with very low SFE and very high surface density (Lancaster et al. 2021a). We have previously argued that with muchreduced momentum injection driven by cooling, winds would be compatible with observed high SFEs in SSCforming clouds; here we directly test this prediction with simulations.
In real clouds, multiple feedback processes operate simultaneously, and several numerical simulations have included treatment of both wind and radiation feedback in some form (Dale et al. 2014;Geen et al. 2021;Wall et al. 2020;Grudić et al. 2021). Here, our goal is instead to isolate the effects of winds independent of other feedback processes, and to systematically explore cloud dynamical evolution and SFH for varying parent cloud properties. This allows us to develop and test simple theoretical models, while providing a baseline for more comprehensive studies.
Our self-gravitating models allow us to investigate another important question. In the literature, two modes of energy loss have been proposed in order to explain weak observed signatures of winds: enhanced radiative cooling (Rosen et al. 2014;Lancaster et al. 2021a) or bulk advection of wind energy out of the cloud -so called "leakage" (Harper-Clark & Murray 2009). In Lancaster et al. (2021b), we showed that interface mixing leading to strong cooling can be very effective, and explains the low pressure of X-ray emitting gas observed in pre-SN systems (see discussion in Lancaster et al. 2021a). Here we are able to investigate the roles of mixing/cooling vs. leakage in removing energy from winds that develop in a cloud with self-consistent star formation. As we shall show, both are important but at different times.
Beyond their potential importance in helping to control star formation, winds are also interesting for the way in which they act to chemically pollute the gas in a molecular cloud, leading to so-called "self-enrichment" of subsequent populations of stars that form within the molecular clouds. The details of this process are important to fields ranging from galactic archaeology, where the process of self-enrichment may undermine the strategy of so-called "chemical tagging" (Freeman & Bland-Hawthorn 2002), to multiple stellar populations with different abundances within globular clusters (GCs), where nearly all proposed models for explaining these populations involve "self-enrichment" (Wünsch et al. 2017;Lochhaas & Thompson 2017;Bailin 2018;Tenorio-Tagle et al. 2019;Bastian & Lardo 2018). To be clear, the winds from massive OB stars explored here do not have drastically different chemical abundances from progenitor gas, especially before the Wolf-Rayet phase, and therefore would not, themselves, cause strong selfenrichment. In tracking the re-integration of wind material in subsequent stars we provide a first look at how this process may work for stellar winds more generally. For this reason we refer to the process explored here as "self-pollution." The plan of this Letter is as follows. In Section 2 we briefly review our numerical methods, explained in full elsewhere. Results of our simulations, and connection to other theoretical work, is presented in Section 3.
Finally, we summarize our conclusions and connect to observations in Section 4.

METHODS
We simulate the collapse and subsequent dispersal by stellar wind feedback of self-gravitating, turbulent gaseous clouds. The simulation framework we employ is built on the 3D magneto-hydrodynamics (MHD) code Athena (Stone et al. 2008), though we do not include magnetic fields in this study. We run all of our simulations on a fixed, Cartesian grid using the Roe approximate Riemann solver (Roe 1981), second order spatial reconstruction, and the unsplit van Leer integrator (Stone & Gardiner 2009).
We include cooling in the gas as implemented in the Athena-TIGRESS code base (Kim & Ostriker 2017). This cooling module assumes a simple cooling coefficient Λ(T ) using the fitting formula of Koyama & Inutsuka (2002, see Kim et al. (2008 for corrections) for T < 10 4.2 K gas and the tabulated results from Sutherland & Dopita (1993) for T > 10 4.2 K gas (assuming collisional ionization equilibrium). We adopt a uniform background heating rate of 2 × 10 −26 erg s −1 H −1 which decreases smoothly to zero in higher temperature gas.
We employ sink/star particles as implemented in Gong & Ostriker (2013), with modifications described in Kim et al. (2020). Sink particles, representing stellar (sub-) clusters, are created when a density threshold (based on the collapse solution of Larson 1969;Penston 1969) is reached, the cell forming the particle is at a minimum of the gravitational potential, and the flow is converging along all grid directions.
We include star formation feedback in the form of energy and momentum representing the collective stellar wind inputs from a coeval population of O and B stars (Vink et al. 2001;Vink & Sander 2021;Leitherer et al. 1992). We use prescriptions for the energy and momentum from Starburst99 (Leitherer et al. 1999(Leitherer et al. , 2014 for two different initial stellar mass functions (IMFs, described further below), both employing the Geneva, non-rotating stellar evolutionary tracks (Ekström et al. 2012). These calculations use the wind prescriptions of Vink et al. (2001), which are thought to over-estimate the mass-loss rates,Ṁ w , by a factor of 2-3 (Smith 2014);ṗ w = √ 2L w M w would also be affected by this. This energy and momentum is injected within a radius r fb /∆x = 3.2 of the star particle using the hybrid kinetic/thermal injection scheme described in Lancaster et al. (2021b), which employs a subcell method for dividing the energy and (vector) momentum more accurately over the injection region (see also Ressler et al. 2020). When mass from a wind is injected to the grid, we also inject an equal mass in a passive scalar which is advected with the grid gas, allowing us to track subsequent mixing of the wind material.
In our simulations, we impose a minimum mass M cut = 10 3 M for the initialization of a stellar wind from a sink particle. That is, for all sink particles with mass M * ,sink < M cut accretion is permitted, but for higher mass particles accretion is turned off and wind injection proportional to M * ,sink is initiated. We do this both to reduce the computational cost of injecting wind material and to represent the correlation of winds with high mass stars, which are relatively rare. In addition, since the most massive stars contribute a disproportionate fraction of the power, the number of "important" wind sources is limited. The mean specific wind luminosity and mass loss rates increase while the variances decrease for more fully-sampled mass functions. For example, clusters of mass 500, 1000, and 5000 M sampled from a Kroupa IMF have respective expectation values of the wind luminosity per stellar mass L w /M * = 4 × 10 32 , 10 33 , and 4 × 10 33 erg/s/M with inter-quartile ranges spanning 1.8, 1.5, and 0.5 dex, respectively. All of our simulations have total mass > 10 4 M , well into the fully-sampled regime, so it is appropriate to assign wind power based on a fully sampled IMF. If we make M cut too low, the simulations become expensive due to having many sources, and they are also unphysical in having the wind power nearly equally divided among an unrealistically large number of sources. If we make M cut too high, we artificially suppress feedback. We therefore choose M cut = 10 3 M as a compromise. For a cluster of mass 3 − 9 × 10 4 M (the range for our simulations) divided among particles of mass 10 3 M , there would be 30-90 active sources. This can be compared to the expectation value 29-85 for the number of stars that would provide 95% of the wind luminosity in a cluster of this mass range that fully samples the Kroupa IMF 1 .
We explore the implications of our prescription and parameter choice, and compare with alternative ways of assigning feedback to star particles, in Appendix A. It is important to note that our simulations are not suited to studying the mass spectrum of either stars or clusters; the initial mass spectrum of sink particles depends on 1 Calculations in the last two paragraphs were performed using the Leitherer et al. (1992) wind prescriptions along with the Choi et al. (2016) isochrones at solar metallicity, without rotation, resulting in small differences between the results of Leitherer et al. (1992). Figure 1. Snapshots from simulations at the time when 75% of the final stellar mass has formed. Models R20 to R2.5 are shown from left to right, with all simulations shown here having the same initial turbulence realization (with differing amplitudes). Top row: the integrated gas surface density, together with projected locations of star particles colored by age with size scaled by mass. The "active" star particles producing winds are shown as stars while inactive star particles are shown as circles. Bottom row: slices of gas temperature (centered on the center of mass of the star particles).
grid resolution, and their growth via accretion is truncated when winds turn on after M = M cut .

Simulation Description
We run a set of simulations of initially uniform density spherical clouds of radius R cl and mass M cl placed at the center of a cubic domain with side length L box = 4R cl . Our standard simulations are run on a uniform grid with L box /∆x = 256, though we explore the effects of resolution in Appendix B. In each simulation the cloud is initially surrounded by gas that is 1000 times less dense and in pressure equilibrium with the cloud, contributing about 1% of the cloud mass to the total simulation domain. We initialize a Gaussian random field for the turbulent velocity, with power spectrum |v k | 2 ∝ k −4 for wave modes 2 ≤ kL box /2π ≤ 64. The normalization of the velocity field is chosen so that the initial kinetic energy of the cloud is equal to its gravitational binding energy (i.e. a virial parameter of α vir = 2).
We set M cl = 10 5 M for all of our simulations, and consider models with four different cloud radii, R cl = 20, 10, 5, and 2.5 pc, denoted R20, R10, R5, and R2.5. For each cloud radius we run five simulations with different random-phase turbulent initial conditions. Our parameter space spans nearly 3 orders of magnitude in initial density, enabling us to investigate the regulation of star formation by stellar winds for conditions ranging from typical Milky Way GMCs (Heyer & Dame 2015;Evans et al. 2021) to the birth clouds of SSCs (Johnson et al. 2015;Oey et al. 2017;Turner et al. 2017;Leroy et al. 2018;Emig et al. 2020). Each simulation is run (at least) until only 5% of the cold gas mass remains on the grid, corresponding to a median duration over the turbulent realizations of 14.8, 6.14, 3.49, and 1.75 Myrs for model R20-R2.5, respectively.
For the purpose of setting stellar wind feedback through the SB99 code, we adopt an IMF consisting of stepped power laws dN/dM ∝ M −α . Our main simulation follow a standard Kroupa IMF with the lowest mass range omitted, that is α = 1.3 for 0.1 < M/M < 0.5 and α = 2.3 for 0.5 < M/M < 100 (Kroupa 2001). We have also experimented with a second, top-heavy IMF which has α = 1.8 for the higher mass range, but is otherwise identical. Except as noted, all models shown adopt a standard IMF.

RESULTS
In Figure 1 we show snapshots of the simulations in gas surface density and gas temperature at the time when 75% of the total stellar mass has formed, corresponding to t − t SF,0 = 5.5, 3.0, 1.5, and 0.6 Myrs for the R20, R10, R5 and R2.5 models, respectively. No- Note-Each model is run with five different Gaussian random-phase turbulent velocity field realizations. The quantities given in the rightmost six columns quote the median and full range over these realizations. ε * ,nfb indicates the final SFE reached in simulations run with the same initial conditionas over the same time without any stellar winds. Figure 2. Evolution of the energy, momentum, and bubble radius from model R20 (top row) and R5 (bottom row) for the same simulations shown in Figure 1.The time axis is shifted to coincide with the formation of the first star particle in each model. The ticks at the top of each panel (vertical thin lines in the upper right panel) indicate certain key times: tpr, when enough momentum to disperse the cloud has been acquired (yellow, see Equation 2); tSF, when 95% of the stellar mass has formed (green); t br , when breakout occurs (blue); and t dest , when the cloud has been effectively destroyed (fuchsia). Left Column: The breakdown of input wind energy "spending" over time, normalized by the instantaneous wind luminosity. Losses due to cooling (dark purple), bulk advection off the grid (red), and the luminosity fraction retained on the grid (teal) are shown.
Middle Column: The total radial momentum in wind-polluted (f w,grid > 10 −4 ) gas for the cold/warm component (blue), the hot component (red), and the equivalent momentum injected by winds for a single source (yellow). The reference radial momentum required to disperse the cloud (pr,esc) is indicated as a horizontal dotted line. Right Column: The effective bubble radius for the simulations, as defined in the text (red lines), with reference predictions for delta-function and linear SFHs (solid and dashed black lines respectively).
tably, cooling is so efficient in the high surface density clouds that even though a very large stellar mass has accumulated and produced winds, there is still no significant volume of hot gas on the grid.

Hot gas losses and cloud evolution
In the high-density regions where star clusters form, the pressure from shocked winds would lead to a momentum input rate exceeding that from radiation feedback unless some of the wind energy is lost (Krumholz et al. 2019;Lancaster et al. 2021a). Two mechanisms for wind energy loss are bulk advection of hot gas out of the dense cloud ("leakage," as advocated by Harper-Clark & Murray (2009)), or turbulent interface mixing followed by cooling (explored in depth in Lancaster et al. 2021a,b). Here we show that both effects are important, but at different times.
In this section we shall consider evolution of energy, momentum, and bubble size, focusing on the R20 and R5 models, which are representative of the trends that we see in the other models.

Energy Evolution
From our simulations, we can measure the total instantaneous mechanical luminosity in winds, L w , the net radiative cooling rate in wind-polluted gas 2 ,Ė cool , and the rate at which energy is advected out of the simulation domain in wind-polluted gas,Ė leak (including contributions from thermal and kinetic energy).
Because we are now considering global clouds, we extend the definition of the fractional energy loss rate from Lancaster et al. (2021a) to include both cooling,Ė cool , and bulk advection out of the domain,Ė leak : We refer to 1 − Θ as the fraction of wind energy that is retained within the cloud.
In the left column of Figure 2 we show the evolution ofĖ cool /L w ,Ė leak /L w , and the retained energy fraction (1 − Θ) over time in models R20 (top row) and R5 (bottom row). Losses transition from being dominated by cooling (at early times) to leakage (at late times). The transition between these two regimes is rapid in the denser models. We define the time at which leakage losses dominate cooling losses as the "breakout time" t br .
The low fraction of retained energy indicated by the left column of Figure 2 is consistent with the picture, as previously demonstrated in Lancaster et al. (2021b), in which turbulently-mediated energy losses are sufficient to explain the low observed hot-gas pressures in wind-driven bubbles (Rosen et al. 2014;Lopez et al. 2011). The measured hot gas pressures in our simulations are broadly consistent with the predictions of Lancaster et al. (2021a). The deviations (at most ∼ 0.5 dex) between the simulations and theory are consistent with moderate momentum enhancement (expected in the theory) or decrement (caused by multiple wind sources, discussed below) along with variations inṗ w /M * .
The fraction of energy that remains after cooling and leakage, 1 − Θ, is only ∼ 1 − 10%, which, during the cooling-dominated period, is consistent with the results of Lancaster et al. (2021a,b). For the denser models, other than the earliest times, 1 − Θ is largest after t br . As we shall show, this is due to an expansion in the volume filled by hot gas.

Momentum Evolution
The middle column of Figure 2 shows the evolution of the gas radial momentum, measured with respect to the center of the simulation domain. We separate into a cold/warm component (T < 2 × 10 4 K), which is mostly cloud gas, and a hot component (T > 2 × 10 4 K), which mostly consists of the expanding wind bubble. For both components, measured momentum includes contributions from gas on the grid as well as wind-polluted (f w,grid > 10 −4 ) gas that has been advected out of the domain. We also show the total outward momentum that has been injected by winds in the simulations for an equivalent single central source ( ṗ w (t)dt). This total equivalent injected momentum can significantly exceed the actual momentum in the gas, as we discuss further below.
For reference we mark the radial momentum needed to disperse the gas remaining after star formation, p r,esc , which we define as where ε * ,f is the final SFE of the simulation and v esc ≡ 2GM cl /R cl is our proxy for the escape speed from the cloud. Figure 2 clearly shows that most of the radial momentum in the cold/warm gas is accumulated before t br . That is, the time t pr when the radial momentum held in cold/warm gas exceeds the reference "escape" value of Equation 2 is prior to t br . For the R5 model, the momentum in the hot gas continues to increase after t br , at a much higher rate than during the period when momentum was mostly stored in cold/warm gas. This likely does not occur in the R20 model due to the decrease ofṗ w /M * with the age of the stellar population.
In the analytic model of Lancaster et al. (2021a) we allowed for the momentum added to the shell to be enhanced above the input wind valueṗ w by a factor α p . By comparing the total momentum input by winds (yellow lines) to the momentum that is actually accumulated in the gas, we can attempt to assess to what degree momentum is accumulated beyond the "Efficiently Cooled" solution (α p > 1). It is important to note that in our current simulations there are additional effects that can change our interpretation of α p when performing the comparison we just described. For example, since stars form in over-dense structures, gravity (and the negative momentum associated with collapse) can act to "crush" expanding bubbles. This could contribute to accumulated momentum lower than the total input value, appearing as α p < 1. At the same time, we are not measuring radial momentum with respect to a fixed source, which can be especially confounding at early times when a "super-bubble" enveloping the multiples sources has not yet formed. In addition, the contribution from background momentum at early times can be non-negligible. The combination of all of these effects could lead to either α p < 1 or α p > 1 at early times depending on the realization.
With these caveats in mind, we can see by comparison of the momenta in Figure 2, that in the R20 model we quickly reach α p ≈ 1. This is consistent with the "Efficiently Cooled" solution with some moderate build up of excess momentum when considering the caveats above.
In the dense R5 model, however, it is apparent that α p 1, which is also true for the R2.5 model. This behaviour can be due both to the effects of gravity described above (which are especially strong in the dense models) and another effect. In the dense clouds, when wind bubbles/shells from individual sources collide with one another at small scales, momentum built up in shells can cancel and the corresponding kinetic energy thermalizes. The swept-up shells in these dense models are dense enough to cool efficiently; with reduced pressure, there is no re-acceleration. This results in a bubble evolution that is technically less than momentumconserving.

Bubble Radius Evolution
In the right column of Figure 2 we show the temporal evolution of the effective radius of the wind-driven bubble, R b , defined as Here V b is the volume of the wind driven bubble, defined as the total volume of cells with temperature T > 2 × 10 6 K or with radial velocity v r > V w /2 as measured with respect to the origin of the simulation domain. Note that in the simulations the bubble is not always a single entity but rather an amalgamation of bubbles from each source which start to expand individually before percolating to become a larger "superbubble." We use a single effective bubble radius here to give the simplest phenomenological comparison to our previous models, which assume a single, central source.
We also sum over all wind inputs in making a prediction for the bubble radius. For reference, in the right column of Figure 2 we show the expected bubble radial evolution for the "Efficiently Cooled" (momentum-driven) solution for two different assumed SFHs: 1) An instantaneous, delta-function SFH, which assumes that the final SFE is reached as soon as the first star forms (solid lines) and 2) a linear SFH given by where t ff,0 is the initial free-fall time of the cloud and the SFE per free-fall time, ε ff = 0.25, is set based on a fit to the early-time SFHs of the simulations (dashed lines). The first method will tend to over-predict the radial size since at early times it overestimates the stellar power driving the wind bubble. In principle, the latter model should work better, since it more closely follows the actual SFH, and indeed Figure 2 shows somewhat better agreement of this model with the numerical results at the time of breakout. Neither model accounts for the effects of gravity (which are most important in the denser clouds), however, or for the varying strength of the stellar winds with age (which are most important in the longer lived, larger clouds).

Star Formation Efficiency
Like other forms of feedback, winds can help to suppress star formation, but their effectiveness can be reduced by energy losses. Having established that wind energy is efficiently lost through both cooling and leakage, it is interesting to investigate the role winds play in suppressing star formation. To do this, we make use of theories that have been proposed in the literature to predict SFEs subject to feedback mechanisms in the form of momentum injection rateṗ. In previous work that accounts for cloud inhomogeneity, this momentum has been conceived of as associated with radiation, but the physical picture of competition between radially outward forces from feedback and inward forces from gravity is quite general, and can be equally well applied to momentum injection from a wind.
Since we evolve our simulations until most of the gas is either driven off the grid or locked in star particles, a straightforward way to assess star formation regulation is to measure the lifetime SFE As a reference value, in Table 1 we also give the SFE reached in simulations run for the same time period but without any feedback from stellar winds, which we denote ε * ,nfb . Comparison between the simulations with and without feedback shows that at the highest surface density, only a few percent of the mass in the cloud is driven out by feedback. In contrast, half of the cloud's mass loss can be attributed to feedback in the lowdensity cloud.
We shall compare the measured SFE in the simulation as a function of cloud parameters to several different theoretical predictions. All of the theoretical models explored here will assume that the wind bubble is in the "Efficiently Cooled" limit in which there is no enhancement of momentum relative to the initial wind injection. That is, we assume α p = 1, even though, as we have seen, it can be either slightly greater or less than one depending on the environment.
In our predictions, we will use the mechanical luminosity of the winds per unit stellar mass and the total radial momentum input rate for an equivalent single central source, here L w is the total mechanical luminosity of all winds collectively, and V w is the wind injection velocity. For our analytic predictions below, we simplify by using the initial value of Ψ w determined by the SB99 code for a fully-sampled Kroupa IMF, which is approximately 10 34 erg/s/M . We note that in detail, the IMFaveraged Ψ w stays relatively constant until an increase due to Wolf-Rayet Stars around t ≈ 3 Myr and subsequent decline beginning at t ≈ 4 Myr (see Figure 3 of Lancaster et al. (2021a)). While our numerical models employ a time-dependent IMF-averaged value of Ψ w for each source, our analytic models adopt this simpler approach of constant Ψ w . Suppression of star formation in clouds occurs by driving mass out of the system. We consider a fluid element of mass M and cross-sectional area A at a distance r from a central cluster. If it absorbs all incident momentumṗ =ṗ w A/(4πr 2 ) from a radially-flowing wind, the Eddington ratio between the wind forceṗ and the gravitational force GM * M/r 2 of the cluster on the fluid element iṡ Here, the gas surface density of the fluid element is Σ = M/A, and we have defined the Eddington surface density Σ Edd ≡ Ψ w /(2πGV w ) as the maximum value for which the feedback wind momentum exceeds gravity; for our adopted parameters this is Σ Edd = 570M pc −2 .

Lancaster et al. (2021a) Prediction
We begin by reviewing the simple prediction of Lancaster et al. (2021a) for the SFE, obtained under the assumption that all star formation occurs instantaneously at the center of a cloud, creating a cluster of mass M * = ε * M cl , and that the rest of the mass (1 − ε * )M cl is swept into a spherical shell at the cloud radius R cl . For the force from the wind to exceed the gravity of the cluster plus the self-gravity of the shell, the conditioṅ must be satisfied. The above can be written as where Σ cl ≡ M cl /πR 2 cl . This inequality may be solved to obtain the minimum ε * given Σ Edd /Σ cl . For low Σ cl , i.e. Σ cl Σ Edd , this yields ε * ≈ Σ cl /(8Σ Edd ) ∝ Σ cl , while for high Σ cl , ε * ≈ 1 − 4Σ Edd /Σ cl (Lancaster et al. 2021a).
The above parallels simple scaling arguments for SFE as limited by radiation pressure, where instead ofṗ w the radiation momentum injection L/c (for L the bolometric luminosity) is used (e.g. Fall et al. 2010;Murray 2011;Raskutti et al. 2016;Li et al. 2019). The assumption of instantaneous star formation is, however, clearly inconsistent with the more realistic simulations carried out here.

Steady Star Formation
An alternative approach allows for the SFH to be extended in time, rather that instantaneous, while still occurring at the center of the cloud. Following the same derivation for bubble expansion as in Lancaster et al. (2021a) (ignoring the effects of gravity and assuming a constant momentum input rate per unit stellar masṡ p w /M * ), while setting M * /M cl = ε * = ε ff t/t ff,0 , the bubble expands as where we have defined the initial free-fall time of the cloud as If we then assume that star formation is immediately halted once the bubble reaches the edge of the cloud, we can obtain an estimate for ε * by setting R b = R cl in the above and solving for t. This gives us the final SFE ε * ,f = [3π 2 ε 2 ff Σ cl /(16Σ Edd )] 1/3 for an assumed ε ff . This model may be accurate in large, low-density clouds where the effects of gravity on the wind bubble expansion are actually quite small. However, there is no consideration for the ability of the swept-up shell to escape the gravity of the system, and if the bubble expands slowly enough this model will predict ε * > 1, which is clearly unphysical.

Thompson & Krumholz (2016) Prediction
The models of Section 3.2.1 and Section 3.2.2 assume the driving source "sees" only a single gas surface density over all solid angles, given by the mean value for the cloud as a whole. In reality, since the clouds in which stars are born are supersonically turbulent, the surface densities viewed either externally or looking outward from the center will have a log-normal distribution (e.g. Raskutti et al. 2017;Kim et al. 2019). Considering momentum injection associated with radiation pressure, Thompson & Krumholz (2016) pointed out that directions where the total line-of-sight surface density is sufficiently low can be "super-Eddington" even if the cloud as a whole is sub-Eddington. Along these lines of sight, gas may be driven out of the cloud. Thompson & Krumholz (2016) denote the massweighted distribution of (logarithmic) surface density x ≡ ln Σ/Σ by p − (x), whereΣ ≡ ΣdΩ/4π is the sky-averaged surface density and the mean of p − (x) is σ 2 ln Σ /2. Integrating over the distribution, Given a distribution of mass, Thompson & Krumholz (2016) posit an instantaneous loss rate of mass ejected from the cloud aṡ where M g (t) and t ff (t) are the instantaneous gas mass and free-fall time of the system 3 , and x crit is determined by the condition for a structure to be super-Eddington. At the same time, star formation is assumed to folloẇ where ε ff is a constant SFE per free-fall time, so that allowing for mass loss to both outflow and star formation (i.e. M g (t) = M g (0) − M * (t) − M cloud,w (t)), the average surface density evolves in time asΣ = 4R cl (t)ρ(t)/3. Thompson & Krumholz (2016) allow the cloud radius to evolve as well, but we keep this fixed. Adapting from the case of radiation feedback to wind feedback, the limit in the integral over the distribution becomes x crit = ln(Σ TK /Σ) where Σ TK = Σ Edd /(1 + M g /M * ) is the Eddington surface density as defined in Thompson & Krumholz (2016), allowing for gravity of both stars and gas.
The above equations can be integrated jointly forward in time until all mass that was originally in gas is either in wind or stars. The final stellar mass then determines the final predicted SFE. This is a function of the initial cloud properties, the adopted Ψ w , and the values assumed for the parameters ε ff and σ ln Σ . Raskutti et al. (2016) also developed a prediction for the SFE of a cloud regulated by radiation pressure feedback that takes in to account the log-normal distribution in surface densities seen by a source particle. The theoretical model is very similar to that of Thompson & Krumholz (2016) except that, instead of integrating ODEs they consider the instantaneous state of a cloud. If some fraction, ε * , of the original cloud mass has been gathered into a central star cluster and the remaining gas is distributed around this source with a log-normal distribution in surface densities, they calculate the fraction of the original cloud that would be super-Eddington as ε of (ε * ) = (1 − ε * )ζ − (x crit (ε * )) ,

Raskutti et al. (2016) Prediction
where the only differences between this x crit and that of Thompson & Krumholz (2016) are (1) a factor of 1/2 in the gas gravity term, and (2) the inclusion of a parameter to account for the possible contraction of the cloud as it evolves (which exists to a different degree in the Thompson & Krumholz (2016) model). Raskutti et al. (2016) argue that star formation proceeds until a SFE of is reached. At this point, the maximum fraction of the original gas satisfies the Eddington condition and will (in their model) be immediately ejected. This corresponds to a minimum prediction for the SFE, since there is still some gas left in the shell at this point that is sub-Eddington. If all of this remaining gas collapses to form more stars then the final SFE will instead be To obtain SFE predictions (upper and lower bounds) a value for the parameter σ ln Σ must be adopted, but no assumption for ε ff is needed.

Comparison to Simulations
In Figure 3 we show the final SFEs for all of our simulations, displaying results based on the five different runs for each parameter set. For comparison we show the predictions of each of the analytic models described in the previous sections.
The assumptions and parameters used for each of these theories are as follows • Lancaster+21a: The only free parameter is the momentum input rate per unit massṗ w /M * = 2Ψ w /V w . We adopt a value 8.6 km/s/Myr, determined using the average over the first 1 Myr using SB99 (Leitherer et al. 1999). We use this value oḟ p w /M * in all models.
• Steady SFR: We adopt ε ff = 0.25, broadly consistent with the SFH in our simulations.
• Thompson & Krumholz : As above we adopt ε ff = 0.25, and for the width of log-normal in surface density we adopt σ ln Σ = 1.5 as a value consistent with many of the simulations over much of the evolution.  Table 1) we show the median SFE as a point, with the range of the SFEs from all runs shown by the vertical error bars. Predictions for the SFE based on several different models are shown, as labeled. For lognormal-based models, solid lines denote σ ln Σ = 1.5 and dashed lines denote σ ln Σ = 0.75 (see text).
• Raskutti+16 : We display the minimum prediction given by Equation 17, also assuming no expansion or contraction of the cloud (x = 1 in their notation). We adopt σ ln Σ = 1.5 as above.
We first note that the agreement between the Lancaster et al. (2021a) prediction and our simulations is best for the highest surface density clouds. This is because the SFH most closely resembles an instantaneous starburst in the high-Σ regime, for which the star formation timescale is much less than the cloud disruption timescale (as we will see in Section 3.3). However, the Lancaster et al. (2021a) prediction significantly underestimates the measured SFE at low surface densities.
The steady star formation prediction does reasonably well at low surface densities, but quite poorly elsewhere. This is because the condition chosen to turn off star formation is arbitrary (when the bubble reaches the edge of the cloud) and does not take into account the effect of gravity in any way. There is also no built-in limit of star formation to prevent ε * from exceeding unity, as it does for Σ cl 10 4 M pc −2 .
The Thompson & Krumholz (2016) model does a much better job at capturing the overall trend in the simulation results, while still slightly overestimating the measured SFE. In principle the model could be made to fit better by adjusting its parameters, but our specific choices for the parameters are motivated by the simulations themselves.
The Raskutti et al. (2016) model, which has no explicit time dependence and assumes a constant PDF in Σ, does reasonably well at high Σ cl and slightly overpredicts ε * at low Σ cl . The latter could be due to a break-down of the single, centralized source or constanṫ p w /M * assumptions (see below).
Since in practice σ ln Σ is not constant in time, we additionally show the predictions of both the Thompson & Krumholz (2016) and Raskutti et al. (2016) models for σ ln Σ = 0.75 in Figure 3 as thin dashed lines. This value is chosen as comparable to the measured PDF variance at the onset of star formation. Adoption of this value leads to a slightly higher (lower) predicted SFE at higher (lower) Σ for both models. However, we caution that there is not a strong physical motivation for preferring this choice. Rather, the comparison provides a feeling for the parameter sensitivity of the two models.
As mentioned above, we have also investigated the effects of top-heavy IMFs in our denser models (R5 and R2.5), motivated by observations of star clusters in dense environments (McCrady et al. 2005;Lu et al. 2013;Hosek et al. 2019). We find that the greater strength of the winds (largerṗ w /M * ) with a top-heavy IMF is noticeable, reducing the final SFEs from 71 to 53% and 85 to 75% in the R5 and R2.5 models respectively.
Finally, as we alluded to above, we note that the approximation ofṗ w ≈ const. breaks down once the Wolf-Rayet stage of stellar evolution occurs around ∼ 3 Myr (see Lancaster et al. (2021a) Figure 3), withṗ w increasing by a factor ∼ 3 and then declining again over the next few Myr. As we shall show in the next section, the two lower-density models have lifetimes that extend into the Wolf-Rayet phase, so one might expect slightly lower ε * than predicted under the model assumption of constantṗ w /M * , as is indeed the case. Because their lifetimes are shorter, the higher Σ cl simulations would not be significantly affected by the increase ofṗ w in the Wolf-Rayet stage.

Timescales
In Figure 4 we show various time-scales measured from the simulations, and how they depend on the surface density of the cloud. All timescales represent the time elapsed between the onset of star formation and the physical event described below, and in Figure 4 are plotted in units of the initial free-fall time in the cloud.
t br : The "breakout time" is when energy losses from the wind become dominated by bulk advection (leakage) of hot gas out of the cloud, rather than turbulent mixing and subsequent cooling (see Section 3.1.1). t SF : The "star formation time" is when 95% of the final stellar mass has been formed in the simulation.
t dest : The "destruction time" is when the gas mass remaining on the grid drops to 5% of the initial cloud mass.
t pr : The "radial momentum time" is when the total radial momentum carried in the cold/warm (T < 2 × 10 4 K) gas exceeds the "escape momentum" defined in Equation 2.
Considering the various timescales together, Figure 4 reveals a common sequence in which winds halt starformation while cooling is still dominating the wind's energy losses and acceleration is still occurring, but the winds eventually break out from the dense cloud gas, and the remaining gas is gradually driven out of the simulation domain. The momentum needed to disrupt the cloud is always accumulated in the gas before advection begins to dominate over cooling losses (t pr < t br ). Interestingly, it is also broadly true that t dest ≈ 2t SF .
The details of evolution nevertheless vary over the range in surface densities probed by our simulations. At the lowest surface densities probed here, the momentum needed to eventually disrupt the cloud is acquired before star formation halts. In this case, star formation is a continuous process that proceeds in parts of molecular clouds while it is simultaneously being halted in other parts. However, at higher surface density, star formation essentially is halted and then winds engage in an extended battle to turn around and expel lower-density gas while also trying to break out of the cloud. In these denser clouds, t br , t pr , and t dest all occur roughly simultaneously. For the highest surface density cloud, as we have shown, very little material is in practice driven away by feedback.
Since the model of Thompson & Krumholz (2016) provides for a full SFH, it gives a prediction for t SF . We include this prediction in Figure 4. The prediction overestimates the star-formation timescale by a factor of ≈ 2 for the σ ln Σ = 0.75 case (dashed line) and by a larger factor for the σ ln Σ = 1.5 value (solid line). The fact that star formation halts much more rapidly than the model predicts indicates that the model assumption of a wind acceleration timescale t ff (t) is not borne out in our simulations. Potentially, other assumptions for the acceleration timescale (e.g. accounting for the distribution of Eddington ratios) could be explored, but a reformulation of this theory is beyond the scope of the present work.

Wind Pollution
One of the central goals of this work is to provide a first investigation of the incorporation of wind material in new stars, and its dependence on cloud environment. As a caveat, we note that the winds investigated here, from massive OB stars, are not very chemically distinct from their parent gas, especially before the onset of the Wolf-Rayet phase (Leitherer et al. 1992), and thus the "pollution" by these winds would not produce compositional signatures. The winds that are thought to selfenrich forming globular clusters presumably would have somewhat different physical properties (e.g. mass-loss rates, wind velocities), and this difference should be kept in mind when interpreting the results we present.
As discussed in Section 2, we use a passive scalar, injected with the wind material, to track how wind material is mixed into the surrounding gas and subsequentlyforming stars. In Figure 5, we present the evolution of the global fraction of mass, in stars (f wind , * ) and in gas (f wind ), that originated from stellar winds. All simulations depicted in Figure 5 have identical random-phase initialization of their turbulence.
From Figure 5, the evolution of f wind in the gas (solid lines) is remarkably similar across the simulations when the evolution is shown in units of the free-fall time. This is consistent with each of the simulations having nearly the same (linear) SFH at early times with ε ff ≈ 0.25, which then implies M w ∝ t 2 (in practice, the slope in the simulations is closer to M w ∝ t 5/2 ). These lines significantly deviate when the wind bubbles break out, first causing f wind to decrease as much wind material is ad-vected off of the grid and then increase as all remaining, non-wind material is dispersed.
For the denser clouds, the fraction of wind material in star particles (dashed lines) very roughly tracks the evolution of f wind in the gas, until star formation slows (dashed lines becoming transparent indicate 95% of star formation is complete). As we have seen in Figure 4, in denser clouds, star formation takes place over more freefall times, resulting here in a more self-polluted stellar population.
In the denser clouds, polluted stars tend to form in gas that is preferentially more polluted than the average gas (points fall above the solid lines), indicating that mixing is less efficient in these clouds. This is somewhat counter-intuitive since one might expect that mixing would be more efficient in the denser clouds where the turbulent Mach numbers are higher, promoting a larger degree of turbulence compared to radial expansion of wind bubbles (see Figure 10 of Lancaster et al. (2021b)). However, as we have seen, when gravity is considered wind bubbles are much more easily contained in the dense clouds. Indeed when we look at the massweighted distribution of f wind in the simulations, the denser clouds are much less well-mixed, explaining the trend we see here.
To the best of our knowledge, this sort of selfconsistent analysis of pollution by the winds from newlyformed stars is the first of its kind that has been shown in the star cluster formation literature. We anticipate that future analyses of this kind, including chemical enrichment from other wind sources, can be used to better understand proposed formation mechanisms for multiple populations in Globular Clusters (Bastian & Lardo 2018;Gratton et al. 2019;Bailin 2018), as well as chemical evolution within normal, redshift-zero star-forming clouds. However, it is important to keep in mind that there are many potentially important physical phenomena which we have not included that may or may not be important in different scenarios (radiative feedback, magnetic fields, dust, varying metallicity).

SUMMARY AND DISCUSSION
The simulations and analyses presented here explore the applicability of the "Efficiently Cooled" stellar wind bubble model (Lancaster et al. 2021a,b) in the more realistic scenario where stars form self-consistently through gravitational collapse, which competes with cloud dispersal by stellar winds. We show the core principle of the "Efficiently Cooled" model holds: enough energy is removed through radiative cooling to make the wind bubbles momentum-driven rather than energy-driven (see Figure 2). Figure 5. Evolution of the global fraction of wind mass f wind in the gas (solid lines), in the stars overall (dashed lines), and in individual star particles (points, where more opaque points indicate more massive star particles). We show results as a function of time in units of the initial free-fall time (top row), and as a function of the fraction of the final stellar mass that has formed (bottom row). For the individual star particles, we show the fraction of their initial mass in wind material as a function of the time when they first formed (or fraction of stellar mass formed when they first formed, in the bottom row). Each panel also includes f wind in gas from other models as thin, transparent lines. We show f wind for the gas (stars) as opaque curves up until t br (tSF), after which more transparent lines are shown. We note that the dashed lines include all star particles, even those with negligible enrichment that do not appear on this plot.
In the absence of losses, stellar winds would dominate other "early feedback" in star-forming clouds. The weak observational signatures of winds has led to suggestions that energy is efficiently removed by either (i) bulk advection of hot gas out of the cloud, or (ii) radiative cooling mediated by turbulent mixing. We show that over the time-scales relevant to star formation (see Figure 4), it is the latter process that is most important for removing energy injected by stellar winds. This conceivably depends upon the exact geometry of the cloud, and it is plausible that for highly porous clouds, bulk advection could dominate energy losses.
We apply theoretical models to make predictions for SFEs in clouds on the basis of "Eddington ratio-like" considerations of competition between gravity and wind forces, taking into account turbulence-driven substructure. We compare predictions of these theoretical models for SFEs to the numerical results from our simulations (see Figure 3). While our simulations include realistically time-varying specific wind momentum input rate, the theoretical models assumeṗ w /M * is constant, for simplicity. Despite this, we find that the models derived from Raskutti et al. (2016) follow the numerical results fairly well, while the Thompson & Krumholz (2016) models slightly over-estimate the final SFEs. Model assumptions regarding (i) constant momentum input rates, (ii) a constant SFE per free-fall time, ε ff , and (iii) a log-normal distribution of surface densities with constant variance are not fully satisfied, however, which may account for some of the deviations between numerical results and models. We also compare the Thompson & Krumholz (2016) prediction of the star formation duration t SF to our numerical results, and find that the former overestimates our measurements of t SF , which we link to their assumption for the timescale for gas expulsion.
We do not attempt to make a direct comparison between our simulated SFE and observed values for different cloud regimes, given that important feedback effects are left out in our numerical models. Nevertheless, it is conceivable that stellar wind feedback is the dominant star formation regulation mechanism in very dense clouds (Levy et al. 2021;Lancaster et al. 2021a), and it is interesting to compare our models to observations in this regime. The SFEs we measure at high surface density are broadly consistent with the range inferred from observations of SSC-forming clouds at the centers of galaxies, which are generally ε * 70% (Leroy et al. 2018;Emig et al. 2020). Due to the range of uncertainties associated with the observations, empirical SFEs are also broadly consistent with the SFEs in our simulations that assume stellar populations with top-heavy IMFs, as might be expected in these clouds (Hosek et al. 2019).
In typical Milky Way (or LMC/SMC) GMCs, similar to our lower density models, momentum-driven stellar wind feedback is not expected to be the dominant star formation regulation mechanism (Lancaster et al. 2021a). However, it is still interesting to quantitatively assess the effects of winds in this regime. Based on the work of Lopez et al. (2011);Rosen et al. (2014) and others, there are large discrepancies between observations of the hot gas pressure produced by winds and the classical Weaver et al. (1977) model, andin Lancaster et al. (2021a,b) we proposed that energy losses due to turbulently-mediated radiative cooling could explain these discrepancies. Inspection of the hot gas pressure in the more-realistic simulations presented here bear out this proposal, where our R20 simulations show hot gas pressures ranging from 10 5 − 10 7 K cm −3 at t 3 Myr, broadly consistent with X-ray observations of nearby massive star-forming regions (Lopez et al. 2011. Finally, using a passive scalar that is injected with the wind material, we investigate the process of wind pollution. The fraction of the gas originating in winds scales roughly as f wind ∝ t 5/2 ∝ M 5/2 * , and for the densest cloud (model R2.5) the mean wind-mass fraction of the stellar population generally follows this up to the time of wind breakout. For the less-dense clouds the mean pollution of the stellar population increases less steeply in time, and individual star particles can have significantly higher or lower f wind than the mean value at the time they form.
The increase of f wind in time, combined with the fact that star-formation takes place over more dynamical times in denser clouds, implies stellar populations with larger fractions of wind material in denser clouds. This investigation opens the door for future direct numerical simulations to explore wind pollution and chemical self-enrichment in cluster-forming clouds in more depth (incorporating additional physical processes and winds that are more chemically enriched), with implications for the study of galactic archaeology and multiple populations in globular clusters (Freeman & Bland-Hawthorn 2002;Bastian & Lardo 2018). ACKNOWLEDGMENTS L.L. thanks Erin Kado-Fong, Todd Thompson, and Madeline Lucey for useful comments and references. We would like to thank the referees for useful comments that helped to improve the paper. This work was partly supported by the National Science Foundation (AARG award AST-1713949) and NASA (ATP grant No. NNX17AG26G). J.-G.K. acknowledges support from the Lyman Spitzer, Jr., Postdoctoral Fellowship at Princeton University. Computational resources were provided by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center at Princeton University.
Software: Athena (Stone et al. 2008;Stone & Gardiner 2009), astropy (Astropy Collaboration et al. 2013, 2018, scipy (Virtanen et al. 2020), numpy (Harris et al. 2020), IPython (Perez & Granger 2007), matplotlib (Hunter 2007), xarray (Hoyer et al. 2017 In this appendix we explore the effects of varying the prescription for the assignment of stellar wind power, L w , to a star particle. All simulations presented in this section have identical conditions to those simulations outlined in columns 2-5 of Table 1 but are run at half the resolution (2∆x as given in column 6 of Table 1).
It has been shown in the literature (e.g. Grudić & Hopkins 2019; Smith 2021) that prescriptions for IMF sampling in imposing feedback can have significant consequences for the key results of a given simulation, despite being a rather arbitrary prescription for physics that is sub-grid to the simulation. To provide some sense of the systematic uncertainty in our results that arises from "model variance" in our simulations, we present here a few numerical experiments covering a small range of different prescriptions for the assignment of feedback to star particles.
The first prescription we use is the same outlined in the text, with M cut = 10 3 M . We also run simulations here with the same prescription except with M cut = 10 2 M . Finally, we include a third prescription based on the Poisson sampling method outlined in Su et al. (2018) and also used in Grudić & Hopkins (2019). This method works by sampling a number of "O-stars" from a Poisson distribution whose mean depends on the mass of the sink particle to which we are assigning feedback. Specifically, for each star particle of mass M * ,sink we sample from a Poisson distribution of mean µ = M * ,sink ∆m , where ∆m = 100 M in our method and is meant to represent the mass at which one might expect to have at least one O-star. We will refer to the number of O-stars sampled from this distribution as N O,sink . Wind luminosity is then assigned to each O-star sampled in this way as The total wind luminosity of a given sink particle is then N O,sink L w,O−star . This method was originally implemented in Su et al. (2018) for galaxy scale simulations where, in general, M * ,sink ∆m. Because of this, it was reasonable to assign N O,sink in an unbiased way after one sampling from a Poisson distribution as outlined above. However, in our simulations, it is quite possible that when a star particle initially forms it has M * ,sink ∆m. In our scenario, when one initially samples from the Poisson distribution it is likely that zero O-stars will be assigned. Since our sinks can grow in mass, we want to allow them eventually to sample from the Poisson distribution again. However, if we allowed this re-sampling every time step, we would essentially be biasing the sampling process, "flipping the coin" over and over again until we got N O,sink > 0 and then halting our sampling. To mitigate this process we sample for N O,sink once when the star particle is formed and then afterwards only every time it has doubled in mass.
In Figure 6 we show a comparison of the results using the different prescriptions described above. For each model we show the history of how some gas forms stars while other gas is expelled from the cloud. In general, the models with M cut = 10 2 M have lower final ε * than when M cut = 10 3 M , and result in more massive outflows earlier in the simulations. This is to be expected, since reducing M cut to 10 2 M from 10 3 M effectively increases the available mass that has active winds, increasing the wind luminosity and helping to drive material out of the cloud.
Interestingly, the Poisson sampling prescription results in very similar behavior to the M cut = 10 2 M simulations, with the outflow and star formation histories only differing significantly for the densest simulations/smallest clouds. From inspecting the wind luminosities and distributions in M * ,sink and N O,sink in these simulations we have determined that this stronger outflow is due to a larger effective wind luminosity caused by more sink particles with M * ,sink < ∆m having N O,sink ≥ 1 than expected from direct sampling. This is due to the above-mentioned sampling problem that was evidently only partially mitigated in our denser simulations by the fix we applied. This could potentially be ameliorated by a stricter sampling prescription, but exploring that possibility is beyond the scope of these tests.
As a comparison for the amount of variation between models, we also show the range of SFHs of the simulations presented in the main text which differ only in the initial turbulent velocity field (but are run at higher resolution than the simulations discussed above). These ranges are shown as shaded regions in Figure 6. For the lowest density clouds Figure 6 shows that the spread in results for varying luminosity assignment choices is similar to the variation seen between models with different turbulence seeds, but this is not true in the higher density simulations. For this reason, the prescription for assigning feedback to sink particles remains an important factor to consider when investigating star-forming cloud simulations.

B. RESOLUTION STUDY
Here, we provide a resolution study of our simulations for the R10 cloud model. The simulations presented here use a single set of random phases to initialize the turbulent velocity spectrum and all of the same numerical techniques as described in Section 2. Here we use the same prescription for assigning wind feedback to sink particles as described in the text, with M cut = 10 3 M . All simulations have the same parameters as outlined in row 2 and columns 2-5 of Table 1. The only parameter varied is the spatial resolution, ∆x, given in column 6 of Table 1. Specifically, we perform a simulation at half the nominal resolution given in Table 1 and a simulation at twice the nominal resolution, but for about half of the full simulation time due to computational cost. Figure 7 presents evolution of several key quantities across the different resolutions. One quantity which we have not shown elsewhere is M >Mcut which we define as the total mass of star particles that are more massive than M cut , Figure 6. A comparison of results from different methods of assigning wind luminosity to star particles. Line style denotes the assignment method (see text). For each model, we show mass of gas on the grid (black), mass of stars (yellow), and mass that has flowed out of the simulation domain (red). The shaded regions show the range of histories given by the higher resolution simulations described in the main body of the text, which differ from one another only in the random phases used to generate the initial turbulent velocity field.
the threshold used to determine whether or not a star particle has a wind. The total stellar mass and Ψ w determine the total wind luminosity (modulo the effects of independently aging star particles, which causes Ψ w to vary mildly from source to source).
The simulations run at low and fiducial resolution (dotted and dashed lines) are quite similar: they have comparable energy loss histories and corresponding radial momentum (t pr ) and breakout (t br ) times, the evolution of their bubble sizes track one another, and their star formation histories are nearly identical. However, there is a very slight increase in the total stellar mass formed going from the low to the fiducial resolution. For the high resolution model, this trend continues and is greater in magnitude, with both a higher star formation rate and higher lifetime SFE than in the fiducial and low-resolution models. In contrast, the histories of the mass in particles that produce winds, M >Mcut (blue lines in bottom right panel), is very similar across all resolutions. This is also mirrored by the similarity in the equivalent injected total wind momentum (yellow lines in upper right panel of Figure 7).
The differences/similarities in M * vs. M >Mcut histories reflect both the reduction in the mean initial sink particle mass at higher resolution, based on the particle creation criteria (see also Grudić et al. (2021) or Haugbølle et al. (2018) for a similar effect), and the prescription for assigning winds to sink particles. This comparison also provides additional evidence of the importance of feedback to limiting star formation. In the highest resolution simulations, Figure 7. A resolution study of key history parameters for model R10. The line style denotes results from simulations with low resolution (dotted), fiducial resolution (dashed), and high resolution (solid). Upper left: wind energy losses and retained fraction (see labels), as in left column of Figure 2. Upper right: radial momentum carried in cool/warm gas and hot gas, along with the net momentum injected by the wind, as in middle column of Figure 2. Lower left: effective radius of the wind driven bubble, as in the right column of Figure 2. Lower right: mass history, showing total gas mass on the grid, gas advected out of the domain, and stellar mass, as in Figure 6. Also shown is the mass of stars above Mcut = 10 3 M , which produces winds. many low-mass sink particles form that never grow sufficiently to have winds associated with them, and therefore do not help to disperse high-density collapsing structures or lower-density surrounding material. Rather, star formation continues until the mass in particles that do have winds is sufficient to remove the remaining cloud material. For the fiducial and low-resolution models, at late stages only ∼ 15 − 20% of the stellar mass is in sink particles that do not provide feedback, but this increases to ∼ 30% in the high-resolution model. This high-resolution model is therefore significantly "underpowered" relative to what it should be for a fully sampled IMF, accounting for its enhanced SFE.
We note that the highest resolution simulation transitions more quickly to the "breakout" phase than the fiducial and low-resolution models, based on the earlier dominance of leakage losses over cooling losses in the upper left panel of Figure 7. This could be due to some combination of the higher-resolution model having less ambient gas due to more total star formation (lower-right of Figure 7) or higher contrasts in the density at high resolution, which increases porosity. The more rapid transition to breakout is also reflected in the more rapid growth of R b for the highest resolution model (lower-left panel of Figure 7).
The differences in results as a function of resolution appear to arise from imposition of a mass cut that is far above the minimum sink mass at high resolution, so that many low mass sinks are "left behind" and never accrete enough ambient gas to reach M cut and generate winds. In contrast, for the fiducial and low resolution models, inspection of late-time mass functions show that almost all particles have grown enough to host winds, consistent with the original design of our wind assignment prescription employing M cut . In future work, it will clearly be important to design a wind-assignment prescription that is more robust to changes in the sink mass function with resolution, while still having an astronomically realistic and numerically practical number of feedback sources. One way to do this is by directly sampling from a mass function to assign individual massive stars (with associated feedback) to the sink particles (e.g. Sormani et al. 2017).