# ASTRAEUS

## X. Indications of a top-heavy initial mass function in highly star-forming galaxies from JWST observations at $z > 10$

Anne Hutter<sup>1,2,*</sup> ◆ Elie R. Cueto<sup>1</sup>, Pratika Dayal<sup>3</sup>, Stefan Gottlöber<sup>4,5</sup>, Maxime Trebitsch<sup>5</sup>, and Gustavo Yepes<sup>6,7</sup>

1. Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark

2. Cosmic Dawn Center (DAWN), Copenhagen, Denmark

3. Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

4. Leibniz-Institut für Astrophysik, An der Sternwarte 16, 14482 Potsdam, Germany

5. Sorbonne Université, Observatoire de Paris, PSL Research University, CNRS, LERMA, 75014 Paris, France

6. Departamento de Fisica Teorica, Modulo 8, Facultad de Ciencias, Universidad Autonoma de Madrid, 28049 Madrid, Spain

7. CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, 28049 Madrid, Spain

Received 2 October 2024 / Accepted 12 December 2024

### ABSTRACT

**Context.** The James Webb Space Telescope (JWST) has uncovered an abundance of $z > 10$ galaxies bright in the ultraviolet (UV), whose presence poses a challenge to traditional theoretical models at high redshifts. Various new models have recently emerged to address this discrepancy by refining their description of star formation.

**Aims.** Here, we investigate whether modifications to the stellar initial mass function (IMF) alone can reproduce the $z > 10$ UV luminosity functions (UV LFs) when the star formation rate is used as a proxy for the fraction of massive stars.

**Methods.** We incorporated an evolving IMF into the ASTRAEUS galaxy evolution and reionisation simulation framework, which becomes increasingly top-heavy as the gas density in a galaxy rises above a given threshold. Our implementation accounts for the IMF’s effects on supernova (SN) feedback and metal enrichment, as well as UV and ionising emissivities.

**Results.** For this evolving IMF model, we find that (i) the maximum UV luminosity enhancement is twice as large in massive galaxies ($\Delta M_{UV} \approx 2.6$) than those where star formation is strongly limited by SN feedback ($\Delta M_{UV} \approx 1.3$); (ii) it successfully reproduces the observed UV LFs at $z = 5-15$; (iii) galaxies with top-heavy IMFs exhibit the highest star formation rates, driven by their location in local density peaks, which facilitates higher gas accretion rates; (iv) the $1\sigma$ variances in the UV luminosity are only slightly higher compared to when assuming a Salpeter IMF, but the $2\sigma$ variances are significantly increased by a factor of $1.4-2$ boosting the abundance of UV-bright galaxies at $z > 10$; and (v) reionisation begins earlier with more extended large ionised regions and fewer smaller ones during its initial stages, although these differences diminish at lower redshifts, leading to a similar end of reionisation at $z \approx 5.6$.

**Key words.** methods: numerical – stars: luminosity function, mass function – galaxies: evolution – galaxies: high-redshift – intergalactic medium – dark ages, reionization, first stars

## 1. Introduction

The first galaxies emerging in our Universe over its first billion years ushered the last major phase transition, the epoch of reionisation (EoR), by emitting energetic photons that gradually ionised the neutral hydrogen in the intergalactic medium (IGM) (Barkana & Loeb 2001; Dayal & Ferrara 2018). Exactly how this phase transition unfolded on global and local scales at $z \approx 20-5$ remains unconstrained, given our incomplete understanding of the population and properties of the first galaxies. Detections of numerous $z \ge 8$ galaxies with the James Webb Space Telescope (JWST) have revealed young, star-forming, compact, metal-poor galaxies whose stellar populations emit harder hydrogen-ionising radiation during the EoR (e.g. Adams et al. 2023; Atek et al. 2024; Bradley et al. 2023; Bunker et al. 2023; Curtis-Lake et al. 2023; Heintz et al. 2023; Rinaldi et al. 2024; Schaerer et al. 2022; Simmonds et al. 2024a; Vanzella et al. 2023). However, these JWST observations also present a puzzle: the detection of a higher abundance of UV-bright galaxies at $z > 10$ than standard theoretical mod-

els predict (e.g. Labbé et al. 2023; Adams et al. 2023, 2024; Arrabal Haro et al. 2023; Atek et al. 2023; Austin et al. 2023; Boylan-Kolchin 2023; Donnan et al. 2024; Harikane et al. 2024; McLeod et al. 2024). This puzzle remains, even when limiting the galaxy candidates to those spectroscopically confirmed (see Harikane et al. 2024).

Various theoretical investigations have been conducted to explain this puzzle, focusing on selection biases and altered physical processes. Some studies suggest that the samples of observed galaxies are biased, due to the stochastic nature of star formation, mostly detecting galaxies undergoing periods of intense star formation, while missing those with strongly suppressed or no star formation at all (e.g. Mason et al. 2023; Mirocha & Furlanetto 2023; Shen et al. 2023; Sun et al. 2023; Kravtsov & Belokurov 2024; Gelli et al. 2024). However, assuming a dark matter (DM) mass-dependent UV luminosity variance, consistent with the FIRE simulation and scaling inversely with the halo’s escape velocity ($\propto M_h^{-1/3}$), has been shown to be insufficient in reproducing the abundance of UV-bright galaxies at $z > 10$ (Gelli et al. 2024). This could indicate that the physical conditions and processes at

<footnote>* Corresponding author; anne.hutter@nbi.ku.dk</footnote>
very high redshifts differ. Indeed, the compactness of these early galaxies implies higher densities of metal-poor gas, where stellar feedback might be less efficient in suppressing subsequent star formation, thereby boosting the star formation efficiency (SFE) (Dekel et al. 2023). Additionally, radiatively driven outflows may eject dust from star formation sites, reducing UV attenuation by dust (Ferrara et al. 2023; Fiore et al. 2023; Mauerhofer & Dayal 2023; Yung et al. 2024; Ziparo et al. 2023), or massive black holes accreting at or slightly above the Eddington rate could contribute to the observed UV luminosity (Pacucci et al. 2022). However, testing whether the SFE is enhanced in dense, metal-poor, star-forming clouds requires radiative hydrodynamical simulations, and even state-of-the-art semi-analytic galaxy evolution models fail to reproduce the abundance of $z > 10$ UV-bright galaxies when neglecting dust attenuation (see e.g. Cueto et al. 2024).

Theoretical models and simulations indicate that the first metal-free stars were massive, with stellar masses of $\ge 60$ M$_\odot$ (Abel et al. 2002; Bromm et al. 2002; Yoshida et al. 2006; Fukushima et al. 2020); however, also see Clark et al. 2011), suggesting that the stellar initial mass function (IMF) in early galaxies was skewed towards higher masses, thereby reducing the overall mass-to-light ratio. Several observations seem to support this idea, with high-redshift galaxies exhibiting higher ionising emissivities (Simmonds et al. 2024a), bluer UV slopes (Cullen et al. 2024), and unusual abundance ratios, such as enhanced [N/O] (e.g. in GN-z11, CEERS-1019, GLASS-150008, JADES-GS-z9; Bunker et al. 2023; Cameron et al. 2024; Isobe et al. 2023; Curti et al. 2024; Senchyna et al. 2024; Topping et al. 2024; Watanabe et al. 2024b), lower [C/O] (e.g. in CEERS-1019; Isobe et al. 2023), and low [O/Fe] abundances (e.g. in GN-z11; Nakane et al. 2024). Enhanced [N/O] abundances, which indicate active CNO cycling in stars (Ekström 2021), are expected in supermassive stars ($\ge 1000$ M$_\odot$, Charbonnel et al. 2023; Nagele & Umeda 2023; Isobe et al. 2023), very massive stars (100–1000 M$_\odot$) (Vink 2023; Curti et al. 2024), rotating massive stars (Nandal et al. 2024), Wolf-Rayet stars in scenarios involving dual star-bursts (Kobayashi & Ferrara 2024), or tidal disruptions. Similarly, lower [C/O] abundances are consistent with contributions from supermassive stars, Wolf-Rayet stars, or tidal disruptions (Isobe et al. 2023; Watanabe et al. 2024a). These observational findings suggest that the stellar populations in the first galaxies formed with an IMF that was more top-heavy than the local Salpeter IMF (Salpeter 1955). First estimates of the associated reduction in the stellar mass-to-UV luminosity ratio have concluded that a top-heavy IMF could resolve the UV-bright galaxy abundance problem at $z > 10$ (Haslbauer et al. 2022; Trinca et al. 2024; Yung et al. 2024).

However, these works omit the IMF dependency of stellar feedback, neglecting the enhanced fraction of stars exploding as supernovae (SNe), which leads to a reduction in subsequent star formation more immediately, while also increasing and altering the metal enrichment of the interstellar gas per stellar mass formed. In Cueto et al. (2024), we addressed this gap using the ASTRAEUS framework, a semi-numerical model that tracks the interdependent evolution of galaxies and reionisation (Hutter et al. 2021; Ucci et al. 2023; Hutter et al. 2023). We adopted the IMF inferred from star-forming cloud simulations in Chon et al. (2022), where the IMF became increasingly top-heavy with rising redshift and decreasing gas-phase metallicity, as a result of the increasing temperature of the cosmic microwave background (CMB) at higher redshifts and the associated decreasing gas cooling efficiency. To reproduce the

observed UV luminosity functions (UV LFs) at $z \lesssim 10$, we adjusted the free model parameters in ASTRAEUS to adopt an overall lower SFE than for the standard assumption of a constant Salpeter IMF. This reduction in SFE roughly corresponded to the decrease caused by enhanced radiative pressure from young stars, as reported by Menon et al. (2024), for cloud surface densities of $\sim 10^3$ M$_\odot$ pc$^{-2}$ at $Z = 0.01$ Z$_\odot$ and a UV luminosity enhancement of 10. It also compensated for the fact that ASTRAEUS did not capture the effects of radiative feedback from massive stars on star formation. Yet, the resulting UV LFs at $z = 5-15$ remained consistent with those produced by models assuming a constant Salpeter IMF. In contrast, Lu et al. (2025) carried out their study without altering the SFE and found that adopting the same top-heavy IMF for all galaxies can effectively match the observed UV LF at $z \lesssim 12$.

While Cueto et al. (2024) assumed top-heavy IMFs for all galaxies, and Lu et al. (2025) applied a universal top-heavy IMF to galaxies in the burst mode of star formation within the GALFORM model, here we consider a scenario where only a fraction of galaxies have a top-heavy IMF, with its degree of top-heaviness varying according to each galaxy's properties. Observations of local star-forming regions and nearby galaxies find the IMF to become more top-heavy in regions of intense star formation and dense gas environments (e.g. Meurer et al. 2009; Watts et al. 2018; Zhang et al. 2018; Gunawardhana et al. 2011; Fontanot et al. 2017, 2018). This trend may also possibly be due to the higher cosmic ray densities in these regions raising gas temperatures (Papadopoulos et al. 2011). For this reason and, in addition, because the radiation and elemental abundances in $z > 5$ galaxies suggest a higher prevalence of massive stars in UV-bright galaxies, here we investigate a scenario where the IMF becomes increasingly top-heavy with rising specific star formation rate (sSFR), referred to as the 'evolving IMF' in this work. Naturally, this IMF dependency also aligns with those found in the galaxy-wide IMF (gwIMF) theory (Kroupa & Weidner 2003; Weidner et al. 2013; Yan et al. 2017; Jeřábková et al. 2018) that relies on empirical relations, such as the relation between the maximum star mass and cluster mass. Assuming that the SFE is unaffected by radiative feedback from massive stars, as suggested for higher cloud densities (Menon et al. 2024), we primarily investigate whether adapting the otherwise assumed Salpeter IMF to a more top-heavy IMF in galaxies with higher sSFR values can reproduce the observed UV LFs at $z = 5-15$, without altering the SFE or SNe wind coupling efficiencies. To further understand the results and their implications, we also explore (i) the extent to which a top-heavy IMF enhances the UV luminosity of galaxies, regardless of whether star formation is suppressed by SN feedback; (ii) the properties and characteristics of galaxies with a top-heavy IMF; and (iii) how reionisation changes when transitioning from a universal Salpeter IMF to our evolving IMF.

