© The Authors 2025

1. Introduction

The emergence of galaxies in the first billion years of the Universe remains a key frontier in the field of galaxy formation. Representing the seeds of all subsequent structure formation, these early sources were responsible for producing and releasing the first heavy elements and dust grains into the interstellar medium (ISM; Maiolino & Mannucci 2019). Crucially, their stars and black holes produced the first hydrogen ionising photons, starting the Epoch of Reionisation (EoR), that marks the last major phase transition of all of the hydrogen in the Universe (for reviews, see e.g. Dayal & Ferrara 2018; Gnedin & Madau 2022). However, the emergence of these sources and their large-scale impacts remain key outstanding issues in the field of physical cosmology.

Our view of the first galaxies has evolved dramatically since the launch of the James Webb Space Telescope (JWST; see reviews by Robertson 2022; Adamo et al. 2024). Its incredible sensitivity has yielded unprecedented insights into galaxy properties well within the EoR, including the redshift evolution of the ultraviolet luminosity function (UV LF) at z ∼ 7 − 20 (e.g. Bouwens et al. 2023; Naidu et al. 2022a; Adams et al. 2024; Donnan et al. 2024; Harikane et al. 2022, 2023, 2025; Leung et al. 2023; McLeod et al. 2024; Pérez-González et al. 2023; Willott et al. 2024; Kokorev et al. 2024; Whitler et al. 2025), the stellar mass function (SMF) at z ∼ 5 − 12 (e.g. Navarro-Carrera et al. 2024; Weibel et al. 2024; Harvey et al. 2025) ISM metallicities at z ∼ 5 − 7 through emission lines (e.g. Curti et al. 2024; Chemerynska et al. 2024), UV spectral slopes at z ∼ 7 − 12 (e.g. Austin et al. 2024; Franco et al. 2024; Kokorev et al. 2024; Rinaldi et al. 2024), Balmer break strengths at z ∼ 6 − 10 (Vikaeus et al. 2024) and even hints on hydrogen ionising photon production efficiencies at z ∼ 5 − 9 (e.g. Simmonds et al. 2024; Llerena et al. 2024; Begley et al. 2025). Crucially, it has also revealed an overabundance of bright sources through the UV LF at z ≳ 11 (e.g. Casey et al. 2024; Donnan et al. 2024; Naidu et al. 2022a; Finkelstein et al. 2024; Robertson et al. 2024). Although a few photometrically selected high redshift galaxies were later found to be low-redshift interlopers with strong nebular emission lines (Naidu et al. 2022b; Arrabal Haro et al. 2023), numerous spectroscopic confirmations by NIRspec have solidified this overabundance of bright sources persisting out to z ∼ 14 (Carniani et al. 2024; Harikane et al. 2024).

The prevalence of such bright sources within the first billion years has led to a number of theoretical explanations, including (i) a decrease in dust attenuation with increasing redshift due to dust ejection in radiation-driven outflows (Ferrara et al. 2023; Fiore et al. 2023; Ziparo et al. 2023; Nakazato & Ferrara 2024), with dust being more dispersed in the host halo (Nikopoulos & Dayal 2024) or being spatially segregated from star-forming regions (Ziparo et al. 2023); (ii) a bias in flux-limited observations that leads to the preferential detection of galaxies with extreme star-bursts or stochastic star formation (Mason et al. 2023; Mirocha & Furlanetto 2023; Shen et al. 2023; Sun et al. 2023; Nikopoulos & Dayal 2024); (iii) accreting black holes contributing to the UV luminosity (Pacucci et al. 2022; Dayal et al. 2025; Hegde et al. 2024); (iv) largely inefficient stellar feedback, which could boost the star formation efficiency (Dekel et al. 2023; Li et al. 2024); and (v) a stellar initial mass function (IMF) that becomes increasingly top-heavy with increasing redshift, resulting in high light-to-mass ratios (Pacucci et al. 2022; Haslbauer et al. 2022; Trinca et al. 2024; Cueto et al. 2024; Hutter et al. 2025; Yung et al. 2024; Jeong et al. 2025; Lu et al. 2025). This last option is extremely promising, given a number of recent observations that seem to hint at a top-heavy IMF at early epochs through enhanced ionising photon production rates (Simmonds et al. 2024), extremely blue UV slopes (Topping et al. 2022; Cullen et al. 2024; Morales et al. 2024), and elevated [N/O] ratios (Bunker et al. 2023; Cameron et al. 2023; Isobe et al. 2023; Curti et al. 2025; Watanabe et al. 2024). Additionally, numerical simulations following the birth of star clusters at high redshifts show that elevated gas densities and low metallicities – conditions typical of the systems detected by JWST – favour top-heavy IMFs (Chon et al. 2022, 2024). Finally, the high values observed for the [O III]_88 μm to [C II]_158 μm ratio at z > 6 (Carniani et al. 2020; Harikane et al. 2020) have been used to infer a low metallicity stellar population with a top-heavy IMF (Katz et al. 2022).

In this work, we explore several of these mechanisms within the semi-analytical DELPHI framework aimed at tracking the assembly of dark matter halos and their baryons in the first billion years (Dayal et al. 2014, 2022; Mauerhofer & Dayal 2023). Compared to the latest version of DELPHI (Mauerhofer & Dayal 2023), referred to as DELPHI23 henceforth, we include the following key enhancements in this work. Firstly, we introduce a stochastic component into the star formation modelling by extracting ISM properties of galaxies from SPHINX²⁰, a state-of-the-art radiation hydrodynamics (RHD) simulation (Rosdahl et al. 2022). The high-resolution, accurate out-of-equilibrium ISM photochemistry and the turbulence-based star formation recipe of this large cosmological simulation provide us with distributions of star formation efficiencies and cold gas fractions as a function of redshift and halo mass. Stochasticity is automatically included in these distributions, since the simulated galaxies’ star formation properties are naturally bursty. Sampling SPHINX²⁰ distributions allows us to introduce random scatter in DELPHI, with some rare galaxies experiencing starbursts. Secondly, we implement a version of the model with an evolving IMF that becomes increasingly top-heavy with decreasing metallicities and increasing redshifts (e.g. Cueto et al. 2024). With these new prescriptions in hand, we aim to analyse the conditions required to reproduce the overabundance of bright galaxies at z > 9, in addition to studying the physical and observable properties of these galaxy populations.

Throughout this work, we assume a ΛCDM model with dark energy, dark matter, and baryonic densities in units of the critical density as Ω_Λ = 0.6889, Ω_m = 0.3097, and Ω_b = 0.04897, respectively. The Hubble constant we adopt is H₀ = 100 h km s⁻¹ Mpc⁻¹, with h = 0.677, a spectral index of n = 0.965, and a normalisation of σ₈ = 0.811 (Planck Collaboration VI 2020). All units in this work are comoving units, and magnitudes are in the standard AB system (Oke & Gunn 1983).

The paper is organised as follows. In Sect. 2, we describe the methodology of DELPHI, including the dark matter merger trees, the derivation of ISM properties from SPHINX²⁰, the modelling of star formation and its associated feedback, the metal and dust enrichment of early systems, and the emerging luminosity for models of varying IMF and star formation efficiency. We undertake a thorough comparison of the properties of these systems against observables – including the redshift evolution of the UV LF and the SMF, dust and metal masses, UV spectral slopes, Balmer break strengths, and the ionising photon production efficiency – in Sect. 3. We then analyse the reionisation history and its key sources for a number of escape fraction scenarios in Sect. 4. We end by summarising our key findings, and discussing caveats and possible improvements in Sect. 5.

2. The theoretical model

The new approach we adopt here combines the DELPHI semi-analytic galaxy formation model with insights on the ISM – in terms of the cold gas fraction and stochastic star formation – from the (very high-resolution) SPHINX²⁰ radiation hydrodynamics simulation (Rosdahl et al. 2022). The models are briefly described here, and readers are referred to previous works describing DELPHI (Dayal et al. 2014, 2022; Mauerhofer & Dayal 2023) and SPHINX²⁰ (Teyssier 2002; Rosdahl et al. 2013; Rosdahl & Teyssier 2015) for complete details. Finally, we use a Kroupa IMF (Kroupa 2001) throughout this work, with a slope of −1.3 (−2.3) between 0.1–0.5 (0.5–100) M_⊙, except for in the evolving IMF model (Sect. 2.6.1).

2.1. Merger trees and gas assembly

We start by building dark matter merger trees, that form the basis of DELPHI, for 600 halos at z = 4.5 with (logarithmically equally spaced) masses between $M_{\mathrm{h}}={10}^{8-14}{\mathrm{M}}_{\odot }$ using a binary merger tree algorithm (Parkinson et al. 2008). In brief, the merger tree for each simulated halo starts at a redshift of z = 4.5 and runs backwards in time up to z = 40.8, with each halo fragmenting into its progenitors. A halo of mass M₀ can either fragment into two halos with masses larger than the resolution mass M_res < M < M₀/2, or lose a fraction of its mass, meaning fragment into progenitor halo(s) below the mass resolution. In the latter case, the mass brought in by unresolved progenitors is denoted as ‘smooth accretion’ from the intergalactic medium (IGM). The merger tree used in this work stores outputs over 44 time-steps equally spaced by 30 Myr with a resolution mass of $M_{\mathrm{res}}={10}^{8.0}{\mathrm{M}}_{\odot }$. Further, each halo at z = 4.5 is assigned a number density by matching to the Sheth-Tormen (Sheth & Tormen 1999) halo mass function (HMF). All the progenitors of a halo are assigned the same number density as their z = 4.5 descendant; we have confirmed the resulting HMFs are in accordance with the Sheth-Tormen HMFs at all higher redshifts.

In terms of baryons, the first halos (‘starting leaves’) of any merger tree are assigned an initial gas mass that is linked to the halo mass through the cosmological baryon-to-dark matter ratio such that $M_{\mathrm{g}}^{\mathrm{i}}=({\mathrm{\Omega }}_b/{\mathrm{\Omega }}_m)M_{\mathrm{h}}$. If any halo has progenitors, its initial gas mass is the sum of the final gas mass brought in by merging progenitors – after star formation and the resulting supernova (SN) feedback (Sect. 2.3) – and the gas mass smoothly accreted from the IGM. For the latter, we reasonably assume smoothly accreted dark matter to drag in a cosmological ratio of baryons.

2.2. Cold gas fractions and star formation efficiencies from the SPHINX²⁰ simulation

In this work we incorporate the cold gas fraction and star formation efficiency (SFE) distributions sampled from the SPHINX²⁰ simulation¹ into the DELPHI semi-analytic model. In brief, SPHINX²⁰ simulates a cubic volume of (20 cMpc)³ using the RAMSES-RT adaptive mesh refinement (AMR) code. With a dark matter particle resolution mass of 2.5 × 10⁵ M_⊙, this simulation resolves halos at the atomic cooling mass of 3 × 10⁷ M_⊙ with 120 dark matter particles. The volume simulated here was selected as the most average amongst 60 cosmological initial conditions in order to minimise cosmic variance biases. The resolution reaches 10 pc in the densest regions of the ISM, sufficient to resolve gas fragmentation into multiple star-forming regions, the escape of ionising radiation and the photochemistry of interstellar gas.

