Star clusters forming in a low metallicity starburst – rapid self-enrichment by (very) massive stars

Stellar winds of massive ( ≳ 9 M ⊙ ) and very massive ( ≳ 100 M ⊙ ) stars may play an important role in the metal-enrichment during the formation of star clusters. With novel high-resolution hydrodynamical griffin-project simulations, we investigate the rapid recycling of stellar wind-material during the formation of massive star clusters up to 𝑀 cluster ∼ 2 × 10 5 M ⊙ in a low-metallicity dwarf galaxy starburst. The simulation realises new stars from a stellar initial mass function (IMF) between 0 . 08 M ⊙ and ∼ 400 M ⊙ and follows stellar winds, radiation and supernova-feedback of single massive stars with evolution tracks. Star clusters form on timescales less than ∼ 5 Myr, and their supernova-material is very inefficiently recycled. Stellar wind-material, however, is trapped in massive clusters resulting in the formation of stars self-enriched in Na, Al, and N within only a few Myr. Wind-enriched (second population) stars can be centrally concentrated in the most massive clusters ( ≳ 10 4 M ⊙ ) and the locked wind-material increases approximately as 𝑀 2cluster . These trends resemble the characteristics of observed second population stars in globular clusters. We fit scaling relations to the log-normal distributed wind-mass fractions and extrapolate to possible globular cluster progenitors of 𝑀 cluster = 10 7 M ⊙ to investigate whether a dominant second population could form. This can only happen if the IMF is well sampled, single massive stars produce at least a factor of a few more enriched winds e.g. through a top-heavy IMF, and a significant fraction of the first population (unenriched) stars is lost during cluster evolution.