This paper is organised as follows. In Section 2, we describe the updates made to the ASTRAEUS model compared to Cueto et al. (2024), focusing on our new parameterisation of the IMF and its scaling with galactic properties to mimic a dependency on the sSFR. Section 3 introduces a simple toy model to explore the maximum possible UV luminosity boost for various halo masses when assuming our top-heavy IMF parameterisation. We then discuss in Section 4 how the free model parameters in our evolving IMF parameterisation affect the UV LFs and we identify the parameter set that best fits the UV LFs at $z = 5-15$. In Sections 5 and 6, we analyse the characteristics of the galaxies with top-heavy IMFs in our model and examine changes in the UV luminosity-halo mass distribution. We briefly discuss the impact of the evolving IMF on reionisation
<table><caption>**Table 1.** Details of the DM-only N-body simulation VSMDPL</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>L</td><td>160 h<sup>-1</sup> cMpc</td><td>Simulation box size</td></tr><tr><td>M<sub>h,min</sub></td><td>10<sup>8.2</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>Minimum halo mass</td></tr></tbody></table>

in Section 7 and the range of the ionising emissivity values of high-redshift galaxies in Section 8. We present our conclusions in Section 9. Throughout this paper, we assume the AB magnitude system (Oke & Gunn 1983) and a ΛCDM universe with the Planck Collaboration Int. XLVI (2016) cosmological parameters: $\Omega_{\Lambda} = 0.692885$, $\Omega_{m} = 0.307115$, $\Omega_{b} = 0.048206$, $H_{0} = 100h = 67.77 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_{s} = 0.96$, and $\sigma_{8} = 0.8228$.

## 2. The model

In this section, we briefly describe the changes we have incorporated in the ASTRAEUS framework to follow an IMF that becomes increasingly top-heavy as the gas content and sSFR in galaxies rise. The ASTRAEUS framework post-processes the VSMDPL (very small multidark planck) DM-only N-body simulation (for details, see Table 1) with an enhanced version of the semi-analytic galaxy evolution model DELPHI (Dayal et al. 2014, 2022) that is coupled to the semi-numerical reionisation scheme CIFOG (Hutter 2018) in a fully self-consistent manner. ASTRAEUS incorporates all the key physical processes that shape the evolution of $z > 5$ galaxies, including gas accretion, mergers, star formation, SN feedback, metal enrichment, dust formation, and radiative feedback from reionisation. It also models reionisation, tracking the time evolution of the spatial distribution of ionised regions. A detailed description of the implementation of these processes can be found in Hutter et al. (2021), Hutter (2018), Ucci et al. (2023), Hutter et al. (2023) and Cueto et al. (2024). We advance the version described in Cueto et al. (2024) by adjusting the parameterisation of the IMF and the ionising and ultraviolet luminosities as we detail in this section. We note that to avoid artificial starbursts of galaxies coming into existence in the simulation (i.e. galaxies without progenitors), particularly since our new IMF parameterisation will depend on the gas-to-DM mass ratio, we alter the assumption on the initial gas mass in a galaxy upon its existence (i.e. when it has no progenitors). Instead of assuming $M_{\text{g}}^{i} = f_{\text{g}} \frac{\Omega_{\text{b}}}{\Omega_{\text{m}}} M_{\text{h}}$, we adopt now $M_{\text{g}}^{i} = 0 M_{\odot}$. Since galaxies without progenitors reside in DM halos with $M_{\text{h}} \lesssim 10^{8.3} M_{\odot}$, their SN feedback limited star formation is highly stochastic (Legrand et al. 2022) oscillating between phases of forming stars and forming no stars; we assume these galaxies to start in a phase of no star formation<fnref n="1" />.

### 2.1. The initial mass function

We adopt an IMF with a Salpeter slope of $\gamma = -2.35$ for stars with masses below a mass, $M_c$, and a log-flat slope of $-1$ for stars with masses above $M_c$, consistent with the findings in Chon et al.

(2022):

$$ \frac{dN}{dM} = \begin{cases} (1 - f_{\text{massive}}) \frac{\gamma+2}{M_c^{\gamma+2} - M_i^{\gamma+2}} M^{\gamma}, & M_i \le M < M_c, \\ f_{\text{massive}} \frac{1}{M_f - M_c} M^{-1}, & M_c \le M \le M_f. \end{cases} \quad (1) $$

Here, $M_c$ represents the stellar mass at which the IMF slope changes and its value depends on the fraction of mass formed in massive stars, $f_{\text{massive}}$. We numerically derive the relation between $f_{\text{massive}}$ and $M_c$ by solving the equation $\int_{M_i}^{M_f} \frac{dN}{dM} M dM = 1 M_{\odot}$. However, in contrast to Cueto et al. (2024), we scale $f_{\text{massive}}$ with the gas content of a galaxy above a redshift-dependent gas density threshold. This scaling and threshold are chosen such that $f_{\text{massive}}$ is effectively proportional to the sSFR of a galaxy above a specified threshold, sSFR$_{thresh}$.

$$ f_{\text{massive}} = k \times \frac{\text{sSFR} - \text{sSFR}_{\text{thresh}}}{\text{sSFR}_{\text{max}} - \text{sSFR}_{\text{thresh}}} \quad (2) $$

In our model, the SFE depends on the galaxy's gravitational potential. While in more massive galaxies, the SFE remains a constant value scaling with the gas free-fall time in halos, $\tau = (4\pi G\rho)^{-1/2} \propto H(z)^{-1} \propto (1+z)^{-3/2}$, at the respective redshift, $z$, it decreases to an efficiency in lower mass galaxies, where galaxies form only as many stars as are necessary to eject all their gas due to SN explosions,

$$ f_{\star}^{\text{ej}}(z) = \frac{v_c^2}{v_c^2 + f_w^{\text{eff}}(z) E_{51} v_z} \left[ 1 - \frac{f_w^{\text{eff}}(z) E_{51} \sum_j v_j M_{\star,j}^{\text{new}}(z_j)}{M_{\star,j}^i(z) v_c^2} \right]. \quad (3) $$

Here, $v_c$ represents the halo's rotational velocity, $E_{51}$ denotes the energy released by a Type II supernova (SNII), and $v_z$ is the IMF-dependent fraction of stellar mass that forms and explodes within the current time step. For the stellar mass formed in the previous time step $j$, $M_{\star,j}^{\text{new}}$, $v_j$ represents the IMF-dependent fraction that explodes in the current time step.

$$ f_w^{\text{eff}}(z) = \frac{f_w}{1 + \left(\frac{\Delta t}{20 \text{ Myr}} - 1\right) \frac{M_{\star,j}^{\text{mer}}(z)}{M_{\star,j}^{\text{g}}(z)}}, \quad (4) $$

describes the fraction of SN energy coupling to the gas and driving winds. Also, $M_{\star,j}^{\text{new}}$ is the stellar mass formed in the current time step with length $\Delta t$, while $M_{\star,j}^{i}$ denotes the gas mass at the beginning of a time step and $M_{\star,j}^{\text{mer}}$ the fraction of gas inherited by the galaxy's progenitors. The SFE is thus given by

$$ f_{\star}^{\text{eff}}(z) = \min \left[ f_{\star}^{\text{ej}}(z), f_{\star} \frac{\tau(z=9)}{\tau(z)} \frac{\Delta t}{20 \text{ Myr}} \right]. \quad (5) $$

We note that although the choice of the normalisation redshift for $f_{\star}$ is arbitrary, we chose $z=9$, as it corresponds to the highest redshifts, where standard IMFs reproduce the observed UV LF. From the above, it is evident that the SFE and thus the sSFR depend on the assumed IMF. For this reason, we focus on the factors that primarily influence the sSFRs in galaxies with efficient star formation, namely, the available gas mass for star formation, $M_{\text{g}}^{i}$, and the redshift, $z^2$,

$$ \text{sSFR} = \frac{\text{SFR}}{M_{\star}} = \frac{f_{\star}^{\text{eff}} M_{\text{g}}^{i}}{\Delta t M_{\star}} \propto \frac{M_{\text{g}}^{i}}{M_{\text{h}}} (1+z)^{3/2}. \quad (6) $$

<footnote n="1">¹ The exact value adopted for $M_{\text{g}}^{i}$ has no significant effect on galaxies with converged star formation histories (SFHs, $M_{\text{h}} \gtrsim 10^{8.6} M_{\odot}$; see Appendix B in Hutter et al. 2021).</footnote>

<footnote n="2">² We note that this approximation becomes inaccurate as star formation is increasingly SN feedback-limited.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. Dependence of $f_{\text{massive}}$ on redshift $z$ and the gas-to-DM mass ratio for the evolving IMF model. The cosmological baryon ratio $\Omega_b/\Omega_m$ corresponds to ~-0.8.</figcaption></figure>

Here, SFR is the star formation rate and $M_\star$ and $M_h$ the stellar and halo masses of the galaxy at $z$. Based on this proportionality, we express the fraction of mass that is bound in massive stars as

$$f_{\text{massive}} = \min \left[ \frac{kx}{(1+z_{\text{max}})^{3/2} - (1+z_{\text{thresh}})^{3/2}} \times H(x), 1 \right] \\ x = \frac{M_g^i \Omega_m}{M_h \Omega_b} (1+z)^{3/2} - (1+z_{\text{thresh}})^{3/2}, \quad (7)$$

with $k$, $z_{\text{thresh}}$ and $z_{\text{max}}$ being free model parameters and $H(x)$ the Heaviside function. Equation (7) shows that $f_{\text{massive}}$ increases towards higher redshifts and higher gas-to-DM mass ratios, resulting in an increasingly top-heavy IMF. However, if $k$ is large, $f_{\text{massive}}$ stalls at 1 for sufficiently high gas-to-DM ratios, corresponding to a log-flat IMF. For $z_{\text{max}} = 20$, $z_{\text{thresh}} = 6$, and $k = 7$, the threshold specific star formation rate, sSFRthresh, is approximately 25 Gyr$^{-1}$, and $f_{\text{massive}}$ approaches 1 as the sSFR increases to around 60 Gyr$^{-1}$. In Fig. 1, we show how the corresponding $f_{\text{massive}}$ values depends on redshift and the gas-to-DM mass ratio of each galaxy. We note that in the remainder of this paper, we refer to an IMF that is top-heavier than the Salpeter IMF as top-heavy.

## 2.2. The ionising and ultraviolet luminosities

As described in Cueto et al. (2024), we used the stellar population synthesis code STARBURST99 (Leitherer et al. 1999) to generate the spectra for the various IMFs with $f_{\text{massive}}$ values ([0, 1]), metallicities ($Z = [0.001, 0.008]$) and ages ($t = [10^5 \text{ yr}, 10^9 \text{ yr}]$). From the simulated spectra, we derive and fit the time evolution of the ionising (number of photons with energies >13.6 eV per second) and UV (at 1500 Å) luminosities for each ($f_{\text{massive}}, Z$) combination. The analytical fitting functions for $N_{\text{ion}}$ [$s^{-1}$] and $L_{\text{UV}}$ [erg $s^{-1}$ Å$^{-1}$] are listed in Appendix A. We note that these fits differ from those presented in Cueto et al. (2024), as we use $f_{\text{massive}}$ and $Z$ as scaling parameters instead of $z$ and $Z$, and we adopt a more complex fitting function for the ultraviolet luminosity (see Appendix A).

By convolving these fitting functions for the time evolution of the ionising and UV luminosities with each galaxy's star

formation history (assuming continuous star formation within each timestep), we derive its intrinsic ionising emissivity, $\dot{Q}$, and intrinsic UV luminosity, $L_{\text{UV}}^{\text{int}}$ (for details see Hutter et al. 2021, 2023). We then modelled the dust as graphite/carbonaceous grains and derive the UV dust attenuation using the galaxy's dust mass, $M_d$, and spatial extension of the gas, $r_g$. The optical depth is calculated as

$$\tau_{\text{UV,c}} = \frac{3\sigma_d}{4as}, \quad (8)$$

