Pre-supernova feedback sets the star cluster mass function to a power law and reduces the cluster formation efficiency

The star cluster initial mass function is observed to have an inverse power law exponent around 2, yet there is no consensus on what determines this distribution, and why some variation is observed in different galaxies. Furthermore, the cluster formation efficiency covers a range of values, particularly when considering different environments. These clusters are often used to empirically constrain star formation and as fundamental units for stellar feedback models. Detailed galaxy models must therefore accurately capture the basic properties of observed clusters to be considered predictive. We use hydrodynamical simulations of a dwarf galaxy as a laboratory to study star cluster formation. We test different combinations of stellar feedback mechanisms, including stellar winds, ionizing radiation, and supernovae. Each feedback mechanism affects the cluster formation efficiency and cluster mass function. Increasing the feedback budget by combining the different types of feedback decreases the cluster formation efficiency by reducing the number of massive clusters. Ionizing radiation is found to be especially influential. This effect depends on the timing of feedback initiation, as shown by comparing early and late feedback. Early feedback occurs from ionizing radiation and stellar winds with onset immediately after a massive star is formed. Late feedback occurs when energy injection only starts after the main-sequence lifetime of the most massive SN progenitor, a timing that is further influenced by the choice of the most massive SN progenitor. Late feedback alone results in a broad, flat mass function, approaching a log-normal shape in the complete absence of feedback. Early feedback, on the other hand, produces a power-law cluster mass function with lower formation efficiency, albeit with a steeper slope than that usually observed.


