Star clusters at cosmic dawn I: the emergence of the initial cluster mass function

Introduction

Stars form in hierarchical structures within molecular clouds, with an environmentally dependent fraction remaining in bound star clusters [104] (Lada & Lada 2003; Bressert et al. 2010; Grasha et al. 2017). The cluster formation efficiency, $\Gamma$, measures the fraction of star formation that ends up in bound clusters and is found to rise with the gas pressure and star formation rate surface density of the host [66] (Goddard et al. 2010; Adamo et al. 2015; Adamo et al. 2020a). The distribution of cluster birth masses determines how stellar feedback is concentrated because radiation and supernovae (SNe) from coeval massive stars can overlap in space and time, changing how their energy couples to the interstellar medium (ISM) (Murray et al. 2010; Gentry et al. 2017; Fielding et al. 2018). It also supplies the initial conditions for the clusters that survive to the present day [67] (Pfeffer et al. 2018; Chen & Gnedin 2024).

Moreover, the youngest and most massive of these clusters are bright enough to be observed individually in the early Universe. Indeed, JWST observations combined with gravitational lensing now allow us to probe their formation. Magnifications of tens to hundreds resolve individual star clusters and compact clumps in galaxies from $z \simeq6$ to beyond $z=10$ (Vanzella et al. 2023; Mowla et al. 2024; Adamo et al. 2024; Messa et al. 2025; Fujimoto et al. 2025; Bradač et al. 2025; Abdulro’uf et al. 2025; Nakane et al. 2025). In the Cosmic Gems arc at $z=9.625$ (Messa et al. 2026), for example, Adamo et al. (2024) find five clusters younger than 50 Myr with masses of $\sim10^{6}\,\mathrm{M}_{\odot}$. These objects appear unusually bright for their host galaxies. In local galaxies, the luminosity of the brightest young cluster rises with the star formation rate (SFR) of its host (Larsen 2002; Weidner et al. 2004; Bastian 2008). At high redshift, Pfeffer et al. (2025) find that the brightest cluster candidates in lensed galaxies are $\approx1$–$1.5$ mag brighter and $\approx0.5$ dex more massive than the median prediction of E-MOSAICS, which draws clusters from a power-law ICMF of slope $-2$ with an environmentally set truncation. The brightest cluster at a given SFR is a useful probe of the upper end of the ICMF (Bastian 2008; Adamo et al. 2015), and these redshifts it is one of the few population statistics that can be measured. Why the first clusters should be so bright is an open question. One possibility, the so-called ‘feedback-free’ scenario, is that early star-forming clouds are dense enough to convert much of their gas into stars before feedback can act (Dekel et al. 2023; Li et al. 2024).

In the local Universe, clouds and young clusters have mass functions close to $\mathrm{d}N/\mathrm{d}M \propto M^{-2}$ [104] (Krumholz et al. 2019; Adamo et al. 2020a), as expected from fragmentation of a scale-free, turbulent ISM (Elmegreen & Efremov 1997; Guszejnov et al. 2018). A direct measurement of the cluster mass function at an earlier epoch is consistent with this form, giving a pure power-law slope of $-1.89^{+0.13}_{-0.12}$ (Claeyssens et al. 2026), whereas local mergers and luminous infrared galaxies have shallower slopes, down to about $-1.4$ (Adamo et al. 2020b; Linden et al. 2021). Similar departures occur in simulations of starbursts in merging dwarfs (Lahén et al. 2019) and in dwarf galaxies before reionization, shortly after a starburst [38, 39]. These shallower slopes can arise if the star formation efficiency, or the fraction of stars that remains bound, rises with cloud density or mass. The most massive clouds then convert more of their gas into bound clusters, shifting the cluster population towards higher masses (Fall et al. 2010; Li et al. 2018; Grudić et al. 2023).

However, the ICMF often departs from a pure power law at its high-mass end. In nearby spiral galaxies it is truncated around $\sim10^{5}\,M_{\odot}$ (Gieles et al. 2006; Adamo et al. 2015), whereas in merging and luminous infrared galaxies the truncation can reach $\sim10^{7}\,M_{\odot}$ (Randriamanakoto et al. 2019; Adamo et al. 2020b; Linden et al. 2021). This environmental dependence suggests that the upper mass scale reflects the gas available for collapse. Analytic models set it by the largest mass that can collapse against galactic shear, or by the fraction of that mass that collapses before stellar feedback disperses the gas (Reina-Campos & Kruijssen 2017), consistent with the observed rise of cluster truncation mass with the host galaxy star formation rate surface density (Johnson et al. 2017). At high redshift, the lensed sample of Claeyssens et al. (2026) does not yet establish a statistically significant truncation. As the sample of lensed star cluster candidates grows, it will enable tests of how both the slope and upper mass scale vary with redshift and environment.

The ICMF depends both on the population of star-forming clouds and on how many clusters each cloud produces, with what stellar masses. Conversion of clouds into clusters depends on gas supply, feedback and the fraction of formed stars that remains bound. Idealised cloud simulations predict how star formation responds to cloud surface density (Fall et al. 2010; Grudić et al. 2018; Fukushima et al. 2020; Fukushima & Yajima 2021; Kim et al. 2018). Clouds in a galaxy, however, are irregular and accrete and merge while forming stars (Li et al. 2017; Han et al. 2022). Their properties also depend on the density threshold used to define them (Quinalha et al. 2026). A cloud can also form several clusters, so its contribution cannot generally be described by a single efficiency applied to a fixed initial mass. Feedback can therefore change both the cloud population and the clusters each cloud produces.

In simulations of an isolated dwarf galaxy, Andersson et al. (2024) found that early feedback, and ionizing radiation in particular, suppresses the most massive clusters and leaves a steep power-law ICMF, whereas SN feedback alone yields a broad and flatter cluster mass function. In the same simulations, Quinalha et al. (2026) found that dense clouds ($n > 100\,\mathrm{cm}^{-3}$) had the same mass function irrespective of feedback, while lower-density complexes ($n > 10\,\mathrm{cm}^{-3}$) became more top-heavy with radiative feedback. They concluded that the ICMF is shaped by feedback within each cloud and, for radiative feedback, by its effect on the temperature of the clouds. Whether these conclusions hold at cosmic dawn, when galaxies are gas-rich, compact, metal-poor and assembling through frequent mergers, has not been tested yet.

Answering this question requires simulations that follow both the cloud population and its conversion into clusters in the same galaxy, with feedback as a controlled variable. In large-scale cosmological simulations, sub-grid models such as E-MOSAICS attach clusters to star particles with a prescribed ICMF and formation efficiency [67] (Pfeiffer et al. 2018), and post-processing models infer cluster populations after the fact (Chen & Gnedin 2024), so neither predicts the ICMF. Simulations that resolve individual stars capture the conversion directly, but their computational cost limits the range of galaxies and environments that can be explored (Lahén et al. 2020; Hislop et al. 2021). Cosmological zooms that produce compact, feedback-free bursts build in the efficiency, either through a delay before any feedback, without radiative transfer (Chen et al. 2026), or by imposing a local star formation efficiency of unity (Pascale et al. 2025).

At intermediate resolution, accreting sink particle models for star cluster formation allow the ICMF to emerge self-consistently from gas supply and stellar feedback on the scale of individual clouds (Li et al. 2017, 2018; Brown & Gnedin 2022; Reina-Campos et al. 2025). We adopt such an approach in this work, and the closest precursor to our model is that described in Li et al. (2017). Their cluster particles grow inside spheres of fixed $10\,\mathrm{pc}$ diameter at a rate set by a prescribed efficiency per free-fall time until their own feedback stops them from doing so. These authors find that such an approach gives rise to a Schechter-like ICMF with a low-mass slope of $\simeq-1.8$ to $-2$ and a cut-off that scales with SFR.

In this paper, we present cosmological radiation-hydrodynamics zoom-in simulations of a halo of $M_{\rm vir}\sim10^{10}\,M_{\odot}$ at $z\geq6$ in which star clusters, represented by sink particles, are the units of star formation and feedback. We build on the sink particle-based star formation model of Kang et al. (2025, 2026), in which rapid gas accretion on to sink particles produces intense radiation that disperses star-forming gas clumps, resulting in bursty star formation. Each star cluster/sink particle accretes the amount of gas flowing into it, with no prescribed efficiency per free-fall time and no fixed formation duration. We run simulations at two different resolutions and compare runs with and without stellar radiation at the coarser resolution. We then link star clusters to their natal clouds, and compare the brightest clusters with those observed in JWST lensed galaxies with a similar SFR. We investigate whether the brightest simulated clusters reach the observed luminosities, how stellar radiation changes the clouds and their cluster populations, and what terminates cluster growth.

§ 2 describes the simulations, the cluster sink particle model and the observational sample. § 3 presents the simulated galaxies, the ICMF and the brightest clusters, and § 4 traces the ICMF back to the natal clouds and the growth of sink particles, compares it to other simulations and observations and discusses the limitations of our model. § 5 summarises our conclusions.

Methods

A detailed description of the physical models used can be found in Kang et al. (2025) and Kang et al. (2026) (referred to hereafter as K25 and K26, respectively), the latter of which also details the initial conditions used for our simulations. We summarise the methods and initial conditions here for convenience, alongside the adaptations made to the model to capture the growth of star clusters.

Table 1 lists the simulation parameters. Our high-resolution run (HiRes) adopts the same spatial resolution and dark-matter particle mass as the sink-particle run of K26, allowing a comparison with the original star formation model. The fiducial run (FID) uses one fewer refinement level than HiRes, increasing the minimum cell size by a factor of two. As the accretion radius, sink exclusion distance and density thresholds change with cell size, this comparison tests the combined resolution-dependence of the star cluster model. Finally, the run NoRad has the same resolution and sink prescription as FID but omits stellar radiation, providing a controlled feedback comparison (§ 4.2).

	FID	NoRad	HiRes	K26
$\Delta x_{\min}$ [pc]	11.4–22.8	11.4–22.8	5.7–11.4	5.7–11.4
$n_{\mathrm{clfind}}$ [cm$^{-3}$]	25	25	100	100
Stellar radiation	✓	✗	✓	✓
Cluster-sink model	✓	✓	✓	✗
$m_{\mathrm{dm}}$ [M$_\odot$]	1600	1600	1600	1600
$m_\star$ [M$_\odot$]	$10^3$	$10^3$	$10^3$	$10^3$
$r_{\mathrm{acc}}$	$\Delta x_{\min}$	$\Delta x_{\min}$	$\Delta x_{\min}$	$\Delta x_{\min}$
Type Ia supernovae	✓	✓	✓	✗
Core-collapse yields	LC18	LC18	LC18	fixed

Numerical framework