where $\sigma_d = \frac{M_d}{\pi r_d^2}$ is the dust surface density, $r_d$ is the dust distribution radius, and $a = 0.03 \ \mu\text{m}$ and $s = 2.25 \ \text{g cm}^{-3}$ are the radius and density of the graphite/carbonaceous dust grains, respectively. We assume the gas and dust to be perfectly mixed, approximating the gas radius as $r_g = 4.5 \lambda r_{\text{vir}} [(1+z)/6]^{1.8}$, where $\lambda$ and $r_{\text{vir}}$ represent the spin parameter and virial radius of the simulated DM halo. The final factor accounts for the redshift evolution of the compactness of galaxies, with gas occupying a larger fraction of halos at higher redshifts$^3$. This assumption is motivated by ALMA observations showing that [CII] sizes for galaxies with fixed UV luminosity remain constant between $z \approx 7$ and 4 (Fujimoto et al. 2020; Fudamoto et al. 2022). Using a slab-like geometry, we derive the UV continuum escape fraction as

$$f_{\text{esc}}^{\text{c}} = \frac{1 - e^{-\tau_{\text{UV,c}}}}{\tau_{\text{UV,c}}}, \quad (9)$$

and calculate the observed UV luminosity of each galaxy as

$$L_{\text{UV}} = f_{\text{esc}}^{\text{c}} L_{\text{UV}}^{\text{int}}. \quad (10)$$

## 2.3. The different IMF scenarios

In the following, we consider two galaxy evolution models that only differ in their assumed initial mass function (IMF): a Salpeter IMF, with a slope of $\gamma = -2.35$ covering stellar masses from 0.1 M$_\odot$ to 100 M$_\odot$, and the evolving IMF described in Section 2.1, which also covers stellar masses from 0.1 M$_\odot$ to 100 M$_\odot$. For both models, we assume the same star formation efficiencies, $f_\star$, and SN wind coupling efficiencies, $f_w$ (see Table 2 for details).

# 3. The UV luminosity boost for top-heavy IMFs

Before analysing how the UV LFs and population of galaxies change for our evolving IMF model, we considered a toy model to estimate the extent to which a more top-heavy IMF or increasing $f_{\text{massive}}$ values can boost the UV luminosity of galaxies with a given halo mass. For this purpose, we consider the SFE (fraction of gas mass forming stars over a given time interval, $f_\star^{\text{eff}}$) when assuming instantaneous SN feedback in the ASTRAEUS framework (see Hutter et al. 2021). Eq. (3) reduces then to

$$f_{\star}^{\text{ej}} = \frac{v_c^2}{v_c^2 + f_w E_{\text{SN}}}, \quad (11)$$

with the rotational velocity of the halo being $v_c = \frac{GM_h}{r_h}(1+z) \propto M_h^{2/3}(1+z)$. We also define $f'_\star = f_\star \frac{\tau(z=9)}{\tau(z)} \frac{\Delta t}{20 \text{ Myr}}$. If we assume

<footnote>$^3$ When this factor is omitted, $f_\star$ needs to be increased to 0.04 to match the UV luminosity functions; while the UV LFs at $z \gtrsim 10$ for the Salpeter IMF underpredict then the observations more significantly, the UV LFs for the evolving IMF are hardly affected (see Sec. 4.2 for explanation).</footnote>
<table><caption>**Table 2.** Best-fit ASTRAEUS model parameters for Salpeter and evolving IMFs.</caption><thead><tr><th>Parameter</th><th>Salpeter</th><th>Evolving</th><th>Description</th></tr></thead><tbody><tr><td>f<sub>*</sub></td><td>0.03</td><td>0.03</td><td>SFE normalisation</td></tr><tr><td>f<sub>w</sub></td><td>0.2</td><td>0.2</td><td>wind coupling efficiency</td></tr><tr><td>IMF</td><td>Salpeter</td><td>Evolving</td><td>IMF</td></tr><tr><td>z<sub>max</sub></td><td>—</td><td>20</td><td>IMF parameter</td></tr><tr><td>z<sub>thresh</sub></td><td>—</td><td>6</td><td>IMF parameter</td></tr><tr><td>k</td><td>—</td><td>7</td><td>IMF parameter</td></tr></tbody></table>

that the SN energy produced by stars forming with a top-heavy IMF ($E_{SN}$) is related to the SN energy produced by stars forming with a Salpeter IMF ($E_{SN,Salpeter}$), such that

$$E_{SN} = \eta E_{SN,Salpeter}, \qquad (12)$$

we yield the relation between the SFEs of a top-heavy and a Salpeter IMF.

$$ \begin{aligned} \frac{f_{\star}^{\text{eff}}}{f_{\star, \text{Salpeter}}^{\text{eff}}} &= \frac{\min\left(\frac{v_c^2}{v_c^2 + f_w E_{\text{SN}}}, f'_{\star}\right)}{\min\left(\frac{v_c^2}{v_c^2 + f_w E_{\text{SN,Salpeter}}}, f'_{\star}\right)} \quad (13) \\ &\approx \begin{cases} 1/\eta & \text{for } \frac{f_w E_{\text{SN,Salpeter}}}{v_c^2} \gtrsim 10, \\ [1/\eta, 1] & \text{for } 10^{-3} \gtrsim \frac{f_w E_{\text{SN,Salpeter}}}{v_c^2} \le 10, \\ 1 & \text{for } \frac{f_w E_{\text{SN,Salpeter}}}{v_c^2} \le 10^{-3}. \end{cases} \end{aligned} $$

In our top-heavy IMF parameterisation, $\eta$ increases from ~1 to ~2.6 (3.4) as $f_{\text{massive}}$ rises, assuming a maximum stellar mass for SNe explosions of $M_{\text{SN}}^{\text{max}} = 50$ (100) M$_{\odot}$. The relations between the SFEs for stars forming according to a top-heavy IMF and those forming according to a Salpeter IMF differ for low-mass (SN feedback-limited) and massive galaxies. To understand these different scalings and their impact on the UV enhancement for top-heavy IMFs, we show in Fig. 2 the halo mass dependence of the SFRs$^4$ (top panel, blue lines) and UV magnitudes (bottom panel, blue lines) of gas-rich galaxies for the Salpeter IMF (black) and increasingly top-heavy IMFs (dark to bright). The green lines in the top panel depict the factor by which stars forming during the past 30 Myr according to a top-heavy IMF are brighter than those forming according to a Salpeter IMF, accounting for the SFE relation in Eq. (13). The assumed values are $f_w = 0.2$ and $f'_{\star} = 0.027$, based on $\Delta t = 11$ Myr and $z = 13$.

In lower mass galaxies ($M_h \lesssim 10^{9.5}$ M$_{\odot}$), where star formation is limited by SN feedback, the SFE decreases for a top-heavy IMF by the same amount that its SN energy increases compared to a Salpeter IMF. As the IMF becomes more top-heavy, the SFE declines because stronger and more instantaneous SN feedback, driven by the presence of more massive stars with shorter lifetimes, suppresses subsequent star formation more quickly (see also Cueto et al. 2024). This can be seen in the top panel of Fig. 2 where the SFR below $M_h \lesssim 10^9$ M$_{\odot}$ decreases with rising $f_{\text{massive}}$ values. At the same time, the resulting UV luminosity per stellar mass that is formed nevertheless increases as $f_{\text{massive}}$ rises (cf. green lines). Yet, most of this UV luminosity enhancement occurs for values of $f_{\text{massive}} \lesssim 0.5$, up to a factor of ~3, corresponding to a shift in $\Delta M_{\text{UV}} \approx -1.2$.

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Halo mass dependence of the top-heavy IMF increase in UV luminosity compared to the Salpeter IMF for a 30 Myr long starburst (top panel, green lines), the SFR of gas-rich galaxies with top-heavy IMFs (top panel, blue lines), and their maximum UV luminosity (bottom panel, blue lines) for $f_{\text{massive}} = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1$ (from dark to light) at $z = 13$. Black lines show the SFR and maximum intrinsic (dust unattenuated) UV luminosity of gas-rich galaxies with a Salpeter IMF. Assumed values are $f'_{\star} = 0.027$ and $f_w = 0.2$.</figcaption></figure>

Any further increase in $f_{\text{massive}}$ enhances the UV luminosity by less than 10%. The reason for this decrease in UV luminosity enhancement is as follows: When $f_{\text{massive}} > 0.525$, the stellar mass at which the IMF slope changes, $M_c$, drops below the minimum stellar mass exploding as SNe, $M_{\text{SN,min}}$; the mass distribution of the stars exploding as SNe is described by a log-flat IMF. Increasing $f_{\text{massive}}$ above 0.525 solely raises the normalisation of the IMF for $M > M_{\text{SN,min}}$, resulting in a more significant increase in the number of SNe. However, the rise in UV luminosity with increasing $f_{\text{massive}}$ remains relatively consistent at both lower and higher $f_{\text{massive}}$ values, as contributions come from stars of all masses. As a result, SN feedback-limited galaxies can only experience a UV luminosity boost up to $\Delta M_{\text{UV}} = -1.3$, even with a fully log-flat IMF. However, due to the stronger SN feedback, a more top-heavy IMF results in more massive galaxies, with their star formation limited by SN feedback. As the break in the maximum intrinsic (i.e. dust unattenuated) UV luminosity, $M_{\text{UV,max}}$, (seen in the bottom panel of Fig. 2) indicates, at $z = 13$, galaxies forming stars according to the Salpeter IMF and expelling all their star-forming gas have halo masses of $M_h \lesssim 10^{8.8}$ M$_{\odot}$.

<footnote>$^4$ SFR = $f_{\star}^{\text{eff}} \frac{\Omega_h}{\Omega_m} M_h$.</footnote>
However, this mass limit increases by ~0.8 dex to $M_h \lesssim 10^{9.6} M_\odot$ for a complete log-flat IMF.

In more massive galaxies ($M_h \gtrsim 10^{9.5} M_\odot$), SN feedback is not strong enough to completely shut down star formation, leaving the SFE unaffected when moving from a Salpeter to a more top-heavy IMF (see the SFR above $M_h \approx 10^{9.5} M_\odot$ in the top panel of Fig. 2). Given there is no decrease in the SFE as the IMF becomes more top-heavy (higher $f_{\text{massive}}$ values), the UV luminosity enhancement due to the presence of more massive stars not only exceeds that in SN feedback-limited galaxies but also does not stall around $f_{\text{massive}} \approx 0.5$; instead, it accelerates. This accelerating rise in UV luminosity enhancement is because replacing multiple lower-mass stars with massive stars increases the net UV luminosity due to the mass-luminosity relation of main-sequence stars where $L \propto M^{3.5}$. A log-flat IMF can boost the UV luminosity of more massive galaxies by a maximum of $\Delta M_{\text{UV}} \approx -2.6$ compared to a Salpeter IMF, which is twice as much as in an SN feedback-limited regime.

# 4. The UV luminosity functions at $z = 5-15$

## 4.1. Model parameter dependence

We found the free model parameters for the Salpeter and evolving IMF models as follows. First, we estimated the SFE ($f_\star$) and SN wind coupling efficiency ($f_w$) that best reproduce the observed UV LFs at $z = 5-10$ for the Salpeter IMF model (see Table 2 for values). Next, we adopted the same $f_\star$ and $f_w$ values for our evolving IMF model, assuming $z_{\text{max}} = 20$, and we found the $z_{\text{thresh}}$ and $k$ values that best reproduce the observed UV LF at $z = 5-15$. The free model parameters $z_{\text{thresh}}$ and $k$ in our evolving IMF model determine how $f_{\text{massive}}$ (which describes the top-heaviness of the IMF) depends on the galaxies' gas content and sSFR values.