Introduction
Star clusters represent an opportunity to constrain star formation and stellar feedback models used to understand galaxies, a notoriously challenging problem due to its physical complexity and wide dynamical range (Somerville & Davé 2015;Naab & Ostriker 2017).Properties of clusters, such as their mass function, have been a subject for observations for a long time in our own Galaxy (see Lada & Lada 2003), in our galactic neighbors the Magellanic Clouds (e.g., Elmegreen & Efremov 1997;Hunter et al. 2003;Larsen 2009) and M31 (e.g., Vansevičius et al. 2009;Johnson et al. 2015Johnson et al. , 2017)), as well as many other galaxies (e.g., Zhang & Fall 1999;Larsen 2009;Cook et al. 2012;Adamo et al. 2015;Chandar et al. 2014Chandar et al. , 2016;;Adamo et al. 2020b,a;Wainer et al. 2022;Cook et al. 2023), providing strong empirical constraints on models.
Typically, the observed mass function for young clusters follows an inverse power law, possibly with a truncation at higher masses (see Portegies Zwart et al. 2010;Krumholz et al. 2019;Adamo et al. 2020b), although the truncation is debated (Mok et al. 2019).The truncation mass varies between different types of galaxies, while the slope α ranges from −1.5 to −2.5 (see figure 10 in Adamo et al. 2020b).The slope α ≈ −2 is somewhat unsurprising, since this is the slope of a scale-free distribution (Guszejnov et al. 2018), which results from gas fragmentation under gravity and turbulence (Elmegreen & Efremov 1997;Elmegreen 2002).The slope of the distribution is similar for molecular clouds (Colombo et al. 2019;Mok et al. 2020;Rosolowsky et al. 2021), indicating that there could be a direct mapping between the two, though this would imply that the star formation efficiency and timescale are independent of the cloud mass (Clark et al. 2007).Nonetheless, understanding the observed variations remains an unsolved problem, particularly their dependence on the environment.
Several numerical studies address the mechanisms that determine the shape of the cluster mass function, primarily focusing on sub-grid modeling of star formation and the effects of prescribed feedback.By treating star formation in units of star clusters in cosmological simulations, Li et al. (2017) found that 100% star formation efficiency per free-fall time (SFE) resulted in a log-normal cluster mass function peaking around 10 4 M ⊙ , while models with lower SFE resulted in a Schechter-like mass function with α ∼ 2 and a cutoff mass with a dependence on the star formation rate (SFR) or SFR surface density (Li et al. 2018).Interestingly, a log-normal cluster mass distribution is often found in simulations of galaxies undergoing major mergers (Renaud et al. 2015;Maji et al. 2017), although exceptions are also reported (Li et al. 2022).Major mergers tend to increase the amount of dense gas, which accelerates star formation (i.e.shorter depletion time, Renaud et al. 2014;Li et al. 2022;Segovia Otero et al. 2022).Similarly to high SFE, this could provide the environment necessary for massive cluster formation (see, e.g., Kravtsov & Gnedin 2005;Parmentier & Gilmore 2007;Renaud et al. 2015).Increasing the feedback strength by boosting the momentum input from supernovae (SNe) by a factor of five (in order to match global SFR with empirical data) resulted in increased variation between different choices of SFE, with higher SFE resulting in mass function extending toward higher masses, shallower slopes (α ≈ 2.5) at their upper end, and a log-normal distribution emerged from high values of SFE (200%, Li et al. 2018).Furthermore, Li et al. (2018) found that reducing the momentum boost factor to three significantly steepened the slope of the mass function for young clusters.It is clear that feedback impacts the cluster mass function; however, there is no systematic study of the impact of different feedback sources in these studies.
Another open question is what fraction of stars emerge in bound clusters.Although most stars originate in groups and sub-clusters that form in gas concentration generated by turbulent motions inside giant molecular clouds (GMCs, see, e.g., Mac Low & Klessen 2004;Gómez & Vázquez-Semadeni 2014;Krumholz et al. 2014) and assemble into clusters (see, e.g., Wall et al. 2019Wall et al. , 2020;;Cournoyer-Cloutier et al. 2023), the vast majority are rapidly dissolved due to gas ejection (see seminal work by Lada & Lada 2003).The cluster formation efficiency (CFE), quantifying the fraction of newly born stars in clusters has been estimated for a large number of observed galaxies (see, e.g., Goddard et al. 2010;Adamo et al. 2011;Annibali et al. 2011;Adamo et al. 2015;Pasquali et al. 2011;Cook et al. 2012;Lim & Lee 2015;Hollyhead et al. 2016;Chandar et al. 2017;Ginsburg & Kruijssen 2018;Messa et al. 2018;Fensch et al. 2019;Adamo et al. 2020a;Chandar et al. 2023;Cook et al. 2023).These observations are typically interpreted as a function of star formation intensity, as measured by star formation rate density Σ SFR to compare different environments.Note that CFE measurements depend on the cluster age (Adamo et al. 2020b), making the interpretation sensitive to observational tracers (Chandar et al. 2017(Chandar et al. , 2023)).For young clusters (< 10 Myr), CFE is typically found to be high (at least a few tens of percent) at Σ SFR > 0.1 M ⊙ yr −1 kpc −2 (see, e.g., Adamo et al. 2011Adamo et al. , 2020a)), while at lower Σ SFR , CFE scatters from a few percent up to values similar to those measured at high Σ SFR (Goddard et al. 2010;Lim & Lee 2015;Chandar et al. 2017).
Theoretical predictions for the CFE have received similar attention to the mass function.Kruijssen (2012) argued that a scaling between CFE and Σ SFR arises due to a higher bound fraction for more intense star formation at higher gas densities, with saturation in CFE of around 70% at the highest Σ SFR .As in Kruijssen (2012), an increasing CFE with respect to Σ SFR was also found in models by Li et al. (2017Li et al. ( , 2018)).In this case, a choice of higher SFE shifted the normalization of the scaling toward higher CFE, best matching observations at SFE ≳ 50%.Surprisingly, Li et al. (2018) found that boosting feedback only marginally changed the normalization of the CFE-Σ SFR relation.Furthermore, they found that at high-Σ SFR , the saturation of CFE was related to the maximum cluster mass rather than tidal effects.
Clustered star formation has received attention for star formation and stellar feedback modeling of galaxies, both for simulations in isolation and cosmological environments (Naab & Ostriker 2017).For example, numerical models of stellar feedback injection are sensitive to clustered SN explosions (Keller et al. 2014;Agertz & Kravtsov 2015;Smith et al. 2021) and runaway stars originating from clusters (Ceverino & Klypin 2009;Andersson et al. 2020;Steinwandel et al. 2022;Andersson et al. 2023).Furthermore, these models can provide predictions for the initial cluster mass function and CFE if measured before the cluster experiences significant dynamical evolution.In simulations resolving individual stars, Hislop et al. (2021) found a variation in the mass function and CFE for bound clusters for different choices of SFE, raising concerns about the ability of star formation recipes with assumed SFE to accurately model cluster formation.In addition, Hislop et al. (2021) found that photoionizing radiation reduced the CFE by almost 40% for 2% SFE.With a similar model assuming 100% SFE, Lahén et al. (2023) found that the CFE is scattered between 1 and 100% between different outputs in a dwarf galaxy with low Σ SFR .This scatter remains unexplained, and to the best of our knowledge, no systematic investigation of how different sources of feedback affect these properties has been conducted, particularly concerning stellar winds.
To this end, we investigate hydrodynamical simulations of dwarf galaxies and systematically compare how the cluster population is affected by different combinations of stellar feedback.We include core-collapse and type Ia SNe, ionizing radiation using a radiative transfer model, and stellar winds.Furthermore, we test the effect of varying the maximum progenitor mass of core-collapse SNe (Janka 2012).
We describe the model and simulation setup in Section 2 and present our results in Section 3. We discuss these results in a broader context and test model assumptions in Section 4. Finally, Section 5 summarizes our work.