In terms of gravity and hydrodynamics, the collisionless dark matter particles and stellar particles are evolved using a particle-mesh solver with cloud-in-cell interpolation (Guillet & Teyssier 2011) and the gas evolution is computed by solving the Euler equations with a second-order Godunov scheme using the HLLC (Harten–Lax–van Leer contact) Riemann solver (Toro et al. 1994). Gas is converted into stars in cells that have a gas density higher than 10 cm⁻³, a locally turbulent Jeans length that is smaller than the cell width and gas that is locally convergent (Rosdahl et al. 2018). When a cell satisfies those conditions, its star formation efficiency per free-fall time is computed as a function of the local virial parameter and turbulence (Federrath & Klessen 2012), and its gas converted stochastically into stellar particles (Rasera & Teyssier 2006) with an initial mass of 400 M_⊙. The feedback implementation of SPHINX²⁰, which is crucial for shaping the cold gas fraction and star formation histories, follows the prescription detailed in Rosdahl et al. (2018) where each SN explosion injects 10⁵¹ erg of energy ‘mechanically’ into the ISM (Kimm et al. 2015). Individual stellar particles undergo several SN explosions from ages 3–50 Myr, with an average of four SN explosions per 100 M_⊙ formed. This SN rate is artificially boosted by a factor of roughly four compared to a Kroupa IMF, which is necessary to suppress star formation and produce realistic high redshift luminosity functions (Rosdahl et al. 2022). The code includes radiative transfer (RT) in two bins of hydrogen-ionising photons using a first-order moment method. RT is an essential ingredient to obtain an accurate non-equilibrium ionisation state and the temperature of hydrogen, crucial for the cold gas fraction considered in this work. The ionising photons are emitted from stellar particles within the simulated galaxies using the BPASSv2.2.1 (Eldridge et al. 2008; Stanway et al. 2016) spectral energy distribution (SED) model.

The most massive halo in SPHINX²⁰ has a virial mass of 4.3 × 10¹¹ M_⊙, 9.5 × 10¹⁰ M_⊙, 2.8 × 10¹⁰ M_⊙ and 4.0 × 10⁹ M_⊙ at z ∼ 5, 7, 9 and 12, respectively. In this paper, we use all the simulated halos from SPHINX²⁰ more massive than 10⁸ M_⊙, corresponding to the lowest halo mass in the merger tree for DELPHI, in eleven snapshots from redshift 5–15. This results in a total of 48 801, 31 338, 15 463 and 3481 SPHINX²⁰ sources at z ∼ 5, 7, 9 and 12, respectively. We note that at any redshift, as expected from hierarchical structure formation, the HMFs from both the DELPHI merger trees and SPHINX²⁰ follow the same shapes, as shown at z ∼ 5 and 9 in Fig. 1. However, given its analytic nature, at any redshift, the DELPHI merger trees contain halos that are up to two orders of magnitude more massive than the most massive SPHINX²⁰ halos.

Fig. 1. Halo mass functions obtained from the SPHINX20 and DELPHI models, shown using blue and green histograms, respectively, as marked. The left and right panels show results at z ∼ 5 and 9, respectively.

In this new version of DELPHI, the first major change is that we split the total gas mass into warm and cold components. We define cold gas as that with a temperature < 1000 K – this forms both the reservoir for star formation and grain growth of dust grains in the ISM (Sect. 2.4). The temperature threshold of 1000 K is chosen to demarcate gas that is cold enough to form stars in some SPHINX²⁰ cells, but not cold enough to enter the regime of very dense molecular clouds that cannot be resolved by the simulation. We then assign a cold gas mass fraction, f_cold, to each halo in DELPHI by sampling halo properties from the SPHINX²⁰ simulation, which introduces a stochastic component in our model. To do so, we sample the probability distribution function (PDFs) for the cold gas fraction in SPHINX²⁰ in stellar mass bins separated by 0.1 dex ranging from 10⁸ M_⊙ up to the largest mass at that redshift. The bins partially overlap so as to avoid discontinuities between neighbouring distributions. The PDFs of this cold gas fraction are plotted in Fig. 2, for a subset of halo mass bins and two redshifts, z ∼ 5 and 9. As seen from this figure, at z ∼ 5, about 70% (58%) of the low-mass halos with M_h ∼ 10⁸ (10^8.4) M_⊙ have no cold gas. This is, in part, because they reside in ionised regions where the UV background from reionisation prevents gas cooling and/or accretion of cold gas (see Katz et al. 2020), and, in part due to internal feedback processes. Only 10 − 20% of such low-mass halos have values of f_cold ∼ 0.2 − 0.4. An increasing fraction of halos start hosting high cold gas fractions with increasing halo mass at z ∼ 5. For example, roughly 50% of the halos with M_h ≳ 10^10.4 M_⊙ show f_cold ∼ 0.2 − 0.45. While roughly the same trends persist at z ∼ 9, a key difference is that low-mass halos (M_h ≲ 10^8.4 M_⊙) contain larger cold gas fractions than at lower redshifts. For example, only 35% of such halos do not contain any cold gas, with about 65% showing values of f_cold ∼ 0.05 − 0.3. This can be explained by the diminution of the UV background compared to that at z ∼ 5. The behaviour at larger masses is roughly the same across these two redshifts, with about 50% of M_h ≳ 10¹⁰ M_⊙ halos showing values of f_cold ∼ 0.2 − 0.4.

Fig. 2. Fraction of halos with a cold gas fraction (fcold) below the value shown on the horizontal axis from the SPHINX20 simulations, for the mass bins marked, at z = 5 (left panel) and z = 9 (right panel). The cold gas fraction is defined as the fraction of gas mass that is below 1000 K, within the virial radius of the halo.

The second key change in our new version of DELPHI is the inclusion of a stochastic star formation parameter, f_*^sphinx. This is the ratio of the stellar mass formed in the last 30 Myr (M_∗[< 30 Myr]) divided by the current mass of cold (< 1000 K) gas; the time-step of 30 Myr is chosen to be coherent with the time-stepping of the underlying DELPHI merger tree. If the cold gas mass of a halo is zero, we define f_*^sphinx = 0, and we cap f_*^sphinx to a maximum of 1. To summarise,

$$\begin{array}{c}\hfill f_{\ast }^{\mathrm{sphinx}}=min(1,\frac{M_{\ast }[<30\mathrm{Myr}]}{M_{\mathrm{g}}^{\mathrm{i}}{f}_{\mathrm{cold}}}).\end{array}$$(1)

The resulting efficiencies, for a number of stellar mass bins, are displayed in Fig. 3. At z ∼ 5, about 90% of M_h ≲ 10^8.4 M_⊙ halos show f_*^sphinx = 0, meaning that either their cold gas fraction is 0 or that they did not form stars in the last 30 Myr. Again, the star formation efficiency increases with an increase in the host halo mass with roughly 50% of M_h ≳ 10^10.4 M_⊙ showing f_*^sphinx ∼  0.1 − 0.4. At z ∼ 9, a larger fraction (∼40%) of low-mass halos with M_h ≲ 10^8.4 M_⊙ show non-zero values of f_*^sphinx, driven by their larger cold gas fractions discussed above. In accordance with z ∼ 5 results, about 50% of the most massive halos with M_h ≳ 10¹⁰ M_⊙ at this redshift show f_*^sphinx ∼ 0.08 − 0.4.

Fig. 3. Fraction of halos with a star formation efficiency value (f*sphinx) below the value shown on the horizontal axis from the SPHINX20 simulations, for the mass bins marked, at z = 5 (left panel) and z = 9 (right panel). The star formation efficiency is defined as the ratio between the mass of stars younger than 30 Myr and the cold gas mass; if a halo has no cold gas mass, we define its star formation efficiency to be f*sphinx = 0.

To assign f_cold and f_*^sphinx values to DELPHI halos, we choose the SPHINX probability distribution at the closest redshift and halo mass bin using a standard Monte Carlo sampling. As shown in Fig. 1 above, DELPHI halos reach higher masses than those in the SPHINX²⁰ simulations given the latter’s relatively small box size; as a result, the most massive DELPHI halos have to be matched with smaller SPHINX²⁰ halos. However, both the f_cold and f_*^sphinx distributions start showing convergence of the PDF for the few most massive SPHINX²⁰ halos. Thus, we expect that matching DELPHI halos with lower mass systems from SPHINX²⁰ does not induce a large error in our results.

2.3. Star formation rates and supernova feedback

The newly formed stellar mass at any time-step can be written as

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\ast }(z)=f_{\ast }^{\mathrm{eff}}\times M_{\mathrm{g}}^{\mathrm{i}},\end{array}$$(2)

where the ‘effective’ star formation efficiency is f_*^eff = min(f_cold ⋅ f_*^sphinx, f_*^ej). Here, $M_{\mathrm{g}}^{\mathrm{i}}$ is the initial gas mass in the halo at the start of the time-step. Further, f_*^ej is the star formation efficiency that produces enough Type II Supernova (SNII) feedback to eject the rest of the gas (see Mauerhofer & Dayal 2023):

$$\begin{array}{c}\hfill f_{\ast }^{\mathrm{ej}}=\frac{v_c^2}{v_c^2+f_w\nu E_{51}},\end{array}$$(3)

where v_c is the halo rotational velocity, ν is the supernova rate, E₅₁ = 10⁵¹ erg is the energy per SNII, and f_w is a free parameter that determines the fraction of SN energy coupling to the gas content. If f_*^ej < f_cold ⋅ f_*^sphinx, a galaxy is ‘feedback limited’; that is, it loses all its gas due to SNII feedback. In this paper, we fix f_w = 4%, in order to match the observed redshift evolution of the UV LF at z ≲ 10, as also detailed in Table 1.

2.4. Dust and metal enrichment

The evolution of the metal (M_Z) and dust contents (M_d) of early galaxies are inextricably coupled through the processes of metal/dust production, astration into further star formation, ejection in outflows, dilution in (metal and dust poor) inflows of gas from the IGM, and dust destruction and growth in the ISM. We assume metals, dust and gas to be perfectly mixed in the ISM. In terms of metal production, we assume all SN and AGB (asymptotic giant branch) stars between 2 − 50 M_⊙ to contribute to the metal yield using the results presented in Kobayashi et al. (2020). The lower limit to this mass is driven by the fact that that stars less than 2 M_⊙ need ∼1.7 Gyr to end their lives, which is larger than the age of the Universe at z > 4.5; stars more massive than 50 M_⊙ are assumed to collapse to black holes without any yield. We consider SNII as the key dust factories, with each SN producing 0.5 M_⊙ of dust (e.g. Leśniewska & Michałowski 2019; Dayal et al. 2022). The coupled evolution of dust and gas-phase metals during a DELPHI time-step is modelled as (Dayal et al. 2022)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{M}_Z}{\mathrm{d}t}& ={\dot{M}}_Z^{\mathrm{pro}}-{\dot{M}}_Z^{\mathrm{eje}}-{\dot{M}}_Z^{\mathrm{ast}}-{\dot{M}}_d^{\mathrm{gro}}+{\dot{M}}_d^{\mathrm{des}}\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{M}_d}{\mathrm{d}t}& ={\dot{M}}_d^{\mathrm{pro}}-{\dot{M}}_d^{\mathrm{eje}}-{\dot{M}}_d^{\mathrm{ast}}+{\dot{M}}_d^{\mathrm{gro}}-{\dot{M}}_d^{\mathrm{des}}.\hfill \end{array}$$(5)

Eqs. (4) and (5) show the time evolution of the metal and dust masses, respectively, considering the processes of production (pro), ejection (eje), astration (ast), the loss of metals into dust grains growing in the cold ISM (gro), and the destruction of dust grains in the warm ISM that add to the metal content (des). We use the newly introduced cold gas fraction to calculate the destruction and growth of dust grains in the warm and cold ISM, respectively; this is a key difference with respect to the DELPHI23 model where we had assumed an equi-partition of the gas into the cold and warm phases.

2.5. Ultraviolet and ionising luminosities