Increasing $z_{\text{thresh}}$ corresponds to decreasing the sSFR threshold (sSFR$_{\text{thresh}}$) above which a galaxy's IMF becomes top-heavy. As a result, fewer galaxies will have a top-heavy IMF because their sSFR values increase with rising redshift due to their decreasing dynamic times, thus increasing SFEs. In the context of the UV LF, this reduction in the fraction of galaxies with top-heavy IMFs shifts the redshift at which the evolving IMF model's UV LF starts exceeding that of the Salpeter IMF model to higher redshifts. For instance, for $z_{\text{thresh}} = 5$, we find the evolving IMF model's UV LF to surpass the Salpeter IMF's UV LF at $z \gtrsim 6$ and for $z_{\text{thresh}} = 6$ at $z \gtrsim 7$ for large $k$ values ($> 5$). However, the exact redshift at which the Salpeter IMF model's UV LF is surpassed depends on the galaxies' actual gas-to-DM mass ratio and the value assumed for $k$, which regulates how top-heavy a galaxy's IMF becomes given its sSFR value or gas-to-DM mass ratio: the larger the value of $k$, the more top-heavy the galaxy's IMF for the same gas-to-DM ratio. In summary, while $z_{\text{thresh}}$ determines how many galaxies exhibit a top-heavy IMF, $k$ regulates how top-heavy their IMFs are.

Raising $k$ leads (on average) to higher $f_{\text{massive}}$ values, more top-heavier IMFs, and a boost of the UV luminosities, shifting the UV LF to higher luminosities. However, as discussed in Section 3, the UV luminosity enhancement stalls for SN feedback-limited galaxies around $f_{\text{massive}} \approx 0.5$ and hardly increases for values beyond. Thus, below a certain UV luminosity threshold, $M_{\text{UV,thresh}}$, increasing $k$ beyond a value that corresponds to the majority of top-heavy galaxies having $f_{\text{massive}} \gtrsim 0.5$ values results in no significant further shift of the UV LF to higher luminosities. Beyond this UV luminosity threshold ($M_{\text{UV}} < M_{\text{UV,thresh}}$), the UV luminosities of top-heavy galaxies

can be further increased by raising $f_{\text{massive}}$ beyond 0.5. However, it is noteworthy that as $f_{\text{massive}}$ increases, the halo mass and UV luminosity thresholds for galaxies to become SN feedback-limited also move to higher values, which is why we do not see abrupt changes. Given the different UV luminosity enhancement of SN feedback-limited and more massive galaxies, raising $k$ above a certain value leads to the bright end of the UV LF being shifted towards higher luminosities than the faint end. Yet, this shift is limited by the maximum shift possible in the UV, $\Delta M_{\text{UV}}$, and the fraction of galaxies with top-heavy IMFs.

The maximum mass above which stars do not explode as energetic SNe but collapse directly into black holes also affects how much the UV luminosities of galaxies with top-heavy IMFs are boosted compared to a Salpeter IMF. Decreasing $M_{\text{SN}}^{\text{max}}$ reduces the SN feedback strength (or number of SNe explosions) fractionally more for a top-heavy than for a Salpeter IMF, since more massive stars directly collapse into black holes. This decrease in SN feedback causes the galaxies with top-heavy IMFs to exhibit higher UV luminosities, shifting the UV LFs to higher luminosities. We find that for $M_{\text{SN}}^{\text{max}} = 100 M_\odot$, the SN feedback is essentially too strong to boost our evolving IMF model's UV LFs to the observational data points. It is only by decreasing $M_{\text{SN}}^{\text{max}}$ to the more realistic value of $\lesssim 50 M_\odot$ that we can achieve a reasonable fit. Given the trends described above, for our evolving IMF, we chose $M_{\text{SN}}^{\text{max}} = 50 M_\odot$, which is supported by the SNe failure rate of $> 30 M_\odot$ stars (e.g. Kobayashi et al. 2020), and find $z_{\text{thresh}} = 6$ and $k = 7$ to best-fit the observed UV LF at $z = 5-15$. These best-fit values have been determined by comparing the simulated dust-attenuated UV LFs to the observational data. We note that due to our coarse griding of $k$ and $z_{\text{thresh}}$ (with $\Delta k = 1$ and $\Delta z_{\text{thresh}} = 1$), the best-fit values have been identified visually. Additionally, we find that both IMF models' stellar mass functions (SMFs) are consistent with the values inferred from observations (Song et al. 2016; Duncan et al. 2014; González et al. 2011).

## 4.2. Best-fit case

Next, we discuss the UV LFs for our best-fit parameter set and the physical processes that drive their shape and evolution at $z > 10$. In Fig. 3, we show the intrinsic (light) and observed (dark) UV luminosity functions (LFs) for the Salpeter IMF (red lines) and evolving IMF (blue lines) models at $z = 6, 8$, and $10-15$.

The UV LFs of both IMF models reflect the physical processes driving their redshift evolution. As redshift decreases, their shift to higher UV luminosities and number densities reflects hierarchical structure formation, with galaxies growing in mass through mergers and accretion and increasingly more low-mass galaxies forming, raising the luminosity of the brightest galaxies (from $M_{\text{UV}} \approx -20$ and $-22$ for the Salpeter and evolving IMF models at $z = 15$ to $-24$ at $z = 6$) and the number of UV-faint galaxies (from $\sim 10^{-2.5} \text{cMpc}^{-3}$ at $z = 15$ to $\sim 10^{-3} \text{cMpc}^{-1.4}$ at $z = 6$ for $M_{\text{UV}} \approx -16$), respectively. Despite their rising numbers, lower mass galaxies not only show bursty star formation (e.g. Legrand et al. 2022) due to SN feedback and radiative feedback from reionisation (leading to a broadening of their UV luminosity values) but are also consumed by mergers in the vicinity of more massive galaxies. Both govern the UV LF's number density and luminosity evolution, causing the flattening of its slope at the faint end and developing a break around $M_{\text{UV}} \approx -18$ towards lower redshifts.

As can be seen from Fig. 3, both IMF models (Salpeter and evolving IMF) are able to reproduce the observed UV LFs at $z \lesssim 10$, capturing both the normalisation and the knee that
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** UV luminosity functions for the Salpeter IMF (red lines) and evolving IMF (blue lines) at $z = 6–15$ including (dark coloured lines) and without (light coloured lines) dust attenuation. Light grey points depict results from observations with the Hubble Space Telescope (Atek et al. 2015, 2018; Bouwens et al. 2021; Castellano et al. 2010; Ishigaki et al. 2018; Livermore et al. 2017; McLeod et al. 2015, 2016; McLure et al. 2013; Oesch et al. 2018; Schenker et al. 2013), while darker grey points show the results from JWST observations (Adams et al. 2024; Bouwens et al. 2023a,b; Donnan et al. 2023a,b, 2024; Finkelstein et al. 2023; McLeod et al. 2024; Harikane et al. 2024; Pérez-González et al. 2023).</figcaption></figure>

develops around $M_{UV} \simeq -18$ at lower redshifts. However, at higher redshifts ($z > 10$), the Salpeter IMF model underpredicts the observed UV LFs, particularly at the bright end ($M_{UV} \le -19$), while the evolving IMF reproduces the observed UV LF at $z = 10–15$.

The reasons for this better fit are as follows. Firstly, towards higher redshifts, a rising fraction of galaxies exhibit top-heavy IMFs, leading to an increasing shift towards higher UV luminosities compared to the Salpeter IMF. For example, at $z = 11$, where only ~0.8% of the simulated galaxies ($M_h \ge 10^{8.3} M_\odot$) have a top-heavy IMF, the UV LF shifts only at the bright end by less than 1 dex towards higher luminosities. However, at $z = 15$, the fraction of galaxies with top-heavy IMFs increases to ~8%, resulting in a shift of the UV LF by about 1 dex. Secondly, the range of IMFs of the galaxies with top-heavy IMFs reaches top-heavy IMFs (i.e. larger $f_{massive}$ values; cf. Eq. (7)), as redshift increases, enabling larger shifts towards brighter UV luminosities even for SN-feedback limited galaxies. At $z = 6$, the UV LF matches that of the Salpeter IMF model, since we assume no galaxies have a top-heavy IMF at this redshift ($f_{massive} = 0$).

A key difference between the UV LFs of the evolving and Salpeter IMF models at $z > 10$ are the shallower slopes observed in the evolving IMF model. This difference arises for two reasons. Firstly, the average gas fraction – and thus the fraction of galaxies with a top-heavy IMF – increases with halo mass (see e.g. Hutter et al. 2021), leading to fewer low-mass and faint galaxies having a top-heavy IMF than at the massive and bright end. Secondly (as shown in Fig. 2 and discussed in Section 3), SN feedback-limited and UV-faint galaxies experience a smaller UV luminosity enhancement for a top-heavy IMF with the same $f_{massive}$ value than more massive and UV brighter galaxies. Both the higher fraction of UV-bright galaxies with a top-heavy IMF

and the more significant UV luminosity enhancement in these galaxies cause the bright end of the UV LF to shift more towards higher UV luminosities than the faint end.

It is also important to note that UV-bright galaxies ($M_{UV} \le -18$) in the evolving IMF model fall into two categories: those that are bright due to their accumulated stellar mass and those that are bright because they form stars according to a top-heavy IMF. The former have stellar and dust masses that are comparable to similarly massive galaxies in the Salpeter IMF model<fnref n="5" />, corresponding to $M_\star \simeq 10^{7.3–8.3} M_\odot$ and $M_d \simeq 10^{4–5} M_\odot$ at $z = 13$, respectively. In contrast, the latter are less massive galaxies ($M_\star \simeq 10^{6.1–7.3} M_\odot$ at $z = 13$) whose UV luminosities are boosted by the higher fraction of massive stars. These galaxies have lower dust masses ($M_d \simeq 10^{2–3.5} M_\odot$ at $z = 13$), which reduces their UV dust attenuation and results in a smaller difference between the intrinsic (bright blue lines) and observed (dark blue lines) UV LFs in Fig. 3.

Finally, we briefly discuss why our results differ from those in Cueto et al. (2024). In Cueto et al. (2024), the $f_{massive}$ values exceed ~0.2 and increase towards higher redshifts and lower gas-phase metallicities. Due to the dependence on gas-phase metallicity, low-mass, SN feedback-limited galaxies generally have larger $f_{massive}$ values than more massive galaxies. When all galaxies form stars according to a top-heavy IMF, the decreasing $f_{massive}$ with increasing galaxy mass and gas-phase metallicity results in similar UV luminosity enhancements for both SN feedback-limited and more massive galaxies (see the top panel of

<footnote n="5">⁵ However, as the fraction of galaxies with top-heavy IMFs increases towards higher redshifts, the stellar-to-halo mass and dust-to-halo mass ratios decrease slightly in massive galaxies, while the dust-to-stellar mass ratio remains unchanged.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** SFR-stellar mass relation for the Salpeter IMF (red) and the evolving IMF (blue) at z = 13. Hashed blue contours show the galaxies with an IMF top-heavier than the Salpeter IMF in the evolving IMF model. Black points show the inferred values from JADES-GS-z14-1, JADES-GS-z13-0, JADES-GS-z12-0, UNCOVER-z13, UNCOVER-z12, and GLASS-z12 (Curtis-Lake et al. 2023; Wang et al. 2023; Bakx et al. 2023). Grey points show the two set of values inferred for JADES-z13-1-LA, where assuming no nebular emission leads to a higher SFR and lower stellar mass (Witstok et al. 2024).</figcaption></figure>

Fig. 2). Consequently, assuming the same $f_\star$ and $f_\text{w}$ values as for the Salpeter IMF would shift the UV LFs to higher luminosities by approximately ~0.5 dex for massive galaxies and ~1 dex for SN feedback-limited galaxies. However, the approach adopted by Cueto et al. (2024) compensates for this shift by decreasing the SFE normalisation by factors of ~0.4 and ~0.7 for massive and SN feedback-limited galaxies, respectively (corresponding to $\Delta M_{\text{UV}} \approx 0.5$ and 1), to match the observed UV LFs at $z = 5-13$.

**5. The galaxy population with top-heavy IMFs**