Numerical method and initial conditions
The suite of galaxies presented are simulated using ramses-rt (Rosdahl et al. 2013;Rosdahl & Teyssier 2015), a multi-group radiative transfer extension of the adaptive mesh refinement and N-body code ramses (Teyssier 2002).ramses has an oct-tree hierarchical grid providing adaptive resolution and solves the fluid equations using a second-order unsplit Godunov method with the HLLC approximate Riemann solver (Toro et al. 1994).To avoid spurious oscillations, we apply a MinMod slope limiter to reconstruct the piecewise linear solution for the Godunov solver (see, e.g., Toro 2009).A monatomic, ideal gas equation of state with an adiabatic index γ = 5/3 is applied to close the set of fluid equations.Radiation in ramses-rt is treated via the momentumbased radiative transfer equations with an M1 closure relation for the Eddington tensor (Levermore 1984).Our simulations account for non-equilibrium gas chemistry and cooling, tracking the ionization fractions of H i, H ii, He ii, and He iii by advecting passive scalars on the grid (see Rosdahl et al. 2013;Rosdahl & Teyssier 2015, for details).Furthermore, the code models the non-equilibrium evolution of molecular hydrogen (H 2 ) following Nickerson et al. (2018).All thermochemistry evolves by a semi-implicit method, described in Rosdahl et al. (2013).
Stars and dark matter are tracked by particles, with a mass resolution of 100 M ⊙ and 10 4 M ⊙ , respectively.For gravity, the Poisson equation is solved on the grid using the multi-grid method (Guillet & Teyssier 2011), with particle mapping via the cloud-in-cell method (Hockney & Eastwood 1988).The refine-ment strategy employs a quasi-Lagrangian technique (aiming for eight particles in each cell) and cell refinement at a mass threshold of 800 M ⊙ .The refinement is limited to 14 levels, with a cell size of 3.6 pc on the maximum refinement level.
The initial conditions mimic an isolated dwarf galaxy using the tool MakeDiscGalaxy (Springel 2005).These include a dark matter halo following a NFW-profile (Navarro et al. 1997) with virial mass M vir = 10 10 M ⊙ , spin parameter λ = 0.04, and concentration parameter c = 15.Embedded within this halo is a gas (M g,disc = 7×10 7 M ⊙ ) and stellar (M s,disc = 10 7 M ⊙ ) disc with an exponential radial density profile with scale length 1.1 kpc.The gas has an initial temperature of 10 4 K and vertical distribution set by hydrostatic equilibrium, while stars have a vertical distribution of a Gaussian with a scale height of 0.7 kpc.The initial metallicity of the gas is set to 0.1 Z ⊙ .These initial conditions were studied using a different model in Andersson et al. (2023), and are highly similar to those studied in Smith et al. (2021).

Star formation and stellar feedback
Star formation proceeds by stochastically forming star particles (each with mass 100 M ⊙ ) in all cells with gas density ρ > 10 3 cm −3 and temperature < 10 4 K at each time step.The stellar mass formed is removed from the gas mass of that cell to ensure mass conservation.In each cell fulfilling the star formation criterion, the number of star particles to form is determined by sampling a Poisson distribution parameterized by m ⋆ per 100 M ⊙ , where the cell volume is V, and the timestep ∆t.A Schmidt-like recipe is used for the local volumetric star formation rate density where ϵ ff = 0.1 is the SFE (see, e.g., Grisdale et al. 2019, for discussion regarding this choice), and t ff = (3π/[32Gρ])1/2 is the local gas free-fall time, where G is the gravitational constant.Each star particle inherits the oxygen and iron abundance of its natal gas to track metallicity M Z = 2.09M O + 1.06M Fe (based on the Solar mixture, see Asplund et al. 2009, for details see, e.g., Kim et al. 2014;Agertz et al. 2021).Elements are tracked by advecting passive scalars on the grid and primarily evolve via enrichment linked to stellar feedback.We use the yield data for stellar winds and core-collapse SNe from NuGrid (Pignatari et al. 2016;Ritter et al. 2018), and for SN Ia from Seitenzahl et al. (2013).
The stellar feedback model employs the method described in Agertz et al. (2020) to model stellar winds and SNe (corecollapse and type Ia).Energy injection from radiation is computed by the radiative transfer solver.Stellar winds and SNe inject the appropriate combinations of mass, metals, momentum, and thermal energy.These quantities are injected equally between eight cells (see Agertz et al. 2013;Agertz & Kravtsov 2015, for details).We assume the initial mass function from Kroupa (2001) for all feedback processes.
Fast stellar winds injected by massive (> 8 M ⊙ ) stars have mass-loss rates computed using functions fitted to starburst99 (Leitherer et al. 1999) evolutionary tracks of single-age stellar populations (see Agertz et al. 2013, for a presentation of these fits).These winds act as a source of momentum with a wind velocity of 1000 km s −1 .Our model also considers mass loss from AGB stars as a source of chemical enrichment (injected without any wind velocity).The AGB mass loss rate calculation uses an IMF-averaged initial-final mass relation (Kalirai et al. 2008) with a mass range 0.5-8 M ⊙ as the AGB progenitor stellar mass.
Core-collapse SNe are stochastically sampled as discrete events using the main sequence lifetimes from Raiteri et al. (1996).For the primary suite of simulations, we assume progenitors in the mass range 8-30 M ⊙ but test a higher limit (100 M ⊙ ) in one simulation.The conservative choice is motivated by the uncertainty of the upper limit due to the direct collapse into a black hole (see Section 4 for a detailed discussion).In practice, we assume that the most massive stars (> 30 M ⊙ ) undergo direct collapse into a black hole without the injection of mass, momentum, or energy unless otherwise stated.In the event of a SN, 10 51 erg of energy is injected into the gas surrounding its location.Furthermore, core-collapse SNe are assumed to release momentum equivalent to 12 M ⊙ at a velocity of 3000 km s −1 , calibrated to starburst99.For explosions where the cooling radius of the SN blast-wave is unresolved by more than 6 cells, the fraction of the energy expected to result in the terminal momentum of the explosion is injected as gas momentum rather than thermal energy following Kim & Ostriker (2015).This model assumes a terminal momentum of 2.95 × 10 5 M ⊙ km s −1 .Mass loss for core-collapse SNe is computed from the fit by Woosley & Weaver (1995).Our model accounts for type Ia SNe via the time-delay distribution following Maoz & Graur (2017), assuming a time delay of 38 Myr and rate of 2.6 × 10 −13 yr −1 M −1 ⊙ to normalize the distribution.For type Ia SNe, we assume a mass loss of 1.4 M ⊙ (Chandrasekhar 1931) and the momentum release equivalent to 10 51 erg.
To incorporate radiation feedback, we apply the evolving spectral energy distribution from Bruzual & Charlot (2003) to each stellar particle, employing six separate photon groups (see Table 1 in Agertz et al. 2020).Briefly described, these groups are selected to incorporate scattered infrared photons (non-thermal), direct radiation pressure from optical and far-ultraviolet, and Lyman-Werner radiation for H 2 dissociation.The remaining three groups account for photo-ionizing radiation around the ionizing energies of H i, He i, and He ii.Production of infrared photons from absorption of optical and ultraviolet photons by dust is included following the methods described in Rosdahl & Teyssier (2015) and Kimm et al. (2017), assuming a dust opacity of 10 Z cm 2 g −1 for infrared photons and 1000 Z cm 2 g −1 for photons of higher energy, where Z is the gas metallicity in units of the solar value Z ⊙ = 0.02 (see Agertz et al. 2020, for details).