We perform cosmological radiation-hydrodynamics simulations with the adaptive mesh refinement code RAMSES-RT (Teyssier 2002; Rosdahl et al. 2013, 2015). The hydrodynamics equations are solved using the Harten–Lax–van Leer–Contact Riemann solver (HLLC) (Toro et al. 1994) with a Courant factor of 0.8, while self-gravity is followed with a multigrid Poisson solver (Guillet & Teyssier 2011). Radiative transfer is evolved using the first-moment method with the M1 closure relation (Levermore 1984) and a global Lax–Friedrichs flux function. To limit computational cost, we adopt a reduced speed of light, $\tilde{c} = 10^{-2}c$, well above the propagation speed of ionization fronts in the star-forming regions of interest.

Radiation is discretised into six photon groups with energy boundaries approximately given by 5.6, 11.2, 13.6, 15.2, 24.6, and 54.4 eV, with the highest-energy bin extending to infinity. These groups are chosen to capture molecular hydrogen dissociation, hydrogen and helium ionization, photoheating, and direct radiation pressure. The radiation field is coupled to a non-equilibrium primordial chemistry network for H$_2$, H i, H ii, He i, He ii, He iii, and free electrons (Rosdahl et al. 2013; Katz et al. 2017). We also include a spatially uniform, redshift-dependent ultraviolet background (Haardt & Madau 2012) when computing photoionization and photoheating rates. It is damped before reionization, by a factor $10^{-20(z/8.5)^6}$, so that it is negligible until $z \approx7$ and it is attenuated in each cell by the self-shielding factor $\exp(-n_{\mathrm{H}}/0.01\ \mathrm{cm}^{-3})$. This correction applies to the background only. Radiation from the stars in the simulation is transported explicitly by the radiative-transfer solver, so its absorption is followed cell by cell.

Gas cooling includes both primordial and metal-line contributions. At temperatures $T \gtrsim10^4$ K, metal-line cooling is computed from tabulated rates as a function of density, temperature, and metallicity using CLOUDY-based tables (Ferland et al. 1998). At lower temperatures, fine-structure cooling follows (Rosen & Bregman 1995). We also include dust following a metallicity-dependent scaling for gas with $T < 10^5$ K based on Rémy-Ruyer et al. (2014).

Cosmological initial conditions and refinement

The simulations use the same cosmological zoom-in initial conditions as the reference model of K26. The initial conditions were generated with GENEric (Stopyra et al. 2021), adopting Planck-compatible cosmological parameters [102], namely $\Omega_{\mathrm{m}} = 0.3158$, $\Omega_{\Lambda} = 0.6842$, $\Omega_{\mathrm{b}} = 0.045$, and $H_0 = 67.32\ \mathrm{km\ s}^{-1}\ \mathrm{Mpc}^{-1}$. A dark-matter-only parent simulation was first run in a $(30\ \mathrm{cMpc})^3$ volume with $256^3$ particles. From this volume, a halo with $M_{\mathrm{vir}} \sim10^{10}\ \mathrm{M}_\odot$ at $z = 6$ was selected for resimulation. The high-resolution region extends to three virial radii around the target halo, within which the dark matter particle mass is $m_{\mathrm{dm}} \simeq1600\ \mathrm{M}_\odot$. The bulk velocity of the zoom region was reduced using the method of Pontzen et al. (2021), giving a halo velocity of approximately 16 km s$^{-1}$ at $z = 6$ and thereby limiting advection-related errors in the grid calculations.

The simulation domain is initially divided into $256^3$ cells and is adaptively refined within the zoom region. A cell is refined if the total mass in dark matter and baryons, divided by the cosmic baryon fraction $\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}$, exceeds $\Omega_{\mathrm{dm}} = 1.24 \times10^4\ \mathrm{M}_\odot$. We also enforce refinement when the thermal Jeans length,

is resolved by fewer than 16 cells, where $c_s$ and $\rho_{\mathrm{gas}}$ are the isothermal sound speed and gas density of the cell. All cells inside the sink particle accretion radius are kept at the maximum refinement level. Additional refinement levels are released at $z = 19$ and $z = 9$, corresponding to $a = 0.05$ and $a = 0.1$, to keep the physical resolution approximately constant. As a result, in FID, the minimum cell size remains in the range $\Delta x_{\min} = 11.4$–22.8 pc, and $\Delta x_{\min} = 5.7$–11.4 pc in the HiRes simulation, over the redshift interval studied.

Haloes are identified using ROCKSTAR-GALAXIES (Behroozi et al. 2013) applied to the dark matter particle distribution. We define the virial radius, $R_{\mathrm{vir}}$, as the radius enclosing the mean density $\Delta_{\mathrm{crit}}(z)$ (Bryan & Norman 1998), and the virial mass $M_{\mathrm{vir}}$ as the mass inside it. At early times, the halo centre is taken to be the dark matter density peak. Once the galaxy contains more than $10^5\ \mathrm{M}_\odot$ in stars (including cluster sink particles), we instead define the centre using a shrinking sphere method on the stellar mass distribution (Power et al. 2003), and recompute $R_{\mathrm{vir}}$ about this centre.

Star formation and cluster sink particles

We employ the Pop III prescription of Kimm et al. (2017). Pop III stars form in maximally refined cells where the gas flow is convergent and the gas metallicity $Z_{\mathrm{gas}} \le10^{-6} Z_\odot$ (with $Z_\odot= 0.0134$ Asplund et al. 2009), once the local turbulent Jeans length falls below the minimum cell size. The masses of individual Pop-III stars are sampled between 10 and $10^3\ \mathrm{M}_\odot$ from a Salpeter-like IMF with an exponential suppression below $100\ \mathrm{M}_\odot$ (Salpeter 1955; Wise et al. 2012).

Seeding, accretion and stellar populations

Our Pop II star formation model is based on the sink particle prescription of K25 and K26, but modified so that, depending on their mass, sink particles represent either individual forming star clusters or temporary reservoirs from which collisionless star particles are spawned. Candidate star-forming gas clumps are identified on the fly with the PHEW algorithm (Bleuler & Teyssier 2014; Bleuler et al. 2015), which segments the gas density field into isodensity contours above a threshold $n_{\mathrm{clfind}}$. We use $n_{\mathrm{clfind}} = 25\ \mathrm{cm}^{-3}$ in FID and NoRad and $100\ \mathrm{cm}^{-3}$ in HiRes, following the $\rho_{\mathrm{th}} \propto\Delta x_{\min}^{-2}$ scaling of Eq. (2) below. A gas density peak is kept as a separate clump only if its density exceeds that of the saddle point connecting it to a neighbouring clump by at least a factor of 2 (the relevance threshold), and is merged with it otherwise. The clouds analysed in § 4 are identified separately in post-processing, without performing this saddle point density splitting step (Section 4.1.1). A sink particle is seeded at the density peak of a clump with an initial mass of $1\,\mathrm{M}_{\odot}$ when all of the following conditions hold: the peak cell is maximally refined, the gas is enriched above the Pop III metallicity threshold, the local gas velocity field is converging, the clump is virialised, and no other sink particle lies within a distance of twice the sink particle accretion radius, $2r_{\mathrm{acc}}$, from the peak ($r_{\mathrm{acc}} = \Delta x_{\min}$ in all our runs). In addition, the peak density must exceed

which corresponds to the density $0.5\Delta x_{\min}$ away from the centre of a collapsing isothermal sphere (Larson 1969; Penston 1969; Gong & Ostriker 2013).

Once formed, sink particles grow by accreting gas according to a flux-based scheme. The accretion rate on to each sink particle is calculated from the net gas flux through the surface of a spherical accretion zone of radius $r_{\mathrm{acc}}$,

where $\boldsymbol{v}_{\mathrm{gas}}$ is the gas velocity, $\boldsymbol{v}_{\mathrm{sink}}$ is the sink particle velocity, and $\hat{\boldsymbol{n}}$ is the outward unit vector normal to the accretion surface, $S$. We use the flux-based rate rather than a Bondi–Hoyle–Lyttleton estimate because the latter would require an ambient density and velocity at infinity, which are not well defined in a turbulent, inhomogeneous ISM. To avoid numerical instabilities, the timestep is limited so that a sink particle cannot double its mass in a single fine step and no cell can lose more than 75 per cent of its gas mass (Bleuler & Teyssier 2014). Pairs of sinks whose accretion regions overlap and that are gravitationally bound are merged, with total mass and momentum conserved.

In the original implementation of K25 and K26, a sink particle acts as a reservoir. Whenever it accretes a mass equivalent to a stellar particle ($1000\,\mathrm{M}_{\odot}$ in K26 and this study, $500\,\mathrm{M}_{\odot}$ in K25) this mass is converted into a collisionless star particle that sources radiative and supernova feedback. A sink particle is not destroyed after the formation of the first $1000\,\mathrm{M}_{\odot}$ star particle, but can continue to accrete and form other star particles of the same mass. We adapt this procedure to model (unresolved) star cluster formation. In our implementation, all the star particles that a sink particle forms remain attached to their parent sink particle rather than becoming freely orbiting collisionless particles. The sink particle therefore represents a forming star cluster that carries one or more embedded stellar populations, each with a unique age and metallicity. We retain the term ’star particle’ for these around $10^{3}\,\mathrm{M}_{\odot}$ populations even though they are technically part of a larger star cluster. Radiation, supernovae, and chemical enrichment from these star particles are injected at the position of the parent sink particle. If two sink particles merge, the star particles associated with both progenitors are assigned to the remnant sink particle. As a result, each sink particle constitutes an unresolved, sub-grid star cluster. We do not attempt to follow the internal structure or bound mass of the cluster here as this will be explored in a companion paper. We refer to sink particles when describing their formation and growth, and to clusters when comparing them with observed star clusters or with star clusters in other simulations.

End of formation and retention

We declare that a sink particle’s star formation episode has ended after its local gas supply has remained dynamically insignificant over a long enough period of time. To estimate the length of this period,

we compute the mean density inside the accretion region,

and the corresponding free-fall time,

We then monitor the gas mass we expect the sink particle to accrete within three local free-fall times, $3M t_{\mathrm{ff}}$, where $M$ is the net mass flux on to the sink particle. If this mass stays below an individual stellar particle mass ($10^{3}\,\mathrm{M}_{\odot}$) for at least 1 Myr, the cluster is considered formed and any further gas accretion on to the sink particle is disabled. The duration of this lull is a numerical choice. It determines when the formation phase of a sink is declared over, whereas the end of accretion itself is usually set by the reversal of the gas inflow (§ 4.2). For all but 0.2 per cent of the sink particles, the accretion rate is already exactly zero throughout this Myr. Accretion therefore ends before the formation phase is declared over, with its termination set by the gas flow rather than this numerical criterion. Once this happens, gas accretion on to the sink particle is prevented from restarting at a later time. The one exception is mergers. If either of two merging sink particles is still accreting, the remnant can accrete as well, so a star cluster whose formation had ended can resume growing if it merges with a forming one. Finally, any remaining non-stellar sink particle mass is converted into a final low-mass star particle when cluster formation ends, and this particle is prevented from producing supernovae if its mass is too small to allow a consistent sampling of the IMF.