To complete the picture, we analysed which galaxies within the evolving IMF model exhibit a top-heavy IMF. Since the IMF’s top-heaviness (defined as the fraction of mass produced in the log-flat part of the IMF) scales with the gas-to-DM mass ratio according to Eq. (7), we find that galaxies with top-heavy IMFs are typically gas-rich and highly star-forming, exceeding a threshold of sSFR $\approx 25 \text{ Gyr}^{-1}$. This is illustrated in Fig. 4, which shows the SFR-stellar mass relations for the evolving IMF (blue contours) and Salpeter IMF (red contours) at $z = 13$, along with the SFR-$M_\star$ combinations for galaxies with top-heavy IMFs ($f_{\text{massive}} > 0$, blue hashed contours). Indeed, the median SFRs of the galaxies with top-heavy IMFs (dashed blue lines) exceed the median SFR of the entire galaxy population (solid blue lines) by ~0.5 dex. The sSFR threshold of $25 \text{ Gyr}^{-1}$ is marked by the grey dotted line.

However, from the blue hashed contours in Fig. 4, we note two deviations from the aforementioned selection criterion for galaxies with top-heavy IMFs, both arising from the limited validity of the assumption that the sSFR scales with the gas-to-DM mass ratio, as assumed in Eq. (7). Firstly, in galaxies with lower stellar masses ($M_\star \lesssim 10^{6.5} M_\odot$), increasingly

lower sSFR values also exhibit top-heavy IMFs; namely, the blue hashed contours extend below the sSFR threshold (grey-dotted line). In these lower-mass galaxies, star formation is increasingly suppressed by SN feedback and, thus, not primarily shaped by the gas accretion rate, resulting in a decrease in the lower sSFR limit$^6$. Secondly, low-mass galaxies with $M_\star \lesssim 10^6 M_\odot$ exhibit no top-heavy IMFs, despite some of them exceeding a sSFR of $25 \text{ Gyr}^{-1}$. Below this mass threshold (corresponding to $M_h \approx 10^{8.7} M_\odot$), the galaxies’ SFHs have not entirely converged: these galaxies have existed for only a few time steps and, thus, they could not build up the stellar mass they would have attained had the merger trees been extended to lower halo masses than $M_h = 10^{8.2} M_\odot$. Both deviations highlight that a better approximation for the sSFR value of SN feedback-limited galaxies would be desirable when determining whether a galaxy exhibits a top-heavy IMF. However, this is challenging (if not impossible) since the SFR depends on the assumed IMF.

Galaxies with high sSFRs are typically gas-rich, sustained by high gas accretion rates that replenish their gas reservoirs and counteract the losses due to SNe. These galaxies are usually found in environments where they are the deepest potential wells or most massive galaxies (see e.g. Legrand et al. 2023). Consequently, low-mass galaxies with top-heavy IMFs are located in slightly above-average dense regions, while their massive counterparts are found in overdense regions, resulting in their median environmental density being slightly above that of the overall galaxy population. Furthermore, due to their higher gas accretion rates continuously diluting the ongoing metal enrichment, galaxies with top-heavy IMFs have the lowest gas-phase metallicities (oxygen-to-hydrogen mass ratio) – despite their top-heavy IMFs producing more oxygen per stellar mass through their SNe explosions.

How does the population of galaxies with top-heavy IMFs evolve? As the matter density in the Universe decreases with cosmic time, the time for matter to collapse and cool in halos of a given mass lengthens, reducing the available cool gas to form stars. Consequently, the sSFRs of galaxies with a given mass decrease over time and the SFR-$M_\star$ relation shifts to lower SFR values. Since our selection of galaxies with top-heavy IMFs effectively corresponds to a cut-off at an sSFR of $25 \text{ Gyr}^{-1}$, fewer galaxies exhibit a top-heavy IMF at lower redshifts; those that do increasingly approach the lower threshold of sSFR$_\text{thresh} \approx 25 \text{ Gyr}^{-1}$, similar to the $f_{\text{massive}}$ values described in Section 4.1. While our simulated galaxies show SFR values up to $130 \text{ Gyr}^{-1}$ at $z = 15$, their SFR only reaches $80 \text{ Gyr}^{-1}$ at $z = 10$. Interestingly, the SFR threshold inferred and required to reproduce the UV LFs at $z = 5-15$ matches the SFR value identified by Ferrara et al. (2023), above which the radiation pressure on dust is sufficient to expel the dust from the star formation sites. Although our predictions for the UV LFs at $z > 10$ also assume negligible dust attenuation, only the evolving IMF model can successfully reproduce the observed number densities.

Finally, we note that assuming more top-heavy IMFs does not significantly alter the simulated SFR-$M_\star$ relation. As explained in Cueto et al. (2024), this insensitivity of the star formation sequence arises from the self-similar mass growth of halos and, consequently, the gas and stellar mass. Making the IMF more top-heavy reduces the SFR and stellar mass, shifting

<footnote>$^6$ In our delayed SN feedback scheme (see Hutter et al. 2021), the specific time delay between star formation and SN explosion is accounted for all stellar masses. Thus, even in short time steps, star formation can be suppressed due to SNe from massive stars formed in the previous time step.</footnote>
the galaxy along the lines of constant sSFR, as illustrated by the black arrow in Fig. 4. This characteristic becomes important when comparing our predictions to observed galaxy properties. Since most observational data analysis pipelines assume Kroupa, Chabrier, or Salpeter IMFs, the stellar masses and SFRs inferred for galaxies with top-heavy IMFs are likely overestimated. These galaxies would then appear to shift along lines of constant sSFR to lower stellar mass and SFR values, which does not affect our criterion for identifying top-heavy IMFs<fnref n="7" />. For this reason and given the complexity of the inference of galaxy properties from observations, we refrained from adjusting the inferred values from the observations to our IMFs in Fig. 4.

**Comparison with observations.** In Fig. 4, we mark with black points where JADES-GS-z14-1, JADES-GS-z13-0 and JADES-GS-z12-0 (Curtis-Lake et al. 2023), UNCOVER-z13 and UNCOVER-z12 (Wang et al. 2023), as well as GLASS-z12 (Bakx et al. 2023) fall within the star formation main sequence of our simulations. Grey points represent the inferred values for JADES-z13-1-LA, both for their default model and an alternative model without nebular emission, with the latter leading to higher SFR and lower stellar mass values. Except for GLASS-z12, the only galaxies with SFR inferred from ALMA observations, JADES-GS-z12-0 is the only galaxy that has an SFR-$M_\star$ combination clearly indicative of a top-heavy in our evolving IMF model. The observed spectrum exhibiting no strong emission lines indeed suggests a high ionising photon production efficiency of $\log_{10}(\xi_{ion}/\text{erg}^{-1} \text{ Hz}) = 25.7$, which is close to the value for typical ISM conditions and normal stellar populations of 25.8, and a red UV slope of $\beta = -1.84$, which is possibly an indication of hot stars and a low-density ISM (Katz et al. 2024). Moreover, its low metallicity value of $\log_{10}(Z/Z_\odot) = -1.4$ is also in agreement with our model predictions. While the SFR-$M_\star$ combinations for JADES-z13-1-LA are inconclusive regarding a top-heavy IMF, its high ionising photon production efficiency ($\log_{10}(\xi_{ion}/\text{erg}^{-1} \text{ Hz}) = 26.5$) and extremely low metallicity ($\log_{10}(Z/Z_\odot) = -2.5$) strongly suggest the presence of such an IMF. Within the uncertainties, JADES-GS-z14-1 and UNCOVER-z13 could also be candidates for containing stellar populations with top-heavy IMFs given their inferred high SFRs of ~2 $M_\odot$/yr, but with metallicities at the upper end of what our simulated galaxies with top-heavy IMFs suggest<fnref n="8" />. However, UNCOVER-z12 and JADES-GS-z13-0 fall not into the range of galaxies with top-heavy IMFs. JADES-GS-z13-0, for instance, has an average UV slope of $\beta = -2.4$ and an ionising photon production efficiency of $\log_{10}(\xi_{ion}/\text{erg}^{-1} \text{ Hz}) = 25.5$, not necessarily pointing towards the extreme conditions associated with a top-heavy IMF.

Finally, we note that at $z = 11$, GN-z11 (Oesch et al. 2016; Bunker et al. 2023) also falls within the SFR-$M_\star$ combinations indicative of a top-heavy IMF in our evolving IMF model. This hypothesis of a top-heavy IMF also aligns with further observational evidence, including the elevated [N/O] abundance, the increased ionising emissivity ($\xi_{ion} = 25.7$) and blue UV slope ($\beta = -2.4$) (see e.g. Cameron et al. 2024; Vink 2023; Senchyna et al. 2024; Nandal et al. 2024).

<figure><img src="image_5.png" /><figcaption>**Fig. 5.** $M_{\text{UV}}-M_h$ relation for the Salpeter IMF (red) and the evolving IMF (blue) at $z = 13$. Hashed blue contours show the galaxies with an IMF top-heavy than the Salpeter IMF in the evolving IMF model.</figcaption></figure>

## 6. The UV luminosity scatter for top-heavy IMFs

In Sect. 3 we discussed how the maximum possible UV luminosity enhancement during the transition from a Salpeter to a top-heavy IMF varies with halo mass. In Fig. 5, we show the relation between UV luminosity ($M_{\text{UV}}$) and halo mass ($M_h$) for our Salpeter IMF (red contours) and evolving IMF model (blue contours) simulations. Consistent with our toy model, we find that simulation-based star formation and gas mass histories also lead to a UV luminosity enhancement of up to ~1 dex for SN feedback-limited galaxies ($M_h \lesssim 10^9 M_\odot$). This enhancement increases towards more massive galaxies, reaching up to ~2.5 dex, as the transition from the SN feedback-limited to constant SFE regime shifts to higher masses with rising $f_{\text{massive}}$ values. Additionally, the evolving IMF model produces by up to ~0.5 dex fainter galaxies at the same halo mass. Stronger and more immediate SN feedback in galaxies with top-heavy IMFs suppresses star formation more quickly than in galaxies with a Salpeter IMF, due to the presence of more massive stars with shorter lifetimes (see also Cueto et al. 2024), making these galaxies fainter in the UV.

Tables 3 and 4 present the median $M_{\text{UV}}$ values along with their $1\sigma$ and $2\sigma$ uncertainties across various halo masses for the Salpeter and evolving IMF models, respectively. While the median $M_{\text{UV}}-M_h$ relations hardly differ between the two IMF models, the evolving IMF model exhibits a larger variance in UV luminosity than the Salpeter IMF model, with $1\sigma_{\text{UV}}$ and $2\sigma_{\text{UV}}$ values being by about a factor 1.1–1.3 and 1.4–2 larger. Although the $\sigma_{\text{UV}}$ values in both IMF models are lower than those inferred from the FIRE simulations ($\sigma_{\text{UV}} \sim 1.6$ to 0.9 from $M_h = 10^{8.5} M_\odot$ to $10^{10.5} M_\odot$; Sun et al. 2023) or those noted as necessary to reproduce the $z > 10$ UV LFs in analytic models ($\sigma_{\text{UV}} \approx 2$ at $z = 13$; Mason et al. 2023; Muñoz et al. 2023; Shen et al. 2023; Kravtsov & Belokurov 2024), the boost in the $2\sigma$ variances when transitioning from the Salpeter to the evolving IMF model are significantly larger and align better with the variances found in the aforementioned analytic models that successfully reproduce the $z > 10$ UV LFs. This indicates that the UV luminosity distribution is stretched more at the extreme

<footnote n="7">⁷ Here, we do not consider potential biases in SFR and stellar mass estimates that may arise from deriving these quantities from different wavelength ranges of galaxy spectra or emission lines.</footnote>

<footnote n="8">⁸ JADES-GS-z14-1 shows no emission lines, indicating either very low metallicity or high ionising escape fractions, and its very blue UV slope ($\beta = -2.7$) points to a young stellar population with low dust attenuation.</footnote>
<table><caption>**Table 3.** UV luminosity median and variance values for the Salpeter IMF.</caption><thead><tr><th>log<sub>10</sub> M<sub>h</sub>/M<sub>☉</sub></th><th>M<sub>UV</sub></th><th>σ<sub>UV</sub></th><th>2σ<sub>UV</sub></th></tr></thead><tbody><tr><td>8.5</td><td>-12.7</td><td>1.22</td><td>5.7</td></tr><tr><td>9</td><td>-14.6</td><td>0.91</td><td>1.71</td></tr><tr><td>9.5</td><td>-16.6</td><td>0.84</td><td>1.59</td></tr><tr><td>10</td><td>-18.1</td><td>0.78</td><td>1.34</td></tr><tr><td>10.5</td><td>-19.7</td><td>0.97</td><td>0.97</td></tr></tbody></table>