Results
Throughout this section, we investigate the properties of young star clusters as a population.Unless otherwise stated, we search the outputs for clusters using a friends-of-friends (FOF) algorithm on the positions of stars younger than 25 Myr with a linking length of 4 pc.The FOF groups are refined by removing stars that are not gravitationally bound 1 to the group and considering the groups with at least five remaining members as clusters.We show how our main result depends on different choices of these parameters in Appendix A.
We present simulations executed from the same initial conditions with different types of feedback activated (see Section 2).SNe and radiation; 5) SNe, stellar winds and radiation (see, e.g., Figure 1 for an example of the labels).The main suite of simulations runs for 1000 Myr with an output cadence of 25 Myr.Furthermore, we omit outputs in the first 300 Myr to avoid spurious effects from the initial conditions.
Figure 1 shows the stellar surface density of the final output in each simulation.In this figure, it is apparent that introducing more feedback decreases the number of visible clusters.This result is persistent throughout the simulation, and, as we will show later, this is primarily the case for more massive star clusters.

Cluster formation efficiency
The upper panel of Figure 2 shows the SFR and cluster formation rate (CFR), while the lower panel shows their ratio (CFE = SFR/CFR) to demonstrate how simulations with more feedback result in galaxies with fewer clusters.We measure both these quantities on the same timescale (25 Myr).For simulations including feedback, the CFE remains roughly constant in each model, albeit with a large short-term variation.With radiation feedback, the mean CFE is ∼ 15%, in stark contrast to the simulations with only SNe and SNe+winds, where the mean CFE is ∼ 65% and ∼ 55%, respectively.Without radiation, stellar winds suppress the CFE compared to only SNe.However, radiation dominates over stellar winds, resulting in similar CFE in the two simulations with radiation and with or without stellar winds.In the case of no feedback, the CFE steadily decreases with time, possibly as a result of the rapid decrease in gas fraction (the mass fraction of cold gas decreases by 40%, in contrast to the simulations with feedback where it remains above 80% throughout the entire simulation).The absence of feedback results in a clumpy morphology as expected for the initially high gas fraction (Li et al. 2005;Renaud et al. 2021).Furthermore, a higher gas fraction can result in the formation of more clusters (Li et al. 2005).
The immediacy of radiation and stellar winds, which begin directly after star formation in our model, implies that these pro- cesses immediately regulate local star formation,2 in contrast to the simulation with only SNe where the onset of feedback after star formation is delayed by ∼ 6 Myr (the lifetime of a 30 M ⊙ progenitor).Note that feedback affects gas on kiloparsec scales; however, complexity in the spatial coupling of the different feed- back sources (e.g.Ohlin et al. 2019) makes it challenging to quantify how clusters affect each other.We leave this to future work.
Although stellar winds also act promptly, they impart less energy and momentum than radiation, explaining the lower CFE in the simulation with SNe+winds+radiation compared to the one with only SNe+winds.We return to this notion in Section 3.3 when we investigate the cluster formation time.
Feedback from stars affects the density structure of the starforming gas, as shown in Figure 3.This figure shows the massweighted gas density distribution for cells inside clusters, calculated by the mean distribution of all clusters in a given simulation.The top axis of this figure show the cell-based depletion time of star-forming gas, defined as gas density divided by star formation rate density, here shown as free-fall time per ϵ ff (see Equation 2).Absent radiation, stellar winds, or feedback from more distant SNe, only local thermal and turbulent forces prevent the gas from collapsing.In this case, dense gas is rapidly consumed (10 Myr at ρ ≈ 10 3.5 cm −3 ) before the first local SNe.In simulations with radiation, we note a sharp drop in gas mass with increasing density above ∼ 1000 cm −3 compared to simulations without radiation.This transition is at the density threshold applied for star formation (see Section 2.1).We also find that stellar winds play a smaller role in shifting the distribution toward lower densities.Notably, the highest densities have depletion times on the order of the timescale for core-collapse SNe, although additional gas is likely accreted from the surrounding material.The increase in ρ⋆ at higher gas densities in simulations without radiation likely plays a role similar to changing the SFE, which has been shown to result in higher CFE (see, e.g., Li et al. 2018;Hislop et al. 2021).Furthermore, this can explain why Hislop et al. (2021) found that including ionizing radiation and thus shifting the density distribution to lower peak values resulted in a ∼ 30% decrease in CFE.We emphasize that the details of the simplified star formation models used in galaxy sim-ulations must be sensitive to numerical resolution, in particular, the choice of star formation density threshold and SFE.