In order to calculate the intrinsic SED and the hydrogen ionising photon production rate for each galaxy, we convolve its star formation history (assuming the stellar mass in any 30 Myr time-step to have formed through constant star formation) with the Starburst99 (SB99) stellar population synthesis code (Leitherer et al. 1999) using the ‘Geneva High mass loss’ tracks (Meynet et al. 1994).

To include the effect of dust, we start by calculating the dust attenuation of UV photons (1450–1550 Å in the rest-frame) as in Mauerhofer & Dayal (2023): we use a single grain size of a = 0.05 μm and a material density of s = 2.25 g cm⁻³ appropriate for graphite/carbonaceous grains (Todini & Ferrara 2001; Nozawa et al. 2003). We then assume this dust to be mixed with gas in a radius given by

$$\begin{array}{c}\hfill r_{\mathrm{gas}}=4.5\lambda r_{\mathrm{vir}}\times \left(\frac{1+z}{6}\right),\end{array}$$(6)

where λ = 0.04 is the spin parameter (e.g. Bullock et al. 2001). As seen, both gas and dust cover an increasing volume fraction of the halo with increasing redshift, as indicated by recent ALMA observations (Fujimoto et al. 2020; Fudamoto et al. 2022). This is used to obtain the dust optical depth to UV photons as τ_c = 3M_d[4πr_gas² as]⁻¹ yielding the observed UV magnitude. We then use the Calzetti extinction curve (Calzetti et al. 2000) convolved with this UV optical depth to account for the dust attenuation of the entire SED. This dust attenuated SED is then used to calculate β slopes (Sect. 3.6) and the Balmer break strength (Sect. 3.7) in order to compare our models to observations.

2.6. Two extreme star formation models

Besides the fiducial model described above, we also explore two alternative models aimed at increasing the UV luminosities of galaxies at very high redshifts. This is motivated by the overabundance of bright sources detected by JWST at z ≳ 11, as discussed in what follows.

2.6.1. A redshift- and metallicity-dependent initial mass function

One way of boosting the UV luminosities of galaxies at z ≳ 10 is by appealing to an evolving IMF (the ‘eIMF’ model) that becomes increasingly top-heavy with redshift or decreasing metallicity as supported by high-resolution star formation simulations (e.g. Chon et al. 2022). Our modified evolving IMF has the following form, inspired by previous results (Chon et al. 2022; Cueto et al. 2024):

$$\begin{array}{c}\hfill \xi (M)\propto \{\begin{array}{cc}(1-f_{\mathrm{massive}})M^{-1.3},\hfill & 0.1{\mathrm{M}}_{\odot }\le M<0.5{\mathrm{M}}_{\odot }\hfill \\ (1-f_{\mathrm{massive}})M^{-2.3},\hfill & 0.5{\mathrm{M}}_{\odot }\le M<M_{\mathrm{c}}\hfill \\ f_{\mathrm{massive}}{M}^{-1},\hfill & M_{\mathrm{c}}\le M\le 100{\mathrm{M}}_{\odot }.\hfill \end{array}\end{array}$$(7)

Here, M_c is the cut-off mass at which the Kroupa IMF transitions to a log-flat component and f_massive is the fraction of stellar mass above M_c. The value of f_massive is given by (Chon et al. 2022):

$$\begin{array}{c}\hfill f_{\mathrm{massive}}=1.07·(1-2^x)+0.04·2.{67}^x·z,\end{array}$$(8)

where z is the redshift and x = 1 + log₁₀(Z/Z_⊙). Here, Z is the metallicity before star formation and Z_⊙ = 0.0134 represents the solar metallicity value (Asplund et al. 2009). M_c is computed by solving the system of equations formed by the two continuity equations of the IMF at 0.5 M_⊙ and at M_c, the normalisation equation requiring ${\displaystyle {\int }_{0.1}^{100}}\xi (m)m\mathrm{d}m=1$, and the definition of f_massive = ∫_Mc¹⁰⁰ξ(m) m dm. We cap f_massive below 0.995 to ensure that M_c is always above 0.5 M_⊙.

The first impact of the eIMF model is on the SNII rate (ν), that impacts both feedback and metal/dust production. As shown in Fig. 4, the SNII rate has a constant value of $\nu =0.011{{\mathrm{M}}_{\odot }}^{-1}$ for the fiducial Kroupa IMF. In the eIMF model, however, the IMF is a function of both redshift and metallicity, which leads to a corresponding change in ν. As shown, the SNII rate increases by about 2.3 times, to $\nu \sim 0.025{{\mathrm{M}}_{\odot }}^{-1}$ for Z ∼ 10⁻⁴ Z_⊙, at all z ∼ 5 − 14. The value of ν decreases with increasing metallicity as expected – at z ∼ 5 (9), the SNII rates in the eIMF and Kroupa IMFs converge for Z ≳ 0.1 (0.25) Z_⊙. However, the SNII rates saturate at a value that is about 1.25 times higher ($\nu \sim 0.014{{\mathrm{M}}_{\odot }}^{-1}$) at z ∼ 14, even for Z ∼ Z_⊙.

Fig. 4. Evolution of the SN rate (ν) as a function of metallicity, at the redshifts marked, for the eIMF model (Sect. 2.6.1). The shape of the IMF in this case is determined by the metallicity and redshift. The purple line shows the supernova rate for the fiducial Kroupa IMF.

We then calculate the impact of this evolving IMF on stellar (dust and gas-phase metals) yields and the total yields (dust, metals and returned gas) using the grids from Kobayashi et al. (2020), the results of which are shown in Fig. 5. For a Kroupa IMF, we find constant dust and metal (total) yield values ∼0.017 (0.147). In the eIMF case, the dust and metals yields at the lowest metallicities are about 0.05 (0.34), independent of redshift. These values decrease with increasing metallicity all at z; however, at any metallicity value, the yields are higher at higher redshifts due to a more top-heavy IMF (see Eq. (8)). At z ∼ 5 (9) the dust and metal yields converge to the Kroupa IMF values at Z ∼ 0.1 (0.25) Z_⊙. At z ∼ 14, the dust and metal yields saturate at values of about 0.028 (0.184) even for Z = Z_⊙.

Fig. 5. Evolution of the stellar metal yields (solid lines) and total yields (metals and gas return fraction – dashed lines) per unit stellar mass as a function of metallicity, for the redshifts marked, for the eIMF model. The dashed purple line shows the metallicity-dependent yield values for a Kroupa IMF used in our other models.

Finally, as expected, the eIMF model also impacts the SEDs for early systems. We compute spectra with SB99 for a list values of M_c distributed between 0.5 and 100 M_⊙ and four values of metallicities between 0.001 and 0.02. For each halo, we interpolate its spectrum using the closest values in our grid of M_c and metallicities. We compute hydrogen-ionising photon production rates, plus luminosities at nine wavelengths between 1300 Å and 2900 Å, and at 3500 Å and 4200 Å, to be able to compute beta slopes and Balmer break strengths (see Sections 3.6 and 3.7). For illustration we show the impact of the eIMF model on the production rate of ionising photons $\left({\dot{N}}_{\mathrm{ion}}^{\mathrm{sf}}\right)$ and on the intrinsic UV luminosity $\left(L_{\mathrm{UV}}^{\mathrm{int}}\right)$ in Fig. 6. At all redshifts, ${\dot{N}}_{\mathrm{ion}}^{\mathrm{sf}}$ $\left(L_{\mathrm{UV}}^{\mathrm{int}}\right)$ for the eIMF model is almost a factor 8 (4) larger than for a Kroupa IMF at Z ∼ 10⁻⁴ Z_⊙. The values slowly decrease with increasing metallicity, until they become equal to that from the Kroupa IMF at Z ∼ 0.15 Z_⊙ (0.5 Z_⊙) at z ∼ 5 (9). For redshifts as high as z ∼ 14, ${\dot{N}}_{\mathrm{ion}}^{\mathrm{sf}}$ $\left(L_{\mathrm{UV}}^{\mathrm{int}}\right)$ saturate at values about 4 (2) times higher than those from the Kroupa IMF, even at solar metallicity.

Fig. 6. As a function of metallicity, we show the ratio between the ionising photon production rates (intrinsic UV luminosities) in the eIMF and fiducial model using solid (dashed) lines at the marked redshifts.

2.6.2. An evolving star formation efficiency

The second model we explore is termed the evolving star formation efficiency model (‘SFE’). Here, both the cold gas fractions and star formation efficiencies evolve as a function of halo mass and redshift. We use f_*(M_h, z) = f_*^sphinx for halos with masses in the range of SPHINX²⁰ halos at the corresponding redshift, as well as for all halos below redshift 7. For more massive halos, at z ≥ 14, we have f_*(M_h, z) = 0.8, and in between those redshifts we do a weighted average: f_*(M_h, z) = (14 − z)/7 ⋅ f_*^sphinx + (z − 7)/7 ⋅ 0.8. The same holds for f_cold(M_h, z), but with a value of 0.5 at z ≥ 14.

We summarise the properties of our fiducial model, our two ‘extreme’ models, as well as the model from Mauerhofer & Dayal (2023) (labelled ‘DELPHI23’) in Table 1. We note that the DELPHI23 model did not model cold gas explicitly, but assumed a fraction of 50%, as mentioned in Sect. 2.4. In the fiducial model presented here, on the other hand, both the cold gas fractions and star formation efficiencies are derived from the SPHINX²⁰ simulations. The only free parameters in the new model are the fraction of SNII energy coupling to gas (f_w) and the gas distribution radius (r_gas) that are chosen so as to simultaneously match to the observed evolution of the UV LF and stellar mass function at z ∼ 5 − 9, and dust observables at z ∼ 5 − 7 from ALMA observations (Béthermin et al. 2020; Bouwens et al. 2022a).

3. Confronting models and observations in the first billion years

We now discuss the calibration of the model parameters against the latest UV LF estimates obtained from the JWST for all of the models considered in this work. We then compare our models against a number of early galaxy observables including the UV luminosity density, the evolving stellar mass function, the dust masses, the mass-metallicity relation, the beta slopes, Balmer break strengths and, the ionising photon production efficiency. In what follows, we compare the results of the new model against those from the last version of the code referred to as DELPHI23 (Mauerhofer & Dayal 2023).

3.1. The redshift evolution of the ultraviolet luminosity function at z ∼5 − 20

We now discuss the redshift evolution of both the intrinsic and observed (dust attenuated) UV LFs. In addition to comparing these to available observations at z ∼ 5 − 14, we show model predictions out to z ∼ 20. Starting at z ∼ 5, we briefly recap that the intrinsic UV LF from the DELPHI23 model over-predicts the bright end of the observed UV LF for M_UV ≲ −21 (−22) at z ∼ 5 (7) as shown in Fig. 7. Adding dust attenuation results in observed UV LFs in accordance with observations at these redshifts. The importance of dust attenuation decreases with increasing redshift in all of the models considered here. This is driven by gas and dust (that are assumed to be perfectly mixed) occupying an increasing volume of the virial radius, that is getting more ‘diffuse’ with redshift, leading to a decrease in the dust optical depth. While dust still plays a (decreasing) role in attenuating the UV luminosity out to z ∼ 12, the intrinsic and observed UV LFs converge at all z ≳ 14, meaning that dust attenuation plays no role at these early epochs. In this model, even the intrinsic UV LF under-predicts the bright end of the observed UV LF, with M_UV ≲ −21, at z ∼ 12 − 14.

Fig. 7. Redshift evolution of the UV LF between z ∼ 5 − 20, as marked. In each panel, we show a comparison of the fiducial model, and the eIMF and eSFE models against the results from DELPHI23 (Mauerhofer & Dayal 2023) for both the intrinsic (dashed lines) and dust attenuated (solid lines) UV LFs. These lines show the means of five runs sampling the PDFs for fcold and f* from SPHINX20. For the sake of clarity we only show the maximum and minimum values associated with these runs for the observed UV LFs through shaded areas. In each panel, as marked, points show observational data (from Adams et al. 2024; Atek et al. 2015, 2018; Bouwens et al. 2021, 2022b, 2023; Bowler et al. 2017; Casey et al. 2024; Castellano et al. 2023; Donnan et al. 2023, 2024; Finkelstein et al. 2024; Harikane et al. 2022, 2023, 2024, 2025; Ishigaki et al. 2018; Leung et al. 2023; McLeod et al. 2024; Naidu et al. 2022a; Napolitano et al. 2025; Pérez-González et al. 2023; Robertson et al. 2024; Willott et al. 2024; Whitler et al. 2025). The last panel shows the UV LF integrated between z ∼ 15 − 20 in order to compare with observations at those redshifts (Kokorev et al. 2024; Whitler et al. 2025).

In the fiducial model presented in this work, as noted, the value of f_*^eff (composed of the cold gas fraction and star formation efficiency) are directly obtained from the SPHINX²⁰ simulations. Compared to the DELPHI23 model, we find f_*^eff values that are a factor 5 (10) lower for M_h ∼ 10¹⁰ M_⊙ halos at z ∼ 5 (17). These lower star formation rates are compensated by invoking a slightly weaker coupling between SNII energy and gas (f_w = 0.04) as compared to the 1.5 times higher value (f_w = 0.06) used in the DELPHI23 model. As a result, in this new fiducial model, galaxies lose less gas at any time-step, both in terms of star formation and outflows. While the lower star formation rates result in the intrinsic UV LF from the fiducial model being underestimated slightly (by about 0.3–0.5 dex for a given UV magnitude) at the bright end (M_UV ≲ −20) at z ≳ 8 compared to that from DELPHI23, self-regulation of star formation and feedback results in the intrinsic UV LFs from both these models converging by z ∼ 5. In this fiducial model, the intrinsic UV LF overestimates the bright end (M_UV ≲ −21) of the observed UV LF at z ≲ 8. Reconciling these requires dust attenuation, the impact of which again decreases (due to a decrease in the dust optical depth) with increasing redshift. Although the bright end of the intrinsic UV LF shows scatter due to the five different random seeds used for f_* and f_cold, the mean UV LF under-predicts the bright end of the UV LF for M_UV ≲ −20.5 at z ∼ 11 − 12, and as faint as M_UV ∼ −19.5 by z ∼ 14: for example, it falls short of observations (Casey et al. 2024; Donnan et al. 2024; Naidu et al. 2022a; Finkelstein et al. 2024; Robertson et al. 2024) by ≳1 dex at z ∼ 11 − 14.

We then appeal to the eIMF model that produces up to four times higher luminosity per unit star formation as compared to the fiducial model. As shown in the same figure, this model yields an intrinsic UV LF that is in very close agreement with that from DELPHI23, that is it overestimates the observations at all z ∼ 5 − 14. Dust is required in order to reconcile this model with observations for M_UV ≲ −21 at z ∼ 5 − 10; by z ∼ 11 the intrinsic and observed UV LFs converge due to a decreasing dust attenuation. Interestingly, this model is able to reproduce all JWST observations, including the bright end of the evolving UV LF between z ∼ 5 − 20 (Casey et al. 2024; Donnan et al. 2024; Finkelstein et al. 2024; Harikane et al. 2025; Robertson et al. 2024; Kokorev et al. 2024; Whitler et al. 2025). The other model designed to boost the UV luminosity of DELPHI halos is the ‘eSFE’ model, where the cold gas fractions and star formation efficiencies of massive halos (above the largest SPHINX²⁰ halo masses) slowly increase with respect to the fiducial model at z > 7. This model also yields an intrinsic UV LF that over-predicts the data at z ∼ 5 − 10 and, including dust, this model successfully reproduces JWST observations of the UV LF, within error bars, at all redshifts z ∼ 5 − 20. From z ∼ 10 to 14, we also see the transition between the mass regime in which f_* is sampled from SPHINX²⁰ and the regime in which it has a redshift-dependent value (see Table 1). This creates a moderate flattening of the UV LF, yielding a slight discrepancy compared to observations at z ∼ 11 − 12 at M_UV = −19. Given the increasing value of f_*^eff, the eSFE model predicts a larger number density of galaxies as compared to the eIMF model for M_UV ∼ −18 to −21 at z ∼ 14 − 17; the eIMF model provides the upper limit to the observed UV LF at all other magnitudes for z ∼ 5 − 21. Indeed, at the highest redshifts of z ∼ 15 − 20, the high SFE values in this model result in the best agreement with observations.

Quantitatively, using the Calzetti extinction curve, we find similar values of dust attenuation E(B–V) of massive galaxies for each model: at z ∼ 5, galaxies above a stellar mass of 10¹⁰ M_⊙ have an average and maximum value of E(B–V) of 0.19 and 0.36, respectively. As expected, given the decreasing impact of dust with increasing redshift, by z ∼ 10, galaxies above a stellar mass of 10⁹ M_⊙ have an average and maximum E(B–V) value of 0.056 and 0.11.

Finally, we note that while the faint ends of the UV LF are relatively similar for all the models studied here, we see a large spread in predictions at the bright end that can be tested with forthcoming JWST observations.

3.2. The redshift evolution of the UV luminosity density at z ∼5 − 21

We now compare the UV luminosity density (ρ_UV) inferred from our models against observationally inferred results in Fig. 8. We start by discussing the dashed lines, that represent the intrinsic UV luminosity density integrating over all galaxies in our models. The fiducial model shows ρ_UV [erg s⁻¹] Hz⁻¹ cMpc⁻³] values that decreases by about four orders of magnitude from ∼10^26.7 at z ∼ 5 to ∼10^22.75 by z ∼ 20. As might be expected from the discussion in Sect. 3.1 above, while the DELPHI23 model shows ρ_UV values higher than the mean fiducial model by about 0.2–0.5 dex at z ∼ 12 − 20, these converge at lower redshifts. Producing the largest amount of light per unit star formation, the eIMF model forms the upper limit to the UV luminosity density at all redshifts with ρ_UV ∼ 10^26.8 (10^23.1) [erg s⁻¹ Hz⁻¹ cMpc⁻³] at z ∼ 5 (20); the eSFE model yields ρ_UV values that lie between the bounds provided by the fiducial and eIMF models. Interestingly, despite the different star formation prescriptions used, all of these models yield ρ_UV values that are different by less than 0.5 dex at any z ∼ 5 − 21.