<table><caption>**Table 4.** UV luminosity median and variance values for the evolving IMF.</caption><thead><tr><th>log<sub>10</sub> M<sub>h</sub>/M<sub>☉</sub></th><th>M<sub>UV</sub></th><th>σ<sub>UV</sub></th><th>2σ<sub>UV</sub></th></tr></thead><tbody><tr><td>8.5</td><td>-12.7</td><td>1.25</td><td>6.2</td></tr><tr><td>9</td><td>-14.6</td><td>0.97</td><td>2.47</td></tr><tr><td>9.5</td><td>-16.4</td><td>0.91</td><td>2.56</td></tr><tr><td>10</td><td>-18.0</td><td>1.00</td><td>2.63</td></tr><tr><td>10.5</td><td>-19.7</td><td>1.94</td><td>1.94</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** $M_{UV}-M_{\star}$ relation for the Salpeter IMF (red) and the evolving IMF (blue) at $z = 13$. Hashed blue contours show the galaxies with an IMF top-heavier than the Salpeter IMF in the evolving IMF model.</figcaption></figure>

values ($>1\sigma$) rather than near the median ($<1\sigma$). This result stems from the assumption that only the most gas-rich galaxies in relatively denser regions exhibit top-heavy IMFs and, thus, enhanced UV luminosities. While this boost in $>1\sigma$ values in clustered locations may enhance the clustering signal, the increased variance may also lead to a decrease in the clustering signal as shown in Gelli et al. (2024). In future works, we will investigate whether the evolving IMF leaves a characteristic imprint in large-scale clustering.

Figure 6 shows the relation between the UV luminosity ($M_{UV}$) and stellar mass ($M_{\star}$). Across all stellar masses where galaxies with top-heavy IMFs exist, we find that the maximum UV enhancement matches the $\Delta M_{UV}$ of ~2.5 dex derived from our toy model. This stands in contrast to the UV luminosity-halo mass relation, where the enhancement increases with rising halo mass. Since the UV luminosity strongly depends on the SFR, this

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Redshift evolution of the global IGM H1 fraction for the Salpeter IMF (red) and evolving IMF (blue) model (bottom) and the ratio between the mass-averaged and volume-averaged IGM H1 fractions (top). Grey points indicate observational constraints from: GRB optical afterglow spectrum analyses (light pentagons; Totani et al. 2006, 2014), quasar sightlines (faint squares; Fan et al. 2006), (medium bright circles; Bosman et al. 2022) Lyman-$\alpha$ LFs (dark upward triangles; Konno et al. 2018; Kashikawa et al. 2011; Ouchi et al. 2010; Ota et al. 2010; Malhotra & Rhoads 2004), Lyman-$\alpha$ emitter clustering (dark downwards triangles; Ouchi et al. 2010), the Lyman-$\alpha$ emitting galaxy fraction (dark leftwards triangles; Pentericci et al. 2011; Schenker et al. 2012; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014; Pentericci et al. 2014), Lyman-$\alpha$ equivalent widths (dark rightwards triangles; Mason et al. 2018, 2019; Bolan et al. 2022), dark pixels (light squares; McGreer et al. 2015), and damping wings (light diamonds; Davies et al. 2018; Greig et al. 2019).</figcaption></figure>

difference arises because the SFR-$M_{\star}$ relation remains largely unchanged when top-heavy IMFs are included, due to the self-similar growth of galaxy mass (see Fig. 4 and Cueto et al. 2024).

## 7. Implications for reionisation

The history and morphology of the reionisation of the neutral hydrogen in the IGM are highly sensitive to the ionising photon output and distribution of the first galaxies. We examine the reionisation process for both IMF models, assuming that the fraction of ionising photons escaping from a galaxy into the IGM scales with the fraction of gas ejected from the galaxy, expressed as $f_{esc} = f_{esc}^0 \min(1, f_{\star}^{\text{eff}}/f_{\star}^{\text{ej}})$ (see e.g. Hutter et al. 2021, 2023), with $f_{esc}^0 = 0.35$. The corresponding reionisation histories, the redshift evolution of the global IGM H1 neutral fraction, $\langle\chi_{HI}\rangle$, are shown in Fig. 7.

Raising the fraction of massive stars increases the ionising emissivity of stellar populations by up to an order of magnitude, similar to the UV luminosity. Consequently, as more galaxies form stars according to a top-heavy IMF at high redshift, the evolving IMF model (blue line) results in an earlier onset of reionisation ($\Delta z \approx 2$ at $\langle\chi_{HI}\rangle \approx 0.99$) compared
to the Salpeter IMF model (red line), characterised by ionised regions growing earlier around the massive galaxies with top-heavy IMFs. At the same global IGM neutral fractions, $\langle\chi_{HI}\rangle$, the ionisation morphology in the evolving IMF model shows more extended large ionised regions and fewer small ones. For example, at $z \simeq 10$, the largest ionised regions in the evolving IMF model extend up to ~10 cMpc, while in the Salpeter IMF model, they barely reach 5 cMpc. As the fraction of galaxies with top-heavy IMFs decreases, the differences in $\langle\chi_{HI}\rangle$ and ionisation morphology diminish: By $z \simeq 7$, the ionisation morphologies become very similar, and the redshift difference at the same $\langle\chi_{HI}\rangle$ values reduces to 0.2. In both IMF models, reionisation completes by $z \simeq 5.6$, consistent with constraints from QSO sightlines (e.g. Bosman et al. 2022). Due to an earlier start of reionisation in the evolving IMF model, the electron optical depth is slightly higher ($\sigma = 0.055$) compared to the Salpeter IMF model ($\sigma = 0.052$), but remains within the limits set by the cosmic microwave background measurements with Planck (Planck Collaboration VI 2020). A more detailed analysis of the differences in ionisation morphology will be presented in future work.

## 8. Observational signatures of a top-heavy IMF

Given the crucial role of the IMF in determining the ionising photon output and, consequently, the reionisation history and morphology, we briefly discuss how the ionising emissivity, $\xi_{ion} = \dot{Q}/L_{UV}^{int}$, differs at $z \gtrsim 10$ in our two IMF models. $\xi_{ion}$ has been inferred from galaxy observations by comparing the luminosity of emission lines (e.g. Hα, OIII) to dust-attenuation-corrected intrinsic UV continua (see e.g. Prieto-Lyon et al. 2023). Since these emission lines are produced by the excitation of hydrogen or oxygen gas from ionising radiation that does not escape the galaxy, the values inferred from observations correspond actually to $\xi_{ion}(1 - f_{esc})$. Therefore, our predictions for $\xi_{ion}$ can be considered as upper limits on the observationally inferred values, albeit only in the case where the assumed dust attenuation of the UV remains below its actual values. Figure 8 shows the relationship between $\xi_{ion}$ and observed UV luminosity ($M_{UV}$) for the Salpeter IMF (red contours) and the evolving IMF (blue contours).

We can see that the maximum $\xi_{ion}$ value is approximately $10^{25.4}$ Hz erg$^{-1}$ when assuming a Salpeter IMF. However, when adopting a top-heavy IMF, ranging up to a log-flat IMF as in the evolving IMF model, $\xi_{ion}$ can reach higher values for UV-bright galaxies, up to around $10^{25.7}$ Hz erg$^{-1}$. The higher abundance of massive stars in a top-heavy IMF shifts stellar emission to higher energies, thereby raising the ratio of H I ionising photons to UV continuum photons (and thus $\xi_{ion}$) for highly star-forming, UV-bright galaxies. Compared to the Salpeter IMF, $\xi_{ion}$ values in the evolving IMF model reach about 2× higher values. Here, it is important to note that our $\xi_{ion}$ values correspond to results derived under the assumption of continuous star formation over the specified time step (~11 Myrs at $z = 13$). If star formation was stochastic on shorter time scales, $\xi_{ion}$ could reach higher values during periods of intense starbursts, presumably up to $10^{26.1}$ Hz erg$^{-1}$ and $10^{25.8}$ Hz erg$^{-1}$ for the evolving and Salpeter IMF, respectively.

Hence, $\xi_{ion}$ values exceeding $10^{25.8}$ Hz erg$^{-1}$ are indicative of star formation driven by a top-heavy IMF. For example, JADES-z12-0 (Curtis-Lake et al. 2023) and JADES-z13-1-LA ($\xi_{ion} \simeq 10^{26.5}$ Hz erg$^{-1}$; Witstok et al. 2024) (for which our evolving IMF model predicts a top-heavy IMF based on their SFR and stellar mass values) show $\xi_{ion}$ values that are consistent

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** $\xi_{ion}-M_{UV}$ relation for the Salpeter IMF (red) and the evolving IMF (blue) at $z = 13$. Solid lines show the median of the corresponding coloured contours. Black points show the inferred values for JADES-GS-z13-0 and JADES-GS-z12-0 (Curtis-Lake et al. 2023), with JADES-GS-z12-0 showing higher $\xi_{ion}$ values that agree more with the evolving IMF than Salpeter IMF model.</figcaption></figure>

with our evolving IMF model, and in the case of JADES-z13-1-LA even exceed them. Observations at slightly lower redshifts indicate that $\xi_{ion}$ tends to increase with redshift (Simmonds et al. 2024a,b), with some galaxies at $z \simeq 7-9$ surpassing the $\xi_{ion}$ threshold for a top-heavy IMF (Endsley et al. 2023; Whitler et al. 2024; Mascia et al. 2024), albeit not reaching the extreme value observed for JADES-z13-1-LA at $z = 13$. This trend is consistent with our evolving IMF model, which predicts galaxies with top-heavy IMFs becoming more prevalent at earlier cosmic times.

## 9. Conclusions

We investigate whether an IMF parameterisation that links the top-heaviness of the IMF to the sSFR of a galaxy, using the gas-to-DM mass ratio as a proxy, can explain the observed abundance of bright $z > 10$ galaxies. Specifically, we consider an IMF that becomes only top-heavy when the sSFR exceeds a certain threshold. To model such a varying IMF, we adopted the IMF shape from Cueto et al. (2024), using a Salpeter IMF below and a log-flat IMF above a stellar mass, $M_c$. Increasing $M_c$ results in a more top-heavy IMF, characterised by a higher fraction of stellar mass produced in massive stars ($M > M_c$), denoted as $f_{massive}$. Using a toy model, we compared how this evolving IMF affects galaxy properties, particularly UV luminosities, and their statistical distributions, relative to the standard assumption of a constant Salpeter IMF. Our findings are as follows:

1. A more top-heavy IMF reduces the SFE of SN feedback-limited galaxies by the factor by which the energy released by SNe is enhanced. Since the SFE of SN feedback-limited galaxies decreases with the rising abundance of massive stars, the UV luminosity enhancement for such galaxies does not exceed a factor of ~3.5, corresponding to a shift of $\Delta M_{UV} = 1.3$. This enhancement level is nearly achieved with a massive star fraction of $f_{massive} \approx 0.5$. Further increasing $f_{massive}$ to 1 results in only an additional 10% enhancement.
However, increasing $f_{\text{massive}}$ shifts the halo mass below which star formation becomes SN feedback-limited to higher masses by ~0.8 dex for a change of $\Delta f_{\text{massive}} = 1$.

2. For more massive galaxies, however, the SFE remains unchanged, leading to a continuous increase in UV luminosity as $f_{\text{massive}}$ rises, reaching a maximum enhancement of ~11, which corresponds to a shift of $\Delta M_{\text{UV}} = 2.6$.

3. Our evolving IMF model successfully reproduces the observed UV LF at $z = 5-15$ due to a growing fraction of galaxies forming stars according to a top-heavy IMF and these IMFs becoming progressively top-heavy towards increasing redshifts. The shallower slopes of the $z > 10$ UV LFs in the evolving IMF model arise from the increasing gas fraction (triggering star formation with top-heavy IMFs) and UV luminosity enhancement with rising halo mass.