INTRODUCTION
The majority of present-day globular clusters (GCs) have been observed to host multiple populations (MPs) of stars, identified through variations in the stellar light-element composition (Gratton et al. 2001;Carretta et al. 2009b;Piotto et al. 2015;Pancino et al. 2017;Marino et al. 2019; see e.g.Renzini et al. 2015;Charbonnel 2016;Bastian & Lardo 2018;Gratton et al. 2019 for comprehensive reviews).The stars in GCs can be divided into normal, field-like stars called the first population (1P) and to anomalous, light-element peculiar second population (2P).The presence of MPs is often surveyed by contrasting certain abundance ratios against each other: for instance, N, Na, Al and Si are typically enhanced when C, O and Mg are depleted.In extreme cases, the abundances of Mg and K are anticorrelated (Cohen et al. 2011;Mucciarelli et al. 2015).Heavy element abundances, on the other hand, show fairly uniform distributions and only small variations across GC stars (Marino et al. 2019), with ironcomplex GCs being the exception (Marino et al. 2015;Johnson et al. 2015).
The spread in the abundance ratios and the fraction of stars in the 2P of GCs correlate with GC mass and anticorrelate with metallicity ★ E-mail: nlahen@mpa-garching.mpg.de(Carretta et al. 2009b(Carretta et al. ,a, 2010;;Milone et al. 2017aMilone et al. , 2020)).Tentative correlation with cluster age has been reported based on the fact that almost all ancient massive clusters host a significant number of 2P stars while young, relatively massive clusters do not (Cabrera-Ziri et al. 2016;Lardo et al. 2017;Martocchia et al. 2017Martocchia et al. , 2018a;;Milone et al. 2020).The radial distribution of the 2P stars within GCs is often centrally concentrated (Kravtsov et al. 2011;Lardo et al. 2011;Simioni et al. 2016;Dalessandro et al. 2019) while some GCs have a centrally concentrated 1P or a radially uniform mix of 1P and 2P stars (Dalessandro et al. 2014;Larsen et al. 2015;Leitinger et al. 2023).The absence of a central concentration in 2P might either be the result of dynamical mixing (Vesperini et al. 2013) or there being no primordial segregation between 1P and 2P.Star clusters also lose mass through gas expulsion and while orbiting their host galaxies (Decressin et al. 2010;Lamers et al. 2010;Khalaj & Baumgardt 2015;Baumgardt et al. 2019;Wang et al. 2016), with complex impact on the MPs depending e.g. on the galactic orbit and the initial radial distribution of the stars (Krause et al. 2016;Milone et al. 2020).
The variety of models that can at least partly explain the abundance and scaling relations of MPs is arguably among the most diverse sets of proposed solutions to an open problem in modern astrophysics (Bastian & Lardo 2018;Gratton et al. 2019).A comparison of the enhanced and depleted elements with the products of nucleosynthesis at increasingly high temperatures (Arnould et al. 1999) indicates such abundance variations may be related to nuclear burning within intermediate mass or massive stars (Denisenkov & Denisenkova 1989;Langer et al. 1993).As the nuclear burning products are released into the surroundings of a massive star, they pollute the local environment.If the ejection process is not energetic enough to expel the material out of a star-forming region, the material may be locked into newly forming stars.This so-called self-enrichment will lead to chemical variations within the population of stars that formed spatially and temporally correlated.Self-enrichment by the first formed stellar population is especially favoured as the origin of the ubiquitous MPs because the peculiar abundance ratios are observed across a range of present-day stellar masses and evolutionary stages (Cannon et al. 1998;Gratton et al. 2001;Yong et al. 2003), which is hard to reconcile only with internal evolution of all those stars.
From a stellar nucleosynthesis point of view, to produce the Na-O anticorrelation observed in the majority of GCs, the Ne-Na chain requires a temperature of tens of MK (Langer et al. 1993;Prantzos et al. 2007).The Al-Mg anticorrelation observed in a limited set of GCs would be obtained through the Mg-Al chain, operating efficiently at temperatures above 70 MK (Prantzos et al. 2007).The less frequently observed Si and K enhancement would require more than 80 MK and 150 MK, respectively (Prantzos et al. 2017).The most popular polluters suggested today include intermediate mass asymptotic giant branch (AGB) stars (Cottrell & Da Costa 1981;Ventura et al. 2001;D'Ercole et al. 2010), fast-rotating massive stars (Norris 2004;Maeder & Meynet 2006;Prantzos & Charbonnel 2006;Decressin et al. 2007), massive binary stars (de Mink et al. 2009;Sills & Glebbeek 2010;Bastian et al. 2013), very massive (Vink 2018) and supergiant stars (Szécsi et al. 2018;Szécsi & Wünsch 2019), supermassive stars (Denissenkov & Hartwick 2014;Gieles et al. 2018;Charbonnel et al. 2023) and stellar mergers (Wang et al. 2020).Tidal disruption events (Kochanek 2016) have also been suggested as a possible source for strong N-enhancements observed in compact high-redshift objects (e.g.Bunker et al. 2023) that may be progenitors to present-day GCs.Many of the proposed scenarios require dilution of the polluting material with the unenriched (i.e.non-polluted) interstellar medium (ISM) either to produce the proper abundance correlations or to alleviate the so-called mass budget problem (Prantzos et al. 2007).The latter describes the discrepancy between the significant, often dominant, fraction of 2P stars per cluster and the fact that stellar populations typically return only a small fraction of their mass for recycling.As pointed out by Renzini et al. (2015) and Bastian & Lardo (2018) among others, none of the proposed scenarios have been able to alone explain all aspects of the ubiquitous MP phenomenon; it may turn out that multiple sources and evolutionary scenarios need to be invoked together (Sills & Glebbeek 2010;Carretta et al. 2012;Carretta 2014;Bastian et al. 2015;Johnson et al. 2019).
Numerical simulations have been used to address a variety of the scenarios that lead to the formation of chemically anomalous stars within GCs.The formation of the 2P polluted by AGB winds in hydrodynamical simulations (D'Ercole et al. 2008;Bekki 2011Bekki , 2019;;Calura et al. 2019;McKenzie & Bekki 2021) is among the most commonly studied scenarios.Semianalytic hydrodynamical models of self-enrichment and gas retainment through stellar winds of lowmetallicity massive stars have been investigated in Szécsi & Wünsch (2019).Such detailed enrichment simulations typically neglect modelling the formation of the 1P and instead impose the polluting component as a single age pre-existing population evolving e.g.within a gaseous background.Self-pollution by massive OB-star winds within a star-forming molecular cloud has been self-consistently modelled by Lancaster et al. (2021c) using a detailed stellar wind prescription but at solar metallicity, incompatible with GC-forming environments, and with no attention to chemical enrichment.Howard et al. (2019) likewise followed the star cluster formation process within giant molecular clouds using sink particles, but painted the stellar enrichment over the cluster stars only in post-processing.Monte Carlo simulations (Vesperini et al. 2021;Hypki et al. 2022) and N-body simulations have been used to probe the dynamical evolution of MPs (Hénault-Brunet et al. 2015;Hong et al. 2017;Tiongco et al. 2019;Lacchin et al. 2024) but again rely on assumptions of the initial distribution of the 1P and 2P stars.So far, the formation of star clusters in the GC mass range, the development of MPs, and their complex variety of correlated abundance ratios have not been investigated in hydrodynamical simulations that realise single stars in an evolving galactic star forming environment.
Here we assess the feasibility of the coeval cluster formation and self-enrichment by following star formation and stellar feedback from galactic to sub-parsec scales throughout a dwarf galaxy starburst.We simulate how the interaction of two dwarf galaxies composed of gas, stars and dark matter leads to compression of dense, high-pressure gas in a starburst phase where the star formation rate (SFR) exceeds the quiescent phase by two orders of magnitude.We follow the formation and evolution of massive star clusters with more than 10 5 stars by applying pre-computed and tabulated stellar-evolutionary models to stellar particles that represent single stars.The stars release chemically enriched material through stellar winds and supernovae, which we trace element-by-element and channel-by-channel throughout the galactic ISM.This allows us to follow when and how a chemically enriched second population of stars forms, including their exact chemical composition.We then compare our theoretical results to observations of present-day GCs, which may have formed in such conditions when the Universe was still young (Vanzella et al. 2023;Adamo et al. 2024).Because MPs are most ubiquitous in GCs with sub-solar metallicity, we concentrate on stellar enrichment processes in a galactic environment with a low metallicity of only 1% of the solar value.
For hydrodynamical simulations we employ the sphgal code (Hu et al. 2014(Hu et al. , 2016(Hu et al. , 2017;;Hu 2019), supplemented with a fully realised stellar initial mass function (IMF; Lahén et al. 2023) and evolutionary tracks for massive stars (Szécsi et al. 2022).The simulations are part of the Galaxy Realizations Including Feedback From INdividual massive stars1 (griffin) project, geared to address outstanding challenges in galaxy formation and evolution (Naab & Ostriker 2017).The methodology has been successfully used to study the formation of stars, star clusters and globular clusters (Lahén et al. 2019(Lahén et al. , 2020a,b;,b;Hislop et al. 2022;Lahén et al. 2023;Steinwandel et al. 2023) and their observational properties (Lahén et al. 2022), as well as the CII-emission and H 2 distribution in the ISM (Bisbas et al. 2022;Szakacs et al. 2022) in low-metallicity dwarf galaxies.
The article is structured as follows.Section 2 introduces the simulation methodology and initial conditions of the hydrodynamical simulations.In Section 3 we discuss the results, starting with general star formation properties of the dwarf galaxy starburst and continuing to detailed investigation of chemical enrichment within the young massive star clusters formed during the starburst.We discuss the implications of our results for the formation of self-enriched MPs in GCs in Section 4. Conclusions are given in Section 5.In the following, we will use the terms enriched or second population (2P) for the chemically peculiar GC stars, and unenriched or first population (1P) for the GC stars that have a chemical composition equivalent to the field stars.

SIMULATIONS
For simulations we use the sphgal code (Hu et al. 2014(Hu et al. , 2016(Hu et al. , 2017;;Hu 2019) that has been modified from gadget-3 (Springel 2005) to introduce improvements to the performance of the smoothed particle hydrodynamics (SPH) method.Modifications to star formation and stellar feedback routines to e.g.account for single massive stars have been outlined in Hu et al. (2014Hu et al. ( , 2016Hu et al. ( , 2017)), Hu (2019), Lahén et al. (2020a) and Lahén et al. (2023).Here we give a brief description of the star formation and stellar feedback models most relevant for the current study and refer the reader for further details to the references listed above.For SPH we adopt 100 neighbours.

Star formation
Non-equilibrium cooling and heating processes in the lowtemperature ISM (< 3 × 10 4 K) are accounted for using a chemical network (Hu et al. 2016) that includes six chemical species (H 2 , H + , H, CO, C + , O) and free electrons.At temperatures exceeding > 3 × 10 4 K, we use the equilibrium metal-dependent cooling rates of Wiersma et al. (2009).
Star formation is implemented as in Lahén et al. (2023) with a few minor modifications.We employ a Jeans mass dependent star formation threshold and sample the gas mass into individual stars after a gravitational collapse time (local dynamical timescale) has elapsed.The local Jeans mass is computed as where  s is the sound speed,  is the gravitational constant, and  is the density of the gas particle.Gas particles for which the Jeans mass crosses below a threshold of half the SPH-kernel mass, 0.5 kernel = 200 M ⊙ , are turned into reservoir particles and decoupled from hydrodynamics.The reservoir particles interact with their surroundings only through gravity, while they remain inert for a duration given by the local dynamical time scale  dyn = (4 ) −1/2 measured for the gas at the time when the particles crossed the star formation threshold.During this inert time, the particles are also decoupled from any nearby stellar feedback: they do not receive metalenriched material once the gravitational collapse has been triggered.After the dynamical time of a reservoir particle has passed, we sample the mass of the particle immediately into single stars randomly according to the Kroupa (2001) IMF between 0.08 M ⊙ and 500 M ⊙ .If the sampled mass exceeds the particle mass by more than 0.08 M ⊙ , we search for neighbouring reservoir particles to account for the missing mass.Compared to Lahén et al. (2023) where we used a fixed search radius of 1 pc, we now use the local Jeans length at the time of conversion as the search radius, given by The Jeans length provides a natural proxy for the size of the region that undergoes gravitational collapse.The radius is allowed to vary between 0.3 pc and 5 pc, limited by the size of the gravitational softening kernel and the size of typical molecular clouds.In case there is not enough mass to satisfy the sampled stellar mass, the last value is discarded and a new random mass is drawn.The sampling and neighbour search is continued until the mass of the particle is successfully sampled, one reservoir particle at a time.As a result, the simulation realises massive (> 9 M ⊙ ) and very massive (> 100 M ⊙ ) stars only in conditions where the local mass reservoir is large enough to satisfy the stellar mass through strict conservation of mass.As shown in Lahén et al. (2023), this results in a realistic relation between the mass of star clusters and the most massive stars they can host, while star formation in the galactic field can only result in relatively low-mass stars.

Stellar feedback
Mass is conserved regardless of particle type throughout star formation and stellar feedback.This is especially important in the case of the chemical elements, as we are interested in the propagation of the elements originating from different enrichment channels.The simulations trace the mass of 13 elements: H, He, N, C, O, Si, Al, Na, Mg, Fe, S, Ca and Ne.In addition to the uniform initial metallicity of 0.016 Z ⊙ , we follow the chemical enrichment through stellar winds of massive stars and AGB stars, as well as core-collapse and pair-instability supernovae (SNe).In the current study we investigate stellar enrichment during the formation of star clusters in a dwarf galaxy starburst mainly over a time scale of a few Myr, therefore we focus on the elements originating from stellar winds of massive stars and supernovae.All stars in the simulation are considered individually while computing the interstellar radiation field through far-ultraviolet radiation (FUV,.The radiation field at the location of each gas particle is attenuated using the treecol algorithm (Clark et al. 2012) along 12 incoming directions, divided evenly around the particle using healpix (Górski & Hivon 2011).In addition, massive stars create H II regions through photoionising radiation propagated using a Strömgren type approximation that iteratively searches for ionization equilibrium.Details of how radiation is implemented in the simulation can be found in Hu et al. (2016) and Hu et al. (2017).
We follow the evolution of single massive and very massive stars until the end of the core helium-burning phase using the Bonn Optimized Stellar Tracks (BoOST, Szécsi et al. 2022).We use the lowmetallicity ( = 0.016 Z ⊙ ) IZw18-grid that contains 1856 slowly rotating stellar models with 608 timestamps each, logarithmically spaced between 9 M ⊙ and 498.4 M ⊙ .As described in (Lahén et al. 2023), each massive star is given an FUV luminosity, a photoionising photon rate, and stellar wind properties such as the outflow rate, chemical composition and velocity, all of which evolve throughout the stellar lifetimes.The choice of applying the BoOST library of massive-stellar wind yields was based on the fact that this library provides a comprehensive list of elements, including Na, Mg and Al, that are under wide-spread scrutiny in GC-studies.Fig. 1 gives examples of how some key abundance ratios on the stellar surface, and hence in the stellar wind, evolve during the stellar lifetime (see Fig. 2 in Lahén et al. 2023 for more examples).The stellar winds of up to ∼ 10 −3 M ⊙ yr −1 are injected radially outward assuming a momentum conserving interaction between the stellar wind and the ISM as described in Lahén et al. (2023, see their Section 2.2).
Once a massive star in our simulation reaches the end of its lifetime (according to Groh et al. 2019), it either dies as a corecollapse SN, a pair-instability SN (PISN) explosion, or in collapse to a black hole.We assume that stars with initial masses between 8 M ⊙ and 40 M ⊙ explode by releasing 10 51 erg of thermal energy 3 and metal-enriched material interpolated using the metallicity dependent yield tables of Chieffi & Limongi (2004).PISNe occur for stars with helium-core masses between 65-133 M ⊙ (Heger & Woosley 2002).In our implementation, we check the helium-core mass of the BoOST stellar models at core-helium exhaustion.In the adopted models at  = 0.016 Z ⊙ , PISN explosions occur in stars with initial mass between 107.2-203.4M ⊙ .As differences between zero and non-zero metallicity PISN models have been shown to be small (Kozyreva et al. 2014;Gilmer et al. 2017), the PISN explosion energies between 5 × 10 51 -9 × 10 52 erg and the chemical yields are interpolated logarithmically according to the zero-metallicity massive star (helium-core) models of Heger & Woosley (2002) 4 .The total mass released in the PISN explosion equals the mass of the progenitor star at helium-exhaustion, with a combined chemical composition given by the helium-core model of Heger & Woosley (2002) and the rest of the star as given by the surface composition of the BoOST model.Because there is no remnant, the stellar particle is removed from the simulation upon the explosion.
Lower-mass stars with initial masses between 0.5-8 M ⊙ are assumed to end their lives in an AGB phase, releasing a single burst of enriched wind-material with yields interpolated from Karakas (2010).In practice, the simulation of 65 Myr, and especially the starburst of < 10 Myr, concern too short timescales to allow for a significant AGB contribution.
The stellar winds and SN injections are implemented using the healpix algorithm to distribute the energy and mass around the stars into 12 equal-sized pixels, each containing 8±2 gas particles.Further details can be found in Hu (2019), Lahén et al. (2023) and references therein.

Stochastic wind output of an idealized stellar population
A few representative examples of stellar tracks selected from the BoOST models are shown here in Fig. 1, in Szécsi et al. (2022, their Figs. 2 and 3) and in Lahén et al. (2023, their Fig. 2).As clustered star formation and feedback are the focus of the current study, we illustrate the output of an idealized stellar population through stellar winds in Fig. 2. We have compiled the averaged wind output of three Monte Carlo cluster populations with total mass of 10 7 M ⊙ composed of 1000, 100 and 10 uniform-age star clusters with masses of 10 4 M ⊙ , 10 5 M ⊙ and 10 6 M ⊙ , respectively.Each cluster has been populated with a stochastically sampled Kroupa (2001) IMF between 0.08-500 M ⊙ .The wind output of stars between 9 M ⊙ and 500 M ⊙ in each single cluster has been integrated to acquire the wind-mass released by the cluster by the time it reaches a stellar age of 2 Myr, 4 Myr and 6 Myr.Fig. 2 shows the mean and mimimum/maximum wind output of the three Monte Carlo sets of clusters (10 4 M ⊙ to 10 6 M ⊙ ) binned by stellar mass and the respective cumulative output with increasing stellar mass, i.e. the contribution of cluster stars up to mass  * .The wind output is shown as normalized both to the total wind output of the averaged cluster sample until cluster age of 25 Myr and to the cluster-mass to indicate the relative contribution of different mass bins.
First we note in Fig. 2 the dominating role of massive stars already very early on in the evolution of a uniform-age cluster.More than 40% of all wind-material has already been released by the time the cluster reaches 2 Myr in age, with the major contribution from very massive (> 100 M ⊙ ) stars.Stars more massive than 50 M ⊙ produce more than 90% of all wind-material in the cluster.
The second aspect emphasized in this simple Monte Carlo experiment is the impact of stochasticity on the output of the star clusters.For lower-mass star clusters (e.g. the 10 4 M ⊙ sample) the variation of stellar masses at the high-mass end produces a significant scatter in the integrated wind output.Depending on the distribution of very massive stars, the maximum wind output of a given stellar mass bin can exceed the mean value of the sample by more than an order of magnitude.The lower limit of zero in the 10 4 M ⊙ and 10 5 M ⊙ samples indicates no stars in those bins in some of the randomly populated clusters and hence a zero contribution to the total wind yield.As can be seen in Fig. 2, the wind output (and, consequently, the total energy output) of massive stars is highly susceptible to stochasticity even in relatively massive star clusters.The star formation and feedback methods introduced here and in Lahén et al. (2023, see also Hirai et al. 2021) have been specifically tailored to capture some of star clusters, computed using the BoOST stellar tracks for massive stars.The wind production is binned (∼ 0.04 dex per bin) according to stellar mass, each cluster has a uniform age, and we show the wind-material produced until cluster age of 2 Myr (yellow), 4 Myr (red) and 6 Myr (blue).Each cluster has been randomly populated with an IMF between 0.08 M ⊙ and 500 M ⊙ , and we sample in total 10 7 M ⊙ of stars by Monte Carlo sampling the 10 4 M ⊙ , 10 5 M ⊙ and 10 6 M ⊙ clusters 1000, 100 and 10 times, respectively, to showcase the effect of stochastic sampling.The solid lines show the mean wind-mass of each Monte Carlo set at a given cluster age, and the shaded regions show the minimum and maximum value in each mass bin to indicate the range of most extreme values with respect to the mean.The top row shows the wind production at a given stellar mass, and the bottom row shows the cumulative distribution of wind-material with increasing stellar mass.The wind-masses are shown normalized to the total produced wind-material by the time the cluster has aged to 25 Myr (left y-axis) and to the mass of the respective cluster (right y-axis).Very massive stars (> 100 M ⊙ ) dominate the mass-loss during the first couple of Myr, 90% of the wind-material is released by > 50 M ⊙ stars already by cluster age of 6 Myr, and the impact of stochasticity on the minimum and maximum wind-mass released per mass bin reduces with increasing cluster-mass.
this stochasticity in a physically motivated manner.The aim is to improve upon traditional galaxy scale simulations that implement stellar feedback either using stellar population particles or assuming always a fully sampled IMF.

Detection of star clusters
In the dwarf starburst simulation, we follow the formation of hundreds of star clusters across the mass range from a few 10 M ⊙ up to more than 10 5 M ⊙ .Star clusters are identified using the structure finding algorithms friend-of-friends and subfind (Springel et al. 2001;Dolag et al. 2009) included in gadget-3.For constructing the friend-of-friends catalogue we use a linking length of 0.2 pc.Once subfind is used to find gravitationally bound substructures in the friend-of-friends catalogue, we consider gravitationally bound star clusters containing at least 50 stars in our final analysis.This corresponds to a typical minimum cluster-mass of ∼ 25 M ⊙ at our adopted input IMF.Unless stated otherwise, we analyse the cluster data as recorded in the final snapshot of the simulation.

Initial conditions
We follow the procedure in Lahén et al. (2020a, and references therein) in setting up the initial conditions for the dwarf galaxy encounter.The initial conditions consist of two identical compact gas rich dwarf galaxies selected from the simulation suite presented in Lahén et al. (2023), set up with methods of Springel (2005).The compact disk configuration is otherwise identical to that used in Hu et al. (2017) and Hislop et al. (2022) except for the initial metallicity which we set to ∼ 0.016 Z ⊙ .The initial dwarf galaxies are chosen right before star formation begins in the isolated configuration of Lahén et al. (2023).The initial disk scale lengths are 0.73 kpc, the disk gas fraction is 66%, and the virial mass and total baryon mass are 2 × 10 7 M ⊙ and 6 × 10 7 M ⊙ , respectively.Gas particles and the old pre-existing disk stars have initially a 4 M ⊙ resolution, dark matter is represented as 6.8 × 10 3 M ⊙ particles, and gravitational softening lengths are 0.1 pc for stars and gas and 62 pc for dark matter.Compared to Lahén et al. (2020a), we reduce the approach velocities of the parabolic encounter by a factor of four to allow a strong tidal interaction between the dwarf galaxies already during the first passage.In addition, the gas distribution is more centrally concentrated in the compact dwarf galaxy models compared to the initial conditions of Lahén et al. (2020a).These modifications shift the starburst observed during and after the second encounter after ∼ 170 Myr of simulation time in Lahén et al. (2020a) to the first passage between the two interacting galaxies already at ∼ 60 Myr.As a result, the galaxies experience an off-nuclear starburst of up to 0.5 M ⊙ yr −1 along the tidal bridge between the galactic disks, ensuring as enhanced star and star cluster formation in as unenriched ISM conditions as possible.We continue the simulation until the most massive star clusters have formed the majority of their stellar mass ∼ 7 Myr into the starburst, corresponding to a total simulation time of ∼ 65 Myr.

Star formation in the dwarf galaxy starburst
A visualization of the off-nuclear starburst and the related gas quantities are given in Fig. 3.The galactic disks can be identified as concentrations of young stars in the top right and bottom left quadrants of the first panel, connected with the star forming tidal bridge.
Similarly vigorous off-nuclear star formation is common in galactic interactions both in observations (Whitmore et al. 1999;Hibbard et al. 2001;Mullan et al. 2011) and simulations (Karl et al. 2010;Teyssier et al. 2010;Li et al. 2022) and qualitatively similar star forming configurations have been reported in interacting dwarf galaxies (Kimbro et al. 2021;Mićić et al. 2023).The tidal bridge is forming the majority of new stars in the simulated galaxies, and the star clusters in the bridge are more than two orders of magnitude more massive than those formed in the identical galaxies run in isolation in Lahén et al. (2023).
Fig. 3 also shows a sequence of star formation within the central ∼ 100 pc of the off-nuclear starburst.The filamentary gas structures form stars and star clusters, some of which accumulate hierarchically into more massive clusters.The formation of the most massive cluster in the simulation with a final bound mass of 2 × 10 5 M ⊙ can be followed ∼ 15 pc above the image centre.Cluster sub-structure and tidal features, resolved with single stars, are visible in and around the massive clusters.Fig. 4 shows the SFR during the dwarf galaxy interaction.The SFR is computed as the sum of stars formed in the past 1 Myr, 5 Myr or 10 Myr, the latter estimating the timescale over which e.g.observed SFRs based on ultraviolet bands tend to smooth the star formation history (Kennicutt & Evans 2012).The system enters into a starburst shortly after the galaxies start interacting for their first encounter.The most extreme SFRs exceed the quiescent phase by more than two orders of magnitude.We identify the starburst as the duration when the SFR is elevated compared to the early quiescent stages of the interaction and the isolated counterpart (see Fig. 9 in Lahén et al. 2023).In the following, we define  0 = 58 Myr as the start of the starburst and centre the time axis at this epoch when discussing the progression of the simulation.The  0 -centreed simulation time is indicated in Fig. 4 as well.
In addition to the SFR, Fig. 4 shows the cluster formation rate (CFR), computed as the sum of all clusters more massive than 100 M ⊙ with mean stellar age younger than 1 Myr, 5 Myr or 10 Myr identified according to Section 2.3, divided by the respective time interval of 1 Myr, 5 Myr or 10 Myr.We also show in Fig. 4   CFR, we observe an enhancement in clustered star formation during the starburst, in agreement with the Γ -SFR surface density correlation found in Lahén et al. (2020a) and various surveys of star forming galaxies (see e.g. the discussion in Adamo et al. 2020).This correlation is in contrast with the isolated dwarf galaxy simulations e.g. in Lahén et al. (2023) and Andersson et al. (2024) as well as recent observations of dwarf galaxies (e.g.Cook et al. 2023;Chandar et al. 2023) where no clear correlation is typically found.The observed dwarf galaxy populations do include dwarf starbursts (e.g.all three galaxies analysed in Chandar et al. 2023) as well as galaxies with very little star formation (e.g.minimum 1-10 Myr Σ SFR in Cook et al. 2023 is less than 10 −3 M ⊙ yr −1 kpc −2 ).It remains unclear whether the undetected Γ-SFR surface density correlation across dwarf galaxy populations is real or the result of observational limitations and low number statistics.
The environmental dependence of the star formation rate is illustrated in Fig. 5.We show the star formation rate surface density Σ SFR as a function of total gas surface density Σ gas in the merging dwarf galaxy system and in one of the progenitor dwarf galaxies run in isolation for reference.The data have been binned in 400 pc pixels and the data shown are averages over pixels that have had star formation within the past 10 Myr.The observational reference data are from various surveys targeting dwarf galaxies and more massive star forming galaxies in the local Universe (Leroy et al. 2008;Elmegreen & Hunter 2015;Roychowdhury et al. 2015) as well as star forming massive galaxies at higher redshifts (Tacconi et al. 2013;Genzel et al. 2010).The isolated dwarf galaxy and the initial stages of the dwarf galaxy interaction match the range of values observed in local systems that host low levels of star formation.As the starburst begins, first the Σ gas increases, followed by increasing amounts of star formation up to 0.1 M ⊙ yr −1 kpc −2 .The most extreme starburst phase The yellow circles show one of the progenitor dwarf galaxies run in isolation over 500 Myr in 50 Myr steps and the squares coloured by time and connected by lines show the progression of the dwarf galaxy merger in 5 Myr steps before the starburst (−48 Myr <  −  0 < 0 Myr) and in 1 Myr steps during the starburst (0 Myr <  −  0 < 7 Myr).The star formation rates are always computed over past 10 Myr and averaged over 400 pc pixels, to match the typical averaging scale in the reference observations.The observed data points are from the THINGS survey (gray dots, Leroy et al. 2008), the LITTLE THINGS survey (orange dots, Elmegreen & Hunter 2015), 5-95 percentile range of the dwarf galaxies in the FIGGS survey (shaded gray, Roychowdhury et al. 2015) and high-z star forming galaxies (up and down pointing triangles, respectively, from Tacconi et al. 2013 andGenzel et al. 2010).The black diagonal line shows the observed fit Σ SFR = 1.6 × 10 −4 Σ 1.4 gas from Kennicutt (1998) that accounts for helium and metals by a correction factor of 1.36.
exceeds the normal star-forming relation indicated in the figure and extends toward the transition between the local star forming galaxies and higher redshift galaxies.Once the starburst subsides, the Σ SFR and Σ gas loop back down toward more quiescent values.

Formation of star clusters
As we already found in Lahén et al. (2020a) using 3D cluster data and in Lahén et al. (2022) using mock observations, the increase in SFR results in the formation of more massive star clusters.Fig. 6 shows the star cluster mass function (CMF = /) throughout the simulation in 10 Myr steps for star clusters with mean stellar age less than 10 Myr.Only the final snapshot, which represents the system after the starburst, includes young clusters more massive than ∼ 10 3 M ⊙ .The high-mass end of the CMF drops off the traditional power-law slope of −2 in the early stages when the SFR is low, similar to the CMFs found in simulated dwarf galaxies in Hislop et al. (2022) and Andersson et al. (2024), while at enhanced SFR during the starburst (Fig. 4) the power-law slope is in agreement with −2 from 10 3 M ⊙ up to more than 10 5 M ⊙ .
Massive star clusters form hierarchically within interconnected filaments along the tidal bridge (see Fig. 3).Fig. 7 shows the individual star formation histories of all clusters more massive than 10 4 M ⊙ .The time axis has been shifted by the mean stellar age of each cluster, so that positive values correspond to the stars formed closer to the end of the simulation.The clusters form rapidly over timescales of less than 10 Myr across the mass range, with the majority of stars forming in a 2-3 Myr timespan.The standard deviation around the mean value spans from 0.7 Myr to 2 Myr.The peak SFR and the width of the star formation profile increase with cluster-mass.Many of the clusters show multiple bursts of star formation.The increasingly high SFRs in massive clusters increases the chance of massive stars to form and the prolonged star formation histories allow the ejected material of those massive stars to be captured in regions of ongoing star formation.On the other hand, the period of star formation is short compared to the lifetime of the progenitors of core-collapse SNe (∼ 5 Myr), thus only PISN ejecta would be able to contribute to self-enrichment on short enough timescale.

Initial stellar mass function
Fig. 8 shows the initial masses of all stars realised in the simulated dwarf galaxy interaction.We also show the initial masses of all stars in clusters more massive than 10 3 M ⊙ , compared to the combined mass function of lower-mass clusters and field stars that are not bound to any cluster.We compare the realized IMF to the IMF shape of  * / * ∝  −1.3 below 0.5 M ⊙ and  * / * ∝  −2.3 above 0.5 M ⊙ that was used to randomly sample the initial masses.In the high-mass end, where we have to group particles to conserve mass, the realized stellar masses are limited by the local sampling reservoir (reservoir particles, Section 2.1).The histograms have been normalized to the total number of stars in each sample to ease comparison.
In contrast to the corresponding isolated compact dwarf galaxy simulation in Lahén et al. (2023), where we had only a handful of stars close to 100 M ⊙ and only one star in the PISN mass range in the entire simulation sample, we now see almost the entire IMF range realised during the starburst.Until the end of the simulation, more than ∼ 16 000 massive stars have formed in the system, out of which 120 had masses in excess of 100 M ⊙ .682 massive stars have already died: there have been 22 PISN explosions and 95 stars have died as direct collapse black holes with final black hole masses in the range  Comparing the initial masses of stars in the field and low-mass clusters to the stars in massive clusters, the mass function in the massive clusters follows better the high-mass slope of  −2.3 * up to tens of solar masses.Very high masses above 200-300 M ⊙ are only present in the most massive clusters that reach more than 10 4 M ⊙ .Conversely, regions corresponding to lower star formation activity exhibit a more pronounced deficit in initially massive and very massive stars.This difference shows the impact of the local mass conservation and the Jeans length dependent mass reservoir in the IMF sampling, which ensure that we allow massive stars to be realised strictly only when there has been enough star forming gas available within the region where stars are being spawned.There is a deficiency of initially very massive stars compared to a fixed power-law IMF across the simulation, therefore a fully sampled IMF is not guaranteed even in the most intense star forming regions that host the most massive clusters of the low-metallicity galaxies.

Enrichment by massive stars
The 13 individual elements (Section 2.2) released through stellar winds and SN explosions are traced throughout the starburst.Fig. 9 shows the integrated total and relative masses of wind and SNmaterial across the starburst.From left to right and top to bottom we show The last two quantities have been computed in 0.1 Myr steps from the wind and SN-material locked in cluster stars at any given moment, divided by the total wind and SN-material produced by massive stars in total.For comparison, we also indicate in Fig. 9 the final values in a single dwarf galaxy analysed in Lahén et al. (2023), measured after 500 Myr of isolated evolution.
As the starburst ramps up, massive stars form in increasing numbers, releasing equivalently increasing amounts of wind-material as they evolve.Meanwhile, the total mass of SN-material is only starting to build up Myrs into the starburst.This is partly because the majority of the most short-lived stars die in a direct collapse without a SN (e.g.87% of > 40 M ⊙ stars that live < 4 Myr), with the PISN progenitors being the exception.SN-material still always dominates the mass returned by stars in the simulation, however by the end of the starburst the integrated winds have narrowed the gap down to only a factor of 4. The wind and SN-material recycled in all new stars is enhanced by two orders of magnitude during the starburst, amounting to approximately 11% of the wind-mass and 6% of the SN-mass produced in total by the time the starburst is over.As indicated in Fig. 9, the value of approximately 10% of recycled feedback material is similar to that accumulated in stars of the isolated compact dwarf galaxy after 500 Myr of evolution in Lahén et al. (2023).The dwarf starburst can therefore recycle wind-material significantly faster than either of the progenitor galaxies evolved in isolation.Winds are more prone to be locked in clustered stars compared to SN-material: by the end of the starburst, 91% of the wind-material recycled in stars resides in clusters, compared to the corresponding value of only 8% for all SN-material.The inefficient recycling of SN-material in clusters is also seen directly in  SN, * ,cluster / SN, from clusters , which shows only 0.6% of the SN-mass produced by cluster stars as recycled in the clusters compared to the value of ∼ 11% for wind-material.
In total, wind and SN-material account for ∼ 10 −4 and a few 10 −4 , respectively, of the total stellar mass formed by the end of the starburst.The difference between the two feedback channels is much smaller than what was found in the case of isolated dwarf galaxies in Lahén et al. (2023), where the SN-material in stars was built up gradually over 500 Myr through cooling and mixing within the galactic ISM.As can be seen in the top left panel of Fig. 9, the spatially clustered, rapidly evolving star formation during the starburst occurs on a shorter time scale than what would be required to cool and recycle the SN-material within the region that is producing it.This is demonstrated in the bottom panels, where the SN-enrichment within star clusters is an order of magnitude less efficient compared to the overall SN recycling.Even though the starburst has both produced and locked in stars a similar amount of SN-material (∼ 10 4 M ⊙ ) compared to the isolated counterpart at 500 Myr, the ten-fold larger stellar mass in the starburst results in a ten-fold smaller contribution of SN-material to the total stellar mass.Meanwhile, because the winds are released continuously with lower injection energies, they are easier to capture in the gas within the same star formation region.The formation of wind-enriched stars is dominated by enrichment from and to the star clusters, as can be seen in the converging lines for all stars and clusters in Fig. 9.We will investigate this selfenrichment scenario within star clusters on a cluster-by-cluster level in the following Sections.

Cluster-averaged wind-mass fraction
To address whether chemical enrichment through stellar winds and SNe may result in self-enrichment during the formation of star clusters, we take a detailed look at the population of bound star clusters formed in the dwarf starburst.We focus on massive star clusters from 10 3 M ⊙ to > 10 5 M ⊙ , which host the majority of the massive and very massive stars (Fig. 8) that are the sources of enrichment in the simulation.In the two following Sections, we consider the properties of the clusters as averaged over bound stars in each cluster at the end of the simulation, as determined in Section 2.3.
We show in Fig. 10 the mass fraction of wind-material in each cluster measured within twice the half-mass radius  1/2 .The correlation of the wind-mass fraction with cluster-mass indicates enhanced selfenrichment via stellar winds in massive star clusters.More massive star clusters form in more dense star-cluster-forming environments, which in our implemented star formation routines (Section 2.1) leads to the formation of more massive stars both in stellar mass and in number.Since the scaling of the wind-mass fraction is linear or superlinear with cluster-mass, it follows that the wind-material recycled in cluster stars scales as  wind, * ,clusters ∝  2 cluster .The more massive clusters therefore produce and/or capture increasing amounts of wind-material per unit cluster-mass.
For comparison, self-enrichment in star-forming molecular clouds of varying initial density has also been investigated by Lancaster et al. (2021c) in high-resolution cloud-scale simulations.Unfortunately the models concerned only solar metallicity, and chemical composition of the winds of their massive stars (up to 100 M ⊙ ) was not a focus.Nevertheless, their higher initial densities resulted in higher efficiency of star formation while the new stars also captured increasing amounts of wind-material, in agreement with our findings  in Fig. 10.Their total wind mass fractions ranged from 10 −4 to a few times 10 −3 in their clusters with masses between 2.8 × 10 4 M ⊙ and 8.5 × 10 4 M ⊙ , indicating a super-linear relation between the windmass fraction and the cluster-mass as well.Given their assumption of a fully sampled Kroupa IMF, they attributed the increased efficiency of wind-capture to two factors: the enhanced trapping of the windbubbles in higher density gas (see the theoretical arguments in Lancaster et al. 2021a andLancaster et al. 2021b), and the longer duration of star formation that can continue for multiple free-fall times in the densest star-forming regions.Krause et al. (2016) have also argued that the compactness of the cluster, defined as  5 = ( cluster /10 5 M ⊙ )/( 1/2 /pc), may govern whether the cluster gas can be retained against expulsion due to stellar winds.Using the definition of cluster compactness from Krause et al. (2016), we recover values between 0.1-4 for clusters more massive than 10 4 M ⊙ with half-mass radii between 0.16-1 pc (see Fig. 15).Krause et al. (2016) report very similar values for young massive star clusters from Bastian et al. (2014), as well as present day GCs that are known to host chemically peculiar stars but should have been more compact in the past.For lower-mass clusters, the compactness increases from 0.001 to 0.4 with increasing mass, similar to those reported for open clusters in Krause et al. (2016) where no multiple populations are found.Wünsch et al. (2017) and Lochhaas & Thompson (2017) have further shown that stellar winds can radiatively cool on a rapid timescale of a couple Myrs if the cluster is compact enough.For us, an additional component is the non-universal high-mass end of the IMF, where the less massive clusters simply form less massive stars.

M
wind, * , cluster /M cluster slope 0.9±0.06M > 10 4 , slope 1.1±0.2 Figure 10.The wind-mass fraction in bound star clusters more massive than 10 3 M ⊙ , measured within 2 1/2 at the end of the simulation.One single cluster with no wind-material has been placed at 10 −7 .The lines show the best-fit power-law relations ∝   cluster to all clusters (dashed) and to clusters more massive than 10 4 M ⊙ , and the power-law slopes are indicated in the legend.Clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig. 7.The wind-mass fraction, sourced mainly through self-enrichment, scales approximately linearly with cluster-mass and hence  wind, * ,clusters ∝  2 cluster .
material or whether the deeper cluster potential is able to capture more wind-material, we inspect in Fig. 11 the wind-mass produced and recycled in the clusters shown in Fig. 10.The relative windproduction is the ratio between the integrated wind-mass produced by all massive stars in the cluster, divided by the cluster mass.The recycled wind-mass fraction has been computed as the ratio between the wind-mass locked in the cluster stars and the wind-mass produced by the clusters.The wind-mass produced by each massive star bound in each cluster has been integrated over the stellar lifetime, assuming that the star was always bound to the cluster in which it resides in the last snapshot.The wind-production may represent lower limits, since massive stars may have escaped and some have already exploded as PISN, leaving no trace of how much wind they produced.Tracing back all stars that contributed to the wind-enrichment would require detailed tracing of all of the clusters, which is beyond the scope of the current study.To minimize the impact of cluster evolution, we analyse the clusters immediately at the end of the starburst, within a few Myr after their formation.If the recycled fraction is above unity, the cluster has been enriched by stars that are not (anymore) bound to the cluster.The large scatter in the both panels is driven by the stochastic IMF.
The relative wind-production scales moderately with cluster-mass ( wind from clusters / cluster ∝  0.3−0.6 cluster ), therefore the total windproduction has a slightly super-linear scaling with cluster mass ( wind from clusters ∝  1.3−1.6  cluster ).This reflects the fact that higher mass clusters are better guaranteed to host a larger number of massive stars, with increasingly high upper mass limit in their IMF.If all wind-material would be captured as efficiently in all clusters, variations in the cluster-to-cluster wind-production alone would produce a too weak scaling in the wind-mass per cluster compared to the one measured in Fig. 10. .Top: the relative wind-mass produced by clusters more massive than 10 3 M ⊙ , measured in the final snapshot as the integrated wind-output of all massive stars bound in a cluster related to the cluster mass.Bottom: the recycled wind-fraction in each cluster in the top panel, computed as the wind-material locked in cluster stars divided by the wind-material produced by the cluster.The lines in both panels show the best-fit power-law relations ∝   cluster to all clusters (dashed) and to clusters more massive than 10 4 M ⊙ , and the power-law slopes are indicated in the legends.Clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig. 7.Even though the majority of the clusters are very young (mean stellar age < 5 Myr), the produced mass-estimate is a lower limit as it excludes escaped stars as well as stars that have exploded as PISN since the PISN events leave no remnants.Due to the lower-limit of total produced mass in the top panel, the retained ratios, especially in the higher cluster-mass range, are upper limits.
The recycled fraction, likewise, scales as  wind, * ,cluster / wind from clusters ∝  0.1−0.7 cluster , with a steeper correlation toward larger cluster masses.If the recycled fraction would be constant, the correlation of wind-mass fraction with cluster mass in Fig. 10 could be expected to emerge directly from variations in the cluster IMF.However, because the recycled fraction has a positive correlation with cluster mass, especially at high cluster masses, there seems to be a contribution of enhanced wind-trapping The SN-mass fraction in clusters more massive than 10 3 M ⊙ , measured within 2 1/2 .Clusters with no feedback material have been placed at 10 −7 and clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig. 7.The errorbars show the linear mean and standard deviation of the binned SN-mass fraction to guide the eye at low cluster-masses where many of the clusters have not been enriched by SNe.The SN-mass fraction remains flat with cluster-mass, reflecting SN-material that is uniformly mixed across the galaxy.
with increasing cluster-mass.The most massive clusters here have the highest concentration  5 , indicating that the most massive clusters should be most resistant against gas expulsion and result in high recycling efficiencies in Fig. 11.The strongly super-linear wind-mass -cluster-mass relation  wind, * ,clusters ∝  2 cluster in Fig. 10 is therefore likely driven by the combination of increasing number of massive stars and more efficient capture of the wind-material.

Cluster-averaged SN-mass fraction
The mass fraction of SN-material in cluster stars is shown in Fig. 12.The SN-fraction remains flat throughout the cluster-population analysed at the end of the simulation.This means that the total SNmaterial locked in cluster-stars scales only linearly with cluster-mass,  SN, * ,cluster ∝  cluster .Since the SN production is only slowly increasing by the time the starburst is already over, the clusters do not have the time to self-enrich by SNe during their formation.This is also demonstrated in Fig. 7, where the duration of the star formation in even the most massive clusters is only 5 Myr or so, barely enough for the first core-collapse SNe to occur.The only possible SN contributor on such short timescales would be the PISNe in clusters massive enough to host very massive stars.PISNe are, however, more energetic and extremely rare in the simulation compared to normal SNe as there is one PISN for every 188 SN events.Instead of local self-enrichment, the flat SN-mass fraction in Fig. 12 illustrates mixing of small amounts of SN-material throughout the galactic ISM: the clusters form during the starburst out of a gas mixture uniformly polluted by traces of SN-material.
The lower-mass clusters that host SN-enriched stars show a large scatter in mean SN-mass.This is partly driven by the adopted SN injection scheme where each SN is on average distributed over a fixed number of ∼ 96 gas particles.The SN-enriched stars occur in clumps of at least 4 M ⊙ whenever a SN-enriched gas particle is sampled into stars, increasing the contribution of one such clump in low-mass star clusters.

Heavy-element enrichment
One of the key features of bona fidé GCs are their fairly uniform heavy element abundances across the stars of a given cluster, indicating no significant SN-enrichment during their formation.Recent spectral studies have however shown that some GCs do in fact host stars that have small spreads e.g. in their iron abundances, especially in the first, chemically normal populations (Lardo et al. 2023;Marino et al. 2023).When there is a spread in [Fe/H], its presence in the first population reflects the initial ISM composition, rather than selfenrichment.This observational interpretation is in agreement with our result of uniform SN-content and little to no SN self-enrichment even within the most massive star clusters simulated here.
Fig. 12 shows that our simulated clusters do, on average, contain a small amount of SN products.To assess the heavy-metal enrichment in our simulated clusters, Fig 13 shows the mean [Fe/H] abundance ratio for cluster stars within 2 1/2 of the centre of each cluster.The abundances have been shifted by the initial galactic ISM composition to highlight the enrichment by heavy elements, or the lack thereof, compared to the initial composition.All the analysed clusters have been enriched in total by less than 0.01 dex in iron compared to the unenriched ISM.The majority of the stars in all clusters have formed out of an ISM composition that does not include any SN products produced during the simulation.The mean [Fe/H] enrichment for most of the low-mass clusters is practically zero even when essentially all of the clusters contain some form of feedback material, because the stellar wind composition is never significantly enriched in heavy elements.
The SN-material is concentrated in up to 1% of the stars in each cluster.These stars can have highly enhanced [Fe/H] values of up to 1.4 dex due to limitations of the SN-injection method that does not include diffusion of metals between particles.Due to these SNenriched stars, the total [Fe/H] standard deviation of the stars in a given cluster is up to 0.07 dex (not shown in Fig. 13 for clarity), even though the mean [Fe/H] is always very low.For comparison, the maximum star-to-star spread is in observed iron-normal GCs Δ[Fe/H] ∼ 0.2-0.3(Lardo et al. 2023;Marino et al. 2023).As noted above, the [Fe/H] spread in our simulated clusters is driven by low levels of SN-enrichment.The total star-to-star [Fe/H] spread from winds alone reaches at most Δ[Fe/H] ∼ 0.01 dex with respect to the unenriched ISM and the 99-percentile of Δ[Fe/H] from windenrichment is within ∼ 0.005 dex in all clusters more massive than 10 3 M ⊙ .

Radial distribution of enriched stars
Fig. 14 shows the azimuthally averaged radial profile of the star-bystar wind-mass fraction in clusters more massive than 10 4 M ⊙ .The profiles show a transition from a flat wind-mass fraction profile at low cluster-masses to a centrally concentrated profile as the cluster-mass approaches 10 5 M ⊙ .The average central wind-mass fractions are up to two orders of magnitude enhanced compared to the outer regions of the clusters.Lancaster et al. (2021c) argued that the efficiency of wind-recycling was essentially determined by the gas-density in the star-forming regions, which governs the ability of winds to be trapped while the star formation continues, possibly for multiple free-fall times.Our centrally enhanced wind-mass profiles may be the result of prolonged star formation histories in the central regions of the most massive clusters as indicated in Fig. 7.However, with our locally limited IMF, the most massive stars are formed in the central regions of the most massive clusters, hence they also deposit their winds in those locations.
Fig. 15 shows a comparison between the radial wind-mass fraction and SN-mass fraction in the 20 most massive clusters from Fig. 14.The wind material that originates from self-enrichment is often centrally concentrated or flat because star formation occurs preferentially in cluster centres during the late stages of cluster assembly and because the central region can capture the ejected wind material.On the other hand, the SN-material that does not have an origin within the clusters themselves is typically either flat or centrally depressed.The most centrally concentrated wind fraction mass profiles also coincide with centrally depressed SN-mass fraction profiles.
In observed GCs, the chemically peculiar stars are often (Dalessandro et al. 2019) but not always (Leitinger et al. 2023) found to be more centrally concentrated than the unenriched stars.Central concentration should be a natural outcome of the self-enrichment scenario if the enriched stars form and release their winds in the deepest parts of the gravitational well of the forming massive cluster.The central radii of a massive star cluster could also be expected to be the only region able to retain the enriched material in the presence of the energetic massive stars that produce the enriched stellar winds.The enrichment properties of our simulated clusters seem to agree with the observations that show a centrally concentrated second population.

Chemical abundances of enriched populations of stars
A detailed look at the wind-material across cluster stars is given in Fig. 16 where we show the distribution of stellar wind-mass fractions  wind in stars within 2 1/2 of each cluster.As reflected in the wind-mass fractions in Fig. 10, the peak of the wind-mass distribution shifts toward higher  wind with increasing cluster-mass.Consequently, the most wind-enriched stars in each cluster also shift towards higher  wind , with the most enriched stars consisting of 10% of wind-material.The width of the distribution increases with cluster-mass, which we will discuss in more detail below.In the massive clusters, typically 10% of the stars have not been enriched by winds at all.The fraction of unenriched stars increases with decreasing cluster-mass, as lower-mass clusters host on average lower-mass stars that evolve slower and release less winds compared to the most massive stars in the more massive clusters.
Next, we inspect the variation of abundances often used to identify MPs in GCs: Na, O, Mg, Al and N. Fig. 17 shows the abundance ratio distribution in [Na/Fe], [O/Fe], [Al/Fe], [Mg/Fe] and [N/Fe] in stars of the most massive clusters (> 10 4 M ⊙ ).The variations in these abundances are dominated in the simulated clusters by windenrichment: stellar winds are the source of at least 95% of the recycled mass in all of these elements.We compute the diluted composition of cluster-stars (left column) as given by the mixture of wind-material and ISM, inherited from the gas out of which the star was born.Since observed GC-stars are often unevolved low-mass stars, we do not here include the internal chemical evolution of the stars themselves.The diluted abundance ratios reflect the chemical composition of each star when they formed, hence any change from the initial pristine value will show the level of chemical pollution in the star-forming environment.The un-diluted values (right column) exclude the ISM component locked in the stars, and thus shows only the composition of the recycled wind-mixture (see Fig. 1).The un-diluted mixture would represent the most extreme chemical variations if the stars were 100% made of such material at formation.Because the enrichment of these elements is driven by the stellar winds, we can safely exclude the SNenriched stars whose chemical contents are artificially enhanced due to the discrete SN injection scheme and the lack of metal diffusion discussed in Sections 2.2 and 3.3.3.This excludes up to 1% of both the stars and the wind-mass locked in those SN-enriched stars in each clusters.
In the right hand column of Fig. 17, the pure wind abundances show a strong anticorrelation between the enhanced [Na/Fe] and depleted [O/Fe], typical for nuclear reactions taking place at high temperatures as shown in Fig. 1, and similar to the abundance ratios as measured in GCs.[N/Fe] is also strongly enhanced.However, when the stellar wind abundances are mixed with the ISM to form the bulk of the stellar mass within the cluster (left column), the spread in abundance ratios narrows down.The diluted [Na/Fe] and [N/Fe] are still enhanced by 0.3-0.5 dex.The enhanced [N/Fe] that persists even after dilution is particularly interesting in relation to recent observations of extreme star forming objects at higher redshifts (Pascale et al. 2023;Bunker et al. 2023).Such N-enhanced systems have been suggested to represent young GCs observed during their formation (Charbonnel et al. 2023;Senchyna et al. 2023;Marques-Chaves et al. 2024).
The [O/Fe] and [Mg/Fe] depletion and the [Al/Fe] enhancement, on the other hand, are almost completely washed out after dilution with the ISM.Even if the stellar wind is strongly enriched in light elements as shown on the right of Fig. 17, the mass recycled in the cluster stars is not enough to produce a significant spread in abundance ratios, and falls short of the extent observed in ancient GCs.The 90% of mass at unenriched ISM composition dilutes the enriched abundances, resulting in only a modest spread in abundances even in the most enriched, massive (> 10 4 M ⊙ ) clusters.This is a manifestation of the so-called mass-budget problem.As a possible Figure 16.The mass fraction of wind-material captured in bound cluster stars, measured within 2 1/2 in clusters more massive than 10 3 M ⊙ .Stars with very little (< 10 −12 ) or no wind-material have been collected in a set of points at 10 −12 .Clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig. 7.The peak of the  wind distribution and the largest  wind value per cluster shift to higher values with increasing cluster-mass.caveat, we note that the stellar wind injection method injects the mass into existing gas particles.Hence, we always end up diluting the wind material by some amount.As a result, the absolute spread in abundance ratios can be viewed as a lower limit.For complementary examples of the upper limit abundance ratios in wind-material released and retained in stellar populations, accounting for radiation and semi-analytic hydrodynamics with no dilution, we refer the reader to Szécsi & Wünsch (2019).Having the possibility of stars made 100% of wind material would increase the number of extremely peculiar stars, possibly even up to the extreme values shown on the right of Fig. 17.Such redistribution of the wind mass would however not change the fact that the clusters return only up to a couple per cent of their mass in stellar winds, that is up to a few 10 3 M ⊙ for a cluster with initial mass of 105 M ⊙ .
The extent of the abundance ratio (anti-)correlations increases with cluster-mass in terms of absolute spread and standard deviation, as shown in the histograms of Fig. 17, following Figs. 10 and 16 where the increasing cluster-mass resulted in increasingly high wind-mass fractions.These correlations with cluster-mass are in qualitative agreement with observed GCs, where the maximum spread and standard deviation of light element abundances increase with the present-day GC mass (Pancino et al. 2017;Nataf et al. 2019).We will discuss the distribution of enriched stars star-by-star in the context of GCs in more detail below.

Implications for GC formation
The simulated star clusters form with up to 2 × 10 5 M ⊙ of stellar mass, corresponding to the lower mass range of evolved GCs.In the simulated clusters, the chemically peculiar stars that differ significantly from the unenriched population in Fig. 17 comprise a minority of the stars.This is in contrast to observed ancient GCs where the peculiar 2P stars start dominating between present-day cluster-masses of a few 10 4 M ⊙ and 10 5 M ⊙ .What would it take to reach a dominant, strongly enriched population of stars, based on the results here so far?Using the nomenclature of observed GCs, we will refer to the enriched population of stars as the second population, 2P, and the unenriched population as the first population, 1P.The defining factor between 1P and 2P is traditionally based on the separation between the populations in the abundance ratio plane, e.g. as a threshold in [Na/Fe] and [O/Fe].As can be seen in Fig. 17, the stars that contain 10% of wind-material are already quite separated from the unenriched composition and almost reach the typical 2P limit of [Na/Fe] ∼ +0.3 dex (Carretta et al. 2009b;Pancino et al. 2017) compared to the unenriched population.Based on the most extreme stars in Fig. 17, we can conclude that the stars must contain at least ∼ 10% of enriched wind-material to be considered part of the 2P.The actual extent of the 2P across the abundance ratio plane of course depends on the detailed combination of abundances and ISM mixture from which the stars form, however we use here the wind-mass fraction as the first order proxy for enrichment when extrapolating our results to higher cluster-masses.
As an illustrative experiment, we fit the star-by-star wind-mass fraction distribution in Fig. 16 with a distribution function and extrapolate the fit parameters to higher star cluster-masses.Fig. 18 shows the data in Fig. 16 for 20 most massive clusters, this time with linear y-axis normalized to the area below each distribution to obtain the probability density function (PDF).Because the massive stars release most of their increasingly enriched winds toward the end of the cluster formation sequence, they result in a skewed distribution toward enriched stars with a tapered cut-off at high values of  wind and a longer tail of low values of  wind .The distributions in Fig. 18 resemble a reversed log-normal or a gamma-distribution that have conveniently a natural lower limit of zero, corresponding to the upper limit of the wind-mass fraction of log 10 (  wind = 100%) = 0 in log-10 scale.After experimenting with fitting the profiles with both distribution functions, the log-normal shape seems to yield more robust results across the cluster population.When we reverse  wind and define  = − log 10 (  wind ) 5 , we can fit the wind-mass fraction PDFs with the log-normal distribution function (Johnson et al. 1994) where  ln and  ln are the mean and standard deviation of the natural logarithm of .The arithmetic mean (expected value) and standard deviation of the log-normal distribution are then given by and The best-fit log-normal distributions of wind-mass fraction locked in our simulated cluster stars are shown in Fig. 18, for the 20 most massive clusters.The x-axis is shown in log-10 scale for simplicity, and the log-normal profiles have been reversed back by multiplying the x-axis by −1.The corresponding values of mean and standard deviation are indicated as well.The  wind profiles show, for the most parts, a good match to a log-normal distribution.In some cases, the distributions deviate in the tails, while some of the profiles have a more peaked distribution than the best-fit log-normal.The final   .Left: the probability density function of the stellar mass composed of wind-material in stars of clusters more massive than 10 4 M ⊙ within 2 ×  1/2 , as in Figs.16 and 18.We only show the most massive clusters for clarity.The best-fit log-normal distributions (Eq.3), using inverted log 10 (  wind ) as the x-axis and zero point at  wind = 100%, are shown with the same colour coding as the simulated cluster data from Fig. 7. Right: the arithmetic mean −  and standard deviation  (Eqs.4 and 5) of the log-normal fits to  wind in clusters more massive than 10 3 M ⊙ , as well as the power-law scaling relations fit to −  and .The lines show the best-fit power-law slopes (∝   cluster , dashed) and the slope upper limits (slope plus bootstrapped standard deviation, solid) , with values indicated in the panels.The black curves in the left panel show log-normal profiles extrapolated according to the best-fit (dashed) and the upper limit (solid) parameter-scaling in the right hand panels for clusters with mass of 10 6 M ⊙ (thin) and 10 7 M ⊙ (thick).The legend indicates the fraction of stars composed of more than 10% of wind-material in the extrapolated distributions.
profile is the convolution between the star formation and enrichment histories in the cluster forming region.Even though the SFRs (Fig. 7) are often gaussian-like, sub-cluster mergers or irregular patterns in the star formation and enrichment histories may cause irregularities in the  wind distribution.
Because the trend between the  wind shape and the cluster-mass is perhaps not obvious in the individual panels of 18, we show in Fig. 19 all the  wind profiles of clusters more massive than 10 4 M ⊙ , together with their best-fit log-normal distributions.As seen already in Fig. 16, the distributions in Fig. 18 tend to shift to higher values of  wind and broaden with increasing cluster-mass.To quantify the scaling of the  wind -distribution with cluster-mass, we also show the correlation of the arithmetic mean (−) and standard deviation () of the log-normal distributions with cluster-mass for all massive clusters in the right hand side of Fig. 19.We fit the log-normal parameters with a power-law relation between − and  against cluster-mass, and estimate the standard deviation of the best-fit power-law slopes by bootstrapping 1000 times.The mean and standard deviation of the  wind -distributions increase with cluster-mass with power-law slopes of 0.82 ± 0.12 and 0.18 ± 0.07, respectively.In the following, we use these power-law scaling relations to extrapolate the distribution of  wind to higher cluster-masses.
Using the power-law fits to the parameters of the log-normal distributions shown in Fig. 19, we first extrapolate the wind-mass fraction distribution to 10 6 M ⊙ and 10 7 M ⊙ clusters.The extrapolated distributions are shown in Fig. 19 with black, with thin lines for 10 6 M ⊙ and thick lines for 10 7 M ⊙ clusters, respectively.The dashed black lines show the distribution extrapolated using the best-fit relation of  and , and the solid lines show the upper limit scenario, extrapolated according to the best-fit plus its standard deviation (as given in the legend of the right hand panels of Fig. 19).The legend in the left hand panel of Fig. 19 indicates the fraction of stars consisting of more than 10% of wind-material, obtained by integrating over the the extrapolated log-normal distributions.These stars are considered as the 2P of the extrapolated clusters, based on the most enriched stars that are on the verge of reaching the traditional 2P limit in the simulated clusters shown in Fig. 17.
Assuming that the stellar wind released and captured in cluster stars scales with cluster-mass as given by the distribution of starby-star wind-mass fractions, the conservative extrapolation would result in only up to ∼ 0.02% and ∼ 2.3% of strongly enriched stars in 10 6 M ⊙ and 10 7 M ⊙ clusters, respectively.The picture changes slightly when we consider the extrapolated distributions using the upper limits of − and : the 10 6 M ⊙ and 10 7 M ⊙ clusters would have ∼ 0.3% and ∼ 11% of strongly enriched stars.We are therefore still not able to produce a cluster with a dominant 2P as observed in many GCs, even after extrapolating the enrichment properties of our simulated cluster to almost two orders of magnitude higher cluster masses.

Evolution of 2P fraction due to GC mass-loss
A ∼ 10% 2P fraction is low compared to the typical values in the Local Group where 10 6 -10 7 M ⊙ GCs show dominant 2P fractions of up to 90%.However, the 2P fractions estimated above describe the initial division between 1P and 2P stars immediately after the clusters have formed.Because the GCs observed in the Local Group are very old and orbit in a galactic potential, they must have lost some stars throughout their evolution in the galactic environment (Lamers et al. 2010;Schaerer & Charbonnel 2011;Baumgardt et al. 2019;Wang et al. 2016;Moreno-Hilario et al. 2024;Lacchin et al. 2024).The fraction of mass lost depends on both the initial properties of the cluster (mass, density profile, IMF) and the galactic tidal field, and estimates for the mass-loss of Milky Way GCs extend from half to more than 90% of the initial GC mass (Webb & Leigh 2015;Baumgardt et al. 2019).Some of the mass-loss is due to stellar evolution, and the rest is caused by stars escaping the cluster.A detailed simulation of massive clusters through Gyrs of evolution is beyond the scope of the current study, especially since such modelling would need to account for the tidal field of the Milky Way in addition to that of our initial dwarf galaxy system.We instead approximate the evolved 2P fraction in the extrapolated clusters by parametrizing the loss of stars over time.We assume that a fraction of the cluster stars are lost throughout the evolution of the clusters, with 1P dominating the mass-loss.This assumption is based on observations of the spatial distribution of the populations and our simulated radial profiles of enriched stars in Fig. 14, wherein the 2P stars are often more centrally concentrated within GCs and therefore less susceptible to be lost to the tidal field.To compute the final extrapolated 2P fraction we use  19 and factoring in the loss of stars due to cluster evolution in a galactic tidal field.The extrapolated relations assume that either 90%, 50% or 0% (from darkest to lightest opacity) of the 1P stars are lost.The dashed lines show the 2P fraction obtained by integrating over the best-fit log-normal wind-mass fraction distributions, and the solid lines show the upper limit using the best-fit parameters plus standard deviation, as illustrated in Fig. 19.The observed Milky Way and Magellanic GCs are from Tables 7, 8 and 9 in Gratton et al. (2019), see text for details.
assuming that the number of 2P stars initially  2P,i does not evolve, i.e.  2P,f =  2P,i , while a fraction  1P of the initial 1P stars  1P,i is lost and results in the final 1P fraction of  1P,f =  1P  1P,i .We assume only two populations, therefore  1P,i = 1 −  2P,i .To gauge the range of  2P,f , we vary  1P from zero to 90%, with  1P = 0 equal to no evolution.For simplicity, we define the final cluster-mass as therefore neglecting mass loss through stellar evolution and dynamical effects that depend on stellar mass.Fig. 20 shows the final fraction of 2P stars in our extrapolated log-normal wind-mass fraction distributions obtained in Fig. 19, assuming  1P = 90%,  1P = 50% or  1P = 0%.We refrain from changing the 2P fraction for simplicity, however since the 1P is always dominant initially, the loss of 2P stars will mostly result in only a minor change in the final 2P fraction.For instance, in the case where 90% of the 1P is lost, the loss of half of the 2P stars would only decrease the final upper limit 2P fraction of a 10 7 M ⊙ cluster from ∼ 55% to ∼ 38%.We compare in Fig. 20 our extrapolated final 2P fractions to GCs in the Milky Way and in the Magellanic clouds collected in Tables 7, 8 and 9 in Gratton et al. (2019).The masses and 2P fractions for clusters in the Milky Way are from Baumgardt et al. (2019) and Milone et al. (2017a); Zennaro et al. (2019), respectively.The masses of SMC/LMC clusters are from Mackey & Gilmore (2003); Glatt et al. (2011) Martocchia et al. (2018b); Hollyhead et al. (2019).As noted earlier, the observed clusters seem to indicate a correlation between the 2P fraction and cluster-mass, although we should stress that there are notable uncertainties in the observed 2P fractions that are at times based on only a handful of stars.
None of the extrapolated models result in 2P dominated clusters when we use the best-fit log-normal parameters from Fig. 19, even after most of the 1P is lost.Only when using the upper limits of the log-normal parameters are we able to reach a dominant 2P population in clusters with initial mass close to 10 7 M ⊙ .Even then, a significant fraction ( 1P ≫ 50%) of the 1P has to be lost.Compared to the observed GCs, the increase in the extrapolated 2P fraction occurs more gradually with cluster-mass, and the highest 2P fractions are still 20-30 percentage units too low compared to the typical 2P fractions in the local GC population with similar present-day masses.The strong increase in the 2P fraction of observed GCs therefore occurs at lower cluster-masses than in any of our extrapolated relations.This again demonstrated the mass-budget problem, in addition to the directly simulated modest chemical spreads in Fig. 17: the selfenriching winds of massive stars self-consistently sampled from a standard Kroupa IMF may not alone produce enough material to result in observationally consistent 2P fractions in a low-metallicity starburst environment.
In addition to dynamical loss of stars, the 2P fraction might evolve also due to changes in stellar evolution induced by the different chemical composition of the 1P and the 2P stars.For the low-mass stars still alive in present-day GCs, the helium content affects the duration of the phases of stellar evolution (Chantereau et al. 2015).Because the chemical pollution from massive stars is produced through central hydrogen burning, the increasing light-element enrichment correlates with increased helium enrichment.Chantereau et al. (2016) argued that the most helium-enriched 2P stars may finish their nuclear burning faster than their helium-normal counterparts, decreasing the present-day 2P fraction in clusters that start with initially extremely helium-enriched stars.The maximum enhancement of star-to-star helium-mass fraction due to stellar winds in our simulated clusters is ∼ 0.02, not enough to drastically change the evolution of the stars.

Avenues to enhance the final 2P fraction
Our extrapolations are based on the simulated clusters that still, even at 10 5 M ⊙ , do not fully realise the input IMF.If all of the winds from a fully realised IMF up to 500 M ⊙ in a very massive cluster would be retained in the 2P, there should potentially be up to a factor of ten more material to be recycled (compare Figs. 2 and 10).A simple exercise of multiplying the wind-material captured in the cluster stars by a factor of two, four and eight results in highest extrapolated 2P fractions of 68%, 79% and 86%, instead of the 55% estimated directly based on the simulated clusters.As an extreme example, we show in the Fig. 21 how the [Na/Fe] abundances and the extrapolated 2P fractions would better agree with the observed distributions if we had ten times more wind-material recycled in the cluster stars.Given ten times more captured wind-material, the extrapolated clusters with initial mass of 10 7 M ⊙ would host a 73% and 88% 2P fraction based on the the best-fit and upper limit log-normal parameters, at final cluster-mass of 2.9 × 10 6 M ⊙ and 4.8 × 10 6 M ⊙ , respectively.In terms of abundance variations, a boost by a factor of ten would increase the largest  wind values to ∼ 80%.Such enhancement in wind-material would shift the most enriched stars up to values in excess of [Na/Fe]∼ 0.6, resulting in a [Na/Fe] spread of 1 dex that is in agreement with observed GCs.Even then, more than 50% of the 1P stars would still have to be lost.the wind-masses captured in cluster stars used to recover the log-normal fits in Fig. 19 have been scaled uniformly up by a factor of ten while keeping the stellar mass the same.Such enhanced wind-output can be expected e.g. from a top-heavy IMF.The colour-coding in the top panel is the same as in Fig. 7 and the observed data points are the same as in the respective Figs. 17 and 20.
To reach a dominant 2P we would therefore require more winds to be captured in the cluster stars.As long as the clusters experience mass-loss dominated by the number of 1P stars, the clusters do not necessarily need to be dominated by 2P stars at birth.An initially higher 2P fraction would help in producing a higher final 2P fraction.Based on Fig. 10, more massive clusters may be able to retain a larger fraction of the released stellar wind, aiding the formation of strongly wind-enriched stars.Dedicated hydrodynamical simulations with similar methodology as presented here would be required to assess whether the initial fraction of 2P stars in more massive clusters would be even more enhanced than what we are able to infer based on clusters up to 2 × 10 5 M ⊙ in Fig. 19.The other option is to enhance the wind production, as long as the added enrichment can be captured within the star forming enrivonment.In addition to actually sampling the IMF in full, enhancing the wind production by a top-heavy IMF (Prantzos & Charbonnel 2006) has been suggested as a particular solution to the mass-budget problem by the semi-analytic computations using similar wind tables to ours presented in Szécsi & Wünsch (2019).Other avenues for enhanced enrichment could be the addition of other sources, such as interacting binaries, rapidly rotating stars or supermassive stars, occurring coevally at such short timescales before the clusters are completely vacuated of star-forming gas.
We show in Fig. 22 examples of wind production by a top-heavy IMF and fast-rotating, chemically homogeneously evolving stars computed using the BoOST stellar tracks at our fiducial metallicity.For the fast-rotating stars, we use the IZw18CHE tracks from Szécsi et al. (2022) where the stars have an initial rotation velocity of 500 km s −1 .The wind production of the two alternative models is shown relative to the fiducial model shown at the bottom right in Fig. 2. The output of the top-heavy IMF is computed by replacing the Kroupa or Salpeter (1955)-like high-mass slope of −2.3 by a relatively extreme −1.3.Such a top-heavy IMF would result in the production of ten times more wind-material, which according to Fig. 20 would in theory be enough to result in a dominant 2P after significant cluster mass-loss ( 1P > 50%).Such enhanced wind production would occur already during the first couple of million years, driven by the increased number of very massive stars.
The integrated stellar wind produced by a population of fastrotating stars in Fig. 22 is also enhanced compared to the fiducial BoOST models, in total by a factor of 2-3.The enhancement increases as the cluster ages due to the increased contribution of massive stars in the 10-100 M ⊙ range.Fast-rotating stars with a few tens of M ⊙ contribute tens of times more winds compared to their slowly rotating counterparts.The capture of the extra wind-material released at later ages may require an extended star formation period, or consecutive episodes of star formation similar to the AGB enrichment scenario.The chemical abundances on the surface of the fast-rotating stars may however result in different abundance patterns (as in Fig. 17), and the stellar winds of such stars are faster, possibly complicating their capture in the cluster centre (see Szécsi et al. 2022).
Lastly, the wind mass-loss rate could be increased by introducing even more massive stars.The line-driven wind mass-loss rate is expected to increase with stellar mass (Vink 2018), thus stellar mass-growth through collisions or further gas-accretion might lead to enhanced production of enriched material.
The enhanced wind-enrichment would result in enhanced heliumenrichment.Multiplying the recycled wind-masses in the stars by a factor of ten, the helium mass-fraction would increase by 0.1-0.2, which might already start impacting stellar evolution.The response of the ISM and the final 2P fraction to changes in the IMF and/or stellar feedback models may therefore be highly non-linear and would require self-consistent testing in a hydrodynamical setting.The scenarios to alleviate the mass-budget problem outlined above offer intriguing avenues for future studies.

Numerical limitations
Because the stellar wind increments released by single massive stars at a rate of up to 10 −3 M ⊙ yr −1 would require a significant improvement to the gas mass resolution to allow their injection as a continuous stream of new particles, the winds are here injected into existing gas particles.For the stellar wind to be retained in the cluster forming region, the wind injection scheme adopted here relies on there being gas particles to be injected with.Once a cluster has expelled all gas within its vicinity, the stellar winds cannot by construction remain in the clusters, nor can they cool down or contribute to further star formation without the presence of gas particles.The results presented here for enrichment within massive star clusters can therefore be viewed as lower limits.If the stellar wind could be retained and cool into star forming conditions without mixing with the ISM, it may be possible to form stars consisting of 100% windmaterial in the later stages of cluster formation.As noted in Fig. 2, the maximum returned mass via winds would still be only a couple per cent.
The protostar phase, especially for low-mass stars, may last for a significant duration up to Myrs (Dunham et al. 2014).Our simulations do not resolve the protostellar phase with gas accretion and jet-like outflows (see e.g.Grudić et al. 2021).To remain conservative, we only inject energy and metal-enriched material into the gaseous component within and surrounding the star-forming regions.Once a gas particle is defined star-forming, we do not allow further enrichment.If we were to consider the inert time as the protostar phase, the metal-content of the star formation reservoir particles (Section 2.1) could be enhanced by allowing them to accrete wind and supernovamaterial.
As indicated in Fig. 2, stochastic sampling of the IMF can also have a strong impact on the integrated feedback of massive stars in clusters as massive as 10 5 M ⊙ .The massive clusters realised in the current study reach stellar masses up to 2 × 10 5 M ⊙ , and remain susceptible to stochasticity.This is best illustrated by concentrating on the most massive clusters in Figs.10-19: the clusters host very massive stars of up to 379 M ⊙ , but their evolution within the cluster formation sequence governs how well their feedback can actually be recycled.As a result, the most massive cluster is less enriched in wind-material compared to the second most massive one.Simulations reaching higher cluster-masses, from 106 M ⊙ upwards, would be required to guarantee the realisation of a wealth of very massive stars spread across the cluster formation sequence.
Regardless of how the massive stars are sampled and how their stellar winds are implemented in the simulations, the stars need to be present within the cluster forming regions during the period when the wind is released to be able to enrich their surroundings.Due to gravitational softening, our simulations do not take into account the dynamical impact of short-range interactions or stellar multiples.Binary interactions are especially relevant for massive stars which are predominantly found in systems of multiple stars (Sana et al. 2012;Chini et al. 2012).Processes including binary stars may result in the ejection of massive stars out of young star clusters through binary dissolution due to exploding companion stars (Blaauw 1961) or through strong gravitational encounters (Poveda et al. 1967).Observed runaway and walkaway stars have indeed been traced back to their cluster origins (Hoogerwerf et al. 2000;de Wit et al. 2005) even at very young cluster ages (e.g.Gvaramadze & Bomans 2008;Stoop et al. 2023).The removal of the massive stars during the early stages of cluster formation may have a direct impact on their role in enrichment and gas dispersal, however direct N-body simulations indicate that not all massive stars are necessarily ejected (Fujii et al. 2022).Molecular cloud scale simulations that model both star formation and stellar two-body dynamics (Fujii et al. 2021;Grudić et al. 2021;Rieder et al. 2022;Cournoyer-Cloutier et al. 2023) have paved the way and we will hopefully in the future be able to study the impact of accurate gravitational dynamics in cluster forming regions within a full galactic environment.

CONCLUSIONS
We have examined star formation and stellar feedback in a lowmetallicity dwarf galaxy starburst, focusing on the stellar selfenrichment during the formation of massive star clusters.The hydrodynamical simulations employ a Jeans mass dependent star formation threshold, a stellar IMF realised with single stars between 0.08 M ⊙ and ∼ 400 M ⊙ , and detailed prescriptions for stellar winds, radiation, and SNe of individual stars.The evolution models for massive (> 9 M ⊙ ) and very massive (> 100 M ⊙ ) stars include the chemical composition of the stellar winds an SNe, which allows us to trace the propagation of individual elements throughout the starburst.
Massive star clusters up to ∼ 2×10 5 M ⊙ form on timescales of less than ∼ 5 Myr.Massive stars evolve on a few-Myr timescale, allowing their light-element enriched winds to be released and recycled in new stars while the clusters are still forming.Following the individual elements released in stellar winds and SNe of massive stars, we showed how the wind-material is on average ten times more efficiently recycled within clusters compared to SN-material.Very massive stars, when present, contribute most of the stellar wind-mass and die either as direct collapse black holes or very energetic PISNe, inhibiting significant self-enrichment by heavy elements in the vicinity of the rapidly forming star clusters.Some of the simulated clusters show small amounts of heavy element enrichment, with the cluster-to-cluster mean [Fe/H] varying by 0.01 dex.Because the measured mass fraction of SN-material locked in cluster stars is independent of the total cluster-mass, the iron enrichment has to occur on galaxy-scale, rather than locally.The iron variations reflect the global mixing of SN-material within the star forming regions throughout the starburst, rather than localized self-enrichment.These results are consistent with detailed spectral observation of GC stars that indicate small but non-zero spreads in iron, especially in the primordial 1P, in iron-normal 6 GCs that were originally thought to have uniform heavy element contents.
We have taken a detailed look at the light-element enrichment through stellar winds in the most massive star clusters that host very massive stars.The stellar wind-material is radially concentrated at the cluster centre in many of the clusters more massive than a few 10 4 M ⊙ , consistent with some of the observed GCs that host a centrally concentrated chemically peculiar 2P.The total mass of recycled wind-material within the simulated clusters increases with the cluster-mass, reflecting the enhanced wind production and possibly easier capture of wind-material during the formation phase of massive star clusters.This is conceptually consistent with observed GCs, where the number fraction and elemental abundance spreads of enriched 2P stars increase with cluster-mass.
The abundance ratios in the pure, un-diluted stellar winds captured in the simulated 2P cluster stars are consistent with the starto-star [Na/Fe]-[O/Fe] and [Al/Fe]-[Mg/Fe] anticorrelation and the enhanced [N/Fe] observed in practically all GCs.After dilution with the star-forming ISM, the spread of the simulated light-element abundance ratios decreases and we run into a manifestation of the mass budget problem with many of the abundance ratios.Even the most enriched stars show barely extreme enough abundance variations in [Na/Fe], [O/Fe], [Al/Fe] and [Mg/Fe] to be considered as observationally selected 2P stars.The largest fraction of stellar mass composed of wind-material reaches 10% per star, resulting in starto-star spread of up to ∼ 0.3 dex in [Na/Fe].[N/Fe] shows the largest spread of up to ∼ 0.5 dex and N-enhancement in particular has been observed in young high-redshift stellar objects that may represent young GCs in formation.
The wind-mass captured in the stars of each simulated cluster is reasonably well fit with a reversed log-normal PDF.Using the bestfit parameters, we obtained scaling relations for the simulated  wind distribution as a function of cluster-mass.We then extrapolated the parameters of the  wind distributions to obtain an observationally motivated 2P fraction (Δ[Na/Fe] > 0.3,  wind > 10%) of clusters that represent the initial mass range of present-day observed GCs.After accounting for significant mass loss (up to 90% of the 1P) via two-body dynamics and tidal effects, we expect winds of evolved massive stars to result in a final 2P fraction of ∼ 55% at most in clusters up to an initial mass of 10 7 M ⊙ .Observed GCs, on the other hand, show 2P fractions of up to 90% and an onset of the increase in 2P fraction already at 10 5 M ⊙ .Such high 2P fractions at so low cluster masses are at first hand inaccessible to our current model because of the mass budget problem.However, we showed that an increase by a factor of ten in enriched material recycled within the simulated clusters could result in an observationally consistent relation between the 2P fraction and final cluster-mass, as well as observationally consistent spreads in star-to-star [Na/Fe] and [O/Fe].Such enhancement in enriched material may be provided by stellar wind output of stars with a top-heavy IMF, or by a combination of various enrichment sources such as interacting binaries, rotating massive stars or supermassive stars.For maximal wind-production, the high-mass end of the IMF would also need to be well populated.In terms of the wind injection method, the recycled wind mass and the most extreme abundance ratios might change by implementing a method that spawns wind-particles instead of relying on existing gas particles that inevitably dilute the chemical composition of the 2P stars.Even then, an enhanced wind-production rate would still be needed if an initially dominant 2P was to be achieved.
Despite the expected mass-budget problem encountered in our study, we can conclude that the observationally consistent, rapidly generated centrally concentrated 2P distribution and increase of recycled wind-mass with cluster-mass support self-enrichment by pre-SN feedback in cluster forming regions as the source of the 2P.

Figure 1 .
Figure 1.The evolution of the abundance ratios [Na/Fe], [O/Fe], [Al/Fe], [Mg/Fe] and [N/Fe] on the surface of the star in four BoOST stellar models with initial masses of 20 M ⊙ , 50 M ⊙ , 150 M ⊙ and 500 M ⊙ .The black dots indicate the zero age main sequence and the abundance ratios get more extreme values as the stars evolve.The models describe the stellar evolution until core-helium exhaustion, for a total of 1.8 Myr, 2.4 Myr, 3.8 Myr and 8.2 Myr in order of increasing mass.The stellar wind that is injected into the interstellar medium surrounding each star is assumed to have the composition of the stellar surface.

Figure 2 .
Figure2.Stellar wind-material produced by 9-500 M ⊙ stars in a randomly sampled stellar population in 10 4 M ⊙ (left), 10 5 M ⊙ (middle) and 10 6 M ⊙ (right) star clusters, computed using the BoOST stellar tracks for massive stars.The wind production is binned (∼ 0.04 dex per bin) according to stellar mass, each cluster has a uniform age, and we show the wind-material produced until cluster age of 2 Myr (yellow), 4 Myr (red) and 6 Myr (blue).Each cluster has been randomly populated with an IMF between 0.08 M ⊙ and 500 M ⊙ , and we sample in total 10 7 M ⊙ of stars by Monte Carlo sampling the 10 4 M ⊙ , 10 5 M ⊙ and 10 6 M ⊙ clusters 1000, 100 and 10 times, respectively, to showcase the effect of stochastic sampling.The solid lines show the mean wind-mass of each Monte Carlo set at a given cluster age, and the shaded regions show the minimum and maximum value in each mass bin to indicate the range of most extreme values with respect to the mean.The top row shows the wind production at a given stellar mass, and the bottom row shows the cumulative distribution of wind-material with increasing stellar mass.The wind-masses are shown normalized to the total produced wind-material by the time the cluster has aged to 25 Myr (left y-axis) and to the mass of the respective cluster (right y-axis).Very massive stars (> 100 M ⊙ ) dominate the mass-loss during the first couple of Myr, 90% of the wind-material is released by > 50 M ⊙ stars already by cluster age of 6 Myr, and the impact of stochasticity on the minimum and maximum wind-mass released per mass bin reduces with increasing cluster-mass.

Figure 3 .
Figure 3. Left half: the surface density of stars formed during the simulation (top left), the gas surface density (top right), the gas temperature (bottom left) and the thermal gas pressure (bottom right) at the end of the off-nuclear dwarf galaxy starburst.The quantities have been computed within ±500 pc of the midplane.The most massive star clusters in the tidal bridge exceed 10 5 M ⊙ .The old stellar disks of the progenitor dwarf galaxies (not shown) coincide with the concentrations of new stars labeled as galaxy 1 and galaxy 2 in the top left image.The image resolution is ∼ 6 pc.Right half: the evolution of the stellar surface density within a 100 pc patch in the central starburst region, shown in 1 Myr steps across the starburst from left to right, top to bottom.The sequence also follows the formation of the most massive star cluster in the simulation (∼ 2 × 10 5 M ⊙ ), circled in each panel.The image resolution is ∼ 0.2 pc.The timestamps in the panels indicate time with respect to  0 , which is the epoch when the starburst begins as defined in Fig. 4.
the cluster formation efficiency Γ = CFR/SFR, computed over the averaging intervals of 5 Myr and 10 Myr.Comparing the SFR and the

Figure 4 .
Figure 4. Top: star formation rate (black solid lines) and the cluster formation rate (red dashed lines,  cluster > 100 M ⊙ ) of the interacting dwarf galaxies.The rates have been averaged in bins of 1 Myr (thin lines), 5 Myr (thicker lines) and 10 Myr (thickest lines).The vertical dotted line indicates the start of the starburst,  0 .Bottom: the cluster formation efficiency CFR/SFR using the 5 and 10 Myr averaged rates in the top panel.

Figure 5 .
Figure5.The relation between star formation rate surface density and gas surface density in the star forming regions of the simulated dwarf galaxies.The yellow circles show one of the progenitor dwarf galaxies run in isolation over 500 Myr in 50 Myr steps and the squares coloured by time and connected by lines show the progression of the dwarf galaxy merger in 5 Myr steps before the starburst (−48 Myr <  −  0 < 0 Myr) and in 1 Myr steps during the starburst (0 Myr <  −  0 < 7 Myr).The star formation rates are always computed over past 10 Myr and averaged over 400 pc pixels, to match the typical averaging scale in the reference observations.The observed data points are from the THINGS survey (gray dots,Leroy et al. 2008), the LITTLE THINGS survey (orange dots,Elmegreen & Hunter 2015), 5-95 percentile range of the dwarf galaxies in the FIGGS survey (shaded gray,Roychowdhury et al. 2015) and high-z star forming galaxies (up and down pointing triangles, respectively, fromTacconi et al. 2013 andGenzel et al. 2010).The black diagonal line shows the observed fit Σ SFR = 1.6 × 10 −4 Σ 1.4 gas from Kennicutt (1998) that accounts for helium and metals by a correction factor of 1.36.

Figure 6 .
Figure 6.The cluster mass function of clusters younger than 10 Myr, starting at simulation time  −  0 = −43 Myr in 10 Myr steps until  −  0 = 7 Myr.The errorbars indicate the Poisson error and the width of the mass bins.The shaded area shows the region below the average cluster-mass limit of ∼ 25 M ⊙ and below the limit of one cluster per mass bin.The line indicates a power-law relation of   /  ∝  −2 cluster .

Figure 7 .
Figure 7.The star formation histories of all clusters more massive than 10 4 M ⊙ , colour-coded by the cluster-mass.Zero age corresponds to the mean stellar age of each cluster and simulation time progresses to the right.All clusters form on a timescale of less than 10 Myr and the majority of the stellar mass forms over a timespan of 2-3 Myr.

•*Figure 8 .
Figure 8.The initial masses of stars realised by the end if the simulation at  −  0 ∼ 7 Myr.All stars are shown in black, stars in clusters more massive than 10 3 M ⊙ in red, and stars in clusters less massive than 10 3 M ⊙ and stars that are not bound to clusters combined in gray.The binned number counts have been normalized to the total number of stars in each sample.The thin solid and dashed lines indicate the input Kroupa IMF slopes of −1.3 and −2.3 used in sampling the stellar masses during the simulation.The shaded regions show the initial mass range where stars die either as direct collapse black holes (light gray) or where the stars explode as PISNe leaving no remnant (dark gray).

Figure 9 .
Figure 9.The integrated wind and SN-enrichment during the dwarf starburst.Solid lines show the values for all stars in the simulation and dashed lines show the equivalent for stars in bound clusters; blue lines are for SN-material and red lines are for wind-material.The light coloured bars on the right of each panel show the equivalent values for a dwarf galaxy run for 500 Myr in isolation inLahén et al. (2023).The labels in the legends are described in more detail in the text.The time axis has been shifted by the start of the starburst  0 as indicated in Fig.4.Top left: total formed stellar mass, bound star cluster mass, produced SN-mass and produced wind-mass.Top right: wind and SN-material recycled in stars and bound star clusters.Bottom left: wind and SN-material recycled in stars and star clusters relative to the total mass in stars, star clusters and the overall recycled mass.Bottom right: wind and SN-material recycled in stars and star clusters relative to the produced wind and SN-material from all stars and cluster clusters.∼ 10% of the wind and SN-material is recycled in new stars ( wind, * ,cluster / wind,tot ∼  SN, * ,cluster / SN,tot ∼ 10%) and the rest (∼ 90%) resides in the gaseous phase.91% of the wind-material that is recycled in new stars remains predominantly in star clusters ( wind, * ,cluster / wind, * = 91%), while only ∼ 8% of the recycled SN-material is in cluster stars ( SN, * ,cluster / SN, * = 7%).
Figure11.Top: the relative wind-mass produced by clusters more massive than 10 3 M ⊙ , measured in the final snapshot as the integrated wind-output of all massive stars bound in a cluster related to the cluster mass.Bottom: the recycled wind-fraction in each cluster in the top panel, computed as the wind-material locked in cluster stars divided by the wind-material produced by the cluster.The lines in both panels show the best-fit power-law relations ∝   cluster to all clusters (dashed) and to clusters more massive than 10 4 M ⊙ , and the power-law slopes are indicated in the legends.Clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig.7.Even though the majority of the clusters are very young (mean stellar age < 5 Myr), the produced mass-estimate is a lower limit as it excludes escaped stars as well as stars that have exploded as PISN since the PISN events leave no remnants.Due to the lower-limit of total produced mass in the top panel, the retained ratios, especially in the higher cluster-mass range, are upper limits.
Figure12.The SN-mass fraction in clusters more massive than 10 3 M ⊙ , measured within 2 1/2 .Clusters with no feedback material have been placed at 10 −7 and clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig.7.The errorbars show the linear mean and standard deviation of the binned SN-mass fraction to guide the eye at low cluster-masses where many of the clusters have not been enriched by SNe.The SN-mass fraction remains flat with cluster-mass, reflecting SN-material that is uniformly mixed across the galaxy.

Figure 13 .
Figure13.The mean abundance ratios of [Fe/H] for stars bound to clusters more massive than 10 3 M ⊙ , measured within 2 1/2 .The material used to compute the ratio for each star is composed of the SN and wind-mass locked in the stars, diluted by the rest of the stellar mass with the initial ISM composition.The mean [Fe/H] have been shifted by the [Fe/H] ratio of the unenriched ISM composition at the beginning of the simulation.Clusters more massive than 10 4 M ⊙ have been colour-coded as in Fig.7.

Figure 14 .
Figure14.The radially binned wind-mass fraction in bound cluster stars in all clusters more massive than 10 4 M ⊙ .The fraction  wind (  ) is computed as the ratio between the wind mass and the stellar mass in bin   .The radii are scaled by the half-mass radii of the clusters and lines are colour-coded by cluster-mass.Clusters up to a few 10 4 M ⊙ have relatively flat wind-mass fraction profiles, while more massive clusters show increasingly centrally concentrated enrichment profiles.

Figure 15 .
Figure15.The radial wind (solid) and SN (dashed) mass fractions of the bound cluster stars in the 20 most massive clusters.The lines use the same colour scale as in Fig.14and the cluster-mass and half-mass radius  1/2 are indicated at the top of each panel.In most of the clusters the radial SN-mass fraction is flat or increases with radius, as opposed to the wind-mass fraction which shows increasingly centrally concentrated profiles toward higher cluster-masses.The wind-mass fraction and the SN-mass fraction are often anti-correlated in clusters that have centrally concentrated wind-mass profiles.
-star wind mass fraction f wind = m wind /m

Figure 17 .
Figure 17.The correlation between various abundance ratios in stars bound to clusters more massive than 10 4 M ⊙ measured within 2 1/2 .The left column shows the diluted abundance ratios measured as the mix of stellar wind and low-metallicity ISM out of which each star formed.The right column shows the abundance ratios for the recycled wind-material only (see Fig. 1), to showcase the most extreme spread that would be obtained if the majority of the mass of a star was composed of the wind-mixture.The histograms show the distribution of stars with each abundance ratio in 0.1 dex bins.The stars and the histograms have been colour-coded by the total mass of the host cluster.The black dot indicates the zero enrichment composition, as in Fig. 1.The observed data include Galactic GC stars from Carretta et al. (2009b, gray crosses) and Pancino et al. (2017, red stars) and LMC GC stars from Mucciarelli et al. (2009, empy circles) in the two uppermost rows, and from Nataf et al. (2019, and references therein) for the bottom row.All observed data have been limited to clusters with [Fe/H] < −0.9, and in the bottom row to stars with signal-to-noise greater than 100.The 0.3 dex bar in the top left panel shows the typical [Na/Fe] threshold for stars considered to be part of the chemically peculiar 2P, which we reach in the most wind-enriched stars .
Figure19.Left: the probability density function of the stellar mass composed of wind-material in stars of clusters more massive than 10 4 M ⊙ within 2 ×  1/2 , as in Figs.16 and 18.We only show the most massive clusters for clarity.The best-fit log-normal distributions (Eq.3), using inverted log 10 (  wind ) as the x-axis and zero point at  wind = 100%, are shown with the same colour coding as the simulated cluster data from Fig.7.Right: the arithmetic mean −  and standard deviation  (Eqs.4 and 5) of the log-normal fits to  wind in clusters more massive than 10 3 M ⊙ , as well as the power-law scaling relations fit to −  and .The lines show the best-fit power-law slopes (∝   cluster , dashed) and the slope upper limits (slope plus bootstrapped standard deviation, solid) , with values indicated in the panels.The black curves in the left panel show log-normal profiles extrapolated according to the best-fit (dashed) and the upper limit (solid) parameter-scaling in the right hand panels for clusters with mass of 10 6 M ⊙ (thin) and 10 7 M ⊙ (thick).The legend indicates the fraction of stars composed of more than 10% of wind-material in the extrapolated distributions.

Figure 20 .
Figure20.The fraction of 2P stars as a function of the final cluster-mass, using the mass-dependent scaling of the log-normal parameters obtained in Fig.19and factoring in the loss of stars due to cluster evolution in a galactic tidal field.The extrapolated relations assume that either 90%, 50% or 0% (from darkest to lightest opacity) of the 1P stars are lost.The dashed lines show the 2P fraction obtained by integrating over the best-fit log-normal wind-mass fraction distributions, and the solid lines show the upper limit using the best-fit parameters plus standard deviation, as illustrated in Fig.19.The observed Milky Way and Magellanic GCs are from Tables 7, 8 and 9 inGratton et al. (2019), see text for details.

Figure 21 .
Figure 21.The stellar [Na/Fe]-[O/Fe] abundance ratios measured at formation equivalent to those shown in the top left panel of Fig. 17(top) and the final 2P fractions equivalent to Fig. 20 (bottom) in an enhanced enrichment scenario:the wind-masses captured in cluster stars used to recover the log-normal fits in Fig.19have been scaled uniformly up by a factor of ten while keeping the stellar mass the same.Such enhanced wind-output can be expected e.g. from a top-heavy IMF.The colour-coding in the top panel is the same as in Fig.7and the observed data points are the same as in the respective Figs. 17 and 20.

M
SN, * /M SN, tot M wind, * /M wind, tot M SN, * , cluster /M SN, from clusters M wind, * , cluster /M wind, from clusters The cumulative wind-mass released by a 10 6 M ⊙ star cluster that either has a top-heavy IMF with high-mass slope −1.3 (solid lines) or consists of fast-rotating massive stars that evolve chemically homogeneously (dashed lines).The stellar winds have been integrated using the appropriate slow and fast-rotating BoOST stellar tracks at 0.016 Z ⊙ metallicity, see text for details.The wind-masses have been normalized to the output of the standard IMF shown in the bottom right panel of Fig.2and the lines are coloured by cluster age up to 25 Myr.The values indicate the mean of a Monte Carlo sample of 10 clusters.