Cluster initial mass function
We show the initial cluster mass function of each simulation in Figure 4, derived by summing all identified clusters, which by assumption have ages < 25 Myr.When all feedback sources are active, the mass function scales as a power-law with a slope α ≈ −3, slightly steeper than the canonical mass function (see, e.g., Portegies Zwart et al. 2010;Krumholz et al. 2019).However, the cluster mass function in dwarf galaxies is often observed to have a steeper slope (Adamo et al. 2020b), although not to the extent found here.Removing sources of feedback shifts the power law section of the mass function toward higher masses, leaving a plateau at lower masses if radiation is turned off.If radiation is included, there is no plateau.Note that our simulations produce only a few massive clusters (particularly for low CFE).We show the robustness of the derived mass function with error bars indicating 16 th and 84 th percentiles, each computed using bootstrapping.
Figure 5 shows the mass function normalized to the number of clusters in a given snapshot, i.e., the average number of young (< 25 Myr) clusters found at any given time.While the number of high-mass clusters differs significantly, there are a similar number of low-mass clusters in all simulations: the average number of clusters per output shows a factor two difference at most.We conclude that the lower value for CFE in simulations with more feedback sources is due to an absence of the most massive clusters.
In the previous section, we noted that the early onset of local star formation regulation from pre-SN feedback affects the formation of star clusters.To test how the mass function is affected by the limitation of the growth time rather than the growth rate, we employ the cluster-finding algorithm with different age limits to the simulation without any feedback.The result is in Figure 6, indicating a truncation at high masses similar to that set by introducing feedback.We show a similar figure for all other simulations in Appendix A. The distribution is steeper at lower mass for age limits < 10 Myr and is most similar to the slopes of the complete feedback model between 2-5 Myr.We find a similar effect in simulations with delayed feedback (only SNe), shown in Appendix A. Since the slope at the upper end of the mass function is similar for all age limits, each cluster in a population must grow exponentially.The rate at which gas is converted into stars depends on the gas mass (see Equation 2); therefore, exponential growth is expected.There is a complexity in this behavior, and we discuss this further in the next Section.

Cluster formation time
Figure 7 shows the time between the formation of the first star in a cluster t 0 and the final time t 100% , and the time it takes each cluster to go from 10% to 90% of its final cluster mass.Note that t 100% − t 0 is more sensitive to outliers and resolution since these points are often the result of a few star particles that formed at a different time compared to the bulk of the cluster members.Without feedback, many clusters grow in mass over long periods, with many cluster masses limited only by the maximum age of 25 Myr for stars considered by our cluster-finding algorithm.The clusters sitting at the upper limit in formation time confirm that the overall mass function is not converged with the age limits applied in this case, due to continued cluster growth (see also  Figure 6).This is not the case for the remaining simulations, and as shown in Appendix A, the mass function has converged for all age limits ≥ 25 Myr.In Figure 7 the lower limit to cluster age at each mass is roughly linearly proportional to the log of the cluster mass.Furthermore, in the simulations with SNe and SNe+winds, all clusters scatter around such a trend.We note that the linear trend limiting the lowest mass cluster has a shallower slope in simulations with radiation, indicating that more rapid growth is necessary to grow to larger masses in these simulations.When including radiation, no single trend is found between the formation time and the cluster mass, indicating a complexity in what determines the growth rate in this case.Similarly, Li et al. (2018) found that the cluster growth rate cannot be described by a simple power law in agreement with predictions for the collapse of turbulent clouds (Murray & Chang 2015;Murray et al. 2017).The top row shows the time between when the first star was born t 0 and when the last star formed t 100% , while the bottom row shows the time between when 10 and 90% of the total cluster mass has assembled.The horizontal dashed line denotes the main sequence lifetime of a 30 M ⊙ star, which is the earliest time a star can deposit energy and momentum in a SN in our simulations.
Many clusters grow past the first possible time for a SN explosion.A stochastic sampling of the IMF for SN explosions implies that not all clusters contain a 30 M ⊙ progenitor, with less massive clusters being more sensitive to this stochasticity.When cluster growth is unlimited, low-mass clusters can grow quickly, thereby eliminating the possibility of clusters with low mass and long formation time in Figure 7.In the case when pre-SNe feedback is active, cluster growth is suppressed by feedback from the time when the first star forms, without any regard for stochasticity. 3This results in low-mass clusters with long formation times, highlighting that pre-SNe feedback is not always strong enough to disrupt the cloud and stop star formation entirely.Many of these low-mass clusters form stars past the time of the first possible SN.