Fig. 8. Redshift evolution of the UV luminosity density at z ∼ 5 − 21. Solid lines show the mean of five runs for each model, using dust attenuated values of ρUV integrating down to observed magnitude limits of MUV ≲ −17. Shaded areas show the maximum and minimum values associated with these runs. Dotted lines show values of the intrinsic luminosity density integrating over all galaxies. As marked, points show observational results (Donnan et al. 2023; McLeod et al. 2024; Finkelstein et al. 2024; Whitler et al. 2025) and the faint pink shaded area corresponds to the 16% and 84% marginal constraints from the JADES origin field (Robertson et al. 2024).

We also compute the ρ_UV including the impact of dust attenuation and integrating only down to observed limits of M_UV ∼ −17, the results of which are shown in the same figure as solid lines. The fiducial model shows ρ_UV values that decrease from 10^26.3 at z ∼ 5 to 10²¹ by z ∼ 20, meaning that galaxies with M_UV ≲ −17 only contribute about 3.5% of the UV luminosity at z ∼ 20, increasing to about 40% by z ∼ 5. Again, the DELPHI23 model yields ρ_UV values that are larger by 0.2–0.5 dex at z ∼ 12 − 20. With its higher production of UV photons, the eIMF model yield upper limits to the UV luminosity density with ρ_UV ∼ 10^26.4 (10²²) at z ∼ 5 (20); while the eSFE model lies very close to the eIMF model at z ≳ 13, at lower redshifts it transitions to the fiducial model as might be expected.

Given the large errors associated with observations (especially Robertson et al. 2024), as of now, all of these models are in reasonable agreement with existing observations at z ∼ 8 − 14. The eSFE and fiducial models are in better agreement with the values inferred at z ∼ 8 − 9 (from Donnan et al. 2023; Finkelstein et al. 2024) and all of the models are in agreement with the data at z ∼ 10 − 11 (Donnan et al. 2023; McLeod et al. 2024; Finkelstein et al. 2024; Whitler et al. 2025), although the eIMF sits slightly above the data. At z ∼ 13 − 15, where the data points from different works show values separated by about 1 dex, the fiducial model is in agreement with the results from Donnan et al. (2023) while the eIMF and eSFE models are in better agreement with the higher values inferred by Finkelstein et al. (2024) and Whitler et al. (2025). As shown, the model differences increase with increasing redshift, with the eIMF (and eSFE) models predicting ρ_UV values an order of magnitude above that from the fiducial model – the slope evolution of the ρ_UV − z relation will therefore be an excellent observable to distinguish these models with forthcoming JWST observations.

3.3. The redshift evolution of the stellar mass function from z ∼5 − 20

We also compare the SMF predicted by our models with those constructed from observations with for example the Hubble Space Telescope (HST) and the JWST, as shown in Fig. 9. At z ∼ 5 − 6, all four of the models studies in this work show SMFs that essentially overlap between M_∗ ∼ 10⁶⁻¹ M_⊙; the DELPHI23 model shows a slightly enhanced high-end tail driven by its large star formation efficiency. As we move to z ≳ 7, these models start showing increasingly larger differences: with its top-heavy IMF, the eIMF model has the largest light-to-mass ratio in addition to the largest SNII fraction per unit mass. As a result, this sets the lower limit to the theoretical SMF evolution. As seen, this model shows a 1.5 (2.5) times lower amplitude as compared to the fiducial model at z ∼ 9 (17). On the other hand, the eSFE model shows the opposite trend. Thanks to its larger effective star formation efficiencies in massive halos above redshift 7, it shows a much larger number of galaxies per unit volume. For example, it shows a factor 10 (100) more galaxies per unit volume at M_* ∼ 10⁹ M_⊙ at z ∼ 9 (14) as compared to the fiducial model.

Fig. 9. Redshift evolution of the SMF at z ∼ 5 − 20, for the different model explored in this work. Solid lines show the mean results from five runs of each model with shaded areas showing the associated maximum and minimum values. The bottom right panel shows the SMF integrated on all outputs between redshifts 15 and 20. As marked, each panel up to z ∼ 12 show observational data, all assuming a Kroupa IMF (from Bhatawdekar et al. 2019; Duncan et al. 2014; Kikuchihara et al. 2020; Navarro-Carrera et al. 2024; Song et al. 2016; Stefanon et al. 2021; Weibel et al. 2024; Harvey et al. 2025).

In terms of observations, the DELPHI23 model is in accordance with the observationally inferred SMF at z ∼ 5 − 9. Within error bars, we find the SMF from the fiducial model to also be in accordance with observations at all z ∼ 5 − 9. It does, however, slightly underestimate the bright end of the SMF at z ∼ 8 at M_∗ ∼ 10^10.25 M_⊙ (Stefanon et al. 2021) and the observations from Kikuchihara et al. (2020) at z ∼ 8 − 9; the latter, however, lie a factor ∼3 above observations from other groups at these redshifts (Song et al. 2016; Bhatawdekar et al. 2019; Stefanon et al. 2021; Navarro-Carrera et al. 2024). Given its increasing star formation efficiency above z ∼ 7, the eSFE model sets the upper limit to the SMF. This globally makes it closer to observations, except for Song et al. (2016) at z ∼ 8 and Stefanon et al. (2021) at z ∼ 9. As the lower mass limit, the eIMF model underestimates the bright end of the observed SMF at all z ∼ 8 − 9 from, for instance, Stefanon et al. (2021) and Weibel et al. (2024), although it is in excellent agreement with the results from Stefanon et al. (2021) at z ∼ 9 at the low-mass end. We caution against a direct comparison of the eIMF model against these observations directly since they all employ a constant Kroupa IMF; assuming an increasingly top-heavy IMF (as in the eIMF) model would naturally result in lower stellar masses being inferred from observations.