After the initial formation phase, if the associated stellar mass satisfies $M_{\mathrm{cl}} \ge10^{4}\,\mathrm{M}_{\odot}$, the sink particle is retained as a long-lived cluster particle, with its star particles attached. If $M_{\mathrm{cl}} < 10^{4}\,\mathrm{M}_{\odot}$, the sink particle is destroyed and its star particles are redistributed as independent collisionless particles, with positions drawn from a Plummer sphere of scale radius $\Delta x_{\min}$ centred on the former sink particle, conserving total mass and momentum. This threshold is a bookkeeping choice that limits the number of long-lived sink particles. It does not remove low-mass objects from the recorded birth catalogue, because every sink particle is followed until cluster formation ends and all birth masses are recorded (§ 3.2).

Stellar feedback

The stellar feedback model follows K26, with the same parameter choices unless otherwise stated. We include radiation from young stellar populations, core-collapse supernovae from Pop II stars, Type Ia supernovae (not included in K26), and mass-dependent explosions of Pop III stars. No additional feedback channel is introduced specifically for star clusters. The main change is geometric and temporal. In the cluster-sink model, each sink particle hosts the star particles it forms or acquires via inter-sink particle mergers, and their feedback is injected at the parent sink position. Concentrating the feedback at one position couples time-resolved sink growth directly to the radiation, mass, metals, momentum and energy returned by its attached stellar populations.