Cluster formation efficiency
We find that feedback introduced immediately after star formation (e.g., stellar winds and radiation) limits the formation of massive star clusters and significantly reduces the cluster formation efficiency.In Figure 8, we show the cluster formation efficiency of each output as a function of the star formation rate surface density Σ SFR .There is notable scatter in the CFE for similar values in Σ SFR in all simulations, as highlighted by the mean value and standard deviation shown with the error bars.Lahén et al. (2023) found a similar scatter in a simulation of a comparable galaxy, albeit with the majority of points sitting at higher efficiency.Note that Lahén et al. (2023) defined clusters with a lower minimum mass (100 M ⊙ ) than us, which by definition may lead to a higher efficiency (our choice does not consider clusters with mass < 500 M ⊙ ).Furthermore, there is a discrepancy between the cluster ages (10 Myr versus our 25 Myr), however, decreasing the star age limit to 10 Myr would not affect our results in a measurable way (see Appendix A).We find a wide range of CFE in multiple snapshots at similar Σ SFR in our simulations with feedback.Furthermore, CFE has roughly the same mean value over almost three orders of magnitude in Σ SFR in SNe and SNe+wind, although the scatter increase towards lower Σ SFR .In the simulation without feedback, we do find a clear trend between different outputs in this space.In this case, the trend has a clear evolution in time (later outputs toward the bottom left) and coincides with a reduced gas fraction over time (our simulations lack significant gas inflow), although causality remains unconfirmed (other factors likely play a role, e.g., morphological evolution).
The red crosses in Figure 8 show observational estimates from the literature (Cook et al. 2023, see also figure caption), and indicate an absence of low CFE at higher Σ SFR .Note that this is debated in the literature, in particular concerning the age limit of the observed clusters (see, e.g., Chandar et al. 2017Chandar et al. , 2023)).Furthermore, these estimates rely on different methods, although all values shown in Figure 8 consider only clusters in the age range 1-10 Myr (see Adamo et al. 2020b;Cook et al. 2023, for discussions about the different methods employed).In our simulations, changing the age limit to 10 Myr does not affect the CFE when including feedback (see Appendix A) since the vast majority of star clusters have stopped growing by this time.Since we only consider one dwarf galaxy in our simulations, we cannot address the question of high CFE and Σ SFR at higher Σ SFR .
Extrapolating the effect of feedback on star cluster formation to more massive galaxies (and thereby higher Σ SFR ) from our models is not trivial, since many properties to which the cluster formation process is sensitive change.For example, more massive systems have higher metal content, which affects gas dynamics through cooling processes (Ferland et al. 1998;Rosen & Bregman 1995) and through the ability for feedback to couple dynamically to the gas, both for radiation (via dust attenuation Rosdahl & Teyssier 2015), and strength of stellar winds (Vink 2011).This is complicated further by other properties that vary between galaxies, such as gas fraction and merger events, which have an impact on the dynamics of the galaxy and the population of clusters (Li et al. 2005;Renaud et al. 2019Renaud et al. , 2021)).

Maximum cluster mass
From Figure 4, it is obvious from our models that increasing the number of active feedback mechanisms reduces the maximum cluster mass.An upper mass limit or a Schechter-like function with an exponential decline above a mass threshold is often reported in the literature (see reviews by Portegies Zwart et al. 2010;Krumholz et al. 2019;Adamo et al. 2020b).In the case of a truncation mass, several works report a value around 10 5 M ⊙ for spiral galaxies (Gieles et al. 2006;Bastian et al. 2012;Adamo et al. 2015;Messa et al. 2018), although notable exceptions are galaxies experiencing mergers (∼ 10 7 M ⊙ , Adamo et al. 2020a) and the M31 galaxy (≲ 10 4 M ⊙ , Johnson et al. 2017).On the other hand, Mok et al. (2019) found many cases where high-mass truncation was not deemed significant when accounting for completeness and errors in the cluster mass estimates, instead reporting a pure power law for several local galaxies.The qualitative shape of the mass functions are different between our model; SNe+winds+radiation is well described by a power law, while no feedback is neither described by a power law nor a Schechter function but instead approaches a log-normal.Extracting a cluster mass function for dwarf galaxies accurately enough to determine the presence of a truncation mass is highly challenging due to the low number of clusters in any one galaxy, and likely requires stacking observations of several galaxies.Nonetheless, we emphasize the importance of considering a range of feedback mechanisms (stellar winds, radiation, and SNe) when formulating a general theoretical framework for star cluster formation in different galaxies.