3.4. The dust-to-stellar mass relations in early systems

We now study the redshift evolution of the dust-to-stellar mass relation, for our four models, at z ∼ 5 − 17, as shown in Fig. 10. In the DELPHI23 model, low-mass galaxies with M_* ≲ 10^8.5 (10⁹) M_⊙ at z ∼ 15 (7) form stars at the ‘ejection’ efficiency, meaning that any episode of star formation ejects the remaining gas and dust, that are assumed to be perfectly mixed, from the system. This results in low-mass systems bringing in little to no gas and dust into larger successors at later times. However, at every step, the processes of star formation and dust are very closely related (Eqs. (4)–(5)), which results in a linear M_d − M_* relation at any redshift with M_d ∼ 10^−3.1M_*.

Fig. 10. Evolution of the dust mass as a function of stellar mass at z ∼ 5 − 17 for the different models studied in this work, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational results from the ALPINE survey (Fudamoto et al. 2020) at z ∼ 5 and from the REBELS survey (Bouwens et al. 2022a) at z ∼ 7, as marked. The thick salmon lines in the first three panels show theoretical results from Popping et al. (2017).

In contrast, the other models used in this work that use the SPHINX²⁰ simulations (fiducial, eIMF and eSFE), have lower star formation efficiencies. This results in a situation where even low-mass galaxies do not eject all their gas and dust and can therefore contribute to enriching larger systems at later times. However, the relation between star formation and the key dust processes ensures an almost linear M_d − M_* relation in all models. The key differences between the fiducial and eIMF, and DELPHI23 models are that: firstly, the dust masses have non-zero values down to the lowest stellar masses, as compared to the DELPHI23 model where much larger systems (∼ M_∗ ∼ 10^8−8.5 M_⊙) are effectively dust-free due to SNII feedback. Secondly, due to lower amounts of gas and dust being lost in outflows, these relations show shallower slopes for the M_d − M_* relation. For example, these new models show dust masses that are about 0.5 dex higher than the DELPHI23 model for M_∗ ∼ 10^9.5 M_⊙ at z ∼ 5 − 7; this difference decreases with increasing redshift such that these results converge (within the dispersions for different random seed runs) for all but the lowest mass systems. Interestingly, the eIMF model results in dust masses close to the fiducial model, even though the metal yields are IMF-dependent. This is because the metal yields change significantly only for extremely low metallicity galaxies (see Fig. 5); as soon as galaxies get metal-enriched, the yields become similar to the fiducial model. As for the eSFE model, at z ∼ 5 − 7, it effectively looks the same as the fiducial model, by construction. However, its behaviour becomes much more complex with increasing redshift. Above z ∼ 12, we notice three regimes in the dust mass – stellar mass relation of the eSFE model. At masses below 10⁷ M_⊙, galaxies are assigned properties from SPHINX²⁰ halos rendering the relations almost identical to the fiducial model. For M_∗ ∼ 10⁷⁻⁹ M_⊙, however, the SFE is boosted compared to the SPHINX²⁰ values; these galaxies are feedback limited, resulting in them being dust-free. Finally, at even larger masses, while the SFE is still boosted above fiducial values, the host halos are massive enough to hold on to a fraction of gas and dust; this model therefore sets the lower limit to the dust-stellar mass relation at z ≳ 12.

In the first two panels of Fig. 10 we also plot the dust and stellar masses based on ALMA observations of FIR emission, from the ALPINE program at z ∼ 5 (Fudamoto et al. 2020) and the REBELS program at z ∼ 7 (Bouwens et al. 2022a). As seen, for the galaxies being observed at z ∼ 5 with M_∗ ∼ 10⁹⁻¹¹ M_⊙, our new models yield dust masses M_d ∼ 10^7−8.5 M_⊙, or M_d ∼ 10^−2.4M_*, which are in excellent accordance with the ALPINE observations (apart from the most massive system that shows a suppressed dust mass). In terms of REBELS systems, with M_∗ ∼ 10^9−10.5 M_⊙, our model yields values of M_d ∼ 10^6.8−7.5 M_⊙ that are again in excellent accordance with observational results. Finally, with their lower star formation rates (and feedback), these new models are also in accordance with the results from the SANTACRUZ semi-analytic model Popping et al. (2017). Overall, given their large error bars and the scatter in the data, all of the new models are equally compatible with currently available observations at z ∼ 5 − 7. It is worth mentioning two key caveats involved in these observations: firstly, most of these sources are detected in a single ALMA band, requiring the assumption of a dust temperature in order to obtain a dust mass (see discussion in Sommovigo et al. 2022). Furthermore, the assumed star formation history can significantly affect the inferred stellar masses (see e.g. Topping et al. 2022). Overall, within dispersions, the new theoretical models explored here do not show any sensible difference in the dust-to-stellar mass relation at z ∼ 5 − 12.

3.5. The mass-metallicity relation and its redshift evolution

We now compare the stellar mass-gas-phase metallicity relation of the four models discussed in this work at z ∼ 5 − 17, as shown in Fig. 11. The metallicity is expressed in standard units in of 12 + log(O/H) where we use a solar metallicity value of 12 + log(O/H)_⊙ = 8.69 (Asplund et al. 2021).

Fig. 11. Gas-phase metallicity as a function of stellar mass at z ∼ 5 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. At z ∼ 5 − 9, orange pentagons show observational results from Curti et al. (2024). At z ∼ 7, grey hexagons show observations of the low-mass end of the mass-metallicity relation from Chemerynska et al. (2024).

As seen in the previous sections, low-mass galaxies (typically M_∗ ≲ 10⁸ M_⊙ at z ∼ 5 − 9) lose all of their baryons due to SNII feedback in the DELPHI23 model. As a result, galaxies below these mass limits are devoid of metals, showing a sharp cut-off in the mass-metallicity relation. For larger masses, the metallicity shows a slight increase with increasing mass – for example, at z ∼ 5 the metallicity value increases from 8.3 to 8.6 as M_* increases from 10^8.5 to 10¹¹ M_⊙; the overall normalisation of metallicity drops by about 0.2 dex between z ∼ 5 and 9. This model shows the largest metallicity values for high stellar masses at z ∼ 5 − 9 given it has the highest values of the star formation efficiency at these redshifts.

In the new models explored in this work (fiducial, eSFE, eIMF), even systems with masses as low as M_* ∼ 10⁶ M_⊙ are not feedback limited and show non-zero metallicities; as expected, the metallicity increases as a function of the stellar mass in all of these models. The fiducial model shows a metallicity value of about 7.0 (8.25) for M_∗ ∼ 10⁷ (10^9.5−11.5) M_⊙ at z ∼ 5. The normalisation of the mass-metallicity relation decreases slightly with redshift – for example, for a typical stellar mass of M_* ∼ 10⁹ M_⊙, the gas-phase metallicity has values of about 8.0, 7.75, and 7.65 at z ∼ 5, 9 and 12.

In the eIMF model, at z ∼ 5 − 7, low-mass systems (M_* ≲ 10⁸ M_⊙) show slightly elevated values of the mean metallicity (by about 0.3 dex) given their top-heavy IMFs; they converge to the fiducial model for more massive systems. Given an increasingly top-heavy IMF (and the associated larger metal yields), this model shows an increase in the mean amplitude of the mass-metallicity relation by about 0.2 (0.4) dex at z ∼ 9 (12), although the slopes remain similar. At z ∼ 12, for example, systems with M_∗ ∼ 10⁹ M_⊙ show metallicity values of about 8.05.

The eSFE model shows a similar behaviour to that discussed in the previous section – at z ∼ 5 − 7, this model yields results in good agreement with the fiducial model. However, given its cold gas fractions and star formation efficiencies that increase with increasing redshifts, high-mass systems (M_∗ ∼ 10⁹ M_⊙) at z ≳ 12 show the maximum metallicity values in this model; indeed, by z ∼ 9, massive systems (M_* ≳ 10⁸ M_⊙) show about 0.5 dex higher metallicities in this model compared to the fiducial one. Quantitatively, for a system with M_∗ ∼ 10⁹ M_⊙, the metallicity value increases from about 8.1 at z ∼ 9 to 8.4 by z ∼ 12 in this model. Further, while galaxies with M_∗ ∼ 10^7.5−8 M_⊙ are completely feedback limited, lower mass galaxies still contain gas due to the low star formation efficiency values from SPHINX²⁰ at z ≳ 12.

We also compare our results with recent observational results from the JWST at z ∼ 5 − 9 (Curti et al. 2024; Chemerynska et al. 2024), as shown in the same figure. At z ∼ 5, within the dispersions associated with the models, the observational values from Curti et al. (2024) cover the parameter space allowed by all of the models roughly equally well. However, possibly due to our assumption of the instantaneous recycling approximation, we are unable to reproduce the high metallicity values (12 + log(O/H) > 8) inferred observationally for sources with M_∗ ∼ 10^7.5−8.25 M_⊙ (see discussion in Sect. 4.1; Ucci et al. 2023). The situation is similar at z ∼ 7 where, within the errors and dispersions associated with the data and theoretical models, respectively, the mass-metallicity relation points from both Curti et al. (2024) and Chemerynska et al. (2024) agree equally well with all of the new models. Overall here too, within dispersions, the new theoretical models explored do not show any sensible differences in the mass-metallicity relation or its redshift evolution at z ≳ 5.

3.6. Beta slopes for early systems

We now show the relation between the observed (dust attenuated) UV spectral slope (β) and the observed UV magnitude, also called the β − M_UV relation, in Fig. 12. The intrinsic spectra for each galaxy is convolved with the Calzetti dust extinction law (Calzetti et al. 2000) using the UV dust attenuation values obtained in Sect. 2.5 – we then obtain the β value by fitting a power-law at nine (equally spaced) wavelengths between 1300 − 2900 Å in the rest-frame of the galaxy. We start by noting that the intrinsic (i.e. dust-unattenuated) UV spectral slope is rather insensitive to the UV magnitude and all four models show values of β ∼ −2.3 − −2.5 at z ∼ 5 − 7 (not shown in the figure, for simplicity) as might be expected for stellar populations that are a few tens of Myr old with metallicity values of a few tenths of the solar value. The observed slopes, on the other hand, become ‘redder’ with an increase in the dust enrichment and/or average stellar ages, meaning that the more massive or luminous a system, the redder is the β slope, for all models. As seen from Sect. 3.4, the DELPHI23 model has the lowest amplitude of the dust mass-stellar mass relation for low-to-intermediate mass objects at z ∼ 5 − 7 driven by its different prescriptions for the gas available for star formation and the star formation efficiency. As a result, low-to-intermediate luminosity galaxies in this model show the bluest slopes with β ∼ −2.3 at M_UV ∼ −20. Given their different assembly histories and overall larger dust masses, the new models in this work (fiducial, eIMF and eSFE) show redder slopes with β ∼ −2 − −2.3 at M_UV ∼ −20; the results from these three models start converging at z ≲ 9 as might be expected. At z ∼ 5 − 9, the most massive systems in the DELPHI23 model show the reddest β slopes as a result of their older stellar populations convolved with high dust masses. At z ≳ 12 the β slopes from all models start converging to values of about −2.4 at M_UV ∼ −20 as a result of younger stellar populations and low dust attenuations. Even for the most massive/luminous systems, the β values from the new models are bluer than −2.