4. In our evolving IMF model, UV-bright galaxies at $z > 10$ fall into two categories: galaxies with the highest stellar and dust masses and galaxies with lower stellar and dust masses whose brightness results from their stars forming according to a top-heavy IMF.

5. Galaxies that form stars according to a top-heavy IMF are typically the most gas-rich and highly star-forming, with sSFR $\ge$ 25 Gyr$^{-1}$. These galaxies are located in environments where they form the deepest potential wells or most massive galaxies, allowing for high gas accretion rates replenishing their gas reservoirs and counteracting any losses attributed to SNe.

6. Adopting our evolving IMF model increases the $1\sigma$ UV luminosity variance (for a given halo mass) only moderately, not sufficient to reproduce the $z > 10$ UV LFs as shown in Gelli et al. (2024), Muñoz et al. (2023), Kravtsov & Belokurov (2024). However, it is the significant increase (by a factor 1.4–2) in the $2\sigma$ UV luminosity variance, driven by the UV luminosity enhancement of the most star-forming galaxies, that reproduces the high abundance of UV-bright $z > 10$ galaxies.

7. Adopting an evolving IMF causes reionisation to start earlier without altering completion time. It also modifies the ionisation morphology during the early stages of reionisation, leading to more extended large ionised regions and fewer small ones.

Along with the work of Cueto et al. (2024), our results highlight how the dependence of the IMF on galaxy properties affects the abundance of UV-bright galaxies at $z > 10$. Specifically, these findings suggest that to reproduce the UV LFs at $z = 5-15$ with a top-heavy IMF, the IMF needs to be more top-heavy for massive than less-massive galaxies. If the IMF is less top-heavy in massive galaxies, the UV luminosity enhancement becomes similar across all galaxy masses, leading to a UV LF slope at $z > 10$ that is steeper than observed. The result of increasing the top-heaviness of the IMF with rising sSFR or gas content is consistent with findings in local observations of star-forming regions and in broad agreement with the IMF trends proposed in the integrated galaxy-wide IMF theory (e.g. Weidner et al. 2013; Jeřábková et al. 2018); namely, it results in a larger fraction of more massive galaxies with a top-heavy IMF, while reproducing the high-redshift UV LFs without necessitating an increase in the SFE or stochasticity beyond what is accounted for by our standard model. Interestingly, according to our model, observed galaxies that classify as those with a top-heavy IMF also display signatures of massive stars, such as elevated ionising emissivities and [N/O] abundances. Furthermore, adopting a top-heavy IMF in galaxies within overdense regions (where the gas accretion rates and SFRs are highest) may help explain the observed clustering bias,

which exceeds simple analytic estimates (e.g. Gelli et al. 2024). However, the greater variance in UV luminosity could offset this effect, which we will explore in future work. In this context, it is interesting to note that simulations suggest that mergers can trigger higher sSFR values or starbursts (e.g. Croton et al. 2016; Patton et al. 2020). Thus, top-heavy IMFs could be predominantly found in merging galaxies, particularly since the galaxy merger rate increases towards higher redshifts (e.g. Gottlöber et al. 2001). We plan to test this hypothesis in subsequent work.

Recently, van Dokkum & Conroy (2024) proposed a concordance IMF for massive galaxies at $z \approx 0-5$, which increasingly becomes both bottom-heavy and top-heavy compared to the Salpeter IMF as the galaxy’s velocity dispersion increases. They suggest this IMF could help explain the bright galaxies at $z > 7$. This IMF exhibits a similar trend to our evolving IMF, where more compact (and likely more star-forming) galaxies have more top-heavy IMFs. However, while the shape of the concordance IMF is similar to that of our evolving IMF in highly star-forming galaxies, its steeper high-mass slope leads to a lower UV boost. Whether this reduced UV boost can be compensated by a higher fraction of galaxies with top-heavy IMFs remains to be tested.

Finally, we note the main limitations of this study. Firstly, we do not adjust the SFE of massive galaxies with a top-heavy IMF. Simulations of star-forming clouds suggest that a higher radiation pressure on dust grains, caused by an increased abundance of massive stars, can reverse the gas accretion and drive gas outward, reducing the SFE by nearly a factor of two at solar metallicity ($Z = Z_\odot$) when the UV emissivity is enhanced by an order of magnitude (Menon et al. 2024). However, this reduction in SFE decreases as the gas becomes more metal-poor (also resulting in a lower dust abundance), reaching values similar to those in our simulated galaxies ($Z \approx 0.01 Z_\odot$); then, it nearly disappears in dense clumps ($\Sigma_{\text{cl}} \approx 3 \times 10^4 M_\odot \text{ pc}^{-2}$). If the SFE of massive galaxies were decreased, it would likely reduce the UV luminosity enhancement in these galaxies. This effect could potentially be counterbalanced by lowering the sSFR threshold above which our model adopts a top-heavy IMF. We plan to explore this in future work.

Secondly, this study explores a specific IMF parameterisation whose shape is motivated by the star-forming cloud simulations by Chon et al. (2022). Alternative parameterisations, such as those with an increasingly shallower slope to represent a more top-heavy IMF, could lead to a different ratio of massive to less massive stars, potentially altering the UV luminosity enhancement. Another scenario, which is increasingly supported by observations suggesting the presence of very massive stars (>100 $M_\odot$) at sub-solar metallicities (e.g. Martins et al. 2023; Meštrić et al. 2023; Wofford et al. 2023; Upadhyaya et al. 2024), posits that the IMF extends to masses larger than the 100 $M_\odot$ assumed here. This extension could also result in a UV luminosity enhancement (Martins & Palacios 2022) and, according to Katz et al. (2024), it would leave more distinct signatures in the galaxy spectra, such as redder UV slopes of $\beta > -2$.

Thirdly, the spectra for our stellar populations do not account for the effects of binary interactions on stellar evolution. However, massive stars, likely to form in dense regions, are expected to often exist in binaries. Nevertheless, as we explored in Hutter et al. (2021), using stellar population synthesis models that include binary interactions, such as BPASS, hardly changes the UV luminosities of galaxies. They only delay the decrease of the ionising photon output following a starburst and could suppress star formation longer due to the associated radiative feedback – which would effectively lead
to an adjustment of ASTRAEUS's SFE and SNe wind coupling parameters.

In summary, our evolving IMF model appears well-aligned with current high-z galaxy observations and observations of local galaxies and star-forming regions. Thus, a top-heavy IMF in strongly star-forming galaxies offers one possible explanation for the high abundance of UV-bright galaxies at z > 10. Yet, alternative explanations also remain plausible, such as feedback-free starbursts, bursty star formation, cosmic variance, contributions from black holes, or a combination of these factors, including different IMF scenarios. An analysis of the UV-luminosity-weighted galaxy clustering bias and gas elemental abundances, such as [N/O] and [C/O], will offer further insights into the applicability of the evolving IMF model.

**Acknowledgements.** We thank the anonymous referee for their careful reading of the manuscript and their comments, which improved the quality of the paper. AH acknowledges support by the VILLUM FONDEN under grant 37459. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. PD acknowledges support from NWO grant 016.VIDI.189.162 ("ODIN") and warmly thanks the European Commission and the University of Groningen's CO-FUND Rosalind Franklin program. GY acknowledges Ministerio de Ciencia e Innovación (Spain) for partial financial support under research grant PID2021-122603NB-C21. The VSMDPL simulation has been performed at LRZ Munich within the project pr87yi. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (www.lrz.de). The CosmoSim database (www.cosmosim.org) provides access to the simulation and the Rockstar data. The database is a service by the Leibniz Institute for Astrophysics Potsdam (AIP). This research made use of matplotlib, a Python library for publication-quality graphics (Hunter 2007); and the Python library numpy (Oliphant 2006).

## References

Abel, T., Bryan, G. L., & Norman, M. L. 2002, *Science*, **295**, 93

Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, **518**, 4755

Adams, N. J., Conselice, C. J., Austin, D., et al. 2024, ApJ, **965**, 169

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, Nature, **622**, 707

Atek, H., Richard, J., Kneib, J.-P., et al. 2015, ApJ, **800**, 18

Atek, H., Richard, J., Kneib, J.-P., & Schaerer, D. 2018, MNRAS, **479**, 5184

Atek, H., Chemerynska, I., Wang, B., et al. 2023, MNRAS, **524**, 5486

Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Nature, **626**, 975

Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJ, **952**, L7

Bakx, T. J. L. C., Zavala, J. A., Mitsuhashi, I., et al. 2023, MNRAS, **519**, 5076

Barkana, R., & Loeb, A. 2001, *Phys. Rep.*, **349**, 125

Bolan, P., Lemaux, B. C., Mason, C., et al. 2022, MNRAS, **517**, 3263

Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, **514**, 55

Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, **162**, 47

Bouwens, R., Illingworth, G., Oesch, P., et al. 2023a, MNRAS, **523**, 1009

Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023b, MNRAS, **523**, 1036

Boylan-Kolchin, M. 2023, *Nature Astronomy*, **7**, 731

Bradley, L. D., Coe, D., Brammer, G., et al. 2023, ApJ, **955**, 13

Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, **564**, 23

Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, **677**, A88

Cameron, A. J., Katz, H., Witten, C., et al. 2024, MNRAS, **534**, 523

Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2012, MNRAS, **427**, 3055

Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2014, MNRAS, **443**, 2831

Castellano, M., Fontana, A., Paris, D., et al. 2010, A&A, **524**, A28

Charbonnel, C., Schaerer, D., Prantzos, N., et al. 2023, A&A, **673**, L7

Chon, S., Ono, H., Omukai, K., & Schneider, R. 2022, MNRAS, **514**, 4639

Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011, ApJ, **727**, 110

Croton, D. J., Stevens, A. R. H., Tonini, C., et al. 2016, ApJS, **222**, 22

Cueto, E. R., Hutter, A., Dayal, P., et al. 2024, A&A, **686**, A138

Cullen, F., McLeod, D. J., McLure, R. J., et al. 2024, MNRAS, **531**, 997

Curti, M., Witstok, J., Jakobsen, P., et al. 2024, A&A, submitted [arXiv:2407.02575]

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, **7**, 622

Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, **864**, 142

Dayal, P., & Ferrara, A. 2018, *Phys. Rep.*, **780**, 1

Dayal, P., Ferrara, A., Dunlop, J. S., & Pacucci, F. 2014, MNRAS, **445**, 2545

Dayal, P., Ferrara, A., Sommovigo, L., et al. 2022, MNRAS, **512**, 989

Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, **523**, 3201

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023a, MNRAS, **518**, 6011

Donnan, C. T., McLeod, D. J., McLure, R. J., et al. 2023b, MNRAS, **520**, 4554

Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, MNRAS, **533**, 3222

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, **444**, 2960

Ekström, S. 2021, *Frontiers in Astronomy and Space Sciences*, **8**, 53

Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, **524**, 2312

Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, **132**, 117

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, **522**, 3986

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, **946**, L13

Fiore, F., Ferrara, A., Bischetti, M., Feruglio, C., & Travascio, A. 2023, ApJ, **943**, L27

Fontanot, F., De Lucia, G., Hirschmann, M., et al. 2017, MNRAS, **464**, 3812

Fontanot, F., De Lucia, G., Xie, L., et al. 2018, MNRAS, **475**, 2467

Fudamoto, Y., Smit, R., Bowler, R. A. A., et al. 2022, ApJ, **934**, 144

Fujimoto, S., Silverman, J. D., Bethermin, M., et al. 2020, ApJ, **900**, 1

Fukushima, H., Hosokawa, T., Chiaki, G., et al. 2020, MNRAS, **497**, 829

Gelli, V., Mason, C., & Hayward, C. C. 2024, ApJ, **975**, 192

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, **735**, L34

Gottlöber, S., Klypin, A., & Kravtsov, A. V. 2001, ApJ, **546**, 223

Greig, P., Mesinger, A., & Bañados, E. 2019, MNRAS, **484**, 5094

Gunawardhana, M. L. P., Hopkins, A. M., Sharp, R. G., et al. 2011, MNRAS, **415**, 1647