Radiative feedback follows Kimm et al. (2017). Each star particle emits photons with age- and metallicity-dependent luminosities in the six radiation groups described above. For Pop II star populations, we use BPASS v2.2.1 stellar population models (Stanway et al. 2016; Stanway & Eldridge 2018), assuming a Kroupa (2001 IMF with lower and upper stellar mass limits of 0.1 and 300 $M_{\odot}$, respectively. Pop III stellar lifetimes and photon yields are taken from Schaerer (2002). Ionizing photons photoionize H i, He i, He ii, and H$_2$ according to the relevant cross sections, while Lyman–Werner photons dissociate H$_2$. Absorbed photons heat the gas through excess photoionization energy and transfer momentum to it through direct radiation pressure.

Core-collapse supernovae are coupled to the gas using the mechanical feedback model of Kimm & Cen (2014); Kimm et al. (2015). All stars with initial masses above 8 $M_{\odot}$ are assumed to explode, corresponding to about one supernova per 100 $M_{\odot}$ of stars formed for our adopted IMF. Each explosion injects $10^{51}$ erg of energy and deposits radial momentum into neighbouring gas cells. The time delays for the supernovae of each star particle are drawn independently from a fit to the cumulative supernova number of a starburst99 population (Kimm et al. 2015). The delay model allows explosions between approximately 3.4 and 44 Myr after star-particle formation. We use its minimum allowed delay of $\sim3.4$ Myr after $t_{\mathrm{first}}$ as a common supernova reference in all three simulations. If the Sedov–Taylor phase of the explosion is unresolved, the injected momentum is set to a terminal value calibrated on resolved supernova remnant calculations, with the terminal momentum depending on the local density and metallicity.

Unlike K26, Pop II core-collapse supernovae return individual chemical species (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) with the mass- and metallicity-dependent yields of the rotating 150 km s$^{-1}$ models of Limongi & Chieffi (2018).

Type Ia supernovae are included as a delayed enrichment channel that couples only thermally to the gas. Delay times follow a $t^{-1}$ distribution between 50 Myr and 13.7 Gyr (Maoz et al. 2012), normalized to $1.3 \times10^{-3}$ events per $M_{\odot}$ of stars formed, and each event injects $10^{51}$ erg of thermal energy and 1.4 $M_{\odot}$ of ejecta enriched in metals with the N100 yields of Seitenzahl et al. (2013).

Pop III stars with masses of 11–20, 20–40 and 140–260 $M_{\odot}$ explode as core-collapse supernovae, hypernovae and pair-instability supernovae respectively, with the same energies and yields as in K26 (Nomoto et al. 2006; Heger & Woosley 2002). Those with other masses collapse directly into black holes.

Cluster photometry and the observational sample

We compute the luminosity of each simulated cluster from the stellar populations it carries, taking into account their ages and metallicities, with the same BPASS v2.2.1 models that drive the on-the-fly radiative transfer in the simulations (§ 2.4). A cluster’s light includes every star it holds, including those brought in by mergers with other clusters. We quote intrinsic, stellar-only rest-frame AB magnitudes (Oke & Gunn 1983) at 1500 Å, $M_{\mathrm{UV}}$, without nebular continuum, emission lines or dust attenuation. For young clusters, the contribution from the nebular continuum is likely non-negligible, and we discuss its potential impact on our comparison to observed clusters in § 4.5.

We compare with the brightest compact sources in ten strongly lensed systems observed by JWST between $z = 4.9$ and 11.5, listed in Table A1. Eight come from Pfeffer et al. (2025), namely RCS0224, the Sunrise arc, the Cosmic Grapes, the three MACS0416 subsystems, the Firefly Sparkle and the Cosmic Gems arc, and we add the Cosmic Spear (Abdurro’uf et al. 2025) and the Misty Moons (Nakane et al. 2025). We consider observed systems for which the brightest source has rest-frame UV photometry and the host galaxy SFR is averaged over a stated time interval, either published or derived in this work from published recombination-line

System	ID	$z$	$M_{1500}$	Age [Myr]	$\log M_{\star,\mathrm{host}}$	SFR$_{100}$	SFR$_{10}$	References
RCS0224	E1	4.88	$-18.94 \pm0.39$	$1$–$20^{e}$	–	–	$8.2 \pm1.4$	1, 2
Sunrise arc	6b	$6.0^{a}$	$-17.10 \pm0.39$	$4^{+2}_{-3}$	$8.5$–$9.3$	–	$6.5 \pm3.5^{b}$	3
Cosmic Grapes	clump 1	6.072	$-18.86 \pm0.17$	$6.9 \pm2.7^{f}$	8.65	2.6	$13.8^{+3.8}_{-0.9}\,{}^{c}$	4
MACS0416	D1	6.144	$-16.26 \pm0.20$	$12^{+0}_{-1}$	$7.5^{d}$	0.35	$0.45^{+0.06}_{-0.04}$	5
MACS0416	T1	6.145	$-15.03 \pm0.25$	$2^{+6}_{-0}$	–	$0.12 \pm0.01$	$0.49^{+0.13}_{-0.10}$	5
MACS0416	UT1	6.145	$-14.73 \pm0.30$	$14^{+57}_{-0}$	–	$0.03 \pm0.01$	–	5
Cosmic Spear	3abc	6.2	$-18.00 \pm0.22$	$5^{+2}_{-1}$	8.64	–	$>12^{+7}_{-1}$	6
Firefly Sparkle	brightest	8.296	$-17.43 \pm0.25$	–	6.8	–	$0.81^{+0.27}_{-0.19}$	7
Cosmic Gems	A.1	9.625	$-16.02 \pm0.62$	$9.2^{+13.5}_{-3.2}$	$>7.4$–$7.8$	$>0.33$	–	8, 9, 10
Misty Moons	ID1 A.1	$11.5^{a}$	$-17.31 \pm0.31$	$42^{+27}_{-19}$	7.57	0.32	$0.60 \pm0.12$	11

fluxes. Appendix § 8 lists the sample and describes how we reduce each observed magnitude to an intrinsic rest-frame value at 1500 Å, as well as how the host galaxy SFRs are estimated.

Results

The simulated galaxies

We first establish the galactic environment in which the clusters form and examine whether the cluster-sink prescription changes the integrated properties and evolution of the galaxy. We compare HiRes, FID and NoRad with K26, which shares our initial conditions and the resolution of HiRes but uses the original non-cluster sink model (§ 2). Table 2 summarises the characteristics of the main galaxy and its cluster population at $z = 6$. We define the central galaxy by a spherical aperture of radius $0.2R_{\mathrm{vir}}$, as in K26, centred as described in § 2.2.

	HiRes	FID	NoRad	K26
$M_{\mathrm{vir}}$ [$10^{10}\ M_{\odot}$]	1.17	1.14	1.18	1.18
$M_{\star}$ [$10^{8}\ M_{\odot}$]	2.23	2.04	3.78	2.15
$Z_{\mathrm{gas}}$ [$Z_{\odot}$]	0.44	0.36	0.64	0.20
$Z_{\star}$ [$Z_{\odot}$]	0.25	0.19	0.51	0.17
SFR$_{100}$ [$M_{\odot}\ \mathrm{yr}^{-1}$]	1.15	0.80	2.80	0.86
$M_{\mathrm{UV}}$ [AB mag]	−18.46	−18.64	−19.60	−18.53
$M_{\mathrm{gas}}$ [$10^{8}\ M_{\odot}$]	2.87	2.55	3.40	3.4
$M_{\mathrm{gas}}/(M_{\mathrm{gas}}+M_{\star})$	0.59	0.58	0.50	0.60
$r_{1/2}$	1.12	0.98	0.92	1.14
$\Sigma_{\mathrm{gas}}$ [$M_{\odot}\ \mathrm{pc}^{-2}$]	36	42	63	41
$N_{\mathrm{cl}}$	1469	385	890	–

HiRes and FID track the stellar mass–halo mass relation of K26 closely (Figure 1), ending within 10 per cent of its stellar mass by $z = 6$. Despite the changes to the cluster, yield and Type Ia models, the stellar mass growth of the galaxy is therefore essentially unchanged, as it is when the spatial resolution of the simulation is increased by a factor 2. Without stellar radiation the galaxy forms more stars at almost every halo mass and ends up with nearly twice the stellar mass of HiRes and FID. All runs lie above the FIRE-2 relation ([83]) and the $z = 6$ relations of empirical models (Behroozi et al. 2019; Moster et al. 2021), by about an order of magnitude in stellar mass at fixed halo mass. This excess indicates a less efficient regulation of galaxy-wide star formation in comparison to these models. Additional feedback channels, including cosmic rays, are possible contributors (Rodríguez Montero et al. 2024), but their effects are not isolated by our experiment.

Image not extracted

HiRes and FID are more metal-rich than K26 at fixed stellar mass, by a factor of $\simeq1.3$–1.5 (Figure 2). The larger metal yield of our enrichment model can account for much of this offset. We adopt the yields of Limongi & Chieffi (2018) for rotating massive stars, whereas K26 use a fixed metal yield (Kang et al. 2026). Our runs thus return about 1.5 times more new metals per unit stellar mass formed. Type Ia supernovae, absent from K26, supply only 2 per cent of the metals returned by $z = 6$ and are unlikely to explain most of the offset. NoRad is more enriched still, because it forms more stars at fixed halo mass. Above $M_\star\approx10^7\,M_\odot$ all runs, including K26, systematically lie within the upper envelope of JWST measurements. We measure the mass-weighted total gas metallicity. Observations give nebular oxygen abundances, converted here assuming solar abundance ratios.

Image not extracted

We estimate SFRs from gas accretion on to sink particles inside the central aperture. Star formation is bursty in all three runs (Figure 3), as commonly found in simulations of low-mass galaxies with strong stellar feedback, where SFRs can easily vary by an order of magnitude over a few Myr (e.g. Sparre et al. 2017). Averaged over 10 Myr, the SFR typically peaks at three to five times its 100 Myr peak value. Such bursts usually coincide with pericentric passages of satellite galaxies and major mergers, after which there are often periods where star formation ceases entirely. The 10 and 100 Myr averages at a given instant can therefore differ widely, which matters when comparing with observations whose SFR tracers probe different time-scales (§ 2.5 and § 3.3). Note that the peak SFR$_{100}$ of our main galaxy lies within the range of the lensed JWST host galaxies at $z \approx6$, although these span a range of masses and environments.

Image not extracted

The main galaxy is gas-rich and compact (Table 2). A few $10^8\,M_\odot$ of gas, i.e. half or more of the baryonic mass inside $0.2R_{\mathrm{vir}}$, lies within a half-mass radius of $\approx1$ kpc, and our star clusters form from metal-poor gas, at a median $\approx0.1\,Z_\odot$. We compare these properties with those of isolated-galaxy simulations in § 4.4.

The initial cluster mass function

Each sink particle represents an unresolved star cluster in formation. We define each sink particle’s birth mass, $M_{\star,\mathrm{birth}}$, as the total stellar

mass formed during the first accretion episode. Mass gained after that episode, through mergers or renewed accretion, is excluded. For the ICMF and subsequent cluster-formation analyses, we select completed first accretion episodes without mergers, requiring the first star particle to form within the central aperture, a measured end of accretion and at least one star particle formed before that time. This gives 1469, 385 and 890 clusters in HiRes, FID and NoRad, respectively. The sample includes low-mass sink particles that are subsequently converted into field stars (§ 2.3). Birth masses measure stars formed, with no correction for unresolved gas-expulsion or dynamical losses, or for the fraction that remains bound. The brightest-cluster comparison instead selects clusters at the time of observation (§ 3.3).

We fit the ICMF above a lower mass threshold $M_{\min}$ with both a power law, $\mathrm{d}N/\mathrm{d}M \propto M^{\beta_{\mathrm{pl}}}$, and an exponentially truncated power law, $\mathrm{d}N/\mathrm{d}M \propto M^{\beta_{\mathrm{tr}}}\exp(-M/M_{\mathrm{c}})$ (Schechter 1976), using the un-binned likelihood of the clusters’ individual masses (Johnston et al. 2017; Messa et al. 2018) and sampling the posterior with EMCEE (Foreman-Mackey et al. 2013) under priors flat in the slopes and in $\log_{10}M_{\mathrm{c}}$, with $10^4 < M_{\mathrm{c}} < 10^{10}\,M_\odot$. We quote posterior medians with $16$th–$84$th percentiles. The fitted slope and cut-off depend on $M_{\min}$, so we scan $M_{\min}$ discretely in each run and pick the lowest value whose Kolmogorov–Smirnov distance is within 3 per cent of the global minimum for the fitted truncated power law. We then test the goodness of fit of both laws with synthetic samples drawn from the fitted ICMFs above $M_{\min}$ (Clauset et al. 2009), and compare one run to another using their own $M_{\min}$ values (see Appendix § 9 for details).

The Schechter fits yield $M_{\mathrm{c}} = 6.6^{+3.8}_{-2.0}$, $6.0^{+3.7}_{-1.9}$ and $6.6^{+4.6}_{-2.2}\times10^6\,M_\odot$ in HiRes, FID and NoRad respectively (Figure 4). The fitted slopes are $\beta_{\mathrm{tr}} = -1.54 \pm0.02$, $-1.34 \pm0.05$ and $-1.60 \pm0.05$. At the highest resolution, the ICMF is therefore shallower with stellar radiation. The Schechter form provides a viable description over the adopted fitting ranges. The pure-power-law fits give $p = 0.015$, $0.0054$ and $0.013$ in HiRes, FID and NoRad, respectively (Appendix § 9).

Image not extracted

The slope of HiRes is steeper than that of FID by $0.19$ and close to that of NoRad, so convergence of the slope is not established. Above a common lower mass threshold of $10^{4}\ {\rm M}_{\odot}$, FID remains shallower than NoRad, although the contrast between these two runs depends on the exact lower fitting limit (Appendix § 9). The fitted high-mass cut-offs of the ICMF overlap, but their uncertainties preclude a precise test of feedback dependence. HiRes forms 156 clusters against 100 in FID above $10^{5}\ {\rm M}_{\odot}$, and 33 against 20 above $10^{6}\ {\rm M}_{\odot}$. We return to this abundance difference in § 4.1.1.

The main-galaxy ICMF is shallower than the slope of about $-2$ typical of young clusters in nearby normal star-forming galaxies ([104] (Krumholz et al. 2019; Adamo et al. 2020a)), and close to the slopes of $-1.4$ to $-1.7$ found in local mergers and luminous infrared galaxies ((Adamo et al. 2020b; Linden et al. 2021)), in all three runs. This contrasts with Andersson et al. (2024), who find that in an isolated dwarf galaxy stellar radiation removes the most massive clusters and leaves a steep power-law ICMF. We discuss this difference in § 4.4, and how stellar radiation acts between clouds and clusters in § 4.1. It is not clear, however, that the ICMF at cosmic dawn should match that of the local Universe, and high-redshift observed cluster populations are sparse. The pure-power-law slope of $-1.89$ measured for young, massive lensed clusters at $1 < z < 10$ ((Claeyssens et al. 2026)) is not directly comparable to the slope below the mass cut-off of our Schechter fits. Their sample is restricted to masses above $2 \times10^{6}\ {\rm M}_{\odot}$, where curvature matters for our fitted $M_{\rm c}$. The local logarithmic slope of $dN/dM$ is $\beta_{\rm tr}-M/M_{\rm c}$, so the exponential factor already contributes about $-0.3$ at their lower value. The observed sample also pools galaxies over a wide redshift range, whereas our primary fit follows one galaxy. A fair quantitative comparison would require the same selection and fitted range. What is better measured at high redshift is the brightest cluster belonging to a galaxy, which samples the upper end of the ICMF, and in the next section we use this statistic to test our simulations.

The brightest clusters and the host star formation rate

The brightest cluster in a galaxy of given SFR is one of the few statistics of the cluster population that can be measured at high redshift ((Bastian 2008; Adamo et al. 2015; Pfeffer et al. 2025)). Figure 5 compares our runs with the observed lensed galaxies listed in Table A1. Every $0.5$ Myr we select the brightest cluster currently inside the central aperture, whatever its age or birth location. Its light includes all attached stellar populations, including those gained through mergers or renewed star formation. We pair it with the SFR of the host galaxy averaged over the preceding 100 or 10 Myr. The tracks give median magnitudes in $0.25$-dex SFR bins containing at least ten $0.5$ Myr instants (i.e. 5 Myr of star formation history across the entire redshift range spanned by the simulations), and therefore weight each SFR by the amount of time the galaxy spends forming stars at this rate. Only 57, 39 and 40 distinct clusters are ever brightest in HiRes, FID and NoRad. Observed SFRs are shown on both panels, as filled symbols for SFR tracers sensitive to $\sim100$ Myr and open symbols for those probing $\sim10$ Myr. Quoted offsets use only the matching tracers. Arrows mark host SFR lower limits. Darker magnitude error bars show photometric and magnification uncertainties, which we do not propagate into the offsets (Appendix § 8).

Image not extracted

To interpret these brightest cluster data, we compute the luminosity of the brightest cluster that we would expect from randomly sampling a given analytic ICMF ((Bastian 2008; Adamo et al. 2015)), which we refer to as the Monte Carlo prediction. In other words, for each SFR bin, we draw clusters between $10^{3}$ and $10^{9}\ {\rm M}_{\odot}$ until their total mass equals the stellar mass formed in the SFR time window, assuming that all stars form in clusters. We assign them ages uniformly over the time window and metallicities drawn from the simulated cluster population, and record the brightest. We use four reference ICMFs, $dN/dM \propto M^{\alpha}$ with $\alpha= -2$ or $-1.5$, either pure or truncated at $M_{\rm c}=7 \times10^{6}\ {\rm M}_{\odot}$, compatible with the fitted Schechter mass cut-offs (§ 3.2). We measure the offset of a sample from each Monte Carlo prediction as the median — either over all the SFRs in a bin for the simulations or over all the systems in that SFR bin for the observations — of the difference between the magnitude of the brightest cluster and its predicted magnitude, and report it in terms of magnitudes above or below the prediction.

Using the nominal $\sim100$ Myr averaged SFR estimates, both the observed and simulated samples lie above the pure $\alpha=-2$ power law ICMF Monte Carlo reference. The brightest clusters of the three runs lie $0.6$–$0.7$ mag above its median prediction, compared with about $1.1$ mag for the six lensed systems with a published SFR$_{100}$. These offsets exclude lens-model systematics and dust-correction uncertainties (see Appendix § 8). The Monte Carlo curves are medians, however, and the brightest cluster of a randomly sampled population is drawn from the sparsely populated upper mass end of the ICMF, so its magnitude varies widely between realizations. For $\alpha= -2$, half of the observed points on the 100 Myr panel and most of those on the 10 Myr panel lie within the 16th–84th percentile range of the brightest cluster, and all of them within the 2.5th–97.5th percentile range. Within that spread the observed points are consistent with an $\alpha=-2$ power law but also, on the 100 Myr panel, with a shallower, truncated ICMF. Such an ICMF predicts a narrower spread, and the same comparison applied to our simulations, whose ICMFs we know, does not recover them consistently. On its own, the brightest cluster at a given SFR therefore cannot clearly distinguish between ICMFs, and the three runs have similar median offsets at this time-scale.

We also note that the simulated tracks end at SFR$_{100}$ values of 1 to 3 $M_{\odot}\,\mathrm{yr}^{-1}$, which are smaller than those of half of the hosts with a measured short-time-scale SFR (left panel). However, when averaged over 10 Myr, the simulated SFRs peak at several times their 100 Myr value (see Figure 3), and the tracks now reach the short-time-scale observations (right panel). The observed hosts move the same way, since every system with an SFR measured on both time-scales has a higher SFR when averaged over 10 Myr than over 100 Myr. This is a signature of bursty star formation, and a galaxy selected to host a bright young cluster is likely to be caught during a burst. On the 10 Myr panel, the simulated brightest clusters of any age lie further above the $\alpha=-2$ median prediction than when the SFR is averaged over 100 Myr, and the tracks lift off the relation at low SFR$_{10}$, whereas the lensed clusters with a measured $\sim$10 Myr SFR scatter about the median curve. Most of the observed clusters are younger than about 15 Myr (see Table A1), so they likely formed during the episode of star formation that the short-time-scale SFR of their host measures, whereas over 100 Myr the same burst is diluted. As a result, observed cluster data points move horizontally towards lower SFRs and end up above the median prediction. As Bastian (2008) sketched for starbursts, a galaxy and its brightest cluster move together along the relation while the galaxy forms clusters. Once the burst ends, the SFR averaged over 10 Myr drops significantly within a few Myr, but the brightest cluster fades more slowly, so the galaxy moves to lower SFR at almost constant cluster brightness, above the predicted relation. Our bursty galaxies spend much of their time in this post-burst state, which is short compared with a 100 Myr averaging window, so the offset is smaller on that longer time-scale. Any comparison with observations therefore depends on how fast bright clusters fade and on the time-scale to which the chosen SFR tracer responds (see § 4.5).

Discussion

We now examine how stellar radiation changes the gas supply and duration of cluster growth, and how these shape the ICMF.

What sets the initial cluster mass function

The ICMF depends on the population of available clouds as well as on the clusters they produce (see § 1). Our runs provide a measurement of both, either with or without radiative feedback. A cloud can form several clusters, so its contribution to the ICMF depends on both their number and their growth. All Pop II stars form through sink particles, but the fraction that remains bound is not modelled (§ 4.6). Our cloud census and natal-cloud sample are also built differently (§ 4.1.2).

The natal clouds

We identify structures in post-processing using phew (Bleuler et al. 2015) above the hydrogen density thresholds used by Quinalha et al. (2026), $n_{\mathrm{H}}>10$, $10^{1.5}$ and $100\ \mathrm{cm}^{-3}$, which are referred to as cloud complexes, big clouds and clouds. We merge all substructure, so each structure is a connected region above its threshold. In the remainder of this paper, ‘clouds’ means structures with $n_{\mathrm{H}}>100\ \mathrm{cm}^{-3}$ unless otherwise stated, and we analyse structures above each density threshold separately.

Figure 6 illustrates the gas structures detected at $z=6$ in the FID and NoRad simulations. Both runs contain a large central complex, while the surrounding gas in NoRad is broken into more structures. The mass functions plotted in Figure 7 sample these structures over the assembly history of the main halo. More specifically, we pool simulation snapshots taken every 10 Myr from $z=12$ to 6, selecting structures whose density peaks lie within $0.2R_{\mathrm{vir}}$. Persistent structures contribute once per snapshot, as in the stacked census of Quinalha et al. (2026). We plot the number of detections per snapshot per dex, including all structures, regardless of how many cells they contain. A block bootstrap over consecutive snapshots accounts for some of the temporal correlation between detections. These mass functions thus

Image not extracted

Image not extracted

describe a time-sampled population, rather than a census of distinct structures.

Figure 7 shows two regimes. Below $10^{6}\,M_{\odot}$ the three runs differ at every density threshold. NoRad holds several times as many complexes per snapshot as FID, and its excess shrinks with increasing density but does not vanish for the dense clouds. HiRes holds more structures than FID at every threshold, most of all among the low-mass dense clouds, which the coarser grids resolve with only a few tens of cells. Above $10^{6}\,M_{\odot}$ the mass functions approach one another, with overlapping bootstrap bands in some bins. There the numbers are small, a few detections per snapshot in the most massive bins and often of the same persistent structure, so we cannot distinguish a feedback effect at the upper end of the mass functions. No single power law describes these upper ends, but slopes fitted above $10^{6}\,M_{\odot}$ are shallower with stellar radiation than without it for the complexes of HiRes and FID and for the dense clouds of FID.

The low-mass counts must be interpreted relative to the resolution scales in Figure 7. At $z=6$, a 32-cell structure at the threshold density has a mass of $4.5\times10^{4}$, $1.4\times10^{5}$ and $4.5\times10^{5}\,M_{\odot}$ for complexes, big clouds and dense clouds respectively, in both FID and NoRad. The corresponding masses in HiRes are one eighth of these values. These scales vary with the physical cell size during the census and are indicative rather than sharp completeness limits. The matched-resolution FID–NoRad comparison isolates the stellar-radiation switch, while comparisons with HiRes also change how well structures are resolved.

Stellar radiation changes the cloud population most strongly below $10^{6}\,M_{\odot}$ and at low density. The finer grid also increases low-mass counts by following fragmentation further before sink particle creation. This contributes to the larger cluster abundance in HiRes, where gas accreted by one sink particle in FID can instead be shared among several. The radiation-on/off comparison at fixed resolution therefore provides the clearest test to isolate the effect of feedback.

In an isolated dwarf, Quinalha et al. (2026) found that complexes become more top-heavy with stellar radiation while the dense cloud mass function keeps its shape. As in Quinalha et al. (2026), we do observe that the feedback contrast weakens at higher density thresholds in our simulations. However, our dense-cloud counts still differ at low mass, while the small number of massive structures prevents a strong test of feedback independence at the upper mass end. These authors’ mass functions are also cut off below $\sim10^{6.4}\,M_{\odot}$, whereas ours extend to $10^{8}\,M_{\odot}$, and the different environments and star formation prescriptions prevent us from attributing the contrast to a unique, well-defined cause.

Where clusters form

We link each sink particle to its birth cloud by tracking its birth position back to the time of the preceding snapshot (typically 2 Myr, and at worst 4 Myr) using a velocity fitted to its early trajectory. If the sink particle’s back-projected birth position lies within an isodensity contour defining a cloud, it is assigned that cloud. Failing that, we assign the cloud that contains the sink particle at the first output after its formation. Clusters born in the same cloud are co-natal, and a cloud that forms at least one cluster is natal. We call a cluster unlinked if neither match succeeds. This can occur when clusters form in gas below the density threshold considered, or when clouds evolve rapidly between successive outputs. Matches use every output and are retained even if the natal cloud lies beyond the $0.2R_{\mathrm{vir}}$ aperture.

Radiative feedback also changes where sink particles form. At $n_{\mathrm{H}} > 100\ \mathrm{cm}^{-3}$, 99, 79 and 43 per cent of the cluster birth sample in HiRes, FID and NoRad are respectively linked to a cloud. Without stellar radiation, more sink particles are therefore not linked to structures selected by this density threshold. The lower natal-cloud matching fraction in NoRad is consistent with sink formation in less dense gas. In NoRad, the colder seeding gas lowers the density required by the sink particle formation criterion, which scales as $c_{\rm s}^{2}/\Delta r_{\rm min}^{2}$ (see equation (2)). The temperature dependence of the seeding criterion provides a second route by which stellar radiation changes cluster formation, alongside its effect on the cloud population.

Converting clouds into clusters

Figure 8 shows the median accretion rate of sink particles, binned by their initial accretion rate $\dot{M}_0$ at $t_{\rm first}$, against the time elapsed since they created their first sub-grid stellar particle, $t - t_{\rm first}$. With stellar radiation turned on, the accretion rate declines from this point and reaches zero before the first supernovae associated with the sink particle go off. Sink particles with higher initial accretion rates generally accrete for longer. In NoRad the accretion rate stays closer to a plateau, across the initial-rate bins, until the first supernovae explode, and then falls abruptly.

Image not extracted

The declining accretion rates motivate an empirical braking law for HiRes and FID, in which the rate of decline increases as the stellar population grows. We first write this law with a general exponent $p$, as:

where $t' = t - t_{\rm first}$, $M_{\star}(0) = m_{\rm p} = 10^{3}\ M_{\odot}$ is the first star particle created, and $\dot{M}_{\star}(0) = \dot{M}_0$. The coupling constant, $K$, absorbs the effects of geometry, gravity and gas supply. We stop the solution at the cluster formation time, $t_{\rm form} \equiv t_{\rm end} - t_{\rm first}$, the first zero of $\dot{M}_{\star}(t')$, where $t_{\rm end}$, should the model accurately describe the simulations, marks the end of the first gas accretion episode onto the sink particle. Here $M_{\star}(t')$ is a smooth approximation to the stellar mass build up in the simulation, which is recorded in the form of discrete stellar particles. Integrating once with respect to time gives:

and the final mass is:

For $p>-1$, the coupling has dimensions $[K] = {\rm mass}^{1-p}\ {\rm time}^{-2}$. At fixed $p$ and with $K$ independent of $\dot{M}_0$, the limit $M_{\star,\rm end}\gg m_{\rm p}$ gives $M_{\star,\rm end}\propto\dot{M}_0^{2/(1+p)}$ and $t_{\rm form}\propto\dot{M}_0^{(1-p)/(1+p)}$. The duration scaling follows by integrating $dt' = dM_{\star}/\dot{M}_{\star}$ from $t' = 0$ to $t' = t_{\rm form}$. Both exponents therefore depend on the adopted braking law. We further retain $p=1/2$ as a simple empirical description of the accretion histories (see Appendix C for a brief discussion of the reasons for this choice), giving the asymptotic scalings $M_{\star,\rm end}\propto\dot{M}_0^{4/3}$ and $t_{\rm form}\propto\dot{M}_0^{1/3}$. The law describes the broad decline of the median accretion tracks in HiRes and FID (see Figure 8), but does not identify which radiative process brakes the inflow. The large-mass approximation fails near the single-particle scale, and the measured duration relation is curved (see § 4.3). We do not fit NoRad, whose sink particles generally accrete until their first supernovae explode, and use FID’s fitted coupling in each initial-rate bin with NoRad’s own initial rate for the reference curves in the bottom panel of Figure 8.

Different sink particles within the same bin stop accreting at different times, so a median track is not the history of a single sink particle. We therefore calibrate $K$ on individual tracks, with $\dot{M}_0$ fixed to the measured rate at $t_{\rm first}$ (see Appendix C for details). Individual-track fits give median couplings of 281 and 260 $M_{\odot}^{1/2}\ {\rm Myr}^{-2}$ in HiRes and FID, with 16th–84th percentile ranges of 172–416 and 154–371 $M_{\odot}^{1/2}\,\mathrm{Myr}^{-2}$ respectively. Similar couplings describe the simulations at the two resolutions. The most massive sink particles scatter more about the model, because fresh gas sometimes reaches them late during their growth (see Appendix C). We therefore regard this as an empirical growth model.

In addition to the cluster formation time, $t_{\mathrm{form}}$, we calculate the accretion-weighted mean time, $\tau_{\mathrm{ave}} = \frac{\int_{t_{\mathrm{first}}}^{t_{\mathrm{end}}}(t-t_{\mathrm{first}})\,\dot{M}\,\mathrm{d}t}{\int_{t_{\mathrm{first}}}^{t_{\mathrm{end}}}\dot{M}\,\mathrm{d}t}$, from the same time origin. Figure 9 shows both time-scales for the cluster birth sample. Cluster median formation times are 0.75, 1.19 and 4.24 Myr, with median accretion-weighted times of 0.25, 0.39 and 1.88 Myr in HiRes, FID and NoRad respectively. Li et al. (2018) observed a similar trend, with mass-weighted formation times shorter than their unweighted counterparts. Most of the growth thus occurs early when simulations account for stellar radiation, and slightly earlier still at the finer resolution. In HiRes and FID, 94.2 and 92.0 per cent of the main-galaxy sink particles, respectively, stop accreting before the minimum supernova delay of $\sim3.4$ Myr after $t_{\mathrm{first}}$, compared to only 20.0 per cent in NoRad. This places the termination of most radiating cluster accretion episodes before their own stellar populations produce a supernova, although supernovae from neighbouring sink particles can already have exploded. With stellar radiation, sink particles also accrete for longer when located in more massive natal clouds. Indeed, the Spearman rank correlation coefficients between $t_{\mathrm{form}}$ and natal-cloud mass are $\rho= 0.55$, 0.65 and $-0.07$ in HiRes, FID and NoRad respectively. More massive clouds end up forming more massive star clusters in every run, but the relation between cloud mass and growth duration changes with stellar radiation, which is critical in mapping the natal-cloud population to the ICMF (§ 4.3).

Image not extracted

Radiative feedback also lowers the stellar mass formed relative to the matched natal-cloud mass (Figure 10). We define

Image not extracted

summing over co-natal sink particles in the main-galaxy birth sample. In HiRes, FID and NoRad the median $\epsilon$ is 0.031, 0.078 and 0.36, with 16th–84th percentile ranges of 0.010–0.082, 0.020–0.19 and 0.10–1.05 respectively. We do not cap this ratio at unity. Values above unity, found for 56 of the 331 natal clouds of NoRad, one of 1090 in HiRes and five of 246 in FID, mean that these sink particles ended up forming more stellar mass than their cloud held in gas when it was matched, consistent with continued gas supply from outside that initial reservoir. $\epsilon$ is therefore not the conversion efficiency of a closed system. Note that $\epsilon$ is only calculated for the cloud-linked sample, which includes a smaller fraction of clusters in NoRad. In every run, $\epsilon$ rises with cloud mass but, as more massive clouds host more sink particles, this does not mean that a single sink particle is more efficient at forming a cluster because it is located in a more massive cloud (see Appendix C). The factor of a few between FID and NoRad at fixed cloud mass is the integrated effect of the shorter accretion episodes of Figure 9, and the factor of two to three between FID and HiRes follows from the lower initial accretion rates of HiRes at fixed cloud mass (see Appendix C), with the same coupling. The natal clouds marked for hosting the brightest young clusters in Figure 10 span the measured range of cloud masses and, with stellar radiation, lie mostly above the median $\epsilon$ of all natal clouds.

Stellar radiation and the end of accretion

The comparison with NoRad shows that stellar radiation significantly shortens these accretion episodes. We compare the expected roles of

photoionization heating and direct radiation pressure using pressure-balance radii (Rosdahl & Teyssier 2015; Kimm et al. 2017; Han et al. 2022), given by:

for photoionization heating and

for direct radiation pressure. In these equations, $Q$ is the stellar ionizing photon rate and $L_{\mathrm{abs}}$ the absorbed luminosity. We adopt an ion temperature $T_{\mathrm{ion}} = 10^{4}$ K and the case-B recombination coefficient $\alpha_{\mathrm{B}} = 2.59 \times10^{-13}\ \mathrm{cm}^{3}\,\mathrm{s}^{-1}$, with $k_{\mathrm{B}}$ the Boltzmann constant and $c$ the speed of light. At each available snapshot between $t_{\mathrm{first}}$ and $t_{\mathrm{end}}$, we calculate individual sink particle radiative outputs by summing over the stellar populations they have already formed, taking into account masses, ages and metallicities, and interpolating from BPASS SEDs. We set $L_{\mathrm{abs}}$ to the six-group luminosity emitted in the simulations as a full-absorption limit.

The definition of the confinement pressure, $P_{\mathrm{conf}}$, requires care at our resolution. Kimm et al. (2017) estimate ambient thermal pressure from neighbouring cells, while Han et al. (2022) include thermal, turbulent and gravitational terms across a resolved cloud core. Our accretion region spans only one finest-cell radius, so it does not generally allow us to separate the ambient cloud from the gas photo-heated by the sink particle. Including that heated gas pressure as a source of confinement would be counting the same thermal feedback as both the driving pressure and the resistance to expansion. Neglecting ambient thermal pressure (see Figure 6 of Han et al. (2022) and discussion at the end of this section for a justification), we therefore use $P_{\mathrm{conf}} = P_{\mathrm{grav}} + P_{\mathrm{ram}}$, with:

We sample the sink particle-centred sphere of radius $r \leq r_{\mathrm{acc}} = \Delta x_{\min}$ on a uniform lattice of spacing $\Delta x_{\min}/4$, assigning to each point the density, $\rho$, velocity, $\mathbf{v}$, and total gravitational acceleration, $\mathbf{g}$, of its enclosing AMR leaf cell. Both terms are averaged over the full selected volume. Here $\mathbf{r}$ denotes position relative to the sink particle and the relative radial velocity is thus $v_{r} = (\mathbf{v} - \mathbf{v}_{\mathrm{sink}}) \cdot\hat{\mathbf{r}}$. We exclude outflowing gas when calculating $P_{\mathrm{ram}}$, while the average remains over the full sampled volume. This volume-averaged inward kinetic stress provides an effective resistance to expansion.

A pressure-balance radius exceeding $r_{\mathrm{acc}}$ indicates that the channel can pressure-balance the whole accretion region under the assumed confinement. For a given $Q$, a larger ionized volume has a lower equilibrium density and therefore lower thermal pressure. The radius $r_{\mathrm{PH}}$ estimates how far photoheated gas could expand before its thermal pressure balances $P_{\mathrm{conf}}$, rather than how much of the inflow is actually ionized. Similarly, $r_{\mathrm{DP}}$ characterises the region over which direct radiation pressure can overcome the confining pressure. When $r_{\mathrm{DP}} < r_{\mathrm{acc}} \leq r_{\mathrm{PH}}$, photoheating alone satisfies this condition across the whole accretion region. If both radii remain below $r_{\mathrm{acc}}$, neither mechanism alone does so. The two mechanisms also operate on different time-scales. Direct radiation pressure transfers momentum upon absorption, while photoheated gas expands over an ionized-gas sound-crossing time. Direct radiation pressure may therefore be particularly effective at rapidly suppressing accretion once $r_{\mathrm{DP}}$ reaches $r_{\mathrm{acc}}$.

Figure 11 shows a mass-dependent change in the expected shutdown mechanism. Below $10^{4}\,\mathrm{M}_{\odot}$, photoheating can typically pressurise the accretion region early and retains the greater reach as direct radiation pressure strengthens. At $10^{4}$–$10^{5}\,\mathrm{M}_{\odot}$, neither channel typically reaches the accretion scale early, but both do later. Above $10^{5}\,\mathrm{M}_{\odot}$, the photoheating radius remains below the accretion scale, while direct radiation pressure typically reaches it late in the formation stage. This favours direct radiation pressure as the accretion shutdown mechanism for the most massive sink particles. The evolution throughout the cluster formation phase mainly shows the accretion scale being crossed while the ordering of the pressure-balance radii persists. The coarse snapshots sample an ensemble of star clusters at different fractions of their formation time, so they do not track individual accretion shutdowns or separate the causal contributions of the two radiative processes to these shutdowns. The measurements presented in Figure 11 only include 570 of 1469 HiRes clusters and 197 of 385 FID clusters, as they require that at least one simulation snapshot falls during the first accretion episode. This requirement favours longer-lived accretion episodes when the formation time of a cluster is shorter than the snapshot spacing.

Image not extracted

For a cold uniform medium, the same argument links the pressurisation to both stellar luminosity and ambient gas density (see Appendix D1 for a derivation of the scalings). Assuming a density of $n_0 = 10^{3}\,\mathrm{cm}^{-3}$ and a temperature of $T_0 = 100\,\mathrm{K}$ for the cold ambient medium, our zero-age BPASS output gives photo-heating and direct radiation pressure threshold masses of about 8–470 and 330–5200 $\mathrm{M}_{\odot}$ respectively across cell sizes of 5.7–22.6 pc. These values are low because this benchmark assumes that the expansion is only opposed by the thermal pressure of the cold ambient gas. In fact, according to this calculation, photoionization heating on its own should be enough to prevent any cluster more massive than $m_{\mathrm{p}} = 10^{3}\,\mathrm{M}_{\odot}$ from forming. However, gravity and radial inflow provide a much stronger confinement of the gas around our sink particles, so these reference masses do not predict a shutdown of their accretion. This illustrates why the thermal pressure of the cold and dense ambient medium alone cannot account for the confinement and supports our approximation of neglecting it when estimating $P_{\mathrm{conf}}$. Note that values of cluster masses below $m_{\mathrm{p}} = 10^{3}\,\mathrm{M}_{\odot}$ also extrapolate our assumption of properly sampling a stellar population beyond its domain of validity.

A further concern is that, unlike in our unresolved clusters, well-resolved cloud cores allow radiation to escape through low-density gas channels, and that opposing momentum input from spatially distributed stellar sources within a forming cluster reduces the net radial gas momentum about the cluster centre ((Raskutti et al. 2016; Kim et al. 2018)). Nevertheless, compact clusters such as the ones formed in our simulations can pass through an early stage dominated by direct radiation pressure before photo-heating takes over as gas is dispersed ((Han et al. 2022)). We also chose to neglect infrared trapping, which would add radiation pressure through dust reprocessing. This is because a large boost is unlikely at our low metallicities, consistent with the small infrared optical depths measured in the simulations of metal-poor cloud collisions of Han et al. (2022). Appendix D (more specifically Figure D1) charts the time evolution of the luminous output of our clusters and the physical state of the gas in their host cell throughout cluster formation in more detail.

Image not extracted

From the cloud to the cluster mass function

We test the connection between gas supply and the ICMF by assigning each main-galaxy sink particle of each run a final mass from equation (8), using its measured initial accretion rate and drawing a value of $K$ from the distribution of couplings obtained by fitting the very same sink particle data. For the NoRad run, we use its measured initial accretion rates and draw $K$ from the fitted distribution in FID to predict a final mass. Figure 13 compares these mass functions with the measured ICMFs and natal-cloud distributions. Because the couplings are calibrated on the same simulations, this tests descriptive consistency rather than independent prediction.

Image not extracted

This procedure reproduces the numbers of clusters with masses above $10^{5}\,\mathrm{M}_{\odot}$ measured in the simulations to within about 20 per cent in HiRes and FID, but underestimates cluster counts at the high-mass end, particularly in HiRes. This is because rapidly fed sink particles can receive fresh gas late in their growth, producing more scatter than the fitted distribution of $K$ captures (see Appendix C). Applying FID’s cluster growth model to NoRad reduces its number of clusters with masses above $10^{5}\,\mathrm{M}_{\odot}$ from 210 to 80. At fixed natal-cloud mass, NoRad’s initial accretion rate on to sink particles are slightly lower than FID’s, so its larger cluster masses arise from longer accretion episodes.

Figure 12 shows duration against initial accretion rate. In the radiation runs, rapidly fed sinks also accrete for longer. Full-range log–log slopes are 0.41 $\pm$ 0.02 in HiRes and 0.35 $\pm$ 0.02 in FID, but the relation is curved. The one-third asymptotic prediction for the $p = 1/2$ empirical model is most useful over $10^{4}$–$10^{5}\,\mathrm{M}_{\odot}\,\mathrm{Myr}^{-1}$. The model curve bends below this scaling at low rates because the finite mass $m_{\mathrm{p}}$ of the first stellar particle cannot be neglected in equation (8). When little mass is subsequently accreted, the first particle already supplies much of the feedback, shortening the remaining growth. The discrete stellar mass scale therefore breaks the model’s scale-free behaviour at low rates. The measured relation also flattens at high rates. In NoRad the duration of the accretion episode is nearly independent of the initial rate and remains close to the first-supernova explosion time delay.

Image not extracted

This difference changes the ICMF slope. If $\mathrm{d}N/\mathrm{d}\dot{M}_0 \propto\dot{M}_0^{\beta_M}$ and $M_{\star} \propto\dot{M}_0^{\kappa}$, then $\beta_{\mathrm{p}} = (\beta_M + 1)/\kappa- 1$. Note that this relation neglects any scatter in cluster growth. Above the adopted mass limits, the fitted mass exponents are $\kappa= 1.29$ and 1.28 in HiRes and FID, compared with 0.83 in NoRad. With stellar radiation, faster feeding is associated with longer growth, so the final mass rises more steeply than linearly with the initial rate and the ICMF becomes shallower. In NoRad, the nearly rate-independent duration is generally limited by the first-supernova delay. Faster feeding then increases mass without also extending the available growth time, giving a mapping that is closer to linear, though still sublinear, thus leading to a steeper ICMF. Radiation changes both the duration of accretion and its dependence on the gas supply. This explains how it can shorten cluster formation overall while producing a shallower mass function. The slopes implied by the empirical growth model agree with the measured power-law slopes to within 0.04 in the radiation runs (see Appendix C). The larger number of low-mass clouds forming sink particles in NoRad also contributes to its steeper ICMF (§ 4.1.1), so growth histories and the cloud population are not cleanly separated.

The growth model does not establish what sets the ICMF upper mass cut-off. Its fitted value is sensitive to the behaviour of a small number of massive clusters. Some rapidly fed sinks stop earlier than predicted and have weaker inflows on larger scales, suggesting that external gas supply also matters at the top end of the cluster mass function. The host galaxy may therefore contribute to setting the upper mass scale of the ICMF, although our fitted mass cut-offs are not accurate enough to rule out an impact from radiative feedback. Observed truncation masses rise with star formation rate surface density ((Johnson et al. 2017; Adamo et al. 2020b)), while analytic models link them to gas surface density, shear and the time available for collapse before feedback (Reina-Campos & Kruijssen 2017).

Testing these dependencies requires a larger sample of galaxies with different masses and environments.

Comparison with other simulations and observed cluster populations

In an isolated dwarf, Andersson et al. (2024) find that early feedback, especially ionizing radiation, suppresses massive clusters and steepens their mass function. Our comparison between FID and NoRad instead shows that the ICMF becomes shallower with stellar radiation. Their clusters assemble from distributed $100\,\mathrm{M}_{\odot}$ star particles formed at a prescribed efficiency per free-fall time, and many low-mass clusters continue forming beyond the first possible supernova. Our sink particles accrete surrounding gas according to how much of this gas flows into them and concentrate all their radiation emission in one spatial point. The differing growth histories and source geometry may contribute to the contrast, alongside the galactic environment. Our galaxy also forms substantially more stellar mass at fixed halo mass than the empirical relations shown in Figure 1. Its shallow ICMF prediction should therefore not be treated as that of the typical galaxy at cosmic dawn.

Li et al. (2017) also obtain an emergent Schechter-like ICMF in cosmological simulations, with a truncation mass at the top end that increases with SFR. Their cluster particles grow within spheres of fixed 10 pc diameter at a prescribed efficiency per free-fall time until feedback interrupts the gas supply. Our growth rate instead follows the inflow through the accretion zone. In both approaches, the mass function emerges from cluster growth, but its relation to the available gas depends on the accretion prescription.

Our truncation masses fall within the range seen in nearby mergers. The HiPEEC sample spans $M_{\mathrm{c}} \sim10^{6}$ to a few $10^{7}\,\mathrm{M}_{\odot}$ (Adamo et al. 2020b), whereas nearby spirals typically truncate near $10^{5}\,\mathrm{M}_{\odot}$ (Gieles et al. 2006; Adamo et al. 2015). The generally steeper ICMFs of these galaxies may reflect different environments or feedback channels omitted in our modelling, such as cosmic rays. Growing lensed samples will also test the slope of the ICMF below the truncation mass at high redshift (Claeyssens et al. 2026) in the near future.

Comparing with JWST lensed clusters

Our brightest clusters reach the UV luminosities of the high-redshift JWST lensed sample at a given nominal SFR averaged over $\sim100$ Myr for the host galaxy (§ 3.3). The similar offsets of the three runs despite their different ICMFs show that this statistic does not uniquely recover the birth mass function. Sampling, burstiness and stellar populations gained after the first episode also affect the brightest source.

The comparison is sensitive to how clusters and host SFRs are measured. Magnification affects both luminosity and SFR, and lens-model systematics are excluded from our quoted offsets. Unresolved neighbours can blend into brighter sources, while detection limits favour hosts with bright clusters (Pfeffer et al. 2025). Our sink particles do not have resolved sizes, so we cannot account for observational size selection. The host SFR tracers also average star formation over different periods. UV light responds to recent star formation history rather than to a fixed 100 Myr window (Flores Velázquez et al. 2021), while recombination lines probe younger stars and SED-based SFRs depend on the detailed assumptions made about the star-formation history (Choe et al. 2026).

Our simulated magnitudes also account for stellar light alone, whereas young observed clusters can have appreciable nebular continuum at 1500 Å (Katz et al. 2025). Assuming that all of a sink particle’s ionizing photons produce nebular emission at $10^{4}\,\mathrm{K}$ gives a brightening of a few tenths of a magnitude in the first few Myr, falling below a tenth after 10 Myr. Escape and dust absorption reduce this upper limit. Correcting the youngest observed sources for nebular light would make their stellar magnitudes fainter and move them towards the $\alpha=-2$ reference power law ICMF.

Monte Carlo predictions derived from specific ICMFs assume constant star formation with all stars forming in clusters. A lower cluster formation efficiency would reduce the expected luminosity of the brightest cluster, and burstiness in the SFR changes the sampled population. With one simulated halo and ten observed lensed systems, the comparison establishes that our model can produce sufficiently bright clusters. It cannot discriminate between different ICMFs or feedback prescriptions.

Limitations and excluded physics

Pre-supernova feedback

We include photoionization heating and direct radiation pressure as pre-supernova feedback, but omit stellar winds (Geen et al. 2020; Lancaster et al. 2021), protostellar jets (Federrath et al. 2014), infrared trapping (Menon et al. 2022) and Lyman-$\alpha$ pressure. The latter can matter in metal-poor gas, depending on its structure and kinematics (Kimm et al. 2018; Nebrin et al. 2025). The geometry and resolution dependence of the direct radiation pressure estimate is discussed in § 4.2.

Radiation has a stronger effect than stellar winds on cluster formation in the dwarf simulations of Andersson et al. (2024), and Quinlahan et al. (2026) find that stellar winds also have little additional effect on their cloud mass function. These comparisons do not determine the importance of these omitted feedback channels in our cosmological galaxy simulations, where they could alter both the cloud population and the local accretion.

Sink particles and a single halo

The cluster seeding threshold depends on gas temperature and resolution through $c_{\mathrm{s}}^{2}/\Delta x_{\min}^{2}$, affecting where sink particles form. Their accretion zone spans one cell radius. The flux scheme reproduces accretion on to a collapsing isothermal sphere and is not strongly sensitive to the chosen accretion radius (Kang et al. 2025), but resolved gas structure still affects local growth and feedback. While sink particles that end up with masses below $10^{4}\,\mathrm{M}_{\odot}$ become field stars after formation and remain in our birth-property analyses, mergers between sink particles occur when they become bound and their accretion zones overlap, potentially restoring accretion on to clusters that have already formed. Thus numerical condition need not imply a physical cluster merger. The accretion rate is usually already zero before our numerical cluster formation shutdown criterion is met (§ 2.3), but disabling subsequent accretion also prevents later growth that might otherwise resume without a merger.

A sink particle is a point mass which carries various stellar populations, so we do not follow internal cluster dynamics, gas-expulsion losses or tidal disruption (Li et al. 2019; Brown & Gnedin 2022). Our birth masses measure stars formed rather than bound mass. Extending the simulations beyond this single halo and to lower redshift will test environmental trends and connect these birth populations to surviving clusters.

Conclusions

We have studied star cluster formation in cosmological zoom-in radiative hydrodynamics simulations of a halo reaching $\sim10^{10}\,M_{\odot}$ at $z = 6$. Sink particles, representing forming star clusters, accrete gas and retain the stellar populations they form without recourse to a prescribed star formation efficiency per free-fall time or fixed accretion duration, so the ICMF naturally emerges from gas supply and stellar feedback in the simulations. We compare simulations of the same halo which include stellar radiation at two different spatial resolutions and a run without stellar radiation at the fiducial resolution. All three include supernova feedback. Our main findings are as follows.

• Schechter fits to the ICMF give truncation masses of about $6 \times10^{6}\,M_{\odot}$, with overlapping posterior distributions and substantial uncertainty. At the fiducial resolution, stellar radiation makes the low-mass power-law slope of the ICMF shallower, from $-1.60$ to $-1.34$. The higher-resolution run, HiRes, yields an intermediate power-law slope of $-1.54$ and produces about 1.6 times as many clusters above $10^{5}\,M_{\odot}$ as its lower-resolution counterpart FID. The low-mass slope of the ICMF and the abundance of star clusters therefore remain sensitive to resolution.

• Stellar radiation reduces the abundance of low-mass clouds, most strongly at low density, and raises the sink particle-seeding density threshold by heating the surrounding ISM gas. The small number of massive clusters (above the truncation mass) produced in the simulations prevents us from reaching robust conclusions regarding the effect of radiative feedback at the upper end of the ICMF.

• In each radiative run, about 92–94 per cent of the analysed main-galaxy clusters stop accreting before the earliest possible supernova from their own stellar populations goes off, with median accretion durations of about a Myr. On average, more rapidly fed clusters accrete for longer, making their final masses rise faster than linearly with the initial gas accretion rate. Without stellar radiation, accretion typically persists beyond the minimum supernova time delay, with a duration only weakly dependent on the initial accretion rate, and the mapping from initial gas accretion rate to final stellar mass is closer to linear. Under the adopted assumptions about pressure confinement and full absorption of stellar luminosity, pressure-balance estimates are consistent with photoionization heating limiting low-mass cluster growth and direct radiation pressure becoming important late in the growth of the most massive clusters. An empirical braking growth law describes much of the mapping from accretion rate to cluster mass in the simulations, but does not fully capture the high-mass tail behaviour or establish the origin of the mass truncation.

• Our brightest simulated clusters reach the UV luminosities of observed JWST lensed clusters hosted by galaxies with similar SFRs. The comparable offsets in cluster brightness relative to fixed ICMF predictions across our runs demonstrate that this statistic alone cannot distinguish the underlying ICMFs. It is affected by sampling, burstiness of star formation histories and potential mergers between clusters.

Despite the limitations of our analysis – omitted stellar feedback channels, cluster birth masses that measure stars formed rather than the fraction remaining bound, and a galaxy sample too limited to test the environmental dependence of the ICMF – our simulations reveal a key link between stellar radiation and the shape of the ICMF at cosmic dawn. Following cluster dynamics and mass loss should allow us to connect these early populations to the fossil record that has survived to the present day. We plan to improve our modelling to address these issues in future work.

Acknowledgements

FST acknowledges support from the Science and Technology Facilities Council (STFC) for a PhD studentship and from Yonsei University for hosting and support during a research visit. TK and CK were supported by the National Research Foundation of Korea (RS-2022-NR070872 and RS-2025-00516961). This work used the DiRAC Data Intensive service (DiLa3) at the University of Leicester, managed by the University of Leicester Research Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Leicester was funded by BEIS, UKRI and STFC capital funding and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure. This work used the DiRAC Memory Intensive service (COSMA8) at Durham University, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Durham was funded by BEIS, UKRI and STFC capital funding, Durham University and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure. The authors thank Jonathan Patterson for maintaining the Glamdring Cluster hosted by the University of Oxford, where part of the data processing was performed. We thank the Astrophysics Data Service and the arXiv preprint repository for providing open-source software and services that were used extensively in this work.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

References

The observed comparison sample

Table A1 lists the lensed systems used in § 2.5. We convert the published photometry of each brightest compact source to an intrinsic rest-frame magnitude at 1500 Å, correcting for magnification where required, for dust using the source’s measured reddening and the Calzetti et al. (2000) law, and for rest wavelength using the measured UV slope. The rest wavelength correction is at most 0.4 mag and the dust correction at most 1.3 mag. Quoted errors include photometric and statistical magnification uncertainties and the wavelength correction, but exclude dust uncertainty and lens-model systematics. The latter span 0.3–2.4 mag and are one-sided towards fainter magnitudes for the Sunrise and Cosmic Gems arcs. The Sunrise magnification is a lower limit, while the Cosmic Gems value we use is the lowest of three lens models. Magnification affects both cluster luminosity and host SFR estimates.

The ∼100 Myr averaged SFRs come from the published UV or SED estimates. For the ∼10 Myr SFR averages we use published recombination-line fluxes and magnifications with the Kennicutt & Evans (2012) Hα calibration for four systems, and the published recent SED-based SFRs otherwise. Incomplete host coverage gives lower SFR limits, including the Cosmic Spear’s single-clump estimate. Host masses and cosmologies are retained as published.

Fitting the initial cluster mass function

Numerical effects may alter the low-mass ICMF, so we restrict the fits to masses above $M_{\min}$. We scan the lower mass fitting limit and pick the lowest $M_{\min}$ whose Kolmogorov–Smirnov distance is within 3 per cent of the minimum for maximum-likelihood Schechter fits (Figure B1). Here $M_{\min}$ is an estimated fitting threshold, rather than a prescribed numerical completeness limit. The adopted mass limits are $3.3 \times10^3$, $9.0 \times10^3$ and $2.3 \times10^4\,M_\odot$ in HiRes, FID and NoRad respectively. Both the power law and the Schechter function are fitted above this Schechter-selected limit in each run.

Image not extracted

We assess goodness of fit with 10 000 synthetic catalogues for each law, adapting the method of Clauset et al. (2009). Each catalogue has the size of the full sample. Each mass is drawn from the maximum-likelihood law above the adopted $M_{\min}$ with probability equal to the measured fraction of clusters above that limit, while masses below the limit are resampled from the measured masses. We then repeat the lower-limit scan and parameter fitting on each catalogue, using the same Schechter-based 3 per cent rule relative to that catalogue’s own minimum. The $p$-value is the fraction of synthetic catalogues whose Kolmogorov–Smirnov distance for the law being tested is at least as large as the measured one. This accounts for the selection of the fitting range.

Table B1 presents the results. At our adopted 5 per cent rejection threshold, the Schechter fits are not rejected in any run, while pure power laws are rejected in all three. These goodness-of-fit tests do not directly measure the relative preference between the two fitted laws.

Run	Schechter	Power law
HiRes	0.072	0.015
FID	0.847	0.0054
NoRad	0.496	0.013

Figure B2 tests fit sensitivity to $M_{\min}$. Power-law fits steepen as the lower mass limit rises, while Schechter fits trade slope against truncation mass. FID remains the run which yields the shallowest slope above a common $10^4\,M_\odot$ mass limit, but the slope contrast with NoRad is smaller than when using separately selected mass limits. In NoRad the fitted truncation mass varies by more than an order of magnitude across the scanned lower mass limits. The posterior intervals quoted in § 3.2 therefore describe uncertainty at the adopted mass limit, rather than this additional sensitivity.

Image not extracted

$^{a}$ Photometric redshift. $^{b}$ Midpoint of the published SED range (3–10). $^{c}$ Derived here from line fluxes or $\xi_{\mathrm{ion}}$ and magnifications.

$^{d}$ The full MACS0416 system. $^{e}$ Range of acceptable age fits. $^{f}$ Compiled by Pfeffer et al. (2025). Ages depend on the assumed star formation history. $>$ denotes incomplete host coverage. Host masses retain the source IMF (Chabrier for the Sunrise arc and Misty Moons, Kroupa otherwise), with the Sunrise range spanning two fits.

References are (1) Messa et al. (2024), (2) Swinbank et al. (2007), (3) Vanzella et al. (2023), (4) Fujimoto et al. (2025), (5) Messa et al. (2025), (6) Abdurro’uf et al. (2025), (7) Mowla et al. (2024), (8) Adamo et al. (2024), (9) Bradley et al. (2025), (10) Messa et al. (2026), (11) Nakane et al. (2025). Photometry for the Cosmic Grapes, MACS0416 and the Firefly Sparkle as compiled by Pfeffer et al. (2025).

The growth model and the cloud-to-cluster mapping

Figure C1 supports the mapping discussed in § 4.3. Natal clouds are matched in every catalogued output, unlike the more widely spaced cloud census. Initial accretion rates rise with cloud mass, while the coupling $K$ remains roughly constant. The apparent high-mass deficit in panel (c) reflects binning by measured mass, which biases towards sink particles with unusually large growth. It disappears when binned by initial rate, although late gas supply broadens the mass distribution beyond the model scatter in $K$.

Image not extracted

More massive clouds also host more sink particles, from one below a cloud mass of a few $10^{4}\,M_{\odot}$ to two or three above $10^{7}\,M_{\odot}$. We predict each sink particle’s mass separately rather than dividing the cloud mass among all sink particles. Above our adopted mass fitting limits for the ICMF, the initial accretion rate slopes are $-1.81$, $-1.68$ and $-1.69$ in HiRes, FID and NoRad respectively. Combined with the measured exponents of the relation between initial accretion rate and final stellar mass, these imply ICMF slopes of $-1.63$, $-1.53$ and $-1.84$, within 0.04 of the measured power-law slopes in the radiation runs and 0.12 steeper in NoRad. These power-law fits use separately selected lower mass limits for the 3 runs and describe the mapping sample rather than the sample used to perform the Schechter fits. The stellar mass already present at the time of initial accretion rate measurement prevents a direct conversion between the fitted duration–rate and mass–rate exponents.

We fit equation (6) with $p = 1/2$ to individual HiRes and FID accretion histories, fixing the initial stellar mass to $10^{3}\,M_{\odot}$ and $\dot{M}_{0}$ to the accretion rate measured when the sink particle reaches this initial stellar mass. The single free parameter of the model, $K$, minimizes squared residuals in $\log\dot{M}$. The mass function mapping draws coupling values from these fits. The prediction bands in Figure 13 represent variation between realizations that draw $K$ from this distribution, rather than uncertainty in the growth law. The predictions in Figure C1 instead use one median $K$ per run, with bands showing the 16th–84th percentiles of the prediction ratios across sink particles. Alternative braking laws, including linear ($p = 1$) mass dependence and saturation of thermal feedback, give similar mass predictions. Allowing the braking exponent, $p$, to vary introduces covariance between $K$ and the initial accretion rate without materially improving the description of the accretion histories. We therefore retain $p = 1/2$ as a simple empirical description, without identifying a unique braking mechanism.

The cell of a forming sink particle

Figure D1 shows rapid heating of the gas cell hosting a forming cluster to $T/\mu= 10^4\,\mathrm{K}$ in HiRes and FID, followed by growth of the nominal Strömgren volume as the ionizing photon rate rises and gas density falls. In NoRad, the host cell remains cold and neutral, as expected. Solid curves describe the overall cluster population rather than individual accretion shutdown histories, and the nominal volume is a photon-budget estimate rather than a resolved ionization front.

Uniform cold-medium benchmark

For the uniform cold-medium benchmark discussed in § 4.2, we set $P_{\mathrm{conf}} = n_0 k_{\mathrm{B}}T_0$ (where $n_0$ and $T_0$ are the density and temperature of the cold ambient medium respectively), $Q = qM_\star$ and $L_{\mathrm{abs}} = f_{\mathrm{abs}}\ell M_\star$ in equations (10) and (11). Equating each radius to $r_{\mathrm{acc}}$ then gives:

Here $q$ and $\ell$ are the source output per unit mass of stars formed and $f_{\mathrm{abs}}$ is the absorbed fraction. The zero-age BPASS stellar populations for a median stellar metallicity of $Z = 0.0018$ give $q \simeq7.8 \times10^{46}\,\mathrm{s}^{-1}\,M_\odot^{-1}$ and $\ell\simeq4.9 \times10^{36}\,\mathrm{erg}\,\mathrm{s}^{-1}\,M_\odot^{-1}$. The threshold masses scale as $r_{\mathrm{acc}}^3n_0^2$ and $r_{\mathrm{acc}}^2n_0/f_{\mathrm{abs}}$ at fixed ambient temperatures. These cold-medium benchmark reference estimates omit gravity and inflow ram pressure as a source of confinement for the gas. A full absorption assumption ($f_{\mathrm{abs}} = 1$) yields the lowest direct-pressure mass threshold.