Fig. 12. Ultraviolet spectral slopes (β), using the Calzetti dust extinction law, as a function of the observed UV magnitude at z ∼ 5 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Grey squares in the top left panel represent observational data from 677 galaxies at z ∼ 5, taken from Bouwens et al. (2014). As marked in other panels, we also show observational estimates of the β slope from recent JWST studies at z ∼ 7 − 17 (Castellano et al. 2023; Austin et al. 2024; Franco et al. 2024; Kokorev et al. 2024; Rinaldi et al. 2024).

Comparing our results to observations, we find the new models to yield a β − M_UV relation that is in excellent agreement with both the slope and amplitude of the (binned) results from Bouwens et al. (2014) at z ∼ 5. Within the current observational error bars, essentially all models agree with results from recent JWST observations at z ∼ 7 − 17 (from Castellano et al. 2023; Austin et al. 2024; Franco et al. 2024; Kokorev et al. 2024; Rinaldi et al. 2024; the confirmation of extremely blue systems (with β ≪ −2.5) would pose a challenge to all of the models discussed here. Further, we note that a caveat that enters these comparisons is that the observed β slopes, derived from photometry, are affected both by emission lines and nebular continuum, both of which are currently not included in our modelling.

3.7. Balmer breaks in the first billion years

We now show the Balmer break strength (the ratio of luminosities at 4200 and 3500 Å) as a function of the observed UV magnitude in Fig. 13. The Balmer break strength depends on the age of stellar populations (with younger populations having weaker breaks), the dust attenuation and the IMF (with more massive stars showing weaker breaks). To compare our models with observations, we compute the Balmer break of all our galaxies assuming the Calzetti attenuation law. At z ∼ 5, the Balmer break strength increases with increasing luminosity in the eIMF model, from ≲1 for low-luminosity systems (M_UV ≳ −18) to ∼1.6 for the brightest systems given the larger fraction of low-mass stars in the former; for the largest systems, the trend is very similar to that shown by the other three models considered. As the IMF becomes increasingly top heavy with increasing redshifts, the Balmer break strength-M_UV relation shows a progressive flattening with increasing redshift. Indeed, the Balmer break strength saturates to a value of about 0.8 − 0.9 at ≳10. Given their similar IMFs, the fiducial and eSFE models show very similar values of the Balmer break strength at z ≲ 8; while the mean values diverge for the most luminous systems, due to the different assembly histories, these differences lie within the dispersions of these two models. Finally, related to its reddest β slopes for luminous systems, the DELPHI23 model shows the largest Balmer break value for these systems. By z ∼ 14, however, the results from the fiducial, eSFE and DELPHI23 models converge to values ∼1.1, about 20% higher than those from the eIMF model.

Fig. 13. Balmer break strength as a function of the observed UV magnitude at z ∼ 6 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational results from Vikaeus et al. (2024).

We also compare our results with recent JWST observations presented in Vikaeus et al. (2024) as shown in the same figure. Given the dispersions both in the observed Balmer break strengths for a given magnitude and the different models used here, as of now, it is hard to infer a particular model being preferred. One might have a slight preference for the eIMF model to explain the very low values of the strength (≲0.9) observed at all z ∼ 6 − 10. We again stress the lack of any nebular continuum in our modelling, which can significantly impact the inferred value of the Balmer Break (Katz et al. 2024). Future datasets will be crucial in building larger statistics to shed light on the stellar populations using this statistic.

3.8. The redshift evolution of the ionising photon production efficiency

A common way of studying the ability of galaxies to contribute to reionisation is through the ionising photon production efficiency, ${\xi }_{\mathrm{ion}}={\dot{N}}_{\mathrm{ion}}^{\mathrm{sf}}/L_{\mathrm{UV}}^{\mathrm{int}}$ where ${\dot{N}}_{\mathrm{ion}}^{\mathrm{sf}}[{\mathrm{s}}^{-1}]$ is the intrinsic production rate of ionising photons from star formation and $L_{\mathrm{UV}}^{\mathrm{int}}[\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}]$ is the intrinsic UV luminosity. This quantity strongly depends on the properties of the underlying stellar population including the IMF and the mean stellar age. We show ξ_ion as a function of UV magnitude for z ∼ 5 − 17 in Fig. 14. As might be expected, the fiducial, eSFE and DELPHI23 models yield very similar results for this quantity given they all use the same IMF. At z ∼ 5 − 9, in these models, we find mean values of log(ξ_ion)∼25.1 − 25.2 [Hz erg⁻¹]. However, producing many more massive stars, galaxies in the eIMF model show a much higher value of ξ_ion: for example, at z ∼ 5 faint galaxies with M_UV ≳ −17 show values as high as log(ξ_ion) = 25.55 [Hz erg⁻¹]. ξ_ion decreases with increasing brightness and metallicities and converges to all the other models by M_UV ∼ −21. With increasing redshifts, even the brightest systems have flatter IMFs resulting in ξ_ion flattening with UV magnitude. Crucially, in the eIMF model galaxies produce 2 − 2.5× as many ionising photons per unit UV luminosity as compared to all the other models at z ∼ 12 − 17.

Fig. 14. Ionising photon production efficiency (ξion) as a function of the observed UV magnitude. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational data as marked (from Simmonds et al. 2024; Llerena et al. 2024; Begley et al. 2025). The crosses at the top right corner of the upper panels indicate the average error on measured magnitudes and ξion as estimated by Simmonds et al. (2024).

We also compare these results with recent JWST observations of the ξ_ion − M_UV relation at z ∼ 5 − 9 (from Begley et al. 2025; Llerena et al. 2024; Simmonds et al. 2024) as shown in Fig. 14. At z ∼ 5 − 9, the observations from Simmonds et al. (2024) show log(ξ_ion) values ranging between 24.5 − 25.8 for systems with M_UV ∼ −16 − −20 and do not show any specific trend as a function of the UV magnitude. This is in accordance with the results from Begley et al. (2025) at z ∼ 7 that range between log(ξ_ion)∼25 − 25.8 for M_UV ∼ −18 − −20 and also do not show any specific trend with the UV magnitude – both of these observations are compatible with all of the four models considered here. On the other hand, at z ∼ 5 − 7, the average results from Llerena et al. (2024) clearly show a decreasing ξ_ion with increasing brightness – for example, their average values of log(ξ_ion) decrease from about 25.6–25 as M_UV decreases from ∼ − 18 to ∼ − 22 at z ∼ 5. Interestingly, this decrease is in excellent accordance with the eIMF model, where low-mass systems produce copious amounts of ionising photons, at z ∼ 5 and reasonably compatible at z ∼ 7. A number of other recent works have also obtained constraints on ξ_ion: Atek et al. (2024) find log(ξ_ion) = 25.8 ± 0.14 for galaxies with M_UV ≳ −16.5 at z ∼ 6 − 7.7, which is a factor four higher than older works (e.g. Robertson et al. 2013) and a factor 1.6 higher than the values from any of our models. Further, Rinaldi et al. (2024) find $log({\xi }_{\mathrm{ion}})=25.{50}_{-0.12}^{+0.10}$ for Hα emitters at z ∼ 7 − 8, which is in good agreement with results from the eIMF model. Finally, Curtis-Lake et al. (2023) derive $log({\xi }_{\mathrm{ion}})=25.{04}_{-0.21}^{+0.17}$ for a source with M_∗ ∼ 10^7.58 M_⊙ at z ∼ 10.38 and $log({\xi }_{\mathrm{ion}})=24.{73}_{-0.15}^{+0.25}$ for another object with with M_∗ ∼ 10^8.67 M_⊙ at z ∼ 11.58. These values are lower than our eIMF model by a factor 3.5 and are in much better agreement with the fiducial and eSFE models. Overall, given the variation in observational results, as of now, they do not favour any specific model presented in this work apart from the observations by Llerena et al. (2024) that seem to show a preference for the eIMF model.

4. The epoch of reionisation in light of JWST data

We show the intrinsic and escaping ionising photon emissivity in all of the models presented in this work before discussing the key reionisation sources, and its progress, in this section.

4.1. The intrinsic and escaping ionising photon emissivities

We start by discussing the intrinsic emissivity of hydrogen ionising photons, weighted over all galaxies, at z ∼ 5 − 21, as shown by the solid lines in Fig. 15. The fiducial and eSFE models show results in good accordance with each other, with the emissivity increasing from about 10^47.7 at z ∼ 21 to 10^51.8 [s⁻¹ Mpc⁻³] by z ∼ 5; given its slightly different star formation efficiency, the DELPHI23 model shows slightly higher values (by about 0.2 dex) at z ≳ 14, converging to the other models at lower redshifts. In contrast, the eIMF model produces a much larger emissivity than the other models at all redshifts, thanks to its flatter IMF slope. Driven by a combination of high redshifts and low metallicities, it results in an emissivity value that is about 4× larger than all the other models at z ≳ 8; the values are a factor ∼1.8 higher by z ∼ 5, driven by low-mass (M_* ≲ 10⁹ M_⊙), low-metallicity systems.

Fig. 15. Evolution of the ionising photon rate density as a function of redshift. In both panels, for reference, the solid lines represent the intrinsic ionising photon density for the different models studied, as marked. Shaded regions highlight the spread from the minimum to maximum values among the five runs with different random seeds. The two panels show the ‘escaping’ emissivity for the two fesc models studied in this work for all galaxies (dashed lines) and those with M* > 109 M⊙ (dotted lines). The left panel shows results for a constant value of fesc needed to match reionisation constraints; the right panel shows results using fesc values from Chisholm et al. (2022) as shown in Eq. (9). In both panels, points show observational constraints (from Becker & Bolton 2013; Bouwens et al. 2015; Becker et al. 2021; Rinaldi et al. 2024).

We also calculate the escaping emissivity using two models for the escape fraction: (i): in the first, we assume a value of f_esc that is constant for all galaxies at all redshifts. This value is obtained by matching to the observed evolution of the volume filling fraction of neutral hydrogen as detailed in what follows; (ii): our second model is driven by results from the Low-redshift Lyman Continuum Survey (LzLCS; PI: Jaskot, HST Project ID: 15626, Flury et al. 2022). From these observations, f_esc depends on the β slope as (Chisholm et al. 2022)

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{chisholm}}=(1.3\pm 0.6)\times {10}^{-4}\times {10}^{(-1.22\pm 0.1)\beta }.\end{array}$$(9)

We compute the f_esc value for each galaxy in each model considered using the β values computed in Sect. 3.6.

A complication is introduced by the fact that the process of reionisation generates a heating UV background (UVB) that can suppress the baryonic content of low-mass halos in ionised regions (e.g. Gnedin 2000; Sobacchi & Mesinger 2013; Hutter et al. 2021). To model this feedback effect, we run our models assuming a complete suppression of the gas mass (and hence star formation) in all halos below M_h = 10⁹ M_⊙ in ionised regions. The total ‘emerging’ ionising emissivity for reionisation is obtained by weighing over the UV-suppressed contribution from low-mass halos in ionised regions as well as that from unsuppressed sources in neutral regions, such that (Dayal et al. 2025)

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{ion}}={\displaystyle \sum {\varphi }_i}(Q_{\mathrm{II}}{\dot{N}}_{\mathrm{ion},\mathrm{II}}^{\mathrm{sf}}(i)f_{\mathrm{esc},\mathrm{II}}(i)+Q_{\mathrm{I}}{\dot{N}}_{\mathrm{ion},\mathrm{I}}^{\mathit{sf}}(i)f_{\mathrm{esc},\mathrm{I}}(i)).\end{array}$$(10)