Maximum SN progenitor mass
In the case of SNe, we limit the injection of mass, momentum, and energy to stars within the mass range 8-30 M ⊙ , assuming direct collapse for more massive stars.The direct collapse to a black hole for M > 30 M ⊙ is predicted by many theoretical models of the core-collapse process (O'Connor & Ott 2011;Ertl et al. 2016;Sukhbold et al. 2016), although there is notable uncertainty, particularly with regards to explosion mechanism (Janka 2012).Our choice of 30 M ⊙ translates to a minimum delay of ∼ 6 Myr between the formation of a star particle and the first release of energy from SNe.Furthermore, there are theoretical predictions that the most massive stars tend to form late during the build-up of the IMF (Haugbølle et al. 2018;Grudić et al. 2023).To investigate how our choice affects our conclusions, we present the mass function for a re-simulated SNe model with an SN mass threshold at 100 M ⊙ (resulting in a 3.3 Myr main sequence lifetime) in Figure 9.In this simulation, we limit the explosions to energy, momentum, and mass (i.e., we do not inject enriched material to maintain a similar gas metallicity).The new simulation ran for 500 Myr.Therefore, we limit the comparison to this time for both simulations presented in Figure 9. Similarly to the simulation with radiation and winds, earlier SNe from the more massive stars reduces the number of the most massive clusters, resulting in a power law-like mass function for the entire range of masses.Note that changing the maximum mass for SN progenitors not only affects the timing of feedback onset but also the energy budget.The sensitivity to choices with regard to uncertainties in SN modeling highlights that uncertainty remains, even for state-of-the-art galaxy models.

Limitations of our model
Our results are derived from numerical simulations of a single galaxy in an isolated environment.This restricts our conclusions to the properties of this specific galaxy, e.g., one initial gas fraction, gravitational background, and galactic dynamical timescale.Furthermore, this limits our range of the turbulence spectrum, Mach numbers, and cloud-to-inter-cloud density contrasts, which are suspected to be key ingredients in the assembly of the star cluster forming sites (e.g.Renaud et al. 2019;Li et al. 2022).Furthermore, our simulations rely on parameter choices that are debated in the literature (e.g., the star formation efficiency per free-fall time).Therefore, we emphasize that our results should be interpreted as relative changes between our simulated models.
Although our model incorporates the most common mechanisms for stellar feedback, it still lacks several processes, e.g., cosmic ray physics, massive rapidly rotating stars (Wolf-Rayet stars), stellar multiples, exotic SNe (e.g., pair-instability SNe), and protostellar jets.Of particular interest are assumptions about the stellar feedback mechanisms that affect the local gas dynamics soon after the star forms (see also Section 4.3).For the star particles, our model assumes a Kroupa IMF (Kroupa 2001) with an upper mass limit of 100 M ⊙ for calibrating the input for the radiative transfer, calculating the mass loss from stellar winds, and stochastically sampling progenitors for SNe events.
Although our model stochastically samples SNe as discrete events, winds, and radiation are calibrated to IMF averages.This can affect both galactic scale properties (Smith 2021), as well as affect the final cluster properties (Lewis et al. 2023).An ideal tool to study these effects are star-by-star models which have recently become common (see, e.g., Hu et al. 2016;Emerick et al. 2018;Andersson et al. 2020;Gutcke et al. 2021;Hirai et al. 2021;Hislop et al. 2021;Andersson et al. 2023;Lahén et al. 2023), and should be applied to this problem in the near future (see Hislop et al. 2021;Lahén et al. 2023, specifically).
Our simulations are initialized from the same initial conditions; however, different feedback mechanisms are active from the start.Since these mechanisms are of a fundamentally distinct nature, differences could accumulate over time and change the environment where the clusters form.As an alternative, the simulations can be run with similar models to some point after which the different feedback mechanisms are applied.This approach is better suited for studying the role different feedback mechanisms play in the cluster formation process in an identical environment.However, a challenge with this approach is the inconsistency in the ISM structure set by a different stellar feedback cycle and thereby different turbulence.Performing these types of tests is beyond the scope of this paper, but see forthcoming work by Ringdahl et al. (in preparation).The cluster formation process in controlled environments can also be tested simulation isolating the environment from initialization on smaller scales (see, e.g., Wall et al. 2019Wall et al. , 2020;;Cournoyer-Cloutier et al. 2021, 2023;Lewis et al. 2023).
Another limitation of our model is unresolved dynamical interactions between particles inside the star clusters because of the lack of star-by-star calculations and systems of multiple stars (e.g., binaries).Furthermore, the variations in the tidal field around the clusters imposed by gas motion are only marginally resolved ( 4 pc).Such variations can affect the cluster evolution (Renaud 2018(Renaud , 2020)).These are important for the dynamical properties of clusters and their evolution, and as such, we limit our study to young clusters and refrain from studying properties where the internal dynamics play a role (e.g., size and velocity dispersion).Due to the computational demand of direct N-body calculations, studying the dynamical evolution of star clusters is typically limited to individual clusters in gas-free simulations (see, e.g., Renaud et al. 2011;Wang et al. 2016), or including gas until it is expelled from an individual cluster (see, e.g., Wall et al. 2019Wall et al. , 2020;;Ballone et al. 2021;Cournoyer-Cloutier et al. 2021;Grudić et al. 2021;Guszejnov et al. 2022;Cournoyer-Cloutier et al. 2023), although preliminary work on the topic was done by Mamikonyan et al. (2017).

Summary
Our work studies how stellar feedback affects the population of young star clusters.To do this, we test different combinations of feedback in hydrodynamical simulations of isolated dwarf galaxies with an initial stellar mass of 10 7 M ⊙ .For the first time, we systematically compare combinations of different feedback sources to investigate the cluster initial mass function and CFE.We find the following: 1. Introducing radiation changes the cluster formation efficiency dramatically (see Figures 2 and 8).The total stellar mass in young clusters relative to the total stellar mass formed is ∼ 0.6 in simulations with only SNe or SNe+winds, while for simulations with additional radiation feedback, the same fraction is ∼ 0.15.2. We find that systematically removing feedback mechanisms transforms the upper end of the initial cluster mass function (see Figure 4).In the simulation with SNe+winds+radiation, the initial cluster mass function follows a power law with a slope around −3.As mechanisms resulting in feedback are removed, the power law shifts to larger masses with a plateau at lower masses.In the absence of feedback, the function resembles a log normal.3. The total number of low mass clusters (M < ∼ 10 3 M ⊙ ) is similar in all simulations (see Figure 5).Therefore, the cluster formation efficiency is determined by the abundance of more massive clusters.4. The timing of feedback initiation is crucial to the formation of star clusters.Feedback with early onset after star formation, such as winds and radiation, prevents cluster growth and explains the lack of massive clusters as compared to simulations without feedback or only SNe. Figure 7 show how the cluster formation time varies between the models with different feedback mechanisms active.5. We employ two tests to check how the cluster formation timescale affects the cluster mass function: (a) by placing progressively lower limits on the age of stars considered to be cluster members, we find that the initial cluster mass function in the simulation without feedback becomes progressively more similar to the simulations with earlier onset of feedback (see Figure 6); (b) by letting massive (up to 100 M ⊙ ) stars deposit mass, momentum, and energy as SNe (relaxing the assumption of direct collapse to black hole for stars with M > 30 M ⊙ ), the initial cluster mass function show only a power law slope (albeit with shallower slope compared to simulation with radiation) even in a simulation with only SNe (see Figure 9).Note that these more massive SN progenitors imply earlier explosions due to the shorter time spent on the main sequence.
In conclusion, the formation of star clusters is sensitive to the timing for the onset of stellar feedback.Therefore, feedback that starts more quickly or slowly results in different CFE and initial cluster mass functions.The immediate regulation of local star formation provided by radiation and stellar winds results in a power-law mass function for young star clusters.Without these sources of pre-SN feedback, star clusters rapidly grow in mass, shifting the entire mass function toward higher masses (with a plateau at lower masses) and high cluster formation efficiencies.

Fig. 1 .
Fig.1.Face on view of the stellar surface density of all simulations.Each panel shows the last snapshot of the simulation (t = 1 Gyr).There is a notable increase in the number of visible star clusters in simulations with less feedback, compare, e.g., SNe and SNe+wind+radiation.

Fig. 2 .
Fig. 2. Top: SFR (filled line) and CFR (dotted line) as a function of time for all simulations.Both these quantities are measured on the same timescale (25 Myr).Bottom: The CFE as a function of time for all simulations.Dashed horizontal lines indicate the fractions of cluster mass to total stellar mass formed after 300 Myr (outside gray shaded region).

Fig. 3 .
Fig.3.Average mass distribution of gas densities inside clusters.The distribution is the mean of the mass-weighted histograms of gas densities inside clusters.The inside of a cluster is the spherical region enclosing all cluster members.The shaded region indicates the standard deviation around the mean value.The dashed vertical line indicates the density threshold we use for star formation.The top horizontal axis shows the ratio between the free-fall time and ϵ ff = 0.1 (see Section 2.1 for more details).

Fig. 4 .
Fig. 4. Cluster initial mass function, normalized to the total number of clusters (defined as objects with ages < 25 Myr) found in each simulation after 300 Myr.The error bars indicate the 16 th and 84 th percentiles computed using bootstrapping.The dashed line shows the slope of scale-free formation models, often associated with the initial cluster mass function (see text for details).

Fig. 5 .
Fig. 5. Histogram showing the cluster mass function, normalized such that the y-axis gives the number of clusters expected when looking at a random simulation snapshot.Notably, all simulations have roughly the same number of low-mass clusters.This figure is identical to Figure 4 except each line is normalized differently.

Fig. 6 .
Fig.6.Cluster mass function for the simulation without feedback applying different age limits to the cluster finding algorithm.The age limit serves as a way to reduce the cluster formation time in a simulation where no stellar feedback is responsible for such limitations.Notably, this shifts the mass function to resemble simulations with feedback.

Fig. 7 .
Fig.7.The formation time of individual clusters as a function of their mass for each simulation.The top row shows the time between when the first star was born t 0 and when the last star formed t 100% , while the bottom row shows the time between when 10 and 90% of the total cluster mass has assembled.The horizontal dashed line denotes the main sequence lifetime of a 30 M ⊙ star, which is the earliest time a star can deposit energy and momentum in a SN in our simulations.

Fig. 9 .
Fig.9.Cluster mass function (cf.Figure4) for the simulation with the fiducial SN feedback model compared to the mass function for a simulation where stars with mass up to 100 M ⊙ can explode as core-collapse SNe (see text for details).Allowing more massive stars to explode as SNe implies an earlier onset of SN feedback due to the shorter time spent on the main sequence.Because the simulation with massive SN progenitors only progressed 500 Myr, we limit the results to clusters formed between 300 and 500 Myr.Note that the earlier onset of SNe results in a mass function without an obvious plateau at lower masses.