Harikane, Y., Inoue, A.K., Ellis, R.S., et al. 2024, ApJ, accepted [arXiv:2406.18352]

Haslbauer, M., Kroupa, P., Zonoozi, A. H., & Haghi, H. 2022, ApJ, **939**, L31

Heintz, K. E., Giménez-Arteaga, C., Fujimoto, S., et al. 2023, ApJ, **944**, L30

Hunter, J. D. 2007, *Computing In Science & Engineering*, **9**, 90

Hutter, A. 2018, MNRAS, **477**, 1549

Hutter, A., Dayal, P., Yepes, G., et al. 2021, MNRAS, **503**, 3698

Hutter, A., Trebitsch, M., Dayal, P., et al. 2027, MNRAS, **524**, 6124

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, **854**, 73

Isobe, Y., Ouchi, M., Tominaga, N., et al. 2023, ApJ, **959**, 100

Jefábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, A&A, **620**, A39

Kashikawa, N., Shimasaku, K., Matsuda, Y., et al. 2011, ApJ, **734**, 119

Katz, H., Cameron, A. J., Saxena, A., et al. 2024, The Open Journal of Astrophysics, submitted [arXiv:2408.03189]

Kobayashi, C., & Ferrara, A. 2024, ApJ, **962**, L6

Kobayashi, C., Karakas, A. I., & Lugaro, M. 2020, ApJ, **900**, 179

Konno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, **70**, S16

Kravtsov, A., & Belokurov, V. 2024, ArXiv e-prints [arXiv:2405.04578]

Kroupa, P., & Weidner, C. 2003, ApJ, **598**, 1076

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, **616**, 266

Legrand, L., Hutter, A., Dayal, P., et al. 2022, MNRAS, **509**, 595

Legrand, L., Dayal, P., Hutter, A., et al. 2023, MNRAS, **519**, 4564

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, **123**, 3

Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, **835**, 113

Lu, S., Frenk, C. S., Bose, S., et al. 2025, MNRAS, **536**, 1018

Malhotra, S., & Rhoads, J. E. 2004, ApJ, **617**, L5

Martins, F., & Palacios, A. 2022, A&A, **659**, A163

Martins, F., Schaerer, D., Marques-Chaves, R., & Upadhyaya, A. 2023, A&A, **678**, A159

Mascia, S., Pentericci, L., Calabrò, A., et al. 2024, A&A, **685**, A3

Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, **856**, 2

Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, **485**, 3947

Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, **521**, 497

Mauerhofer, V., & Dayal, P. 2023, MNRAS, **526**, 2196

McGreer, I. D., Mesinger, A., & D'Odorico, V. 2015, MNRAS, **447**, 499

McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, **450**, 3032

McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, **459**, 3812

McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, **527**, 5004

McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, **432**, 2696

Menon, S. H., Lancaster, L., Burkhart, B., et al. 2024, ApJ, **967**, L28

Meurer, G. R., Wong, O. I., Kim, J. H., et al. 2009, ApJ, **695**, 765

Meštrić, V., Vanzella, E., Upadhyaya, A., et al. 2023, A&A, **673**, A50

Mirocha, J., & Furlanetto, S. R. 2023, MNRAS, **519**, 843

Muñoz, J. B., Mirocha, J., Furlanetto, S., & Sabti, N. 2023, MNRAS, **526**, L47

Nagele, C., & Umeda, H. 2023, ApJ, **949**, L16

Nakane, M., Ouchi, M., Nakajima, K., et al. 2024, ApJ, **976**, 122

Nandal, D., Sibony, Y., & Tsiatsiou, S. 2024, A&A, **688**, A142

Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, **819**, 129

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanov, M. 2018, ApJ, **855**, 105

Oke, J. B., & Gunn, J. E. 1983, ApJ, **266**, 713

Oliphant, T. 2006, NumPy: A guide to NumPy (USA: Trelgol Publishing)

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, **744**, 83

Ota, K., Iye, M., Kashikawa, N., et al. 2010, ApJ, **722**, 803
Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869

Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6

Papadopoulos, P. P., Thi, W.-F., Miniati, F., & Viti, S. 2011, MNRAS, 414, 1705

Patton, D. R., Wilson, K. D., Metrow, C. J., et al. 2020, MNRAS, 494, 4969

Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, ApJ, 743, 132

Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, ApJ, 793, 113

Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJ, 951, L1

Planck Collaboration VI. 2020, A&A, 641, A6

Planck Collaboration Int. XLVI. 2016, A&A, 596, A107

Prieto-Lyon, G., Strait, V., Mason, C. A., et al. 2023, A&A, 672, A186

Rinaldi, P., Caputi, K. I., Iani, E., et al. 2024, ApJ, 969, 12

Salpeter, E. E. 1955, ApJ, 121, 161

Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4

Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179

Schenker, M. A., Robertson, B. E., Ellis, R. S., et al. 2013, ApJ, 768, 196

Senchyna, P., Plat, A., Stark, D. P., et al. 2024, ApJ, 966, 92

Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254

Simmonds, C., Tacchella, S., Hainline, K., et al. 2024a, MNRAS, 527, 6139

Simmonds, C., Tacchella, S., Hainline, K., et al. 2024b, MNRAS, 535, 2998

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023, ApJ, 955, L35

Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, ApJ, submitted [arXiv:2407.19009]

Totani, T., Kawai, N., Kosugi, G., et al. 2006, PASJ, 58, 485

Totani, T., Aoki, K., Hattori, T., et al. 2014, PASJ, 66, 63

Treu, T., Trenti, M., Stiavelli, M., Auger, M. W., & Bradley, L. D. 2012, ApJ, 747, 27

Trinca, A., Schneider, R., Valiante, R., et al. 2024, MNRAS, 529, 3563

Ucci, G., Dayal, P., Hutter, A., et al. 2023, MNRAS, 518, 3557

Upadhyaya, A., Marques-Chaves, R., Schaerer, D., et al. 2024, A&A, 686, A185

van Dokkum, P., & Conroy, C. 2024, ApJ, 973, L32

Vanzella, E., Loiacono, F., Bergamini, P., et al. 2023, A&A, 678, A173

Vink, J. S. 2023, A&A, 679, L9

Wang, B., Fujimoto, S., Labbé, I., et al. 2023, ApJ, 957, L34

Watanabe, K., Ouchi, M., Isobe, Y., & Tominaga, N. 2024a, EAS2024, 1243

Watanabe, K., Ouchi, M., Nakajima, K., et al. 2024b, ApJ, 962, 50

Watts, A. B., Meurer, G. R., Lagos, C. D. P., et al. 2018, MNRAS, 477, 5554

Weidner, C., Kroupa, P., Pflamm-Altenburg, J., & Vazdekis, A. 2013, MNRAS, 436, 3309

Whitler, L., Stark, D. P., Endsley, R., et al. 2024, MNRAS, 529, 855

Witstok, J., Jakobsen, P., Maiolino, R., et al. 2024, ArXiv e-prints [arXiv:2408.16608]

Wofford, A., Sixtos, A., Charlot, S., et al. 2023, MNRAS, 523, 3949

Yan, Z., Jerabkova, T., & Kroupa, P. 2017, A&A, 607, A126

Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2024, MNRAS, 527, 5929

Zhang, Z.-Y., Romano, D., Ivison, R. J., Papadopoulos, P. P., & Matteucci, F. 2018, Nature, 558, 260

Ziparo, F., Ferrara, A., Sommovigo, L., & Kohandel, M. 2023, MNRAS, 520, 2445
# Appendix A: Time evolution of emissivities

For an IMF described by Eq.1, we first derive fitting functions to model the time evolution of ionising and ultraviolet emissivities for a 1M☉ starburst, considering all available combinations of $f_{massive}$ and Z values ($f_{massive} \in [0, 1]$, Z ∈ [0.001, 0.008]). For metallicity values that are lower or larger than the covered interval, the respective values for the lower or upper limit are assumed.

## A.1. Ionising emissivity

We approximate the time evolution of the number of photons above an energy of 13.6 eV per second for a starburst of 1 M☉ with the following fitting function:

$$ \frac{N_{\text{ion}}}{\text{s}^{-1}} = \begin{cases} N_0, & t \le t_{\text{break}}, \\ N_0 (t/t_{\text{break}})^{-\gamma}, & t > t_{\text{break}}, \end{cases} \qquad (A.1) $$

where the normalisation $N_0$, the break $t_{break}$, and slope $\gamma$ are given by

$$ \lg\left(\frac{t_{\text{break}}}{\text{Myr}}\right) = \begin{cases} \max[0.49(1-23.37Z) + 2.64(f_{\text{massive}} - 0.209)^2, 0.477], & f_{\text{massive}} < 0.209 \\ \max[0.47(1-17.15Z) + 0.06(1 - e^{-7.5(f_{\text{massive}}-0.212)}), 0.477], & f_{\text{massive}} \ge 0.209 \end{cases} \qquad (A.2) $$

$$ \gamma = \begin{cases} -5.07(1 + 2.16Z) - 0.49(1 + 8.96Z), & f_{\text{massive}} \ge 0.7 - 0.5Z \\ -5.07(1 + 2.16Z) + 0.49(1 + 8.96Z) \\ \times \cos[-5.15f_{\text{massive}} - 6.16(1 + 2048Z)], & f_{\text{massive}} < 0.7 - 0.5Z \end{cases} \qquad (A.3) $$

$$ \lg N_0 = 3.4 \times 10^{-9}(1 + 0.62Z) + 53.66(1 - 0.32Z) \\ \times [f_{\text{massive}} + 0.051(1 + 227.6Z)]^{0.00766(1+50.01Z)}. \qquad (A.4) $$

## A.2. Ultraviolet luminosity

We approximate the time evolution of the UV luminosity for a starburst of 1 M☉ with the following fitting function

$$ \lg\left(\frac{L_{\text{UV}}}{\text{erg s}^{-1}\text{Å}^{-1}}\right) = \begin{cases} \lg N_0, & t \le 1.26 \text{ Myr}, \\ \lg N_0 + 26.08 \lg N_1 \lg\left(\frac{t}{1.26 \text{ Myr}}\right), & 1.26 \text{ Myr} < t \le 3.5 \text{ Myr}, \\ \lg N_0 + \lg N_2 + \gamma \lg\left(\frac{t}{3.5 \text{ Myr}}\right), & t > 3.5 \text{ Myr}, \end{cases} \qquad (A.5) $$

where $\gamma$, $N_0$, $N_1$ and $N_2$ are given by

$$ \gamma = \begin{cases} -1.920 + 0.123e^{-2\times10^3(Z-0.001)} - 0.564, & f_{\text{massive}} \le 0.7, \\ -1.920 + 0.123e^{-2\times10^3(Z-0.001)} + 0.564\cos(-4.155f_{\text{massive}} - 0.385), & f_{\text{massive}} > 0.7, \end{cases} \qquad (A.6) $$

$$ \lg N_0 = 36.0(1 + 0.1Z) + 3.82(1 + 1.1Z)(f_{\text{massive}} + 0.13)^{0.146} - 6, \qquad (A.7) $$

$$ \lg N_2 = \begin{cases} \max[-0.072(1 + 98.64Z) + 4.92(f_{\text{massive}} - 0.252)^2, 0], & f_{\text{massive}} < 0.252, \\ 0, & 0.252 \le f_{\text{massive}} \le 0.452, \\ \max[0.023(1 + -53691Z^2) + 0.06(1 - e^{-7.5(f_{\text{massive}}-0.452)}), 0], & f_{\text{massive}} > 0.452, \end{cases} \qquad (A.8) $$

$$ \lg N_1 = \begin{cases} 0.064 [1 + 8166(Z - 0.0035)^2] e^{-3.29 [f_{\text{massive}} - 0.5(1+9.89(Z-0.006)^2)]^2}, & f_{\text{massive}} \le 0.283, \\ 0.044 [1 + 9311(Z - 0.0041)^2] + 0.009 [1 + 24847(Z - 0.002)^2], & f_{\text{massive}} > 0.283. \end{cases} \qquad (A.9) $$