Here, Q_II and Q_I represent the volume filling fractions of ionised and neutral hydrogen, respectively. Further, ϕ_i is the number density associated with any halo and the first and second terms on the right-hand side show the contribution from suppressed sources in ionised regions and that from unsuppressed sources in neutral regions, respectively.

The evolution of Q_II can be calculated as (Madau 2017)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{Q}_{\mathrm{II}}}{\mathrm{d}t}=\frac{{\dot{\mathrm{n}}}_{\mathrm{ion}}}{⟨n_{\mathrm{H}}⟩(1+\langle {\kappa }_{{\nu }_L}^{\mathrm{LLS}}\rangle /\langle {\kappa }_{{\nu }_L}^{\mathrm{IGM}}\rangle )}-\frac{Q_{\mathrm{II}}}{{\overline{t}}_{\mathrm{rec}}}·\end{array}$$(11)

Here, ⟨n_H⟩ is the mean hydrogen density of the Universe. The terms ⟨κ_νL^LLS⟩ and ⟨κ_νL^IGM⟩ are the volume averaged absorption probabilities per unit length due to Lyman limit systems and the uniform IGM, respectively, and are given by Equations (11) and (13) of Madau (2017). Finally, ${\overline{t}}_{\mathrm{rec}}$ is the effective recombination timescale in the IGM, and is equal to [(1+χ)⟨n_H⟩α₀C_R]⁻¹, where χ = 0.083 accounts for the presence of electrons from He II, α₀ is the case-A recombination coefficient assuming T = 10^4.3 K and C_R = 2.9[(1+z)/6]^−1.1 is the clumping factor (Shull et al. 2012). We solve Equations (10) and (11) jointly using a backwards differentiation formula method with SCIPY’s OdeSolver², as in Trebitsch et al. (2022).

We show the resulting escaping emissivities of ionising photons as dashed lines also in Fig. 15 for both f_esc models considered in this work. As discussed in what follows, in the constant f_esc scenario, matching to the redshift evolution of Q_I requires similar average values of f_esc ∼ 15.2% and 15.4% for the fiducial and eSFE models, respectively. With its larger star formation efficiency value, the DELPHI23 model requires a slightly smaller value of f_esc ∼ 13.1%. However, as discussed above, at almost all z ∼ 6 − 21, galaxies in the eIMF model produce about three times as many ionising photons per unit volume as compared to the fiducial model; compensating for this requires a much lower value of f_esc ∼ 5.7% in this model in order to match the observed constraints. With these values, we obtain essentially similar evolutions of the escaping emissivities in all four models, which match observed emissivities (from Becker & Bolton 2013; Bouwens et al. 2015; Becker et al. 2021; Rinaldi et al. 2024). Further, low-mass systems (M_* ≲ 10⁹ M_⊙) provide ∼85% of ionising photons down to the midpoint of reionisation at z ∼ 7, dominating the reionisation budget. The contribution of higher mass systems increases with decreasing redshift from about 15% at z ∼ 7 to about 30% by the end of reionisation, as shown in Fig. 15. In the model considering f_esc values using the Chisholm et al. (2022) relations, the escaping emissivities again increase with decreasing redshift. While the fiducial, eSFE and DELPHI23 models show very similar results, the overall higher ξ_ion values in the eIMF model lead to a three times higher escaping emissivity at any redshift. Given their larger dust contents (i.e. redder β slopes), in this case too, high-mass systems (M_* > 10⁹ M_⊙) have a sub-dominant contribution to the escaping emissivity and contribute at most 15% to the process at any z ∼ 5 − 21. Hence, in our model, low-mass galaxies are the key drivers of the reionisation process, even accounting for the impact of reionisation feedback.

4.2. The key sources and progress of reionisation

We then explore the cumulative contribution of galaxies of different stellar masses to the escaping emissivity of ionising photons in Fig. 16. For the sake of clarity, we omit the DELPHI23 model and limit to showing results for a constant f_esc value; we note that the results remain qualitatively unchanged when using the f_esc relation from Eq. (9). At z ∼ 10, in the fiducial and eSFE models, 50% of ionising photons come from low-mass systems with M_∗ ≲ 10^6.75 M_⊙ and M_* ≲ 10⁷ M_⊙ respectively; 90% of escaping ionising photons can be traced to systems with M_∗ ≲ 10⁸ (10^8.5) M_⊙ in the fiducial and eSFE models, respectively. With its larger production rate of ionising photons per unit stellar mass, a larger fraction of the emissivity comes from lower-mass galaxies in the eIMF model: for example, at z ∼ 10, 50% (90%) of the total ionising photons come from systems with M_∗ ≲ 10^6.3 (10^7.8) M_⊙. By z ∼ 7, the fiducial and eSFE model are almost indistinguishable with ∼50%(90%) of the emissivity being provided by galaxies with M_∗ ≲ 10⁸ (10^9.4) M_⊙. In the eIMF model, by z ∼ 7, the emissivity again shifts to slightly lower-mass systems: ∼50%(90%) of the emissivity being provided by galaxies with M_∗ ≲ 10^7.6 (10^9.1) M_⊙. To summarise, in any of the models discussed here, galaxies with M_* ≲ 10⁹ M_⊙ provide the bulk (∼75%) of the escaping ionising emissivity for reionisation in agreement with our previous works (e.g. Dayal et al. 2025).

Fig. 16. Cumulative emissivity of ‘escaping’ ionising photons as a function of stellar mass for the new models explored in this work (fiducial, eIMF and eSFE), as marked. Solid and dashed lines show results at z ∼ 7 and 10, respectively. For the sake of clarity, we limit results to our constant fesc models, matched to reionisation constraints (see Sect. 4).

We then discuss the progress of reionisation through the evolution of Q_I in Fig. 17 for both a constant f_esc (left panel) and the f_esc − β relation detailed in Eq. (9) (right panel). In the constant f_esc model, that has been tuned to reproduce a reionisation history matching observational constraints including on the tail end of the process (Jung et al. 2020; Gaikwad et al. 2023; Nakane et al. 2024; Umeda et al. 2024), reionisation starts at z ∼ 13, is 50% complete by z ∼ 7 and finishes by z ∼ 5.5. In the second model, that uses $f_{\mathrm{esc}}^{\mathrm{chisholm}}$, while the fiducial, eSFE and DELPHI23 models match constraints at z ≳ 7, they predict a later end of the reionisation process (by Δz ∼ 0.2) compared to z ≲ 6 data (Gaikwad et al. 2023); we explore different threshold masses for reionisation feedback ranging between 10^8.5−10^9.5 M_⊙, which do not affect our results sensibly. Finally, with its much higher escaping emissivity, the eIMF model leads to a very fast progress of reionisation. In this case, reionisation is 50% over at z ∼ 8.5, compared to at z ∼ 7 for the other models, and finishes between z ∼ 5.8 − 6.8 depending on the exact halo mass below which sources are suppressed due to reionisation feedback. We conclude that the escape fraction prescription of Chisholm et al. (2022), based on UV β slopes, is not suitable to produce realistic reionisation histories within our framework. A fixed global escape fraction for all halos better reproduce observations. It is possible that the relation between β slopes and f_esc evolve with redshift, which would make the low redshift fits of Chisholm et al. (2022) invalid during the epoch of Reionisation. Another possibility is that the lack of nebular continuum and the simple dust geometry of our model does not produce fully realistic β slopes.

Fig. 17. Volume filling fraction of neutral hydrogen as a function of redshift. In the left panel, we carry out calculations assuming a constant global escape fraction (fesc), as marked, given by the best fit to a selection of observational data (Jung et al. 2020; Gaikwad et al. 2023; Nakane et al. 2024; Umeda et al. 2024). Shaded areas represents the range of values obtained with the five runs using different random seeds. In the right panel, we show results using the escape fractions as in Chisholm et al. (2022), using our beta slopes as presented in Sect. 3.6. Solid, dashed, and dotted lines show results for the different halo masses (109, 109.5 and 108.5 M⊙) below which the baryonic content is completely suppressed due to the UVB. For clarity, shaded area are omitted in this panel, but would be similar as in the left one.

5. Summary and discussion

We use the DELPHI semi-analytic framework to study the galaxy populations being detected by the JWST in the first billion years of the Universe. Crucially, we extract ISM properties – namely the cold gas fraction and star formation efficiency – from SPHINX²⁰, a state-of-the-art radiation hydrodynamics simulation (Rosdahl et al. 2022). Both the cold gas fraction (f_cold) and star formation efficiency (f_*^sphinx) values are mass- and redshift-dependent. For example, at z ∼ 9 while a bulk (∼80%) of low-mass halos (M_h ≲ 10^8.4 M_⊙) show suppressed values of both (f_cold≲0.18; f_*^sphinx≲0.05), both quantities increase with halo mass such that half of the most massive systems (with M_h ≳ 10¹⁰ M_⊙) show f_cold ∼ 0.2 − 0.4 and f_*^sphinx ∼ 0.05 − 0.4. This is a key first step in implementing physically motivated and tested prescriptions into semi-analytic models. In addition to inducing stochasticity in star formation, the fraction of gas in cold and warm phases is crucial for the dust calculations (ISM growth and dust destruction) included in our model that are base-lined against the latest ALMA results.

While the augmentation with stochastic star formation allows us to reproduce the redshift evolution of the UV LF up to z ∼ 10, this fiducial model under-predicts the bright end of the z > 10 redshift Universe as probed recently by JWST. We explore two different approaches to match this observed overabundance of bright sources in the first billion years: in the eIMF model, the gas-phase metallicity- and redshift-dependent IMF becomes increasingly top-heavy at higher redshifts and lower metallicities as in previous works (e.g. Cueto et al. 2024). In the eSFE model, on the other hand, both the cold gas fractions and the star formation efficiencies of massive halos increase with increasing redshift; we maintain a Kroupa IMF as in the fiducial model. Both these models reproduce the ‘dust attenuated’ (observed) UV LF from z ∼ 5 out to z ∼ 15 − 20. Besides the UV LF, our models predict a wide range of observables at z ∼ 5 − 20 including the stellar mass function, the dust mass-stellar mass relations, mass-metallicity relations, the UV spectral slopes (β), the Balmer break strength and the ionising photon production efficiency. Our key findings are:

• Dust attenuation plays an important role in determining the bright end (M_UV ≲ −21) of the UV LF at z ≲ 11 in all of the models studied in this work. The impact of dust attenuation decreases with increasing redshift – this is due to dust and gas being more dispersed (covering a larger fractional volume of the host halo) with increasing redshift, leading to a decrease in the dust optical depth. Dust plays no role in attenuating the UV luminosity at z ≳ 11 − 12.

• The varying mass-to-light ratios in the different models explored in this work yield stellar mass functions that differ significantly, specially at the massive end at high redshifts; within observational dispersions, all models are in reasonable agreement with available data. For example, the number density associated with stellar mass of M_* ∼ 10⁹ M_⊙ at z ∼ 12 varies by more than 2.5 orders of magnitude between the upper limit set by the eSFE model and the lower limit set by the eIMF model. A conjunction of the UV LF and SMF will be crucial in demarcating models predicting the correct mass-to-light ratios.

• The mass-metallicity relation is in place as early as z ∼ 17 in all models. The fiducial model shows a metallicity value of about 2% (36%) Z_⊙ for M_∗ ∼ 10⁷ (10^9.5−11.5) M_⊙ at z ∼ 5. The eIMF model shows an increase in the mean amplitude of the mass-metallicity relation by about 0.2 (0.4) dex at z ∼ 9 (12) compared to the fiducial model, although the slopes remain similar. At z ≳ 12, the eSFE model shows the largest metallicity values.

• We find values of β ≲ −2.3 for systems fainter than UV magnitudes of about −20. Massive systems at z ∼ 5 reach β slopes as red as ∼ − 0.75, but become increasingly bluer with increasing redshift.

• Galaxies in the eIMF model show the steepest Balmer break strength-M_UV relation at z ∼ 6, increasing from from ∼0.9 to 1.6 as M_UV increases from −16 to −22. All the other models find relatively flat slopes with Balmer break strengths of about 1.4 at z ∼ 5.

• Galaxies in the eIMF model show elevated values of ξ_ion: for example, at z ∼ 5 galaxies with M_UV ≳ −17 show values as high as log(ξ_ion) = 25.55 [Hz erg⁻¹] compared to values of about 25.2 in all other models; systems also saturate to log(ξ_ion) = 25.6 [Hz erg⁻¹] at z ≳ 12, a factor 2–2.5 times higher than those in all other models. These spectral signatures will be crucial in differentiating between models as more data comes in.

• In all scenarios compatible with reionisation observables, star formation in low-mass galaxies (M_* ≲ 10⁹ M_⊙) provide about 85% of the ‘escaping’ emissivity for reionisation down to the midpoint at z ∼ 7, irrespective of the model considered; higher mass systems contribute at most 30% by the end stages of reionisation. Indeed, at z ∼ 10, 90% of the cumulative escaping ionising photons can be traced to systems with M_∗ ≲ 10^8−8.5 M_⊙ in the different models considered here. The (constant) values of f_esc range between 5.7% for the eIMF model to ∼15% for the eSFE and fiducial models.

Given their dispersions, existing datasets are, as of now, equally compatible with the fiducial, eIMF and eSFE models. Further constraints on mass-to-light ratios (through joint analyses of the UV LF and SMF), spectral signatures (such as Balmer break strengths and the ionising photon production efficiency) and stronger hints on f_esc (e.g. Choustikov et al. 2024; Jaskot et al. 2024; Jung et al. 2024; Kerutt et al. 2024) will be crucial in differentiating between these classes of models.

Finally, our models can be further improved in the future to better reflect on the complex and diverse nature of high redshift galaxies. We have seen, for example, that our models do not recover extremely metal rich, low-mass systems, β slopes bluer than −2.5 or Balmer break strengths weaker than 0.8. These could be driven by, for example, stochasticity in the star formation process where extremely young stars (< 5 Myr) dominate the spectrum or the presence of very massive stars (up to 400 M_⊙) boosting the luminosity of star-forming galaxies (Schaerer et al. 2025). Our metal prescriptions use the instantaneous recycling approximation that could be responsible for the lack of low-mass high-metallicity systems (Sect. 4.1; Ucci et al. 2023). Further we assume perfect mixing of gas and dust and assume a simple dust geometry to calculate the total optical depth, ignoring both clumping of dust and the spatial segregation between star forming regions and dust (Inami et al. 2022). We aim at including these crucial effects into our models in future works to keep pace with the incredible datasets available in the first billion years.

Acknowledgments

VM and PD acknowledge support from the NWO grant 016.VIDI.189.162 (‘ODIN’). PD warmly thanks the European Commission’s and University of Groningen’s CO-FUND Rosalind Franklin program. MGH has been supported by STFC consolidated grants ST/N000927/1 and ST/S000623/1. TK was supported by the National Research Foundation of Korea (RS-2022-NR070872 and RS-2025-00516961) and by the Yonsei Fellowship, funded by Lee Youn Jae. The authors thank S. Gazagnes, A. Hutter, J. Kerutt, C. Leiterer, and M. Trebitsch for their helpful comments and insightful discussions.

References

All Tables

All Figures

	Fig. 1. Halo mass functions obtained from the SPHINX20 and DELPHI models, shown using blue and green histograms, respectively, as marked. The left and right panels show results at z ∼ 5 and 9, respectively.
In the text

	Fig. 2. Fraction of halos with a cold gas fraction (fcold) below the value shown on the horizontal axis from the SPHINX20 simulations, for the mass bins marked, at z = 5 (left panel) and z = 9 (right panel). The cold gas fraction is defined as the fraction of gas mass that is below 1000 K, within the virial radius of the halo.
In the text

Fig. 3. Fraction of halos with a star formation efficiency value (f*sphinx) below the value shown on the horizontal axis from the SPHINX20 simulations, for the mass bins marked, at z = 5 (left panel) and z = 9 (right panel). The star formation efficiency is defined as the ratio between the mass of stars younger than 30 Myr and the cold gas mass; if a halo has no cold gas mass, we define its star formation efficiency to be f*sphinx = 0.

In the text

	Fig. 4. Evolution of the SN rate (ν) as a function of metallicity, at the redshifts marked, for the eIMF model (Sect. 2.6.1). The shape of the IMF in this case is determined by the metallicity and redshift. The purple line shows the supernova rate for the fiducial Kroupa IMF.
In the text

	Fig. 5. Evolution of the stellar metal yields (solid lines) and total yields (metals and gas return fraction – dashed lines) per unit stellar mass as a function of metallicity, for the redshifts marked, for the eIMF model. The dashed purple line shows the metallicity-dependent yield values for a Kroupa IMF used in our other models.
In the text

	Fig. 6. As a function of metallicity, we show the ratio between the ionising photon production rates (intrinsic UV luminosities) in the eIMF and fiducial model using solid (dashed) lines at the marked redshifts.
In the text

Fig. 7. Redshift evolution of the UV LF between z ∼ 5 − 20, as marked. In each panel, we show a comparison of the fiducial model, and the eIMF and eSFE models against the results from DELPHI23 (Mauerhofer & Dayal 2023) for both the intrinsic (dashed lines) and dust attenuated (solid lines) UV LFs. These lines show the means of five runs sampling the PDFs for fcold and f* from SPHINX20. For the sake of clarity we only show the maximum and minimum values associated with these runs for the observed UV LFs through shaded areas. In each panel, as marked, points show observational data (from Adams et al. 2024; Atek et al. 2015, 2018; Bouwens et al. 2021, 2022b, 2023; Bowler et al. 2017; Casey et al. 2024; Castellano et al. 2023; Donnan et al. 2023, 2024; Finkelstein et al. 2024; Harikane et al. 2022, 2023, 2024, 2025; Ishigaki et al. 2018; Leung et al. 2023; McLeod et al. 2024; Naidu et al. 2022a; Napolitano et al. 2025; Pérez-González et al. 2023; Robertson et al. 2024; Willott et al. 2024; Whitler et al. 2025). The last panel shows the UV LF integrated between z ∼ 15 − 20 in order to compare with observations at those redshifts (Kokorev et al. 2024; Whitler et al. 2025).

In the text

Fig. 8. Redshift evolution of the UV luminosity density at z ∼ 5 − 21. Solid lines show the mean of five runs for each model, using dust attenuated values of ρUV integrating down to observed magnitude limits of MUV ≲ −17. Shaded areas show the maximum and minimum values associated with these runs. Dotted lines show values of the intrinsic luminosity density integrating over all galaxies. As marked, points show observational results (Donnan et al. 2023; McLeod et al. 2024; Finkelstein et al. 2024; Whitler et al. 2025) and the faint pink shaded area corresponds to the 16% and 84% marginal constraints from the JADES origin field (Robertson et al. 2024).

In the text

Fig. 9. Redshift evolution of the SMF at z ∼ 5 − 20, for the different model explored in this work. Solid lines show the mean results from five runs of each model with shaded areas showing the associated maximum and minimum values. The bottom right panel shows the SMF integrated on all outputs between redshifts 15 and 20. As marked, each panel up to z ∼ 12 show observational data, all assuming a Kroupa IMF (from Bhatawdekar et al. 2019; Duncan et al. 2014; Kikuchihara et al. 2020; Navarro-Carrera et al. 2024; Song et al. 2016; Stefanon et al. 2021; Weibel et al. 2024; Harvey et al. 2025).

In the text

Fig. 10. Evolution of the dust mass as a function of stellar mass at z ∼ 5 − 17 for the different models studied in this work, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational results from the ALPINE survey (Fudamoto et al. 2020) at z ∼ 5 and from the REBELS survey (Bouwens et al. 2022a) at z ∼ 7, as marked. The thick salmon lines in the first three panels show theoretical results from Popping et al. (2017).

In the text

Fig. 11. Gas-phase metallicity as a function of stellar mass at z ∼ 5 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. At z ∼ 5 − 9, orange pentagons show observational results from Curti et al. (2024). At z ∼ 7, grey hexagons show observations of the low-mass end of the mass-metallicity relation from Chemerynska et al. (2024).

In the text

Fig. 12. Ultraviolet spectral slopes (β), using the Calzetti dust extinction law, as a function of the observed UV magnitude at z ∼ 5 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Grey squares in the top left panel represent observational data from 677 galaxies at z ∼ 5, taken from Bouwens et al. (2014). As marked in other panels, we also show observational estimates of the β slope from recent JWST studies at z ∼ 7 − 17 (Castellano et al. 2023; Austin et al. 2024; Franco et al. 2024; Kokorev et al. 2024; Rinaldi et al. 2024).

In the text

	Fig. 13. Balmer break strength as a function of the observed UV magnitude at z ∼ 6 − 17, as marked. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational results from Vikaeus et al. (2024).
In the text

Fig. 14. Ionising photon production efficiency (ξion) as a function of the observed UV magnitude. Lines represent the mean of the five runs with different seeds for each model; shaded areas represent the range between the minima and maxima shown by the 16th and 84th percentiles, respectively. Points show observational data as marked (from Simmonds et al. 2024; Llerena et al. 2024; Begley et al. 2025). The crosses at the top right corner of the upper panels indicate the average error on measured magnitudes and ξion as estimated by Simmonds et al. (2024).

In the text

Fig. 15. Evolution of the ionising photon rate density as a function of redshift. In both panels, for reference, the solid lines represent the intrinsic ionising photon density for the different models studied, as marked. Shaded regions highlight the spread from the minimum to maximum values among the five runs with different random seeds. The two panels show the ‘escaping’ emissivity for the two fesc models studied in this work for all galaxies (dashed lines) and those with M* > 109 M⊙ (dotted lines). The left panel shows results for a constant value of fesc needed to match reionisation constraints; the right panel shows results using fesc values from Chisholm et al. (2022) as shown in Eq. (9). In both panels, points show observational constraints (from Becker & Bolton 2013; Bouwens et al. 2015; Becker et al. 2021; Rinaldi et al. 2024).

In the text

	Fig. 16. Cumulative emissivity of ‘escaping’ ionising photons as a function of stellar mass for the new models explored in this work (fiducial, eIMF and eSFE), as marked. Solid and dashed lines show results at z ∼ 7 and 10, respectively. For the sake of clarity, we limit results to our constant fesc models, matched to reionisation constraints (see Sect. 4).
In the text

Fig. 17. Volume filling fraction of neutral hydrogen as a function of redshift. In the left panel, we carry out calculations assuming a constant global escape fraction (fesc), as marked, given by the best fit to a selection of observational data (Jung et al. 2020; Gaikwad et al. 2023; Nakane et al. 2024; Umeda et al. 2024). Shaded areas represents the range of values obtained with the five runs using different random seeds. In the right panel, we show results using the escape fractions as in Chisholm et al. (2022), using our beta slopes as presented in Sect. 3.6. Solid, dashed, and dotted lines show results for the different halo masses (109, 109.5 and 108.5 M⊙) below which the baryonic content is completely suppressed due to the UVB. For clarity, shaded area are omitted in this panel, but would be similar as in the left one.

In the text